Changes to eSRRF 3D

setup.py

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

src/include/_c_sr_radial_gradient_convergence.c

@@ -140,7 +140,7 @@ float _c_get_bound_value(float* im, int slices, int rows, int cols, int s, int r
 
 float _c_calculate_rgc3D(int xM, int yM, int sliceM, float* imIntGx, float* imIntGy, float* imIntGz, int colsM, int rowsM, int slicesM, int magnification_xy, int magnification_z, float voxel_ratio, float fwhm, float fwhm_z, float tSO, float tSO_z, float tSS, float tSS_z, float sensitivity) {
 
-    float vx, vy, vz, Gx, Gy, Gz, dx, dy, dz, dz_real, distance, distance_xy, distance_z, distanceWeight, distanceWeight_xy, distanceWeight_z, GdotR, Dk;
+    float vx, vy, vz, Gx, Gy, Gz, dx, dy, dz, distance, distance_xy, distance_z, distanceWeight, distanceWeight_xy, distanceWeight_z, GdotR, Dk;
 
     float xc = (float)(xM) / magnification_xy;
     float yc = (float)(yM) / magnification_xy;
@@ -171,11 +171,9 @@ float _c_calculate_rgc3D(int xM, int yM, int sliceM, float* imIntGx, float* imIn
                         if (0 < vz && vz <= (slicesM/magnification_z) - 1) {
                             dx = vx - xc;
                             dy = vy - yc;
-                            dz = vz - zc; 
-                            dz_real = dz * voxel_ratio;
-                            distance = sqrt(dx * dx + dy * dy + dz_real * dz_real);
+                            distance_z = vz - zc;
+                            distance = sqrt(dx * dx + dy * dy + distance_z * distance_z);
                             distance_xy = sqrt(dx * dx + dy * dy);
-                            distance_z = dz_real;
 
                             if (distance != 0 && distance_xy <= tSO && distance_z <= tSO_z) {
                                 int linear_index = (int)(vz * magnification_z) * rowsM * colsM +
@@ -189,10 +187,10 @@ float _c_calculate_rgc3D(int xM, int yM, int sliceM, float* imIntGx, float* imIn
                                 distanceWeight = _c_calculate_dw3D(distance, distance_xy, distance_z, tSS, tSS_z);
 
                                 distanceWeightSum += distanceWeight;
-                                GdotR = Gx*dx + Gy*dy + Gz*dz_real;
+                                GdotR = Gx*dx + Gy*dy + Gz*distance_z;
 
                                 if (GdotR < 0) {
-                                    Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, dz_real, distance);
+                                    Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, distance_z, distance);
                                     RGC += (Dk * distanceWeight);
                                 }
                             }

src/liquid_benchmarks/_le_esrrf3d/eSRRF3D.yml

@@ -1,48 +1,53 @@
 opencl:
-  ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'':
-    2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'':
-    1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})'
-  : - 12363612.0
-    - 0.534501292015193
-    - 0.4527143339801114
-    - 0.4526337919814978
+  ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'':
+    5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'':
+    1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})'
+  : - 50399999.99999999
+    - 6.163972874986939
+    - 6.691725915996358
+    - 5.600578582962044
+  ? '([''shape(1000, 30, 29, 30)''], {''magnification_xy'': 2, ''magnification_z'':
+    2, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'':
+    1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})'
+  : - 657719999.9999999
+    - 51.01166029192973
 threaded:
-  ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'':
-    2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'':
-    1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})'
-  : - 12363612.0
-    - 1.9790987920132466
-    - 1.9240770420001354
-    - 1.895047333004186
+  ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'':
+    5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'':
+    1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})'
+  : - 50399999.99999999
+    - 19.27102233399637
+    - 18.66210575005971
+    - 18.99253195791971
 threaded_dynamic:
-  ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'':
-    2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'':
-    1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})'
-  : - 12363612.0
-    - 1.7660593329928815
-    - 1.719053625012748
-    - 1.7240022089972626
+  ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'':
+    5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'':
+    1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})'
+  : - 50399999.99999999
+    - 17.991310499957763
+    - 18.140572750009596
+    - 18.62613470805809
 threaded_guided:
-  ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'':
-    2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'':
-    1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})'
-  : - 12363612.0
-    - 1.781112333002966
-    - 1.8617710830003489
-    - 2.0126889999955893
+  ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'':
+    5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'':
+    1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})'
+  : - 50399999.99999999
+    - 18.974687915993854
+    - 18.85600141703617
+    - 18.480262708035298
 threaded_static:
-  ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'':
-    2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'':
-    1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})'
-  : - 12363612.0
-    - 1.8515684169833548
-    - 1.8730584580043796
-    - 1.8657334169838578
+  ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'':
+    5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'':
+    1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})'
+  : - 50399999.99999999
+    - 19.723863166058436
+    - 19.532742790994234
+    - 19.2554613329703
 unthreaded:
-  ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'':
-    2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'':
-    1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})'
-  : - 12363612.0
-    - 4.863078458001837
-    - 4.9081600419885945
-    - 4.920950291998452
+  ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'':
+    5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'':
+    1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})'
+  : - 50399999.99999999
+    - 35.281799208023585
+    - 35.78980629099533
+    - 35.422467958997004

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

@@ -34,25 +34,28 @@ class eSRRF3D(LiquidEngine):
         self._designation = "eSRRF_3D"
         super().__init__(clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose)
 
-    def run(self, image, magnification_xy: int = 2, magnification_z: int = 2, radius: float = 1.5, radius_z: float = 1.5, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True, run_type=None):
+    def run(self, image, magnification_xy: int = 2, magnification_z: int = 2, radius: float = 1.5, PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True, run_type=None):
         # TODO: complete and check _run inputs, need to complete variables?
+        radius_z = radius * PSF_ratio / voxel_ratio
         if len(image.shape) == 3:
             image = image.reshape((1, image.shape[0], image.shape[1], image.shape[2]))
         if len(image.shape) != 4:
             print("Warning:image must either be 3D or 4D. If 3D, it will be reshaped to 4D.")
             return None
-        if radius * 2 > (image.shape[2]) / 2 or radius * 2 > (image.shape[3] / 2):
-            print("Warning: Radius is too big for the image. Half the radius must be smaller than both half the number of columns and half number of rows of the image.")
+        if radius > (image.shape[2]) / 2 or radius * 2 > (image.shape[3] / 2):
+            print("Warning: Radius is too big for the image. The radius must be smaller than half the number of columns and half number of rows of the image.")
             return None
-        if radius_z * 2 > image.shape[1] / 2:
-            print("Warning: Radius_z is too big for the image. Half the radius_z must be smaller than half of number of Z planes.")
+        if radius_z > image.shape[1] / 2:
+            print("Warning: Radius_z is too big for the image. The radius_z must be smaller than half of number of Z planes. Radius_z is automatically calculated from radius * PSF_ratio / voxel_ratio.")
             return None
         if image.dtype != np.float32:
             image = image.astype(np.float32)
 
         return self._run(image, magnification_xy=magnification_xy, magnification_z=magnification_z, radius=radius, radius_z=radius_z, voxel_ratio=voxel_ratio, sensitivity=sensitivity, mode=mode, doIntensityWeighting=doIntensityWeighting, run_type=run_type)
 
-    def benchmark(self, image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, radius_z: float = 1.5, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True):
+    def benchmark(self, image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True):
+
+        radius_z = radius * PSF_ratio / voxel_ratio
         if image.dtype != np.float32:
             image = image.astype(np.float32)
         if len(image.shape) == 4:
@@ -70,7 +73,6 @@ class eSRRF3D(LiquidEngine):
         % endif
         @cython
         """
-
         time_start = time.time()
         # calculate all constants
         cdef float sigma = radius / 2.355

src/nanopyx/core/transform/_le_esrrf3d.pyx

@@ -32,25 +32,28 @@ class eSRRF3D(LiquidEngine):
         self._designation = "eSRRF_3D"
         super().__init__(clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose)
 
-    def run(self, image, magnification_xy: int = 2, magnification_z: int = 2, radius: float = 1.5, radius_z: float = 1.5, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True, run_type=None):
+    def run(self, image, magnification_xy: int = 2, magnification_z: int = 2, radius: float = 1.5, PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True, run_type=None):
         # TODO: complete and check _run inputs, need to complete variables?
+        radius_z = radius * PSF_ratio / voxel_ratio
         if len(image.shape) == 3:
             image = image.reshape((1, image.shape[0], image.shape[1], image.shape[2]))
         if len(image.shape) != 4:
             print("Warning:image must either be 3D or 4D. If 3D, it will be reshaped to 4D.")
             return None
-        if radius * 2 > (image.shape[2]) / 2 or radius * 2 > (image.shape[3] / 2):
-            print("Warning: Radius is too big for the image. Half the radius must be smaller than both half the number of columns and half number of rows of the image.")
+        if radius > (image.shape[2]) / 2 or radius * 2 > (image.shape[3] / 2):
+            print("Warning: Radius is too big for the image. The radius must be smaller than half the number of columns and half number of rows of the image.")
             return None
-        if radius_z * 2 > image.shape[1] / 2:
-            print("Warning: Radius_z is too big for the image. Half the radius_z must be smaller than half of number of Z planes.")
+        if radius_z > image.shape[1] / 2:
+            print("Warning: Radius_z is too big for the image. The radius_z must be smaller than half of number of Z planes. Radius_z is automatically calculated from radius * PSF_ratio / voxel_ratio.")
             return None
         if image.dtype != np.float32:
             image = image.astype(np.float32)
 
         return self._run(image, magnification_xy=magnification_xy, magnification_z=magnification_z, radius=radius, radius_z=radius_z, voxel_ratio=voxel_ratio, sensitivity=sensitivity, mode=mode, doIntensityWeighting=doIntensityWeighting, run_type=run_type)
 
-    def benchmark(self, image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, radius_z: float = 1.5, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True):
+    def benchmark(self, image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True):
+
+        radius_z = radius * PSF_ratio / voxel_ratio
         if image.dtype != np.float32:
             image = image.astype(np.float32)
         if len(image.shape) == 4:
@@ -65,7 +68,6 @@ class eSRRF3D(LiquidEngine):
         @threaded
         @cython
         """
-
         time_start = time.time()
         # calculate all constants
         cdef float sigma = radius / 2.355
@@ -142,7 +144,6 @@ class eSRRF3D(LiquidEngine):
         @threaded
         @cython
         """
-
         time_start = time.time()
         # calculate all constants
         cdef float sigma = radius / 2.355
@@ -219,7 +220,6 @@ class eSRRF3D(LiquidEngine):
         @threaded
         @cython
         """
-
         time_start = time.time()
         # calculate all constants
         cdef float sigma = radius / 2.355
@@ -296,7 +296,6 @@ class eSRRF3D(LiquidEngine):
         @threaded
         @cython
         """
-
         time_start = time.time()
         # calculate all constants
         cdef float sigma = radius / 2.355
@@ -372,7 +371,6 @@ class eSRRF3D(LiquidEngine):
         @cpu
         @cython
         """
-
         time_start = time.time()
         # calculate all constants
         cdef float sigma = radius / 2.355

src/nanopyx/core/transform/_le_esrrf3d_.cl

@@ -239,11 +239,9 @@ float _c_calculate_rgc3D(int xM, int yM, int sliceM, __global float* imIntGx, __
                         if (0 < vz && vz <= (slicesM/magnification_z) - 1) {
                             dx = vx - xc;
                             dy = vy - yc;
-                            dz = vz - zc; 
-                            dz_real = dz * voxel_ratio;
-                            distance = sqrt(dx * dx + dy * dy + dz_real * dz_real);
+                            distance_z = vz - zc; 
+                            distance = sqrt(dx * dx + dy * dy + distance_z * distance_z);
                             distance_xy = sqrt(dx * dx + dy * dy);
-                            distance_z = dz_real;
 
                             if (distance != 0 && distance_xy <= tSO && distance_z <= tSO_z) {
                                 int linear_index = (int)(vz * magnification_z) * rowsM * colsM +
@@ -264,10 +262,10 @@ float _c_calculate_rgc3D(int xM, int yM, int sliceM, __global float* imIntGx, __
                                 distanceWeightSum = distanceWeightSum + distanceWeight;
                                 // distanceWeightSum_xy += distanceWeight_xy;
                                 // distanceWeightSum_z += distanceWeight_z;
-                                GdotR = Gx*dx + Gy*dy + Gz*dz_real;
+                                GdotR = Gx*dx + Gy*dy + Gz*distance_z;
 
                                 if (GdotR < 0) {
-                                    Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, dz_real, distance);
+                                    Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, distance_z, distance);
                                     RGC = RGC + (Dk * distanceWeight);
                                 }
                             }