src/Objects/DataManip/mrcImage/src/lmrcImageAverage.c

   1 /*
   2 # lmrcImageAverage : $Revision$
   3 # $Date$
   4 # Created by $Author$
   5 # Usage : lmrcImageAverage
   6 # Attention
   7 #   $Loccker$
   8 #       $State$
   9 #
  10 */
  11 /* $Log$ */
  12
  13 #define DEBUG
  14 #include "genUtil.h"
  15 #include "Memory.h"
  16 #include "./lmrcImageAverage.h"
  17
  18 void
  19 lmrcImageAverage(mrcImage* out, mrcImage* in, int n, int mode)
  20 {
  21         int i;
  22         long num=0;
  23
  24         for(i=0; i<n; i++) {
  25                 lmrcImageAdd(out, &(in[i]), &num);
  26         }
  27         if(num==n) {
  28                 lmrcImageDevidedByReal(out, (double)num);
  29         } else {
  30                 fprintf(stderr, "Something wrong: n %d != num %d\n", n, num);
  31         }
  32 }
  33
  34 void
  35 lmrcImageAverageWithWeight(mrcImage* out, mrcImage* outW, mrcImage* in, mrcImage* inW, int n, int mode)
  36 {
  37         int i;
  38         mrcImage* tmpIn;
  39         mrcImage* tmpOut;
  40         mrcImageParaTypeReal x, y, z, X, Y, Z;
  41         double re, im, reOut, imOut, w, data;
  42
  43         DEBUGPRINT("Start lmrcImageAverageWithWeight\n");
  44         /* Allocation */
  45         outW->Header = inW[0].Header;
  46         mrcInit(outW, NULL); /* All pixel data are zero */
  47
  48         /* preTreatment: Fit Mode of Image to that of Weight */
  49         if(IsFT(&in[0], "in AverageWithWeight", 0)) {
  50                 if(IsFT(&inW[0], "in AverageWithWeight", 0)) {
  51                         /*      in:  FT
  52                                 inW: FT */
  53                         DEBUGPRINT("FTvsFT\n");
  54                         tmpIn  = in;
  55                         out->Header = in[0].Header;
  56                         mrcInit(out, NULL);
  57                         tmpOut = out;
  58                 } else {
  59                         /*      in:  FT
  60                                 inW: Image */
  61                         DEBUGPRINT("FTvsImg\n");
  62                         tmpIn = (mrcImage*)memoryAllocate(sizeof(mrcImage)*n, "in lmrcImageAverageWithWeight");
  63                         for(i=0; i<n; i++) {
  64                                 lmrcImageFFT(&tmpIn[i], &in[i], 0);
  65                         }
  66                         tmpOut = (mrcImage*)memoryAllocate(sizeof(mrcImage), "in lmrcImageAverageWithWeight");
  67                         tmpOut->Header = in[0].Header;
  68                         mrcInit(tmpOut, NULL);
  69                         out->Header = in[0].Header;
  70                 }
  71         } else {
  72                 if(IsFT(&inW[0], "in AverageWithWeight", 0)) {
  73                         /*      in:  Image
  74                                 inW: FT */
  75                         DEBUGPRINT("ImgvsFT\n");
  76                         tmpIn = (mrcImage*)memoryAllocate(sizeof(mrcImage)*n, "in lmrcImageAverageWithWeight");
  77                         for(i=0; i<n; i++) {
  78                                 lmrcImageFFT(&tmpIn[i], &in[i], 0);
  79                         }
  80                         tmpOut = (mrcImage*)memoryAllocate(sizeof(mrcImage)*1, "in lmrcImageAverageWithWeight");
  81                         tmpOut->Header = inW[0].Header;
  82                         tmpOut->HeaderMode = mrcComplexFloatFT;
  83                         mrcInit(tmpOut, NULL);
  84                 } else {
  85                         /*      in:  Image
  86                                 inW: Image */
  87                         DEBUGPRINT("Imgvs\n");
  88                         tmpIn  = in;
  89                         out->Header = in[0].Header;
  90                         mrcInit(out, NULL);
  91                         tmpOut = out;
  92                 }
  93         }
  94
  95         /* Averaging */
  96         if(IsFT(&tmpIn[0], "in AverageWithWeight", 0)) {
  97                 DEBUGPRINT("Weight: FT\n");
  98                 /* Sum */
  99                 for(i=0; i<n; i++) {
 100                         DEBUGPRINT4("%d %d %d %d\n", i, tmpIn[i].HeaderN.x, tmpIn[i].HeaderN.y, tmpIn[i].HeaderN.z);
 101                         for(Z=-tmpIn[i].HeaderN.z/2.0; Z< tmpIn[i].HeaderN.z/2.0; Z++) {
 102                         for(Y=-tmpIn[i].HeaderN.y/2.0; Y< tmpIn[i].HeaderN.y/2.0; Y++) {
 103                         for(X=0;                       X<=tmpIn[i].HeaderN.x/2.0; X++) {
 104                                 mrcPixelDataGet(&tmpIn[i], X, Y, Z, &re,    mrcPixelRePart, mrcPixelHowNearest);
 105                                 mrcPixelDataGet(&tmpIn[i], X, Y, Z, &im,    mrcPixelImPart, mrcPixelHowNearest);
 106                                 mrcPixelDataGet(&inW[i],   X, Y, Z, &w,     mrcPixelRePart, mrcPixelHowNearest);
 107
 108                                 mrcPixelDataGet(tmpOut,    X, Y, Z, &reOut, mrcPixelRePart, mrcPixelHowNearest);
 109                                 mrcPixelDataGet(tmpOut,    X, Y, Z, &imOut, mrcPixelImPart, mrcPixelHowNearest);
 110                                 mrcPixelDataGet(outW,      X, Y, Z, &data,  mrcPixelRePart, mrcPixelHowNearest);
 111
 112                         //      DEBUGPRINT6("%f %f %f %f %f %f\n", re, im, w, reOut, imOut, data);
 113
 114                                 mrcPixelDataSet(tmpOut,    X, Y, Z, reOut+re*w, mrcPixelRePart);
 115                                 mrcPixelDataSet(tmpOut,    X, Y, Z, imOut+im*w, mrcPixelImPart);
 116                                 mrcPixelDataSet(outW,      X, Y, Z, data +w,    mrcPixelRePart);
 117                         }
 118                         }
 119                         }
 120                 }
 121
 122                 /* Average */
 123                 DEBUGPRINT("Start Average");
 124                 for(Z=-tmpOut->HeaderN.z/2.0; Z< tmpOut->HeaderN.z/2.0; Z++) {
 125                 for(Y=-tmpOut->HeaderN.y/2.0; Y< tmpOut->HeaderN.y/2.0; Y++) {
 126                 for(X=0;                      X<=tmpOut->HeaderN.x/2.0; X++) {
 127                         mrcPixelDataGet(tmpOut,    X, Y, Z, &reOut, mrcPixelRePart, mrcPixelHowNearest);
 128                         mrcPixelDataGet(tmpOut,    X, Y, Z, &imOut, mrcPixelImPart, mrcPixelHowNearest);
 129                         mrcPixelDataGet(outW,      X, Y, Z, &w,     mrcPixelRePart, mrcPixelHowNearest);
 130                         //DEBUGPRINT6("(%f,%f) w %f at (%f %f %f)\n", reOut, imOut, w, X, Y, Z);
 131                         if(0<w) {
 132                                 mrcPixelDataSet(tmpOut,    X, Y, Z, reOut/w, mrcPixelRePart);
 133                                 mrcPixelDataSet(tmpOut,    X, Y, Z, imOut/w, mrcPixelImPart);
 134                         } else {
 135                                 mrcPixelDataSet(tmpOut,    X, Y, Z, 0.0, mrcPixelRePart);
 136                                 mrcPixelDataSet(tmpOut,    X, Y, Z, 0.0, mrcPixelImPart);
 137                         }
 138                 }
 139                 }
 140                 }
 141                 DEBUGPRINT("End Average");
 142         } else {
 143                 DEBUGPRINT("Weight: Img\n");
 144                 for(i=0; i<n; i++) {
 145                         for(z=0; z< tmpIn[i].HeaderN.z; z++) {
 146                         for(y=0; y< tmpIn[i].HeaderN.y; y++) {
 147                         for(x=0; x< tmpIn[i].HeaderN.x; x++) {
 148                                 mrcPixelDataGet(&tmpIn[i], x, y, z, &re,    mrcPixelRePart, mrcPixelHowNearest);
 149                                 mrcPixelDataGet(&inW[i],   x, y, z, &w,     mrcPixelRePart, mrcPixelHowNearest);
 150                                 mrcPixelDataGet(tmpOut,    x, y, z, &reOut, mrcPixelRePart, mrcPixelHowNearest);
 151                                 mrcPixelDataGet(outW,      x, y, z, &data,  mrcPixelRePart, mrcPixelHowNearest);
 152
 153                                 mrcPixelDataSet(tmpOut,    x, y, z, reOut+re*w, mrcPixelRePart);
 154                                 mrcPixelDataSet(outW,      x, y, z, data +w,    mrcPixelRePart);
 155                         }
 156                         }
 157                         }
 158                 }
 159
 160                 for(z=0; z< tmpOut->HeaderN.z; z++) {
 161                 for(y=0; y< tmpOut->HeaderN.y; y++) {
 162                 for(x=0; x< tmpOut->HeaderN.x; x++) {
 163                         mrcPixelDataGet(outW,      x, y, z, &w,     mrcPixelRePart, mrcPixelHowNearest);
 164                         mrcPixelDataGet(tmpOut,    x, y, z, &reOut, mrcPixelRePart, mrcPixelHowNearest);
 165                         if(0<w) {
 166                                 mrcPixelDataSet(tmpOut, x, y, z, reOut/w, mrcPixelRePart);
 167                         } else {
 168                                 mrcPixelDataSet(tmpOut, x, y, z, 0.0,     mrcPixelRePart);
 169                         }
 170                 }
 171                 }
 172                 }
 173         }
 174
 175         /* Post Treatment */
 176         if(IsFT(&in[0], "in AverageWithWeight", 0)) {
 177                 if(IsFT(&inW[0], "in AverageWithWeight", 0)) {
 178                         /*      in:  FT
 179                                 inW: FT */
 180                 } else {
 181                         /*      in:  FT
 182                                 inW: Image */
 183                         for(i=0; i<n; i++) {
 184                                 mrcImageFree(&tmpIn[i], 0);
 185                         }
 186                         memoryFree(tmpIn);
 187                         lmrcImageFFT(out, tmpOut, 0);
 188                         mrcImageFree(tmpOut, 0);
 189                         memoryFree(tmpOut);
 190                 }
 191         } else {
 192                 if(IsFT(&inW[0], "in AverageWithWeight", 0)) {
 193                         /*      in:  Image
 194                                 inW: FT */
 195                         DEBUGPRINT("Start FREE")
 196                         for(i=0; i<n; i++) {
 197                                 mrcImageFree(&tmpIn[i], 0);
 198                         }
 199                         memoryFree(tmpIn);
 200                         DEBUGPRINT("Start FFT")
 201                         DEBUGPRINT1("out->HeaderMode: %d\n", out->HeaderMode);
 202                         DEBUGPRINT1("tmpOut->HeaderMode: %d\n", tmpOut->HeaderMode);
 203                         lmrcImageFFT(out, tmpOut, 0);
 204
 205                         DEBUGPRINT("End FFT")
 206                         mrcImageFree(tmpOut, 0);
 207                         memoryFree(tmpOut);
 208                 } else {
 209                         /*      in:  Image
 210                                 inW: Image */
 211                 }
 212         }
 213 }