2 # lmrcImageAverage : $Revision$
5 # Usage : lmrcImageAverage
16 #include "./lmrcImageAverage.h"
19 lmrcImageAverage(mrcImage* out, mrcImage* in, int n, int mode)
25 lmrcImageAdd(out, &(in[i]), &num);
28 lmrcImageDevidedByReal(out, (double)num);
30 fprintf(stderr, "Something wrong: n %d != num %d\n", n, num);
35 lmrcImageAverageWithWeight(mrcImage* out, mrcImage* outW, mrcImage* in, mrcImage* inW, int n, int mode)
40 mrcImageParaTypeReal x, y, z, X, Y, Z;
41 double re, im, reOut, imOut, w, data;
43 DEBUGPRINT("Start lmrcImageAverageWithWeight\n");
45 outW->Header = inW[0].Header;
46 mrcInit(outW, NULL); /* All pixel data are zero */
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)) {
53 DEBUGPRINT("FTvsFT\n");
55 out->Header = in[0].Header;
61 DEBUGPRINT("FTvsImg\n");
62 tmpIn = (mrcImage*)memoryAllocate(sizeof(mrcImage)*n, "in lmrcImageAverageWithWeight");
64 lmrcImageFFT(&tmpIn[i], &in[i], 0);
66 tmpOut = (mrcImage*)memoryAllocate(sizeof(mrcImage), "in lmrcImageAverageWithWeight");
67 tmpOut->Header = in[0].Header;
68 mrcInit(tmpOut, NULL);
69 out->Header = in[0].Header;
72 if(IsFT(&inW[0], "in AverageWithWeight", 0)) {
75 DEBUGPRINT("ImgvsFT\n");
76 tmpIn = (mrcImage*)memoryAllocate(sizeof(mrcImage)*n, "in lmrcImageAverageWithWeight");
78 lmrcImageFFT(&tmpIn[i], &in[i], 0);
80 tmpOut = (mrcImage*)memoryAllocate(sizeof(mrcImage)*1, "in lmrcImageAverageWithWeight");
81 tmpOut->Header = inW[0].Header;
82 tmpOut->HeaderMode = mrcComplexFloatFT;
83 mrcInit(tmpOut, NULL);
87 DEBUGPRINT("Imgvs\n");
89 out->Header = in[0].Header;
96 if(IsFT(&tmpIn[0], "in AverageWithWeight", 0)) {
97 DEBUGPRINT("Weight: FT\n");
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);
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);
112 // DEBUGPRINT6("%f %f %f %f %f %f\n", re, im, w, reOut, imOut, data);
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);
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);
132 mrcPixelDataSet(tmpOut, X, Y, Z, reOut/w, mrcPixelRePart);
133 mrcPixelDataSet(tmpOut, X, Y, Z, imOut/w, mrcPixelImPart);
135 mrcPixelDataSet(tmpOut, X, Y, Z, 0.0, mrcPixelRePart);
136 mrcPixelDataSet(tmpOut, X, Y, Z, 0.0, mrcPixelImPart);
141 DEBUGPRINT("End Average");
143 DEBUGPRINT("Weight: Img\n");
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);
153 mrcPixelDataSet(tmpOut, x, y, z, reOut+re*w, mrcPixelRePart);
154 mrcPixelDataSet(outW, x, y, z, data +w, mrcPixelRePart);
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);
166 mrcPixelDataSet(tmpOut, x, y, z, reOut/w, mrcPixelRePart);
168 mrcPixelDataSet(tmpOut, x, y, z, 0.0, mrcPixelRePart);
176 if(IsFT(&in[0], "in AverageWithWeight", 0)) {
177 if(IsFT(&inW[0], "in AverageWithWeight", 0)) {
184 mrcImageFree(&tmpIn[i], 0);
187 lmrcImageFFT(out, tmpOut, 0);
188 mrcImageFree(tmpOut, 0);
192 if(IsFT(&inW[0], "in AverageWithWeight", 0)) {
195 DEBUGPRINT("Start FREE")
197 mrcImageFree(&tmpIn[i], 0);
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);
205 DEBUGPRINT("End FFT")
206 mrcImageFree(tmpOut, 0);