From f2f2ffd16cbffc9314410502f192824c7a55053a Mon Sep 17 00:00:00 2001 From: Bruno Manuel Santos Saraiva Date: Thu, 11 Sep 2025 11:20:42 +0100 Subject: [PATCH 1/2] reverting esrrf implementation to be the same as NanoJ and added _new files with the new versions (not available in methods) --- .../_c_sr_radial_gradient_convergence.c | 54 +--- .../_c_sr_radial_gradient_convergence.h | 2 +- .../_c_sr_radial_gradient_convergence_new.c | 228 ++++++++++++++ .../nanopyx.core.transform._le_esrrf.pyx | 124 +++----- .../nanopyx.core.transform._le_esrrf_new.pyx | 228 ++++++++++++++ ...nsform._le_radial_gradient_convergence.pyx | 46 ++- src/nanopyx/core/transform/_le_esrrf.pyx | 172 ++++------- src/nanopyx/core/transform/_le_esrrf_new.pyx | 290 ++++++++++++++++++ .../_le_radial_gradient_convergence.cl | 62 ++-- .../_le_radial_gradient_convergence.pyx | 76 +++-- .../_le_radial_gradient_convergence_new.cl | 146 +++++++++ src/nanopyx/methods/esrrf/eSRRF_workflow.py | 2 - 12 files changed, 1093 insertions(+), 337 deletions(-) create mode 100644 src/include/_c_sr_radial_gradient_convergence_new.c create mode 100644 src/mako_templates/nanopyx.core.transform._le_esrrf_new.pyx create mode 100644 src/nanopyx/core/transform/_le_esrrf_new.pyx create mode 100644 src/nanopyx/core/transform/_le_radial_gradient_convergence_new.cl diff --git a/src/include/_c_sr_radial_gradient_convergence.c b/src/include/_c_sr_radial_gradient_convergence.c index e7b9cee9..b9ccffa0 100644 --- a/src/include/_c_sr_radial_gradient_convergence.c +++ b/src/include/_c_sr_radial_gradient_convergence.c @@ -14,64 +14,36 @@ double _c_calculate_dk(float Gx, float Gy, float dx, float dy, float distance) { return Dk; } -void _rotate_vector(float* Gx, float* Gy, float angle) { - float cos_angle = cos(angle); - float sin_angle = sin(angle); - - float original_Gx = *Gx; - float original_Gy = *Gy; - - *Gx = original_Gx * cos_angle - original_Gy * sin_angle; - *Gy = original_Gx * sin_angle + original_Gy * cos_angle; -} - -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, float offset, float xyoffset, float angle) { +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) { float vx, vy, Gx, Gy, dx, dy, distance, distanceWeight, GdotR, Dk; - float correct_vx, correct_vy; - float xc = (float)xM / magnification + offset; // offset in non-magnified space - float yc = (float)yM / magnification + offset; + float xc = (xM + 0.5) / magnification; + float yc = (yM + 0.5) / magnification; float RGC = 0; float distanceWeightSum = 0; - int _start = -(int)(fwhm); - int _end = (int)(fwhm + 1); + int _start = -(int)(Gx_Gy_MAGNIFICATION * fwhm); + int _end = (int)(Gx_Gy_MAGNIFICATION * fwhm + 1); for (int j = _start; j < _end; j++) { - vy = yc + j; + vy = (int)(Gx_Gy_MAGNIFICATION * yc) + j; + vy /= Gx_Gy_MAGNIFICATION; - if (0 < vy && vy <= rowsM/magnification - 1) { + if (0 < vy && vy <= rowsM - 1) { for (int i = _start; i < _end; i++) { - vx = xc + i; + vx = (int)(Gx_Gy_MAGNIFICATION * xc) + i; + vx /= Gx_Gy_MAGNIFICATION; - if (0 < vx && vx <= colsM/magnification - 1) { + if (0 < vx && vx <= colsM - 1) { dx = vx - xc; dy = vy - yc; distance = sqrt(dx * dx + dy * dy); if (distance != 0 && distance <= tSO) { - - correct_vx = vx+xyoffset; - correct_vy = vy+xyoffset; - - if (correct_vx +#include + +double _c_calculate_dw(double distance, double tSS) { + return pow((distance * exp((-distance * distance) / tSS)), 4); +} + +double _c_calculate_dk(float Gx, float Gy, float dx, float dy, float distance) { + float Dk = fabs(Gy * dx - Gx * dy) / sqrt(Gx * Gx + Gy * Gy); + if (isnan(Dk)) { + Dk = distance; + } + Dk = 1 - Dk / distance; + return Dk; +} + +void _rotate_vector(float* Gx, float* Gy, float angle) { + float cos_angle = cos(angle); + float sin_angle = sin(angle); + + float original_Gx = *Gx; + float original_Gy = *Gy; + + *Gx = original_Gx * cos_angle - original_Gy * sin_angle; + *Gy = original_Gx * sin_angle + original_Gy * cos_angle; +} + +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, float offset, float xyoffset, float angle) { + + float vx, vy, Gx, Gy, dx, dy, distance, distanceWeight, GdotR, Dk; + float correct_vx, correct_vy; + + float xc = (float)xM / magnification + offset; // offset in non-magnified space + float yc = (float)yM / magnification + offset; + + float RGC = 0; + float distanceWeightSum = 0; + + int _start = -(int)(fwhm); + int _end = (int)(fwhm + 1); + + for (int j = _start; j < _end; j++) { + vy = yc + j; + + if (0 < vy && vy <= rowsM/magnification - 1) { + for (int i = _start; i < _end; i++) { + vx = xc + i; + + if (0 < vx && vx <= colsM/magnification - 1) { + dx = vx - xc; + dy = vy - yc; + distance = sqrt(dx * dx + dy * dy); + + if (distance != 0 && distance <= tSO) { + + correct_vx = vx+xyoffset; + correct_vy = vy+xyoffset; + + if (correct_vx= 0 && sensitivity > 1) { + RGC = pow(RGC, sensitivity); + } else if (RGC < 0) { + RGC = 0; + } + + 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 = sqrt((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 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, distance, distance_xy, distance_z, distanceWeight, distanceWeight_xy, distanceWeight_z, GdotR, Dk; + + float xc = (float)(xM) / magnification_xy; + float yc = (float)(yM) / magnification_xy; + float zc = (float)(sliceM) / magnification_z; + + float RGC = 0; + float distanceWeightSum = 0; + float distanceWeightSum_xy = 0; + float distanceWeightSum_z = 0; + + int _start = -(int)(fwhm); + int _end = (int)(fwhm + 1); + + int _start_z = -(int)(fwhm_z); + int _end_z = (int)(fwhm_z + 1); + + for (int j = _start; j <= _end; j++) { + vy = yc + j; + + if (0 < vy && vy <= (rowsM/magnification_xy) - 1) { + for (int i = _start; i <= _end; i++) { + vx = xc + i; + + if (0 < vx && vx <= (colsM/magnification_xy) - 1) { + for (int k = _start_z; k <= _end_z; k++) { + vz = zc + k; + + if (0 < vz && vz <= (slicesM/magnification_z) - 1) { + dx = vx - xc; + dy = vy - yc; + distance_z = vz - zc; + distance = sqrt(dx * dx + dy * dy + distance_z * distance_z); + distance_xy = sqrt(dx * dx + dy * dy); + + if (distance != 0 && distance_xy <= tSO && distance_z <= tSO_z) { + int linear_index = (int)(vz * magnification_z) * rowsM * colsM + + (int)(magnification_xy * vy) * colsM + + (int)(magnification_xy * vx); + + Gx = imIntGx[linear_index]; + Gy = imIntGy[linear_index]; + Gz = imIntGz[linear_index]; + + distanceWeight = _c_calculate_dw3D(distance, distance_xy, distance_z, tSS, tSS_z); + + distanceWeightSum += distanceWeight; + GdotR = Gx*dx + Gy*dy + Gz*distance_z; + + if (GdotR < 0) { + Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, distance_z, distance); + RGC += (Dk * distanceWeight); + } + } + } + } + } + } + } + } + + RGC /= distanceWeightSum; + + if (RGC >= 0 && sensitivity > 1) { + RGC = pow(RGC, sensitivity); + } else if (RGC < 0) { + 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; +} diff --git a/src/mako_templates/nanopyx.core.transform._le_esrrf.pyx b/src/mako_templates/nanopyx.core.transform._le_esrrf.pyx index 1fc91a8c..66b9f0c1 100644 --- a/src/mako_templates/nanopyx.core.transform._le_esrrf.pyx +++ b/src/mako_templates/nanopyx.core.transform._le_esrrf.pyx @@ -14,7 +14,6 @@ from .__interpolation_tools__ import check_image, value2array from ...__liquid_engine__ import LiquidEngine from ...__opencl__ import cl, cl_array, _fastest_device -from ._le_convolution import Convolution from ._le_interpolation_catmull_rom import ShiftAndMagnify from ._le_roberts_cross_gradients import GradientRobertsCross from ._le_radial_gradient_convergence import RadialGradientConvergence @@ -29,15 +28,15 @@ class eSRRF(LiquidEngine): self._designation = "eSRRF_ST" super().__init__(clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose) - def run(self, image, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True, run_type=None): + def run(self, image, magnification: int = 5, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True, run_type=None): image = check_image(image) - return self._run(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=run_type) + return self._run(image, magnification=magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=run_type) - def benchmark(self, image, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True): + def benchmark(self, image, magnification: int = 5, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True): image = check_image(image) - return super().benchmark(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) + return super().benchmark(image, magnification=magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) - def _run_opencl(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, mem_div=1): + def _run_opencl(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, mem_div=1): """ @gpu """ @@ -51,52 +50,41 @@ class eSRRF(LiquidEngine): output_shape = (image.shape[0], int(image.shape[1]*magnification), int(image.shape[2]*magnification)) - # needs input image, 2 conv kernels, first interpolation output, roberts cross output, magnified gradients and output image + # needs input image, first interpolation output, roberts cross output, magnified gradients and output image - total_memory = (3*image[0, :, :].nbytes) + (2*np.zeros((2,2), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1], output_shape[2]), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1]*grad_magnification, output_shape[2]*grad_magnification), dtype=np.float32).nbytes) + 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) - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) - - offset = 0 - xyoffset = -0.5 - angle = np.pi/4 - 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_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty_like(kernely).nbytes, dc, max_slices)) - row_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty_like(kernelx).nbytes, dc, max_slices)) - cl.enqueue_copy(cl_queue, col_kernel_cl, kernely).wait() - cl.enqueue_copy(cl_queue, row_kernel_cl, kernelx).wait() 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*grad_magnification, image.shape[2]*magnification*grad_magnification), 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*grad_magnification, image.shape[2]*magnification*grad_magnification), 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 - conv_code = self._get_cl_code("_le_convolution.cl", device['DP']) - conv_prg = cl.Program(cl_ctx, conv_code).build(options=["-cl-mad-enable -cl-fast-relaxed-math"]) - conv_knl = conv_prg.conv2d_2 + 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) - lowest_row = margin # TODO discuss edges calculation - highest_row = output_shape[1] - margin - lowest_col = margin - highest_col = output_shape[2] - margin + 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: @@ -106,7 +94,7 @@ class eSRRF(LiquidEngine): cr_knl(cl_queue, (n_slices, int(image.shape[1]*magnification), int(image.shape[2]*magnification)), - None, + None, input_cl, magnified_cl, np.float32(0), @@ -114,41 +102,34 @@ class eSRRF(LiquidEngine): np.float32(magnification), np.float32(magnification)).wait() - conv_knl(cl_queue, - (n_slices, image.shape[1], image.shape[2]), - None, + rc_knl(cl_queue, + (n_slices,), + None, input_cl, col_gradients_cl, - col_kernel_cl, - np.int32(2)).wait() - - conv_knl(cl_queue, - (n_slices, image.shape[1], image.shape[2]), - None, - input_cl, row_gradients_cl, - row_kernel_cl, - np.int32(2)).wait() + np.int32(image.shape[1]), + np.int32(image.shape[2])).wait() cr_knl(cl_queue, - (n_slices, int(image.shape[1]*magnification*grad_magnification), int(image.shape[2]*magnification*grad_magnification)), + (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*grad_magnification), - np.float32(magnification*grad_magnification)).wait() + np.float32(magnification*2), + np.float32(magnification*2)).wait() cr_knl(cl_queue, - (n_slices, int(image.shape[1]*magnification*grad_magnification), int(image.shape[2]*magnification*grad_magnification)), + (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*grad_magnification), - np.float32(magnification*grad_magnification)).wait() + np.float32(magnification*2), + np.float32(magnification*2)).wait() rgc_knl(cl_queue, (n_slices, highest_row - lowest_row, highest_col - lowest_col), @@ -157,18 +138,15 @@ class eSRRF(LiquidEngine): row_magnified_gradients_cl, magnified_cl, output_cl, - np.int32(output_shape[2]), # Ensure correct dimensions - np.int32(output_shape[1]), + np.int32(image.shape[2]*magnification), + np.int32(image.shape[1]*magnification), np.int32(magnification), - np.float32(grad_magnification), + np.float32(2), np.float32(radius), - np.float32(2 * (radius / 2.355) + 1), # Match sigma calculation - np.float32(2 * (radius / 2.355) ** 2), + np.float32(2 * (radius / 2.355) + 1), + np.float32(2 * (radius / 2.355) * (radius / 2.355)), np.float32(sensitivity), - np.int32(doIntensityWeighting), - np.float32(offset), - np.float32(xyoffset), - np.float32(angle)).wait() + np.int32(doIntensityWeighting)).wait() cl.enqueue_copy(cl_queue, output_image[i:i+n_slices,:,:], output_cl).wait() @@ -181,7 +159,7 @@ class eSRRF(LiquidEngine): return output_image % for sch in schedulers: - def _run_${sch}(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_${sch}(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -189,40 +167,32 @@ class eSRRF(LiquidEngine): """ runtype = "${sch}".capitalize() crsm = ShiftAndMagnify(verbose=False) - conv = Convolution(verbose=False) + rbc = GradientRobertsCross(verbose=False) rgc = RadialGradientConvergence(verbose=False) - - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype) - gradient_col = conv.run(image, kernely, run_type=runtype) - gradient_row = conv.run(image, kernelx, run_type=runtype) - gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, offset=0, xyoffset=-0.5, angle=np.pi/4, run_type=runtype) + 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) + 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 % endfor - def _run_unthreaded(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_unthreaded(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @cython """ runtype = "Unthreaded" crsm = ShiftAndMagnify(verbose=False) - conv = Convolution(verbose=False) + rbc = GradientRobertsCross(verbose=False) rgc = RadialGradientConvergence(verbose=False) - - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype) - gradient_col = conv.run(image, kernely, run_type=runtype) - gradient_row = conv.run(image, kernelx, run_type=runtype) - gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, offset=0, xyoffset=-0.5, angle=np.pi/4, run_type=runtype) + 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) + 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 \ No newline at end of file diff --git a/src/mako_templates/nanopyx.core.transform._le_esrrf_new.pyx b/src/mako_templates/nanopyx.core.transform._le_esrrf_new.pyx new file mode 100644 index 00000000..c0dc776e --- /dev/null +++ b/src/mako_templates/nanopyx.core.transform._le_esrrf_new.pyx @@ -0,0 +1,228 @@ +<%! +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 + +import numpy as np + +cimport numpy as np + +from cython.parallel import parallel, prange + +from libc.math cimport cos, sin + +from .__interpolation_tools__ import check_image, value2array +from ...__liquid_engine__ import LiquidEngine +from ...__opencl__ import cl, cl_array, _fastest_device + +from ._le_convolution import Convolution +from ._le_interpolation_catmull_rom import ShiftAndMagnify +from ._le_roberts_cross_gradients import GradientRobertsCross +from ._le_radial_gradient_convergence import RadialGradientConvergence + + +class eSRRF(LiquidEngine): + """ + eSRRF using the NanoPyx Liquid Engine and running as a single task. + """ + + def __init__(self, clear_benchmarks=False, testing=False, verbose=True): + self._designation = "eSRRF_ST" + super().__init__(clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose) + + def run(self, image, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True, run_type=None): + image = check_image(image) + return self._run(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=run_type) + + def benchmark(self, image, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True): + image = check_image(image) + return super().benchmark(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) + + def _run_opencl(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, mem_div=1): + """ + @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, 2 conv kernels, first interpolation output, roberts cross output, magnified gradients and output image + + total_memory = (3*image[0, :, :].nbytes) + (2*np.zeros((2,2), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1], output_shape[2]), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1]*grad_magnification, output_shape[2]*grad_magnification), 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) + + kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) + kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) + + offset = 0.5 + xyoffset = -0.5 + angle = np.pi/4 + + 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_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty_like(kernely).nbytes, dc, max_slices)) + row_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty_like(kernelx).nbytes, dc, max_slices)) + cl.enqueue_copy(cl_queue, col_kernel_cl, kernely).wait() + cl.enqueue_copy(cl_queue, row_kernel_cl, kernelx).wait() + 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*grad_magnification, image.shape[2]*magnification*grad_magnification), 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*grad_magnification, image.shape[2]*magnification*grad_magnification), 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 + + conv_code = self._get_cl_code("_le_convolution.cl", device['DP']) + conv_prg = cl.Program(cl_ctx, conv_code).build(options=["-cl-mad-enable -cl-fast-relaxed-math"]) + conv_knl = conv_prg.conv2d_2 + + 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) + lowest_row = margin # TODO discuss edges calculation + highest_row = output_shape[1] - margin + lowest_col = margin + 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() + + conv_knl(cl_queue, + (n_slices, image.shape[1], image.shape[2]), + None, + input_cl, + col_gradients_cl, + col_kernel_cl, + np.int32(2)).wait() + + conv_knl(cl_queue, + (n_slices, image.shape[1], image.shape[2]), + None, + input_cl, + row_gradients_cl, + row_kernel_cl, + np.int32(2)).wait() + + cr_knl(cl_queue, + (n_slices, int(image.shape[1]*magnification*grad_magnification), int(image.shape[2]*magnification*grad_magnification)), + None, + col_gradients_cl, + col_magnified_gradients_cl, + np.float32(0), + np.float32(0), + np.float32(magnification*grad_magnification), + np.float32(magnification*grad_magnification)).wait() + + cr_knl(cl_queue, + (n_slices, int(image.shape[1]*magnification*grad_magnification), int(image.shape[2]*magnification*grad_magnification)), + None, + row_gradients_cl, + row_magnified_gradients_cl, + np.float32(0), + np.float32(0), + np.float32(magnification*grad_magnification), + np.float32(magnification*grad_magnification)).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(output_shape[2]), # Ensure correct dimensions + np.int32(output_shape[1]), + np.int32(magnification), + np.float32(grad_magnification), + np.float32(radius), + np.float32(2 * (radius / 2.355) + 1), # Match sigma calculation + np.float32(2 * (radius / 2.355) ** 2), + np.float32(sensitivity), + np.int32(doIntensityWeighting), + np.float32(offset), + np.float32(xyoffset), + np.float32(angle)).wait() + + cl.enqueue_copy(cl_queue, output_image[i:i+n_slices,:,:], output_cl).wait() + + if i+n_slices(gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION) @@ -67,14 +67,14 @@ class RadialGradientConvergence(LiquidEngine): for rM in range(margin, rowsM - margin): for cM in range(margin, colsM - margin): if _doIntensityWeighting: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) * image_interp[f, rM, cM] + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) * image_interp[f, rM, cM] else: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) return np.asarray(rgc_map,dtype=np.float32) % for sch in schedulers: - def _run_${sch}(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, grad_magnification=1, radius=1.5, sensitivity=1, doIntensityWeighting=True, offset: float =0, xyoffset: float =0, angle: float =0): + def _run_${sch}(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -82,13 +82,13 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius - cdef int margin = int(fwhm * magnification) cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = grad_magnification + cdef float Gx_Gy_MAGNIFICATION = 2.0 cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting + cdef int margin = int(fwhm * magnification) cdef int nFrames = gradient_col_interp.shape[0] cdef int rowsM = (gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION) @@ -106,14 +106,14 @@ class RadialGradientConvergence(LiquidEngine): % endif for cM in range(margin, colsM - margin): if _doIntensityWeighting: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) * image_interp[f, rM, cM] + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) * image_interp[f, rM, cM] else: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) return np.asarray(rgc_map,dtype=np.float32) % endfor - def _run_opencl(self, gradient_col_interp, gradient_row_interp, image_interp, magnification=5, grad_magnification=1, radius=1.5, sensitivity=1, doIntensityWeighting=True, offset=0, xyoffset=0, angle=0, device=None, int mem_div=1): + def _run_opencl(self, gradient_col_interp, gradient_row_interp, image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, int mem_div=1): """ @gpu """ @@ -131,16 +131,13 @@ class RadialGradientConvergence(LiquidEngine): # Parameters cdef float sigma = radius / 2.355 cdef float fwhm = radius - cdef int margin = int(fwhm*magnification) cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = grad_magnification + cdef float Gx_Gy_MAGNIFICATION = 2.0 cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting - cdef float _offset = offset - cdef float _xyoffset = xyoffset - cdef float _angle = angle + cdef int margin = int(fwhm * magnification) # Sizes cdef int nFrames = gradient_col_interp.shape[0] @@ -148,7 +145,7 @@ class RadialGradientConvergence(LiquidEngine): cdef int cols_interpolated = (gradient_row_interp.shape[2] / Gx_Gy_MAGNIFICATION) # Grid size of the global work space - lowest_row = margin # TODO discuss edges calculation + lowest_row = margin highest_row = rows_interpolated - margin lowest_col = margin highest_col = cols_interpolated - margin @@ -200,10 +197,7 @@ class RadialGradientConvergence(LiquidEngine): np.float32(tSO), np.float32(tSS), np.float32(_sensitivity), - np.int32(_doIntensityWeighting), - np.float32(_offset), - np.float32(_xyoffset), - np.float32(_angle)).wait() + np.int32(_doIntensityWeighting)).wait() # Copy output cl.enqueue_copy(cl_queue, rgc_map[i:i+n_slices,:,:], rgc_map_out).wait() diff --git a/src/nanopyx/core/transform/_le_esrrf.pyx b/src/nanopyx/core/transform/_le_esrrf.pyx index ec511196..b1e4fbb9 100644 --- a/src/nanopyx/core/transform/_le_esrrf.pyx +++ b/src/nanopyx/core/transform/_le_esrrf.pyx @@ -12,7 +12,6 @@ from .__interpolation_tools__ import check_image, value2array from ...__liquid_engine__ import LiquidEngine from ...__opencl__ import cl, cl_array, _fastest_device -from ._le_convolution import Convolution from ._le_interpolation_catmull_rom import ShiftAndMagnify from ._le_roberts_cross_gradients import GradientRobertsCross from ._le_radial_gradient_convergence import RadialGradientConvergence @@ -27,15 +26,15 @@ class eSRRF(LiquidEngine): self._designation = "eSRRF_ST" super().__init__(clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose) - def run(self, image, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True, run_type=None): + def run(self, image, magnification: int = 5, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True, run_type=None): image = check_image(image) - return self._run(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=run_type) + return self._run(image, magnification=magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=run_type) - def benchmark(self, image, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True): + def benchmark(self, image, magnification: int = 5, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True): image = check_image(image) - return super().benchmark(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) + return super().benchmark(image, magnification=magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) - def _run_opencl(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, mem_div=1): + def _run_opencl(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, mem_div=1): """ @gpu """ @@ -49,52 +48,41 @@ class eSRRF(LiquidEngine): output_shape = (image.shape[0], int(image.shape[1]*magnification), int(image.shape[2]*magnification)) - # needs input image, 2 conv kernels, first interpolation output, roberts cross output, magnified gradients and output image + # needs input image, first interpolation output, roberts cross output, magnified gradients and output image - total_memory = (3*image[0, :, :].nbytes) + (2*np.zeros((2,2), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1], output_shape[2]), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1]*grad_magnification, output_shape[2]*grad_magnification), dtype=np.float32).nbytes) + 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) - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) - - offset = 0 - xyoffset = -0.5 - angle = np.pi/4 - 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_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty_like(kernely).nbytes, dc, max_slices)) - row_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty_like(kernelx).nbytes, dc, max_slices)) - cl.enqueue_copy(cl_queue, col_kernel_cl, kernely).wait() - cl.enqueue_copy(cl_queue, row_kernel_cl, kernelx).wait() 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*grad_magnification, image.shape[2]*magnification*grad_magnification), 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*grad_magnification, image.shape[2]*magnification*grad_magnification), 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 - conv_code = self._get_cl_code("_le_convolution.cl", device['DP']) - conv_prg = cl.Program(cl_ctx, conv_code).build(options=["-cl-mad-enable -cl-fast-relaxed-math"]) - conv_knl = conv_prg.conv2d_2 + 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) - lowest_row = margin # TODO discuss edges calculation - highest_row = output_shape[1] - margin - lowest_col = margin - highest_col = output_shape[2] - margin + 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: @@ -104,7 +92,7 @@ class eSRRF(LiquidEngine): cr_knl(cl_queue, (n_slices, int(image.shape[1]*magnification), int(image.shape[2]*magnification)), - None, + None, input_cl, magnified_cl, np.float32(0), @@ -112,41 +100,34 @@ class eSRRF(LiquidEngine): np.float32(magnification), np.float32(magnification)).wait() - conv_knl(cl_queue, - (n_slices, image.shape[1], image.shape[2]), - None, + rc_knl(cl_queue, + (n_slices,), + None, input_cl, col_gradients_cl, - col_kernel_cl, - np.int32(2)).wait() - - conv_knl(cl_queue, - (n_slices, image.shape[1], image.shape[2]), - None, - input_cl, row_gradients_cl, - row_kernel_cl, - np.int32(2)).wait() + np.int32(image.shape[1]), + np.int32(image.shape[2])).wait() cr_knl(cl_queue, - (n_slices, int(image.shape[1]*magnification*grad_magnification), int(image.shape[2]*magnification*grad_magnification)), + (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*grad_magnification), - np.float32(magnification*grad_magnification)).wait() + np.float32(magnification*2), + np.float32(magnification*2)).wait() cr_knl(cl_queue, - (n_slices, int(image.shape[1]*magnification*grad_magnification), int(image.shape[2]*magnification*grad_magnification)), + (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*grad_magnification), - np.float32(magnification*grad_magnification)).wait() + np.float32(magnification*2), + np.float32(magnification*2)).wait() rgc_knl(cl_queue, (n_slices, highest_row - lowest_row, highest_col - lowest_col), @@ -155,18 +136,15 @@ class eSRRF(LiquidEngine): row_magnified_gradients_cl, magnified_cl, output_cl, - np.int32(output_shape[2]), # Ensure correct dimensions - np.int32(output_shape[1]), + np.int32(image.shape[2]*magnification), + np.int32(image.shape[1]*magnification), np.int32(magnification), - np.float32(grad_magnification), + np.float32(2), np.float32(radius), - np.float32(2 * (radius / 2.355) + 1), # Match sigma calculation - np.float32(2 * (radius / 2.355) ** 2), + np.float32(2 * (radius / 2.355) + 1), + np.float32(2 * (radius / 2.355) * (radius / 2.355)), np.float32(sensitivity), - np.int32(doIntensityWeighting), - np.float32(offset), - np.float32(xyoffset), - np.float32(angle)).wait() + np.int32(doIntensityWeighting)).wait() cl.enqueue_copy(cl_queue, output_image[i:i+n_slices,:,:], output_cl).wait() @@ -178,7 +156,7 @@ class eSRRF(LiquidEngine): return output_image - def _run_threaded(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_threaded(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -186,21 +164,17 @@ class eSRRF(LiquidEngine): """ runtype = "threaded".capitalize() crsm = ShiftAndMagnify(verbose=False) - conv = Convolution(verbose=False) + rbc = GradientRobertsCross(verbose=False) rgc = RadialGradientConvergence(verbose=False) - - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype) - gradient_col = conv.run(image, kernely, run_type=runtype) - gradient_row = conv.run(image, kernelx, run_type=runtype) - gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, offset=0, xyoffset=-0.5, angle=np.pi/4, run_type=runtype) + 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) + 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 - def _run_threaded_guided(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_threaded_guided(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -208,21 +182,17 @@ class eSRRF(LiquidEngine): """ runtype = "threaded_guided".capitalize() crsm = ShiftAndMagnify(verbose=False) - conv = Convolution(verbose=False) + rbc = GradientRobertsCross(verbose=False) rgc = RadialGradientConvergence(verbose=False) - - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype) - gradient_col = conv.run(image, kernely, run_type=runtype) - gradient_row = conv.run(image, kernelx, run_type=runtype) - gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, offset=0, xyoffset=-0.5, angle=np.pi/4, run_type=runtype) + 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) + 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 - def _run_threaded_dynamic(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_threaded_dynamic(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -230,21 +200,17 @@ class eSRRF(LiquidEngine): """ runtype = "threaded_dynamic".capitalize() crsm = ShiftAndMagnify(verbose=False) - conv = Convolution(verbose=False) + rbc = GradientRobertsCross(verbose=False) rgc = RadialGradientConvergence(verbose=False) - - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype) - gradient_col = conv.run(image, kernely, run_type=runtype) - gradient_row = conv.run(image, kernelx, run_type=runtype) - gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, offset=0, xyoffset=-0.5, angle=np.pi/4, run_type=runtype) + 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) + 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 - def _run_threaded_static(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_threaded_static(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -252,39 +218,31 @@ class eSRRF(LiquidEngine): """ runtype = "threaded_static".capitalize() crsm = ShiftAndMagnify(verbose=False) - conv = Convolution(verbose=False) + rbc = GradientRobertsCross(verbose=False) rgc = RadialGradientConvergence(verbose=False) - - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype) - gradient_col = conv.run(image, kernely, run_type=runtype) - gradient_row = conv.run(image, kernelx, run_type=runtype) - gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, offset=0, xyoffset=-0.5, angle=np.pi/4, run_type=runtype) + 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) + 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 - def _run_unthreaded(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_unthreaded(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @cython """ runtype = "Unthreaded" crsm = ShiftAndMagnify(verbose=False) - conv = Convolution(verbose=False) + rbc = GradientRobertsCross(verbose=False) rgc = RadialGradientConvergence(verbose=False) - - kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) - kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) magnified_image = crsm.run(image, 0, 0, magnification, magnification, run_type=runtype) - gradient_col = conv.run(image, kernely, run_type=runtype) - gradient_row = conv.run(image, kernelx, run_type=runtype) - gradient_col_interp = crsm.run(gradient_col, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - gradient_row_interp = crsm.run(gradient_row, 0, 0, magnification*grad_magnification, magnification*grad_magnification, run_type=runtype) - radial_gradients = rgc.run(gradient_col_interp, gradient_row_interp, magnified_image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, offset=0, xyoffset=-0.5, angle=np.pi/4, run_type=runtype) + 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) + 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 \ No newline at end of file diff --git a/src/nanopyx/core/transform/_le_esrrf_new.pyx b/src/nanopyx/core/transform/_le_esrrf_new.pyx new file mode 100644 index 00000000..bd83f2e8 --- /dev/null +++ b/src/nanopyx/core/transform/_le_esrrf_new.pyx @@ -0,0 +1,290 @@ +# cython: infer_types=True, wraparound=False, nonecheck=False, boundscheck=False, cdivision=True, language_level=3, profile=False, autogen_pxd=False + +import numpy as np + +cimport numpy as np + +from cython.parallel import parallel, prange + +from libc.math cimport cos, sin + +from .__interpolation_tools__ import check_image, value2array +from ...__liquid_engine__ import LiquidEngine +from ...__opencl__ import cl, cl_array, _fastest_device + +from ._le_convolution import Convolution +from ._le_interpolation_catmull_rom import ShiftAndMagnify +from ._le_roberts_cross_gradients import GradientRobertsCross +from ._le_radial_gradient_convergence import RadialGradientConvergence + + +class eSRRF(LiquidEngine): + """ + eSRRF using the NanoPyx Liquid Engine and running as a single task. + """ + + def __init__(self, clear_benchmarks=False, testing=False, verbose=True): + self._designation = "eSRRF_ST" + super().__init__(clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose) + + def run(self, image, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True, run_type=None): + image = check_image(image) + return self._run(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=run_type) + + def benchmark(self, image, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True): + image = check_image(image) + return super().benchmark(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) + + def _run_opencl(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, mem_div=1): + """ + @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, 2 conv kernels, first interpolation output, roberts cross output, magnified gradients and output image + + total_memory = (3*image[0, :, :].nbytes) + (2*np.zeros((2,2), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1], output_shape[2]), dtype=np.float32).nbytes) + (2*np.zeros((1, output_shape[1]*grad_magnification, output_shape[2]*grad_magnification), 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) + + kernelx = np.array([[0, -1], [1, 0]]).astype(np.float32) + kernely = np.array([[-1, 0], [0, 1]]).astype(np.float32) + + offset = 0.5 + xyoffset = -0.5 + angle = np.pi/4 + + 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_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty_like(kernely).nbytes, dc, max_slices)) + row_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty_like(kernelx).nbytes, dc, max_slices)) + cl.enqueue_copy(cl_queue, col_kernel_cl, kernely).wait() + cl.enqueue_copy(cl_queue, row_kernel_cl, kernelx).wait() + 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*grad_magnification, image.shape[2]*magnification*grad_magnification), 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*grad_magnification, image.shape[2]*magnification*grad_magnification), 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 + + conv_code = self._get_cl_code("_le_convolution.cl", device['DP']) + conv_prg = cl.Program(cl_ctx, conv_code).build(options=["-cl-mad-enable -cl-fast-relaxed-math"]) + conv_knl = conv_prg.conv2d_2 + + 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) + lowest_row = margin # TODO discuss edges calculation + highest_row = output_shape[1] - margin + lowest_col = margin + 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() + + conv_knl(cl_queue, + (n_slices, image.shape[1], image.shape[2]), + None, + input_cl, + col_gradients_cl, + col_kernel_cl, + np.int32(2)).wait() + + conv_knl(cl_queue, + (n_slices, image.shape[1], image.shape[2]), + None, + input_cl, + row_gradients_cl, + row_kernel_cl, + np.int32(2)).wait() + + cr_knl(cl_queue, + (n_slices, int(image.shape[1]*magnification*grad_magnification), int(image.shape[2]*magnification*grad_magnification)), + None, + col_gradients_cl, + col_magnified_gradients_cl, + np.float32(0), + np.float32(0), + np.float32(magnification*grad_magnification), + np.float32(magnification*grad_magnification)).wait() + + cr_knl(cl_queue, + (n_slices, int(image.shape[1]*magnification*grad_magnification), int(image.shape[2]*magnification*grad_magnification)), + None, + row_gradients_cl, + row_magnified_gradients_cl, + np.float32(0), + np.float32(0), + np.float32(magnification*grad_magnification), + np.float32(magnification*grad_magnification)).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(output_shape[2]), # Ensure correct dimensions + np.int32(output_shape[1]), + np.int32(magnification), + np.float32(grad_magnification), + np.float32(radius), + np.float32(2 * (radius / 2.355) + 1), # Match sigma calculation + np.float32(2 * (radius / 2.355) ** 2), + np.float32(sensitivity), + np.int32(doIntensityWeighting), + np.float32(offset), + np.float32(xyoffset), + np.float32(angle)).wait() + + cl.enqueue_copy(cl_queue, output_image[i:i+n_slices,:,:], output_cl).wait() + + if i+n_slices(gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION) @@ -65,13 +65,13 @@ class RadialGradientConvergence(LiquidEngine): for rM in range(margin, rowsM - margin): for cM in range(margin, colsM - margin): if _doIntensityWeighting: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) * image_interp[f, rM, cM] + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) * image_interp[f, rM, cM] else: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) return np.asarray(rgc_map,dtype=np.float32) - def _run_threaded(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, grad_magnification=1, radius=1.5, sensitivity=1, doIntensityWeighting=True, offset: float =0, xyoffset: float =0, angle: float =0): + def _run_threaded(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -79,13 +79,13 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius - cdef int margin = int(fwhm * magnification) cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = grad_magnification + cdef float Gx_Gy_MAGNIFICATION = 2.0 cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting + cdef int margin = int(fwhm * magnification) cdef int nFrames = gradient_col_interp.shape[0] cdef int rowsM = (gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION) @@ -99,11 +99,11 @@ class RadialGradientConvergence(LiquidEngine): for rM in prange(margin, rowsM - margin): for cM in range(margin, colsM - margin): if _doIntensityWeighting: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) * image_interp[f, rM, cM] + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) * image_interp[f, rM, cM] else: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) return np.asarray(rgc_map,dtype=np.float32) - def _run_threaded_guided(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, grad_magnification=1, radius=1.5, sensitivity=1, doIntensityWeighting=True, offset: float =0, xyoffset: float =0, angle: float =0): + def _run_threaded_guided(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -111,13 +111,13 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius - cdef int margin = int(fwhm * magnification) cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = grad_magnification + cdef float Gx_Gy_MAGNIFICATION = 2.0 cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting + cdef int margin = int(fwhm * magnification) cdef int nFrames = gradient_col_interp.shape[0] cdef int rowsM = (gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION) @@ -131,11 +131,11 @@ class RadialGradientConvergence(LiquidEngine): for rM in prange(margin, rowsM - margin, schedule="guided"): for cM in range(margin, colsM - margin): if _doIntensityWeighting: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) * image_interp[f, rM, cM] + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) * image_interp[f, rM, cM] else: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) return np.asarray(rgc_map,dtype=np.float32) - def _run_threaded_dynamic(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, grad_magnification=1, radius=1.5, sensitivity=1, doIntensityWeighting=True, offset: float =0, xyoffset: float =0, angle: float =0): + def _run_threaded_dynamic(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -143,13 +143,13 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius - cdef int margin = int(fwhm * magnification) cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = grad_magnification + cdef float Gx_Gy_MAGNIFICATION = 2.0 cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting + cdef int margin = int(fwhm * magnification) cdef int nFrames = gradient_col_interp.shape[0] cdef int rowsM = (gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION) @@ -163,11 +163,11 @@ class RadialGradientConvergence(LiquidEngine): for rM in prange(margin, rowsM - margin, schedule="dynamic"): for cM in range(margin, colsM - margin): if _doIntensityWeighting: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) * image_interp[f, rM, cM] + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) * image_interp[f, rM, cM] else: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) return np.asarray(rgc_map,dtype=np.float32) - def _run_threaded_static(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, grad_magnification=1, radius=1.5, sensitivity=1, doIntensityWeighting=True, offset: float =0, xyoffset: float =0, angle: float =0): + def _run_threaded_static(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -175,13 +175,13 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius - cdef int margin = int(fwhm * magnification) cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = grad_magnification + cdef float Gx_Gy_MAGNIFICATION = 2.0 cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting + cdef int margin = int(fwhm * magnification) cdef int nFrames = gradient_col_interp.shape[0] cdef int rowsM = (gradient_row_interp.shape[1] / Gx_Gy_MAGNIFICATION) @@ -195,13 +195,13 @@ class RadialGradientConvergence(LiquidEngine): for rM in prange(margin, rowsM - margin, schedule="static"): for cM in range(margin, colsM - margin): if _doIntensityWeighting: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) * image_interp[f, rM, cM] + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) * image_interp[f, rM, cM] else: - rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity, offset, xyoffset, angle) + rgc_map[f, rM, cM] = _c_calculate_rgc(cM, rM, &gradient_col_interp[f,0,0], &gradient_row_interp[f,0,0], colsM, rowsM, _magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, _sensitivity) return np.asarray(rgc_map,dtype=np.float32) - def _run_opencl(self, gradient_col_interp, gradient_row_interp, image_interp, magnification=5, grad_magnification=1, radius=1.5, sensitivity=1, doIntensityWeighting=True, offset=0, xyoffset=0, angle=0, device=None, int mem_div=1): + def _run_opencl(self, gradient_col_interp, gradient_row_interp, image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, int mem_div=1): """ @gpu """ @@ -219,16 +219,13 @@ class RadialGradientConvergence(LiquidEngine): # Parameters cdef float sigma = radius / 2.355 cdef float fwhm = radius - cdef int margin = int(fwhm*magnification) cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = grad_magnification + cdef float Gx_Gy_MAGNIFICATION = 2.0 cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting - cdef float _offset = offset - cdef float _xyoffset = xyoffset - cdef float _angle = angle + cdef int margin = int(fwhm * magnification) # Sizes cdef int nFrames = gradient_col_interp.shape[0] @@ -236,7 +233,7 @@ class RadialGradientConvergence(LiquidEngine): cdef int cols_interpolated = (gradient_row_interp.shape[2] / Gx_Gy_MAGNIFICATION) # Grid size of the global work space - lowest_row = margin # TODO discuss edges calculation + lowest_row = margin highest_row = rows_interpolated - margin lowest_col = margin highest_col = cols_interpolated - margin @@ -288,10 +285,7 @@ class RadialGradientConvergence(LiquidEngine): np.float32(tSO), np.float32(tSS), np.float32(_sensitivity), - np.int32(_doIntensityWeighting), - np.float32(_offset), - np.float32(_xyoffset), - np.float32(_angle)).wait() + np.int32(_doIntensityWeighting)).wait() # Copy output cl.enqueue_copy(cl_queue, rgc_map[i:i+n_slices,:,:], rgc_map_out).wait() diff --git a/src/nanopyx/core/transform/_le_radial_gradient_convergence_new.cl b/src/nanopyx/core/transform/_le_radial_gradient_convergence_new.cl new file mode 100644 index 00000000..c1caea2e --- /dev/null +++ b/src/nanopyx/core/transform/_le_radial_gradient_convergence_new.cl @@ -0,0 +1,146 @@ +float _c_calculate_rgc(int xM, int yM, __global float* imIntGx, __global float* imIntGy, int colsM, int rowsM, int magnification, float Gx_Gy_MAGNIFICATION, float fwhm, float tSO, float tSS, float sensitivity, float offset, float xyoffset, float angle); +double _c_calculate_dw(double distance, double tSS); +double _c_calculate_dk(float Gx, float Gy, float dx, float dy, float distance); + +float2 _rotation_matrix(float2 point, float angle); + +// RGC takes as input the interpolated intensity gradients in the x and y directions + +// calculate distance weight +double _c_calculate_dw(double distance, double tSS) { + return pow((distance * exp((-distance * distance) / tSS)), 4); +} + +// calculate degree of convergence +double _c_calculate_dk(float Gx, float Gy, float dx, float dy, float distance) { + float Dk = fabs(Gy * dx - Gx * dy) / sqrt(Gx * Gx + Gy * Gy); + if (isnan(Dk)) { + Dk = distance; + } + Dk = 1 - Dk / distance; + return Dk; +} + +float2 _rotate_vector(float Gx, float Gy, float angle) { + float cos_angle = cos(angle); + float sin_angle = sin(angle); + float rotated_Gx = Gx * cos_angle - Gy * sin_angle; + float rotated_Gy = Gx * sin_angle + Gy * cos_angle; + return (float2)(rotated_Gx, rotated_Gy); +} + +float2 _rotation_matrix(float2 point, float angle){ + + //xcos - ysin + //xsin + ycos + + float x = point.x; + float y = point.y; + + return (float2)(x*cos(angle) - y*sin(angle), x*sin(angle) + y*cos(angle)); +} + +// calculate radial gradient convergence for a single subpixel + +float _c_calculate_rgc(int xM, int yM, __global float* imIntGx, __global float* imIntGy, int colsM, int rowsM, int magnification, float Gx_Gy_MAGNIFICATION, float fwhm, float tSO, float tSS, float sensitivity, float offset, float xyoffset, float angle) { + + float vx, vy, Gx, Gy, dx, dy, distance, distanceWeight, GdotR, Dk; + float2 correctedv; + float2 correctedd; + float correct_vx, correct_vy; + + float xc = (float)xM / magnification + offset; // offset in non-magnified space + float yc = (float)yM / magnification + offset; + + float RGC = 0; + float distanceWeightSum = 0; + + int _start = -(int)(fwhm); + int _end = (int)(fwhm + 1); + + for (int j = _start; j < _end; j++) { + vy = yc + j; + + if (0 < vy && vy <= (rowsM/magnification) - 1) { + for (int i = _start; i < _end; i++) { + vx = xc + i; + + if (0 < vx && vx <= (colsM/magnification) - 1) { + dx = vx - xc; + dy = vy - yc; + distance = sqrt(dx * dx + dy * dy); + + if (distance != 0 && distance <= tSO) { + + correct_vx = vx+xyoffset; + correct_vy = vy+xyoffset; + + if (correct_vx= 0 && sensitivity > 1) { + RGC = pow(RGC, sensitivity); + } else if (RGC < 0) { + RGC = 0; + } + + return RGC; +} + + + __kernel void calculate_rgc(__global float* imIntGx, __global float* imIntGy, __global float* imInt, __global float* image_out, int nCols, int nRows, int magnification, float Gx_Gy_MAGNIFICATION, float fwhm, float tSO, float tSS, float sensitivity, int doIntensityWeighting, float offset, float xyoffset, float angle) { + + // Index of the current pixel + int f = get_global_id(0); + int row = get_global_id(1); + int col = get_global_id(2); + + // Output image dimensons + int nPixels_out = nRows * nCols; + + // gradient image dimensions + int nPixels_grad = (int) (nRows*Gx_Gy_MAGNIFICATION * nCols*Gx_Gy_MAGNIFICATION); + + row = row + (int)(fwhm*magnification); + col = col + (int)(fwhm*magnification); + + if (doIntensityWeighting == 1) { + image_out[f * nPixels_out + row * nCols + col] = _c_calculate_rgc(col, row, &imIntGx[f * nPixels_grad], &imIntGy[f * nPixels_grad], nCols, nRows, magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, sensitivity, offset, xyoffset, angle) * imInt[f * nPixels_out + row * nCols + col]; + } + else { + image_out[f * nPixels_out + row * nCols + col] = _c_calculate_rgc(col, row, &imIntGx[f * nPixels_grad], &imIntGy[f * nPixels_grad], nCols, nRows, magnification, Gx_Gy_MAGNIFICATION, fwhm, tSO, tSS, sensitivity, offset, xyoffset, angle); + } + } diff --git a/src/nanopyx/methods/esrrf/eSRRF_workflow.py b/src/nanopyx/methods/esrrf/eSRRF_workflow.py index 51de6986..325e839b 100644 --- a/src/nanopyx/methods/esrrf/eSRRF_workflow.py +++ b/src/nanopyx/methods/esrrf/eSRRF_workflow.py @@ -12,7 +12,6 @@ def eSRRF( image, magnification: int = 5, - grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1, frames_per_timepoint: int = 0, @@ -84,7 +83,6 @@ def eSRRF( ), { "magnification": magnification, - "grad_magnification": grad_magnification, "radius": radius, "sensitivity": sensitivity, "doIntensityWeighting": doIntensityWeighting, From bf53b8221f5a0d86336b42eb01e6a91854fe2984 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ant=C3=B3nio=20Brito?= <50997716+antmsbrito@users.noreply.github.com> Date: Thu, 11 Sep 2025 11:33:39 +0100 Subject: [PATCH 2/2] Version bump to 1.3.0 --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 41ee8237..6c87a225 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ ] EXTRA_LING_ARGS = [] -VERSION = "1.2.3" # sets version number for whole package +VERSION = "1.3.0" # sets version number for whole package def run_command(command: str) -> str: