OSDN Git Service

mrcFFTExpression modified for 3dimension
authorhonda <honda.yasuhisa624@mail.kyutech.jp>
Wed, 16 Jun 2021 11:36:48 +0000 (20:36 +0900)
committerhonda <honda.yasuhisa624@mail.kyutech.jp>
Wed, 16 Jun 2021 11:36:48 +0000 (20:36 +0900)
src/Objects/DataManip/mrcImage/src/lmrcFFTExpression.c

index 7453f3c..23e9a55 100755 (executable)
@@ -91,16 +91,18 @@ lmrcFFTExpressionOffset(mrcImage* fft)
 void
 lmrcFFTAmplitude(mrcImage* img, mrcImage* fft)
 {
-    long ix, iy, offset;
+    long ix, iy, iz, offset;
     double data, re, im;
 
        offset = lmrcFFTExpressionOffset(fft);
-    for(ix=0; ix<img->HeaderN.x; ix++) {
+    for(iz=0; iz<img->HeaderN.z; iz++){
         for(iy=0; iy<img->HeaderN.y; iy++) {
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &re, mrcPixelRePart, mrcPixelHowNearest);
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &im, mrcPixelImPart, mrcPixelHowNearest);
-            data = sqrt(re*re + im*im);
-            mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, 0.0, data, mrcPixelRePart);
+            for(ix=0; ix<img->HeaderN.x; ix++) {
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &re, mrcPixelRePart, mrcPixelHowNearest);
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &im, mrcPixelImPart, mrcPixelHowNearest);
+                data = sqrt(re*re + im*im);
+                mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, (mrcImageParaTypeReal)iz, data, mrcPixelRePart);
+            }
         }
     }
 }
@@ -108,22 +110,24 @@ lmrcFFTAmplitude(mrcImage* img, mrcImage* fft)
 void 
 lmrcFFTLogAmplitude(mrcImage* img, mrcImage* fft ,double th)
 {
-    long ix, iy;
+    long ix, iy, iz;
     double data, re, im;
        long offset;
 
        offset = lmrcFFTExpressionOffset(fft);
-    for(ix=0; ix<img->HeaderN.x; ix++) {
+    for(iz=0; iz<img->HeaderN.z; iz++){
         for(iy=0; iy<img->HeaderN.y; iy++) {
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &re, mrcPixelRePart, mrcPixelHowNearest);
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &im, mrcPixelImPart, mrcPixelHowNearest);
-           data=re*re + im*im;
-           if (data < th){
-             data = log10(th)/2.0;
-           } else {
-             data = log10(data)/2.0;
-           }
-            mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, 0.0, data, mrcPixelRePart);
+            for(ix=0; ix<img->HeaderN.x; ix++) {
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &re, mrcPixelRePart, mrcPixelHowNearest);
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &im, mrcPixelImPart, mrcPixelHowNearest);
+                data = sqrt(re*re + im*im);
+                if (data < th){
+                    data = log10(th)/2.0;
+                } else {
+                    data = log10(data)/2.0;
+                }
+                mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, (mrcImageParaTypeReal)iz, data, mrcPixelRePart);
+            }
         }
     }
 }
@@ -131,17 +135,19 @@ lmrcFFTLogAmplitude(mrcImage* img, mrcImage* fft ,double th)
 void 
 lmrcFFTPower(mrcImage* img, mrcImage* fft)
 {
-    long ix, iy;
+    long ix, iy, iz;
     double data, re, im;
        long offset;
 
        offset = lmrcFFTExpressionOffset(fft);
-    for(ix=0; ix<img->HeaderN.x; ix++) {
+    for(iz=0; iz<img->HeaderN.z; iz++){
         for(iy=0; iy<img->HeaderN.y; iy++) {
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &re, mrcPixelRePart, mrcPixelHowNearest);
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &im, mrcPixelImPart, mrcPixelHowNearest);
-            data = re*re + im*im;
-            mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, 0.0, data, mrcPixelRePart);
+            for(ix=0; ix<img->HeaderN.x; ix++) {
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &re, mrcPixelRePart, mrcPixelHowNearest);
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &im, mrcPixelImPart, mrcPixelHowNearest);
+                data = re*re + im*im;
+                mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, (mrcImageParaTypeReal)iz, data, mrcPixelRePart);
+            }
         }
     }
 }
@@ -149,37 +155,41 @@ lmrcFFTPower(mrcImage* img, mrcImage* fft)
 void 
 lmrcFFTLogPower(mrcImage* img, mrcImage* fft ,double th)
 {
-    long ix, iy;
+    long ix, iy, iz;
     double data, re, im;
        long offset;
 
        offset = lmrcFFTExpressionOffset(fft);
-    for(ix=0; ix<img->HeaderN.x; ix++) {
+    for(iz=0; iz<img->HeaderN.z; iz++){
         for(iy=0; iy<img->HeaderN.y; iy++) {
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &re, mrcPixelRePart, mrcPixelHowNearest);
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &im, mrcPixelImPart, mrcPixelHowNearest);
-           data = re*re + im*im;
-           if (data < th){
-             data = log10(th);
-           } else {
-             data = log10(re*re + im*im);
-           }       
-            mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, 0.0, data, mrcPixelRePart);
+            for(ix=0; ix<img->HeaderN.x; ix++) {
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &re, mrcPixelRePart, mrcPixelHowNearest);
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &im, mrcPixelImPart, mrcPixelHowNearest);
+                data = re*re + im*im;
+                if (data < th){
+                    data = log10(th);
+                } else {
+                    data = log10(re*re + im*im);
+                }
+                mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, (mrcImageParaTypeReal)iz, data, mrcPixelRePart);
+            }
         }
     }
 }
 void 
 lmrcFFTPhase(mrcImage* img, mrcImage* fft)
 {
-    long ix, iy;
+    long ix, iy, iz;
     double data, re, im;
        long offset;
 
        offset = lmrcFFTExpressionOffset(fft);
-    for(ix=0; ix<img->HeaderN.x; ix++) {
+    for(iz=0; iz<img->HeaderN.z; iz++){
         for(iy=0; iy<img->HeaderN.y; iy++) {
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &data, mrcPixelPhase, mrcPixelHowNearest);
-               mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, 0.0, data, mrcPixelRePart);
+            for(ix=0; ix<img->HeaderN.x; ix++) {
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &data, mrcPixelRePart, mrcPixelHowNearest);
+                mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, (mrcImageParaTypeReal)iz, data, mrcPixelRePart);
+            }
         }
     }
 }
@@ -187,16 +197,18 @@ lmrcFFTPhase(mrcImage* img, mrcImage* fft)
 void
 lmrcFFTReal(mrcImage* img, mrcImage* fft)
 {
-    long ix, iy;
+    long ix, iy, iz;
     double data, re, im;
        long offset;
 
        offset = lmrcFFTExpressionOffset(fft);
-    for(ix=0; ix<img->HeaderN.x; ix++) {
+    for(iz=0; iz<img->HeaderN.z; iz++){
         for(iy=0; iy<img->HeaderN.y; iy++) {
-            mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, 0.0, &re, mrcPixelRePart, mrcPixelHowNearest);
-            data = re;
-            mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, 0.0, data, mrcPixelRePart);
+            for(ix=0; ix<img->HeaderN.x; ix++) {
+                mrcPixelDataGet(fft, (mrcImageParaTypeReal)ix-img->HeaderN.x/2, (mrcImageParaTypeReal)iy-img->HeaderN.y/2+offset, (mrcImageParaTypeReal)iz-img->HeaderN.z/2, &re, mrcPixelRePart, mrcPixelHowNearest);
+                data = re;
+                mrcPixelDataSet(img, (mrcImageParaTypeReal)ix, (mrcImageParaTypeReal)iy, (mrcImageParaTypeReal)iz, data, mrcPixelRePart);
+            }
         }
     }
 }