Implement FHT interpolations for eSRRF and Error Map

setup.py

@@ -22,7 +22,7 @@
 ]
 EXTRA_LING_ARGS = []
 
-VERSION = "1.3.0"  # sets version number for whole package
+VERSION = "1.3.1"  # sets version number for whole package
 
 
 def run_command(command: str) -> str:

src/include/_c_gradients.c

@@ -24,20 +24,26 @@ void _c_gradient_radiality(float* image, float* imGc, float* imGr, int rows,
   }
 }
 
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
 // as in REF:
 // https://github.com/HenriquesLab/NanoJ-eSRRF/blob/785c71b3bd508c938f63bb780cba47b0f1a5b2a7/resources/liveSRRF.cl
 // under calculateGradient_2point
-void _c_gradient_2point(float* image, float* imGc, float* imGr, int rows,
+void _c_gradient_two_point(float* image, float* imGc, float* imGr, int rows,
                      int cols) {
-  int c0, r0, c1, r1;
-  for (int j = 1; j < rows; j++) {
-    r1 = j * cols;
-    r0 = (j - 1) * cols;
-    for (int i = 1; i < cols; i++) {
-      c1 = i;
-      c0 = i - 1;
-      imGc[r1 + i] = image[r1 + c1] - image[r1 + c0];
-      imGr[r1 + i] = image[r1 + c1] - image[r0 + c1];
+  int c1, r1;
+  int c0, r0;
+  int offset;
+
+  for (r1 = 0; r1 < rows; r1++) {
+      for (c1 = 0; c1 < cols; c1++) {
+        offset = r1 * cols + c1;
+        r0 = MAX(r1 - 1, 0);
+        c0 = MAX(c1 - 1, 0);
+
+        imGc[offset] = image[offset] - image[r1 * cols + c0];
+        imGr[offset] = image[offset] - image[r0 * cols + c1];
+
     }
   }
 }

src/include/_c_gradients.h

@@ -4,7 +4,7 @@
 void _c_gradient_radiality(float* image, float* imGc, float* imGr, int rows,
                            int cols);
 
-void _c_gradient_2point(float* image, float* imGc, float* imGr, int rows,
+void _c_gradient_two_point(float* image, float* imGc, float* imGr, int rows,
                      int cols);
 
 void _c_gradient_roberts_cross(float* image, float* imGc, float* imGr, int rows, int cols);

src/mako_templates/nanopyx.core.transform._le_esrrf.pyx

@@ -1,6 +1,6 @@
 <%!
 schedulers = ['threaded','threaded_guided','threaded_dynamic','threaded_static']
-%># cython: infer_types=True, wraparound=False, nonecheck=False, boundscheck=False, cdivision=True, language_level=3, profile=False, autogen_pxd=False
+%># cython: infer_types=True, wraparound=False, nonecheck=False, boundscheck=False, cdivision=True, language_level=3, profile=False, autogen_pxd=True
 
 import numpy as np
 
@@ -17,6 +17,7 @@ from ...__opencl__ import cl, cl_array, _fastest_device
 from ._le_interpolation_catmull_rom import ShiftAndMagnify
 from ._le_roberts_cross_gradients import GradientRobertsCross
 from ._le_radial_gradient_convergence import RadialGradientConvergence
+from ._interpolation import fht_space_interpolation as dht_interpolation
 
 
 class eSRRF(LiquidEngine):
@@ -40,123 +41,18 @@ class eSRRF(LiquidEngine):
         """
         @gpu
         """
-        if device is None:
-            device = _fastest_device
-
-        # TODO doIntensityWeighting is irrelevant on gpu2
-        cl_ctx = cl.Context([device['device']])
-        dc = device['device']
-        cl_queue = cl.CommandQueue(cl_ctx)
-
-        output_shape = (image.shape[0], int(image.shape[1]*magnification), int(image.shape[2]*magnification))
-
-        # needs input image, first interpolation output, roberts cross output, magnified gradients and output image
-
-        total_memory = (3*image[0, :, :].nbytes) + (2*np.zeros((1, output_shape[1], output_shape[2]), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1]*2, output_shape[2]*2), dtype=np.float32).nbytes)
-        output_image = np.zeros(output_shape, dtype=np.float32)
-
-        max_slices = int((dc.global_mem_size // total_memory)/mem_div)
-        max_slices = self._check_max_slices(image, max_slices)
-
-        mf = cl.mem_flags
-        input_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(image[0:max_slices, :, :].nbytes, dc, max_slices))
-        output_cl = cl.Buffer(cl_ctx, mf.WRITE_ONLY, self._check_max_buffer_size(np.empty((max_slices, output_shape[1], output_shape[2]), dtype=np.float32).nbytes, dc, max_slices))
-        magnified_cl = cl.Buffer(cl_ctx, mf.READ_WRITE, self._check_max_buffer_size(np.empty((max_slices, output_shape[1], output_shape[2]), dtype=np.float32).nbytes, dc, max_slices))
-        col_gradients_cl = cl.Buffer(cl_ctx, mf.READ_WRITE, self._check_max_buffer_size(np.empty((max_slices, image.shape[1], image.shape[2]), dtype=np.float32).nbytes, dc, max_slices))
-        row_gradients_cl = cl.Buffer(cl_ctx, mf.READ_WRITE, self._check_max_buffer_size(np.empty((max_slices, image.shape[1], image.shape[2]), dtype=np.float32).nbytes, dc, max_slices))
-        col_magnified_gradients_cl = cl.Buffer(cl_ctx, mf.READ_WRITE, self._check_max_buffer_size(np.empty((max_slices, image.shape[1]*magnification*2, image.shape[2]*magnification*2), dtype=np.float32).nbytes, dc, max_slices))
-        row_magnified_gradients_cl = cl.Buffer(cl_ctx, mf.READ_WRITE, self._check_max_buffer_size(np.empty((max_slices, image.shape[1]*magnification*2, image.shape[2]*magnification*2), dtype=np.float32).nbytes, dc, max_slices))
-        cl.enqueue_copy(cl_queue, input_cl, image[0:max_slices,:,:]).wait()
-
-        cr_code = self._get_cl_code("_le_interpolation_catmull_rom_.cl", device['DP'])
-        cr_prg = cl.Program(cl_ctx, cr_code).build(options=["-cl-mad-enable -cl-fast-relaxed-math"])
-        cr_knl = cr_prg.shiftAndMagnify
-
-        rc_code = self._get_cl_code("_le_roberts_cross_gradients.cl", device['DP'])
-        rc_prg = cl.Program(cl_ctx, rc_code).build(options=["-cl-mad-enable -cl-fast-relaxed-math"])
-        rc_knl = rc_prg.gradient_roberts_cross
-
-        rgc_code = self._get_cl_code("_le_radial_gradient_convergence.cl", device['DP'])
-        rgc_prg = cl.Program(cl_ctx, rgc_code).build(options=["-cl-mad-enable -cl-fast-relaxed-math"])
-        rgc_knl = rgc_prg.calculate_rgc
-
-        margin = int(radius * magnification)
-        cdef int lowest_row = margin
-        cdef int highest_row = output_shape[1] - margin
-        cdef int lowest_col = margin
-        cdef int highest_col = output_shape[2] - margin
-
-        for i in range(0, image.shape[0], max_slices):
-            if image.shape[0] - i >= max_slices:
-                n_slices = max_slices
-            else:
-                n_slices = image.shape[0] - i
-
-            cr_knl(cl_queue,
-                (n_slices, int(image.shape[1]*magnification), int(image.shape[2]*magnification)),
-                None,
-                input_cl,
-                magnified_cl,
-                np.float32(0),
-                np.float32(0),
-                np.float32(magnification),
-                np.float32(magnification)).wait()
-
-            rc_knl(cl_queue,
-                (n_slices,),
-                None,
-                input_cl,
-                col_gradients_cl,
-                row_gradients_cl,
-                np.int32(image.shape[1]),
-                np.int32(image.shape[2])).wait()
-
-            cr_knl(cl_queue,
-                (n_slices, int(image.shape[1]*magnification*2), int(image.shape[2]*magnification*2)), 
-                None, 
-                col_gradients_cl,
-                col_magnified_gradients_cl,
-                np.float32(0),
-                np.float32(0),
-                np.float32(magnification*2),
-                np.float32(magnification*2)).wait()
-
-            cr_knl(cl_queue,
-                (n_slices, int(image.shape[1]*magnification*2), int(image.shape[2]*magnification*2)), 
-                None, 
-                row_gradients_cl,
-                row_magnified_gradients_cl,
-                np.float32(0),
-                np.float32(0),
-                np.float32(magnification*2),
-                np.float32(magnification*2)).wait()
-
-            rgc_knl(cl_queue,
-                (n_slices, highest_row - lowest_row, highest_col - lowest_col),
-                None,
-                col_magnified_gradients_cl,
-                row_magnified_gradients_cl,
-                magnified_cl,
-                output_cl,
-                np.int32(image.shape[2]*magnification),
-                np.int32(image.shape[1]*magnification),
-                np.int32(magnification),
-                np.float32(2),
-                np.float32(radius),
-                np.float32(2 * (radius / 2.355) + 1),
-                np.float32(2 * (radius / 2.355) * (radius / 2.355)),
-                np.float32(sensitivity),
-                np.int32(doIntensityWeighting)).wait()
-
-            cl.enqueue_copy(cl_queue, output_image[i:i+n_slices,:,:], output_cl).wait()
-
-            if i+n_slices<image.shape[0]:
-                cl.enqueue_copy(cl_queue, input_cl, image[i+n_slices:i+2*n_slices,:,:]).wait()
-
-            cl_queue.finish()
-
-
-        return output_image
+        runtype = "OpenCL"
+        crsm = ShiftAndMagnify(verbose=False)
+        rbc = GradientRobertsCross(verbose=False)
+        rgc = RadialGradientConvergence(verbose=False)
+
+        magnified_image = dht_interpolation(image, magnification)
+        gradient_col, gradient_row = rbc.run(image, run_type=runtype)
+        gradient_col_interp = dht_interpolation(gradient_col, magnification*2)
+        gradient_row_interp = dht_interpolation(gradient_row, magnification*2)
+        radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=runtype)
+
+        return radial_gradients
 
     % for sch in schedulers:
     def _run_${sch}(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True):
@@ -169,11 +65,11 @@ class eSRRF(LiquidEngine):
         crsm = ShiftAndMagnify(verbose=False)
         rbc = GradientRobertsCross(verbose=False)
         rgc = RadialGradientConvergence(verbose=False)
-        
-        magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype)
+
+        magnified_image = dht_interpolation(image, magnification)
         gradient_col, gradient_row = rbc.run(image, run_type=runtype)
-        gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*2, magnification*2, run_type=runtype)
-        gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*2, magnification*2, run_type=runtype)
+        gradient_col_interp = dht_interpolation(gradient_col, magnification*2)
+        gradient_row_interp = dht_interpolation(gradient_row, magnification*2)
         radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=runtype)
 
         return radial_gradients
@@ -185,14 +81,13 @@ class eSRRF(LiquidEngine):
         @cython
         """
         runtype = "Unthreaded"
-        crsm = ShiftAndMagnify(verbose=False)
         rbc = GradientRobertsCross(verbose=False)
         rgc = RadialGradientConvergence(verbose=False)
-        
-        magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype)
+
+        magnified_image = dht_interpolation(image, magnification)
         gradient_col, gradient_row = rbc.run(image, run_type=runtype)
-        gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*2, magnification*2, run_type=runtype)
-        gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*2, magnification*2, run_type=runtype)
+        gradient_col_interp = dht_interpolation(gradient_col, magnification*2)
+        gradient_row_interp = dht_interpolation(gradient_row, magnification*2)
         radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=runtype)
 
         return radial_gradients

src/mako_templates/nanopyx.core.transform._le_radial_gradient_convergence.pyx

@@ -56,8 +56,8 @@ class RadialGradientConvergence(LiquidEngine):
         cdef int margin = int(fwhm * magnification)
 
         cdef int nFrames = gradient_col_interp.shape[0]
-        cdef int rowsM = <int>(gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION)
-        cdef int colsM = <int>(gradient_row_interp.shape[2] / Gx_Gy_MAGNIFICATION)
+        cdef int rowsM = image_interp.shape[1]
+        cdef int colsM = image_interp.shape[2]
 
         cdef float [:,:,:] rgc_map = np.zeros((nFrames, rowsM, colsM), dtype=np.float32)
 
@@ -91,8 +91,8 @@ class RadialGradientConvergence(LiquidEngine):
         cdef int margin = int(fwhm * magnification)
 
         cdef int nFrames = gradient_col_interp.shape[0]
-        cdef int rowsM = <int>(gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION)
-        cdef int colsM = <int>(gradient_row_interp.shape[2] / Gx_Gy_MAGNIFICATION)
+        cdef int rowsM = image_interp.shape[1]
+        cdef int colsM = image_interp.shape[2]
 
         cdef float [:,:,:] rgc_map = np.zeros((nFrames, rowsM, colsM), dtype=np.float32)
 
@@ -141,8 +141,8 @@ class RadialGradientConvergence(LiquidEngine):
 
         # Sizes
         cdef int nFrames = gradient_col_interp.shape[0]
-        cdef int rows_interpolated = <int>(gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION)
-        cdef int cols_interpolated = <int>(gradient_row_interp.shape[2] / Gx_Gy_MAGNIFICATION)
+        cdef int rows_interpolated = image_interp.shape[1]
+        cdef int cols_interpolated = image_interp.shape[2]
 
         # Grid size of the global work space
         lowest_row = margin 

src/mako_templates/nanopyx.core.transform._le_roberts_cross_gradients.pyx

@@ -10,6 +10,8 @@ from ...__liquid_engine__ import LiquidEngine
 from cython.parallel import prange
 from .__interpolation_tools__ import check_image
 
+from scipy.stats import pearsonr as pearson_correlation
+
 cdef extern from "_c_gradients.h":
     void _c_gradient_roberts_cross(float* pixels, float* GxArray, float* GyArray, int w, int h) nogil
 
@@ -28,6 +30,25 @@ class GradientRobertsCross(LiquidEngine):
     def benchmark(self, image):
         image = check_image(image)
         return super().benchmark(image)
+
+    def _compare_runs(self, output_1, output_2):
+        """@public"""
+        if isinstance(output_1, tuple):
+            output_1 = np.concatenate((output_1[0], output_1[1]), axis=0)
+        if isinstance(output_2, tuple):
+            output_2 = np.concatenate((output_2[0], output_2[1]), axis=0)
+        if output_1.ndim > 2:
+            pcc = 0
+            for i in range(output_1.shape[0]):
+                pcc += pearson_correlation(output_1[i, :, :].flatten(), output_2[i, :, :].flatten()).statistic
+            pcc /= output_1.shape[0]
+        else:
+            pcc = pearson_correlation(output_1.flatten(), output_2.flatten()).statistic
+
+        if pcc > 0.8:
+            return True
+        else:
+            return False
 
     def _run_unthreaded(self, float[:,:,:] image):
         """