eSRRF3D merge

notebooks/SRMetrics.ipynb

@@ -667,7 +667,11 @@
    ]
   }
  ],
- "metadata": {},
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  }
+ },
  "nbformat": 4,
  "nbformat_minor": 5
 }

src/include/_c_gradients.c

@@ -1,3 +1,4 @@
+#include <math.h>
 
 void _c_gradient_radiality(float* image, float* imGc, float* imGr, int rows,
                           int cols) {
@@ -66,23 +67,27 @@ void _c_gradient_roberts_cross(float* image, float* imGc, float* imGr, int rows,
 
 // as in https://www.nature.com/articles/s41592-022-01669-y#MOESM1
 // 3D Gradient calculation
-void _c_gradient_3d(float* image, float* imGc, float* imGr, float* imGs, int slice,
+void _c_gradient_3d(float* image, float* imGc, float* imGr, float* imGs, int slices,
                  int rows, int cols) {
   float ip0, ip1, ip2, ip3, ip4, ip5, ip6, ip7;
 
-  int z_i, y_i, x_i;
+  int z_i, y_i, x_i, z_1, y_1, x_1;
 
-  for (z_i = 0; z_i < slice - 1; z_i++) {
-    for (y_i = 0; y_i < rows - 1; y_i++) {
-      for (x_i = 0; x_i < cols - 1; x_i++) {
+  for (z_i = 0; z_i < slices; z_i++) {
+    for (y_i = 0; y_i < rows; y_i++) {
+      for (x_i = 0; x_i < cols; x_i++) {
+
+        z_1 = z_i >= slices - 1 ? slices - 1 : z_i + 1;
+        y_1 = y_i >= rows - 1 ? rows - 1 : y_i + 1;
+        x_1 = x_i >= cols - 1 ? cols - 1 : x_i + 1;
         ip0 = image[z_i * rows * cols + y_i * cols + x_i];
-        ip1 = image[z_i * rows * cols + y_i * cols + x_i + 1];
-        ip2 = image[z_i * rows * cols + (y_i + 1) * cols + x_i];
-        ip3 = image[z_i * rows * cols + (y_i + 1) * cols + x_i + 1];
-        ip4 = image[(z_i + 1) * rows * cols + y_i * cols + x_i];
-        ip5 = image[(z_i + 1) * rows * cols + y_i * cols + x_i + 1];
-        ip6 = image[(z_i + 1) * rows * cols + (y_i + 1) * cols + x_i];
-        ip7 = image[(z_i + 1) * rows * cols + (y_i + 1) * cols + x_i + 1];
+        ip1 = image[z_i * rows * cols + y_i * cols + x_1];
+        ip2 = image[z_i * rows * cols + y_1 * cols + x_i];
+        ip3 = image[z_i * rows * cols + y_1 * cols + x_1];
+        ip4 = image[z_1 * rows * cols + y_i * cols + x_i];
+        ip5 = image[z_1 * rows * cols + y_i * cols + x_1];
+        ip6 = image[z_1 * rows * cols + y_1 * cols + x_i];
+        ip7 = image[z_1 * rows * cols + y_1 * cols + x_1];
         imGc[z_i * rows* cols + y_i * cols + x_i] =
             (ip1 + ip3 + ip5 + ip7 - ip0 - ip2 - ip4 - ip6) / 4;
         imGr[z_i * rows* cols + y_i * cols + x_i] =
@@ -93,3 +98,21 @@ void _c_gradient_3d(float* image, float* imGc, float* imGr, float* imGs, int sli
     }
   }
 }
+
+void _c_gradient_2_point_3d(float* image, float* imGc, float* imGr, float* imGs, int slices, int rows, int cols) {
+  int s_i, y_i, x_i, s0, y0, x0;
+
+  for (s_i=0; s_i<slices; s_i++) {
+    for (y_i=0; y_i<rows; y_i++) {
+      for (x_i=0; x_i<cols; x_i++) {
+        s0 = s_i > 0 ? s_i - 1 : 0;
+        y0 = y_i > 0 ? y_i - 1 : 0;
+        x0 = x_i > 0 ? x_i - 1 : 0;
+        imGc[s_i*rows*cols + y_i*cols + x_i] = image[s_i*rows*cols + y_i*cols + x_i] - image[s_i*rows*cols + y_i*cols + x0];
+        imGr[s_i*rows*cols + y_i*cols + x_i] = image[s_i*rows*cols + y_i*cols + x_i] - image[s_i*rows*cols + y0*cols + x_i];
+        imGs[s_i*rows*cols + y_i*cols + x_i] = image[s_i*rows*cols + y_i*cols + x_i] - image[s0*rows*cols + y_i*cols + x_i];
+
+      }
+    }
+  }
+}

src/include/_c_gradients.h

@@ -9,7 +9,9 @@ void _c_gradient_2point(float* image, float* imGc, float* imGr, int rows,
 
 void _c_gradient_roberts_cross(float* image, float* imGc, float* imGr, int rows, int cols);
 
-void _c_gradient_3d(float* image, float* imGc, float* imGr, float* imGs, int slice,
+void _c_gradient_3d(float* image, float* imGc, float* imGr, float* imGs, int slices,
                  int rows, int cols);
 
+void _c_gradient_2_point_3d(float* image, float* imGc, float* imGr, float* imGs, int slices, int rows, int cols);
+
 #endif

src/include/_c_sr_radial_gradient_convergence.c

@@ -1,4 +1,5 @@
 #include <math.h>
+#include <stdio.h>
 
 double _c_calculate_dw(double distance, double tSS) {
   return pow((distance * exp((-distance * distance) / tSS)), 4);
@@ -68,3 +69,136 @@ float _c_calculate_rgc(int xM, int yM, float* imIntGx, float* imIntGy, int colsM
 
     return RGC;
 }
+
+double _c_calculate_dw3D_isotropic(double distance, double tSS) {
+  return pow((distance * exp(-(distance * distance) / tSS)), 4);
+}
+
+double _c_calculate_dw3D(double distance, double distance_xy, double distance_z,double tSS, double tSS_z) {
+  float D_weight_xy, D_weight_z, D_weight;
+  D_weight_xy = (exp((-distance_xy * distance_xy) / tSS));
+  D_weight_z = (exp((-distance_z * distance_z) / tSS_z));
+  D_weight = distance * (D_weight_xy * D_weight_z);
+  return pow(D_weight, 4);
+}
+
+double _c_calculate_dw_xy(double distance_xy, double tSS) {
+  return pow((distance_xy * exp((-distance_xy * distance_xy) / tSS)), 4);
+}
+
+double _c_calculate_dw_z(double distance_z, double tSS_z) {
+  return pow((distance_z * exp((-distance_z * distance_z) / tSS_z)), 4);
+}
+
+double _c_calculate_dk3D(float Gx, float Gy, float Gz, float dx, float dy, float dz, float distance) {
+  float Dk = sqrtf((Gy * dz - Gz * dy) * (Gy * dz - Gz * dy) + (Gz * dx - Gx * dz) * (Gz * dx - Gx * dz) + (Gx * dy - Gy * dx) * (Gx * dy - Gy * dx)) / sqrt(Gx * Gx + Gy * Gy + Gz * Gz);
+  if (isnan(Dk)) {
+    Dk = distance;
+  }
+  Dk = 1 - Dk / distance;
+  return Dk;
+}
+
+float _c_get_bound_value(float* im, int slices, int rows, int cols, int s, int r, int c){
+    int _s = s > 0 ? s : 0;
+    _s = _s < slices - 1 ? _s : slices - 1;
+    int _r = r > 0 ? r : 0;
+    _r = _r < rows - 1 ? _r : rows - 1;
+    int _c = c > 0 ? c : 0;
+    _c = _c < cols - 1 ? _c : cols - 1;
+
+    return im[_s * rows * cols + _r * cols + _c];
+}
+
+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 ratio_px, float Gx_Gy_MAGNIFICATION, float Gz_MAGNIFICATION, 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 xc = (xM + 0.5) / magnification_xy;
+    float yc = (yM + 0.5) / magnification_xy;
+    float zc = (sliceM + 0.5) / magnification_z;
+
+    float RGC = 0;
+    float distanceWeightSum = 0;
+    float distanceWeightSum_xy = 0;
+    float distanceWeightSum_z = 0;
+
+    int _start = -(int)(Gx_Gy_MAGNIFICATION * fwhm);
+    int _end = (int)(Gx_Gy_MAGNIFICATION * fwhm + 1);
+
+    int _start_z = -(int)(Gz_MAGNIFICATION * fwhm_z);
+    int _end_z = (int)(Gz_MAGNIFICATION * fwhm_z + 1);
+
+    for (int j = _start; j <= _end; j++) {
+        vy = ((float) ((int) (Gx_Gy_MAGNIFICATION*yc)) + j)/(float) Gx_Gy_MAGNIFICATION;
+
+        if (0 < vy && vy < rowsM - 1) {
+            for (int i = _start; i <= _end; i++) {
+                vx = ((float) ((int) (Gx_Gy_MAGNIFICATION*xc)) + i)/(float) Gx_Gy_MAGNIFICATION;
+
+                if (0 < vx && vx < colsM - 1) {
+                    for (int k = _start_z; k <= _end_z; k++) {
+                        vz = ((float) ((int) (zc)) + k);
+
+                        if (0 < vz && vz < slicesM - 1) {
+                            dx = vx - xc;
+                            dy = vy - yc;
+                            dz = vz - zc; 
+                            dz_real = dz * ratio_px; // This has been already divided by magnification_z
+                            distance = sqrt(dx * dx + dy * dy + dz_real * dz_real);
+                            distance_xy = sqrt(dx * dx + dy * dy);
+                            distance_z = dz_real;
+
+                            if (distance != 0 && distance_xy <= tSO && distance_z <= tSO_z) {
+                                Gx = _c_get_bound_value(imIntGx, slicesM*Gz_MAGNIFICATION, rowsM*Gx_Gy_MAGNIFICATION, colsM*Gx_Gy_MAGNIFICATION, vz * magnification_z * Gz_MAGNIFICATION, magnification_xy * Gx_Gy_MAGNIFICATION * vy, magnification_xy * Gx_Gy_MAGNIFICATION * vx);
+                                Gy = _c_get_bound_value(imIntGy, slicesM*Gz_MAGNIFICATION, rowsM*Gx_Gy_MAGNIFICATION, colsM*Gx_Gy_MAGNIFICATION, vz * magnification_z * Gz_MAGNIFICATION, magnification_xy * Gx_Gy_MAGNIFICATION * vy, magnification_xy * Gx_Gy_MAGNIFICATION * vx);
+                                Gz = _c_get_bound_value(imIntGz, slicesM*Gz_MAGNIFICATION, rowsM*Gx_Gy_MAGNIFICATION, colsM*Gx_Gy_MAGNIFICATION, vz * magnification_z * Gz_MAGNIFICATION, magnification_xy * Gx_Gy_MAGNIFICATION * vy, magnification_xy * Gx_Gy_MAGNIFICATION * vx);
+
+                                // distanceWeight = _c_calculate_dw3D_isotropic(distance, tSS);
+                                distanceWeight = _c_calculate_dw3D(distance, distance_xy, distance_z, tSS, tSS_z);
+
+                                // distanceWeight_xy = _c_calculate_dw_xy(distance_xy, tSS); 
+                                // distanceWeight_z = _c_calculate_dw_z(distance_z, tSS_z);
+
+                                distanceWeightSum += distanceWeight;
+                                // distanceWeightSum_xy += distanceWeight_xy;
+                                // distanceWeightSum_z += distanceWeight_z;
+                                GdotR = Gx*dx + Gy*dy + Gz*dz_real;
+
+                                if (GdotR < 0) {
+                                    Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, dz_real, distance);
+                                    RGC += (Dk * distanceWeight);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    RGC /= distanceWeightSum;
+
+    if (RGC >= 0) {
+        RGC = pow(RGC, sensitivity);
+    } else {
+        RGC = 0;
+    }
+
+    return RGC;
+}
+
+float _c_calculate_dk_3d(float dx, float dy, float dz, float Gx, float Gy, float Gz, float Gmag, float distance) {
+    float G1 = dy*Gz - dz*Gy;
+    float G2 = dz*Gx - dx*Gz;
+    float G3 = dx*Gy - dy*Gx;
+    float cross_product = sqrt(G1*G1 + G2*G2 + G3*G3)/Gmag;
+
+    if (isnan(cross_product)) {
+        cross_product = distance;
+    }
+
+    cross_product = 1 - cross_product / distance;
+
+    return cross_product;
+}

src/include/_c_sr_radial_gradient_convergence.h

@@ -9,4 +9,18 @@ double _c_calculate_dk(float Gx, float Gy, float dx, float dy, float distance);
 
 float _c_calculate_rgc(int xM, int yM, float* imIntGx, float* imIntGy, int colsM, int rowsM, int magnification, float Gx_Gy_MAGNIFICATION, float fwhm, float tSO, float tSS, float sensitivity);
 
+double _c_calculate_dw3D_isotropic(double distance, double tSS);
+
+double _c_calculate_dw3D(double distance, double distance_xy, double distance_z, double tSS, double tSS_z);
+
+double _c_calculate_dw_xy(double distance_xy, double tSS);
+
+double _c_calculate_dw_z(double distance_z, double tSS_z);
+
+double _c_calculate_dk3D(float Gx, float Gy, float Gz, float dx, float dy, float dz, float distance);
+
+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 ratio_px, float Gx_Gy_MAGNIFICATION, float Gz_MAGNIFICATION, float fwhm, float fwhm_z, float tSO, float tSO_z, float tSS, float tSS_z, float sensitivity);
+
+float _c_calculate_rgc_3d(int cM, int rM, int sM, float* imIntGx, float* imIntGy, float* imIntGz, int colsM, int rowsM, int slicesM, int magnification_xy, int magnification_z, float Gx_Gy_MAGNIFICATION, float fwhm, float fwhm_z, float tSO, float tSS, float sensitivity);
+
 #endif

src/nanopyx/core/transform/_interpolation.pyx

@@ -5,6 +5,8 @@ cimport numpy as np
 
 import cython
 
+from ._le_interpolation_catmull_rom import ShiftAndMagnify as ShiftMagnify_CR
+
 
 cdef extern from "_c_interpolation_catmull_rom.h":
     float _c_cr_interpolate "_c_interpolate" (float *image, float row, float col, int rows, int cols) 
@@ -18,3 +20,37 @@ def cr_interpolate(img_stack, row, col):
     img_stack = np.array(img_stack, dtype=np.float32)
     return _cr_interpolate(img_stack, row, col)
 
+def interpolate_3d(image, magnification_xy: int = 5, magnification_z: int = 5):
+    interpolator = ShiftMagnify_CR()
+
+    xy_interpolated = interpolator.run(image, 0, 0, magnification_xy, magnification_xy)
+
+    yz_interpolated = np.transpose(interpolator.run(np.transpose(xy_interpolated, axes=[1, 0, 2]).copy(), 0, 0, magnification_z, 1), axes=[1, 0, 2]).copy()
+    xz_interpolated = np.transpose(interpolator.run(np.transpose(xy_interpolated, axes=[2, 1, 0]).copy(), 0, 0, 1, magnification_z), axes=[2,1,0]).copy()
+
+    return np.mean(np.array([yz_interpolated, xz_interpolated]), axis=0)
+
+def linear_interpolation_1D_z(image, magnification):
+    # linear interpolator that takes a 3D image and do linear interpolation in a 1D column in the chosen axis
+
+    image_interpolated = np.zeros((image.shape[0] * magnification, image.shape[1], image.shape[2]), dtype=np.float32)
+
+    z_coords = np.linspace(0, image.shape[0], image.shape[0])
+    new_z_coords = np.linspace(0, image.shape[0], image.shape[0] * magnification)
+
+    for r in range(image.shape[1]):
+        for c in range(image.shape[2]):
+            slc = image[:, r, c]
+            image_interpolated[:, r, c] = np.interp(new_z_coords, z_coords, slc)
+
+
+    return image_interpolated
+
+
+def interpolate_3d_zlinear(image, magnification_xy: int = 5, magnification_z: int = 5):
+    interpolator_xy = ShiftMagnify_CR()
+
+    xy_interpolated = interpolator_xy.run(np.ascontiguousarray(image), 0, 0, magnification_xy, magnification_xy)
+    z_interpolated = linear_interpolation_1D_z(xy_interpolated, magnification_z)
+
+    return z_interpolated