src/Objects/DataManip/mrcImage/src/lmrcImageMontageCreate.c

   1 /*
   2 # %M% %Y% %I%
   3 # The latest update : %G% at %U%
   4 #
   5 #%Z% lmrcImageMontageCreate ver %I%
   6 #%Z% Created by
   7 #%Z%
   8 #%Z% Usage : lmrcImageMontageCreate
   9 #%Z% Attention
  10 #%Z%
  11 */
  12 #undef DEBUG
  13 #include "genUtil.h"
  14 #include "Memory.h"
  15 #include "lmrcImageMontage.h"
  16 #include "lmrcImageEdgeAverage.h"
  17
  18 static char __sccs_id[] = "%Z%lmrcImageMontageCreate ver%I%; Date:%D% %Z%";
  19
  20 int
  21 lmrcImageMontageCreate(mrcImage* out, mrcImage* in, int imageNum, lmrcImageMontageInfo linfo, int mode)
  22 {
  23         int i;
  24         mrcImage tmp;
  25         mrcImageParaTypeReal x,   y,  z;
  26         mrcImageParaTypeReal sx, sy, sz;
  27         double data;
  28         double* avg;
  29         double avgall;
  30         double  orgValue;
  31         lmrcImageEdgeAverageInfo* llinfo;
  32         lmrcImageEdgeAverageInfo  llinfotmp;
  33         mrcImageParaTypeInteger icenter = imageNum/2; /* Center */
  34         mrcImageParaTypeInteger icenterx;
  35         mrcImageParaTypeInteger icentery;
  36
  37         icenterx = icenter%linfo.nx;
  38         icentery = icenter/linfo.nx;
  39
  40         out->Header    = in[0].Header;
  41         out->HeaderN.x = linfo.PadInfo.Width *linfo.nx;
  42         out->HeaderN.y = linfo.PadInfo.Height*linfo.ny;
  43         out->HeaderN.z = 1;
  44         mrcInit(out, NULL);
  45
  46         for(x=0; x<out->HeaderN.x; x++) {
  47                 for(y=0; y<out->HeaderN.y; y++) {
  48                         for(z=0; z<out->HeaderN.z; z++) {
  49                                 mrcPixelDataSet(out, x, y, z, linfo.PadInfo.Value, mrcPixelRePart);
  50                         }
  51                 }
  52         }
  53
  54         if(linfo.flagEdgeAverage) {
  55                 avg    = (double*)memoryAllocate(sizeof(double)*imageNum, "in lmrcImageMontageCreate");
  56                 llinfo = (lmrcImageEdgeAverageInfo*)memoryAllocate(sizeof(lmrcImageEdgeAverageInfo)*imageNum, "in lmrcImageMontageCreate");
  57                 for(i=0; i<imageNum; i++) {
  58                         llinfo[i].width  = linfo.EdgeAverageWindow;
  59                         llinfo[i].devide = 2;
  60                         DEBUGPRINT1("#### Edge %d ####\n", i);
  61                         lmrcImageEdgeAverage(&(in[i]), &(llinfo[i]), 0);
  62                 }
  63                 switch(imageNum) {
  64                         case 9: {
  65                     /*
  66                                Edge 2
  67                                            d0 d1
  68                                -----
  69                            d1 |     | d1
  70                     Edge 3    |     |    Edge 1
  71                            d0 |     | d0
  72                                -----
  73                                            d0 d1
  74                                Edge 0
  75
  76                     shotID
  77
  78                    6    7    8
  79
  80                                    3    4    5
  81
  82                    0    1    2
  83
  84                                 llinfo[shotID].avg[Edge][Devide]
  85
  86                     */
  87                                 avg[6] = sqrt((llinfo[6].avg[0][1]/(llinfo[3].avg[2][1]/(llinfo[3].avg[1][0]/llinfo[4].avg[3][0])))
  88                                              *(llinfo[6].avg[1][0]/(llinfo[7].avg[3][0]/(llinfo[7].avg[0][0]/llinfo[4].avg[2][0]))));
  89
  90                                 avg[7] = (llinfo[7].avg[0][0]/llinfo[4].avg[2][0]);
  91
  92                                 avg[8] = sqrt((llinfo[8].avg[0][0]/(llinfo[5].avg[2][0]/(llinfo[5].avg[3][0]/llinfo[4].avg[1][0])))
  93                                             * (llinfo[8].avg[3][0]/(llinfo[7].avg[1][0]/(llinfo[7].avg[0][0]/llinfo[4].avg[2][0]))));
  94
  95                                 avg[3] = (llinfo[3].avg[1][0]/llinfo[4].avg[3][0]);
  96
  97                                 avg[4] = 1;
  98
  99                                 avg[5] = (llinfo[5].avg[3][0]/llinfo[4].avg[1][0]);
 100
 101
 102                                 avg[0] = sqrt((llinfo[0].avg[2][1]/(llinfo[3].avg[0][1]/(llinfo[3].avg[1][0]/llinfo[4].avg[3][0])))
 103                              *(llinfo[0].avg[1][1]/(llinfo[1].avg[3][1]/(llinfo[1].avg[2][0]/llinfo[4].avg[0][0]))));
 104
 105                                 avg[1] = (llinfo[1].avg[2][0]/llinfo[4].avg[0][0]);
 106
 107                                 avg[2] = sqrt((llinfo[2].avg[2][0]/(llinfo[5].avg[0][0]/(llinfo[5].avg[3][0]/llinfo[4].avg[1][0])))
 108                                                      *(llinfo[2].avg[3][1]/(llinfo[1].avg[1][1]/(llinfo[1].avg[2][0]/llinfo[4].avg[0][0]))));
 109
 110                                 DEBUGPRINT3("avg[7] %f = [7]avg[0][0] %f / [4].avg[2][0] %f\n", avg[7], llinfo[7].avg[0][0], llinfo[4].avg[2][0]);
 111                                 break;
 112                         }
 113                         default: {
 114                                 fprintf(stderr, "Not supported imageNum for Edge-Averaging: %d\n", imageNum);
 115                                 exit(EXIT_FAILURE);
 116                         }
 117                 }
 118
 119         }
 120         if(linfo.flagNoAverage) {
 121                 avgall = 0;
 122                 for(i=0; i<imageNum; i++) {
 123                         avgall += in[i].HeaderAMean;
 124                 }
 125                 avgall /= imageNum;
 126         }
 127         orgValue = linfo.PadInfo.Value;
 128
 129         for(i=0; i<imageNum; i++) {
 130                 mrcStatDataSet(&(in[i]), 0);
 131                 if(in[i].HeaderN.z>1) {
 132                         fprintf(stderr, "Not supported: Nz=%f\n", in[i].HeaderN.z);
 133                         exit(EXIT_FAILURE);
 134                 }
 135
 136                 if(linfo.flagEdgeAverage) {
 137                         linfo.PadInfo.Value = orgValue*(in[i].HeaderAMean/in[icenter].HeaderAMean)/avg[i];
 138
 139                         DEBUGPRINT2("##### Edge %d center %d####\n", i, icenter);
 140                         DEBUGPRINT2("ImageAverage: center %f  current %f \n", in[icenter].HeaderAMean, in[i].HeaderAMean);
 141                         DEBUGPRINT2("PadInfo : %f -> %f \n", orgValue, linfo.PadInfo.Value);
 142                         DEBUGPRINT2("Edge: %f == %f \n", llinfo[icenter].avg[2][0]/in[icenter].HeaderAMean*orgValue, llinfo[i].avg[0][0]/in[i].HeaderAMean*linfo.PadInfo.Value);
 143
 144                 } else if(linfo.flagNoAverage) {
 145                         DEBUGPRINT("No Avearge");
 146                         linfo.PadInfo.Value = orgValue - (avgall - in[i].HeaderAMean);
 147                 } else {
 148                         /* No Action */;
 149                 }
 150                 lmrcImagePad(&tmp, &(in[i]), &(linfo.PadInfo), linfo.PadMode);
 151
 152 #ifdef DEBUG
 153                 DEBUGPRINT1("##### Edge %d ####\n", i);
 154                 llinfotmp = llinfo[i];
 155                 lmrcImageEdgeAverage(&tmp, &llinfotmp, 0);
 156                 mrcStatDataSet(&tmp, 0);
 157                 DEBUGPRINT2("ImageAverage: original %f  current %f \n", in[i].HeaderAMean, tmp.HeaderAMean);
 158 #endif
 159
 160                 sx = (i%linfo.nx)*tmp.HeaderN.x;
 161                 sy = (i/linfo.nx)*tmp.HeaderN.y;
 162                 DEBUGPRINT3("s(%f %f) %d\n", sx, sy, i);
 163                 for(x=0; x<tmp.HeaderN.x; x++) {
 164                         for(y=0; y<tmp.HeaderN.y; y++) {
 165                                 for(z=0; z<tmp.HeaderN.z; z++) {
 166                                         mrcPixelDataGet(&tmp,    x,    y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 167                                         if(linfo.flagMaxValue) {
 168                                                 if(linfo.MaxValue < data) {
 169                                                         if(linfo.flagValueAssignedToMax) {
 170                                                                 data = linfo.ValueAssignedToMax;
 171                                                         } else {
 172                                                                 data = linfo.PadInfo.Value;
 173                                                         }
 174                                                 }
 175                                         }
 176                                         mrcPixelDataSet( out, sx+x, sy+y, z,  data, mrcPixelRePart);
 177                                 }
 178                         }
 179                 }
 180                 mrcImageFree(&tmp, "main");
 181         }
 182 }
 183