Merge pull request #5 from 3dem/fix_aberration_at_nyquist

Make convention for Nyquist components consistent.

src/ctf.h

@@ -302,6 +302,11 @@ class CTF
 
 	/// Generate (Fourier-space, i.e. FFTW format) image with all CTF values.
 	/// The dimensions of the result array should have been set correctly already
+	//  FFTW format means the Nyquist component is on the positive side.
+	//   e.g. for N = 6,  kx = [0, 1, 2, 3], ky = [0, 1, 2, 3, -2, -1]
+	//  Unfortunately, codes in jaz/ use different convention.
+	//   e.g. for N - 6, kx = [0, 1, 2, 3], ky = [0, 1, 2, -3, -2, -1]
+	//  TODO: FIXME: Thus, the returned values at Nyquist are WRONG!!!
 	void getFftwImage(MultidimArray < RFLOAT > &result, int orixdim, int oriydim, RFLOAT angpix,
 	                  bool do_abs = false, bool do_only_flip_phases = false, bool do_intact_until_first_peak = false,
 	                  bool do_damping = true, bool do_ctf_padding = false, bool do_intact_after_first_peak = false) const;

src/fftw.h

@@ -78,6 +78,11 @@
  * physical definition. It also defines 'kp', 'ip' and 'jp', which are the logical coordinates
  * It also works for 1D or 2D FFTW transforms
  *
+ * FFTW format means, for N = 6,
+ *  kx = [0, 1, 2, 3], XSIZE = 4 = N/2 + 1
+ *  ky = [0, 1, 2, 3, -2, -1], YSIZE = 6
+ *  kz = [0, 1, 2, 3, -2, -1], ZSIZE = 6
+ *
  * @code
  * FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM(V)
  * {

src/jaz/single_particle/ctf/tilt_helper.cpp

@@ -321,7 +321,7 @@ std::vector<Image<RFLOAT> > TiltHelper::computeOddZernike(
 		for (int x = 0; x < sh; x++)
 		{
 			const double xx0 = x/as;
-			const double yy0 = y < sh? y/as : (y-s)/as;
+			const double yy0 = y < sh-1? y/as : (y-s)/as;
 
 			const double xx = mag(0,0) * xx0 + mag(0,1) * yy0;
 			const double yy = mag(1,0) * xx0 + mag(1,1) * yy0;
@@ -515,7 +515,7 @@ std::vector<Image<RFLOAT> > TiltHelper::computeEvenZernike(
 		for (int x = 0; x < sh; x++)
 		{
 			const double xx0 = x/as;
-			const double yy0 = y < sh? y/as : (y-s)/as;
+			const double yy0 = y < sh-1? y/as : (y-s)/as;
 
 			const double xx = mag(0,0) * xx0 + mag(0,1) * yy0;
 			const double yy = mag(1,0) * xx0 + mag(1,1) * yy0;

src/jaz/single_particle/ctf/tilt_helper.h

@@ -160,6 +160,7 @@ class TiltHelper
                 const std::vector<double>& coeffs,
                 Image<RFLOAT>* fit);
 
+		// Nyquist X is positive, Y is negative (non-FFTW!!)
 		static std::vector<Image<RFLOAT>> computeOddZernike(
 				int s, double angpix, const Matrix2D<RFLOAT>& mag, int n_max);
 
@@ -199,6 +200,7 @@ class TiltHelper
                 const std::vector<double>& coeffs,
                 Image<RFLOAT>* fit);
 
+		// Nyquist X is positive, Y is negative (non-FFTW!!)
 		static std::vector<Image<RFLOAT>> computeEvenZernike(
 				int s,
 				double angpix, const Matrix2D<RFLOAT>& mag,

src/jaz/single_particle/fftw_helper.cpp

@@ -25,25 +25,40 @@ void FftwHelper::decenterUnflip2D(const MultidimArray<RFLOAT> &src, MultidimArra
     const long int w = src.xdim;
     const long int h = src.ydim;
 
-    dest.reshape(h, 2*(w - 1));
+    if (h != 2 * (w - 1))
+    {
+        std::cerr << "w = " << w << " h = " << h << std::endl;
+        REPORT_ERROR("FftwHelper::decenterUnflip2D: illegal input size");
+    }
+    dest.reshape(h, h);
 
-    const long int yc = dest.ydim/2;
+    const long int origin = w - 1;
+    const long int nyquist = origin;
 
     for (long int y = 0; y < dest.ydim; y++)
     for (long int x = 0; x < dest.xdim; x++)
     {
-        long int xs = x - w;
+        // Logical coordinates of destination
+        long int xl = x - origin;
+        long int yl = y - origin;
+        int sign = 1;
 
-        if (xs < 0)
-        {
-            long int ys = (y + yc - 1) % dest.ydim;
-            DIRECT_A2D_ELEM(dest, y, x) = -DIRECT_A2D_ELEM(src, h-ys-1, -xs-1);
-        }
-        else
+        if (xl < 0)
         {
-            long int ys = (y + yc) % dest.ydim;
-            DIRECT_A2D_ELEM(dest, y, x) =  DIRECT_A2D_ELEM(src, ys, xs);
+            sign = -1;
+            xl = -xl;
+            yl = -yl;
         }
+
+	// Cannot trust Nyquist
+        if (xl == nyquist)
+            xl--;
+        if (yl == nyquist)
+            yl--;
+        if (xl == -nyquist)
+            yl++;
+
+        DIRECT_A2D_ELEM(dest, y, x) = sign * DIRECT_A2D_ELEM(src, (yl + h) % h, xl); 
     }
 }
 
@@ -52,24 +67,37 @@ void FftwHelper::decenterDouble2D(const MultidimArray<RFLOAT> &src, MultidimArra
     const long int w = src.xdim;
     const long int h = src.ydim;
 
-    dest.reshape(h, 2*(w - 1));
+    if (h != 2 * (w - 1))
+    {
+        std::cerr << "w = " << w << " h = " << h << std::endl;
+        REPORT_ERROR("FftwHelper::decenterDouble2D: illegal input size");
+    }
+    dest.reshape(h, h);
 
-    const long int yc = dest.ydim/2;
+    const long int origin = w - 1;
+    const long int nyquist = origin;
 
     for (long int y = 0; y < dest.ydim; y++)
     for (long int x = 0; x < dest.xdim; x++)
     {
-        long int xs = x - w;
+        // Logical coordinates of destination
+        long int xl = x - origin;
+        long int yl = y - origin;
 
-        if (xs < 0)
-        {
-            long int ys = (y + yc - 1) % dest.ydim;
-            DIRECT_A2D_ELEM(dest, y, x) = DIRECT_A2D_ELEM(src, h-ys-1, -xs-1);
-        }
-        else
+        if (xl < 0)
         {
-            long int ys = (y + yc) % dest.ydim;
-            DIRECT_A2D_ELEM(dest, y, x) =  DIRECT_A2D_ELEM(src, ys, xs+1);
+            xl = -xl;
+            yl = -yl;
         }
+
+	// Cannot trust Nyquist
+        if (xl == nyquist)
+            xl--;
+        if (yl == nyquist)
+            yl--;
+        if (xl == -nyquist)
+            yl++;
+
+        DIRECT_A2D_ELEM(dest, y, x) = DIRECT_A2D_ELEM(src, (yl + h) % h, xl); 
     }
 }

src/jaz/single_particle/fftw_helper.h

@@ -39,6 +39,21 @@ class FftwHelper
         template <typename T>
         static void recenterFull(const MultidimArray<T> &src, MultidimArray<T> &dest);
 
+	/*
+ 	    decenterUnflip2D and decenterDouble2D uses Jasenko's convention of
+	    taking the X Nyquist at positive and the Y Nyquist negative.
+	    For N = 6, the input is treated as kx = [0, 1, 2, 3] and ky = [0, 1, 2, -3, -2, -1].
+
+	    The output obeys the XMIPP convention, so both X and Y are [-3, -2, -1, 0, 1, 2].
+  
+	    We have two problems:
+	        - Outside the jaz folder, the Y Nyquist is positive. So ky = [0, 1, 2, 3, -2, -1].
+	          See comments in fftw.h and ctf.h.
+	        - For aberration functions, f(kx = 3) is not always the same as f(kx = -3), because it is not periodic.
+
+	    Because the Nyquist component is unreliable due to these issues, values from neighbouring row/column
+            are taken instead.
+  	*/
         static void decenterUnflip2D(const MultidimArray<RFLOAT> &src, MultidimArray<RFLOAT> &dest);
         static void decenterDouble2D(const MultidimArray<RFLOAT> &src, MultidimArray<RFLOAT> &dest);
 };

src/jaz/single_particle/obs_model.cpp

@@ -1113,7 +1113,7 @@ const BufferedImage<RFLOAT>& ObservationModel::getMtfImage(int optGroup, int s)
 			for (int x = 0; x < sh; x++)
 			{
 				const double xx = x/as;  // logical X-coordinate in 1/A
-				const double yy = y < sh? y/as : (y-s)/as; // logical Y-coordinate in 1/A
+				const double yy = y < sh-1? y/as : (y-s)/as; // logical Y-coordinate in 1/A
 
 				RFLOAT res = sqrt(xx*xx + yy*yy); // get resolution in 1/Ang
 				int i_0 = FLOOR(res / res_per_elem);

src/jaz/single_particle/obs_model.h

@@ -106,6 +106,7 @@ class ObservationModel
 				bool do_multiply_instead = false, bool do_correct_average_mtf = true);
 
 		// 2D image with the MTF (cached)
+		// Nyquist X is positive, Y is negative (non-FFTW!!)
 		const BufferedImage<RFLOAT>& getMtfImage(int optGroup, int s);
 
 		// 2D image with the average MTF (cached)
@@ -123,11 +124,12 @@ class ObservationModel
 				long int particle, MultidimArray<Complex>& obsImage,
 				bool do_modulate_instead = false);
 
-
 		// effect of antisymmetric aberration (cached)
+		// Nyquist X is positive, Y is negative (non-FFTW!!)
 		const BufferedImage<Complex>& getPhaseCorrection(int optGroup, int s);
 
 		// effect of symmetric aberration (cached)
+		// Nyquist X is positive, Y is negative (non-FFTW!!)
 		const BufferedImage<RFLOAT>& getGammaOffset(int optGroup, int s);
 
 		Matrix2D<RFLOAT> applyAnisoMag(Matrix2D<RFLOAT> A3D, int opticsGroup);