diff --git a/src/include/_c_sr_radial_gradient_convergence.c b/src/include/_c_sr_radial_gradient_convergence.c index d66a04ff..5f6850c3 100644 --- a/src/include/_c_sr_radial_gradient_convergence.c +++ b/src/include/_c_sr_radial_gradient_convergence.c @@ -14,36 +14,47 @@ double _c_calculate_dk(float Gx, float Gy, float dx, float dy, float distance) { return Dk; } -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) { +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 xc = (xM + 0.5) / magnification; - float yc = (yM + 0.5) / magnification; + 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)(Gx_Gy_MAGNIFICATION * fwhm); - int _end = (int)(Gx_Gy_MAGNIFICATION * fwhm + 1); + int _start = -(int)(2 * fwhm); + int _end = (int)(2 * fwhm + 1); for (int j = _start; j < _end; j++) { - vy = (int)(Gx_Gy_MAGNIFICATION * yc) + j; - vy /= Gx_Gy_MAGNIFICATION; + vy = yc + j; - if (0 < vy && vy <= rowsM - 1) { + if (0 < vy && vy <= rowsM/magnification - 1) { for (int i = _start; i < _end; i++) { - vx = (int)(Gx_Gy_MAGNIFICATION * xc) + i; - vx /= Gx_Gy_MAGNIFICATION; + vx = xc + i; - if (0 < vx && vx <= colsM - 1) { + if (0 < vx && vx <= colsM/magnification - 1) { dx = vx - xc; dy = vy - yc; distance = sqrt(dx * dx + dy * dy); if (distance != 0 && distance <= tSO) { - Gx = imIntGx[(int)(vy * magnification * Gx_Gy_MAGNIFICATION * colsM * Gx_Gy_MAGNIFICATION) + (int)(vx * magnification * Gx_Gy_MAGNIFICATION)]; - Gy = imIntGy[(int)(vy * magnification * Gx_Gy_MAGNIFICATION * colsM * Gx_Gy_MAGNIFICATION) + (int)(vx * magnification * Gx_Gy_MAGNIFICATION)]; + Gx = imIntGx[(int)((vy+xyoffset) * magnification * Gx_Gy_MAGNIFICATION * colsM * Gx_Gy_MAGNIFICATION) + (int)((vx+xyoffset) * magnification * Gx_Gy_MAGNIFICATION)]; + Gy = imIntGy[(int)((vy+xyoffset) * magnification * Gx_Gy_MAGNIFICATION * colsM * Gx_Gy_MAGNIFICATION) + (int)((vx+xyoffset) * magnification * Gx_Gy_MAGNIFICATION)]; + + _rotate_vector(&Gx, &Gy, angle); distanceWeight = _c_calculate_dw(distance, tSS); distanceWeightSum += distanceWeight; diff --git a/src/include/_c_sr_radial_gradient_convergence.h b/src/include/_c_sr_radial_gradient_convergence.h index b9793678..60d48153 100644 --- a/src/include/_c_sr_radial_gradient_convergence.h +++ b/src/include/_c_sr_radial_gradient_convergence.h @@ -7,7 +7,9 @@ double _c_calculate_dw(double distance, double tSS); 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); +void _rotate_vector(float* Gx, float* Gy, 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 offset, float xyoffset, float angle); double _c_calculate_dw3D_isotropic(double distance, double tSS); @@ -21,6 +23,7 @@ double _c_calculate_dk3D(float Gx, float Gy, float Gz, float dx, float dy, float 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 diff --git a/src/liquid_benchmarks/_le_convolution/Convolution.yml b/src/liquid_benchmarks/_le_convolution/Convolution.yml index 1c3a22a0..83a19941 100644 --- a/src/liquid_benchmarks/_le_convolution/Convolution.yml +++ b/src/liquid_benchmarks/_le_convolution/Convolution.yml @@ -1,5 +1,5 @@ Cuda: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.0004861999999974387 - 0.0004765000000048758 @@ -22,7 +22,7 @@ Cuda: - 0.00047459999998977764 - 0.0004846000000213735 Dask: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.00937319999999886 - 0.00935059999999055 @@ -45,7 +45,7 @@ Dask: - 0.009444899999976997 - 0.008983199999988756 Numba: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.005503999999987741 - 0.005404399999989096 @@ -68,7 +68,7 @@ Numba: - 0.006126800000004096 - 0.00559069999999906 OpenCL: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.09695740000000796 - 0.0791657999999984 @@ -91,7 +91,7 @@ OpenCL: - 0.07781260000001566 - 0.07502910000002316 Python: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 2.4462038000000064 - 2.4630003999999985 @@ -114,7 +114,7 @@ Python: - 2.3908779999999865 - 2.423667199999983 Threaded: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.000824199999996722 - 0.0008181999999976597 @@ -137,7 +137,7 @@ Threaded: - 0.0007410999999990509 - 0.00074050000000625 Threaded_dynamic: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.0004450999999932037 - 0.00044110000000330274 @@ -160,7 +160,7 @@ Threaded_dynamic: - 0.00042849999999816646 - 0.00042849999999816646 Threaded_guided: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.00045450000000357704 - 0.0004828000000003385 @@ -183,7 +183,7 @@ Threaded_guided: - 0.0004787000000021635 - 0.000453700000008439 Threaded_static: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.0004264999999890051 - 0.0006580000000013797 @@ -206,7 +206,7 @@ Threaded_static: - 0.00047909999997841624 - 0.0004892999999981384 Transonic: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.005411799999990308 - 0.00537510000000907 @@ -229,7 +229,7 @@ Transonic: - 0.005397600000009106 - 0.005529200000012224 Unthreaded: - (['shape(100, 100)', 'shape(23, 23)'], {}): + (['shape(1, 100, 100)', 'shape(23, 23)'], {}): - 5290000 - 0.006530499999996664 - 0.0064456000000063796 diff --git a/src/mako_templates/nanopyx.core.analysis._le_channel_registration.pyx b/src/mako_templates/nanopyx.core.analysis._le_channel_registration.pyx index ec8db173..2146207a 100644 --- a/src/mako_templates/nanopyx.core.analysis._le_channel_registration.pyx +++ b/src/mako_templates/nanopyx.core.analysis._le_channel_registration.pyx @@ -85,10 +85,9 @@ def _gaussian_filter(inpimg, _runtype, sigma): knl1 = np.exp(-0.5 / (sigma*sigma) * x1 ** 2) knl1 = knl1 / knl1.sum() knl1 = knl1.astype(np.float32) - img1 = conv.run(inpimg,knl1.reshape((len(x1), 1)),run_type=_runtype) + img1 = np.expand_dims(img1, axis=0) img2 = conv.run(img1,knl1.reshape((1, len(x1))),run_type=_runtype) - return img2 class ChannelRegistrationEstimator(LiquidEngine): diff --git a/src/mako_templates/nanopyx.core.transform._le_convolution.pyx b/src/mako_templates/nanopyx.core.transform._le_convolution.pyx index 3974ce96..dc56d34c 100644 --- a/src/mako_templates/nanopyx.core.transform._le_convolution.pyx +++ b/src/mako_templates/nanopyx.core.transform._le_convolution.pyx @@ -13,7 +13,7 @@ from cython.parallel import parallel, prange from libc.math cimport cos, sin from .__interpolation_tools__ import check_image, value2array -from .convolution import check_array, convolution2D_cuda, convolution2D_dask, convolution2D_numba, convolution2D_python, convolution2D_transonic +from .convolution import convolution2D_cuda, convolution2D_dask, convolution2D_numba, convolution2D_python, convolution2D_transonic from ...__liquid_engine__ import LiquidEngine from ...__opencl__ import cl, cl_array, _fastest_device @@ -29,14 +29,15 @@ class Convolution(LiquidEngine): clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose) def run(self, image, kernel, run_type=None): - image = check_array(image) + image = check_image(image) return self._run(image, kernel, run_type=run_type) def benchmark(self, image, kernel): + image = check_image(image) return super().benchmark(image, kernel) % for sch in schedulers: - def _run_${sch}(self, float[:,:] image, float[:,:] kernel): + def _run_${sch}(self, float[:,:,:] image, float[:,:] kernel): """ @cpu % if sch!='unthreaded': @@ -44,8 +45,9 @@ class Convolution(LiquidEngine): % endif @cython """ - cdef int nRows = image.shape[0] - cdef int nCols = image.shape[1] + cdef int nFrames = image.shape[0] + cdef int nRows = image.shape[1] + cdef int nCols = image.shape[2] cdef int nRows_kernel = kernel.shape[0] cdef int nCols_kernel = kernel.shape[1] @@ -61,29 +63,30 @@ class Convolution(LiquidEngine): cdef float acc = 0 - conv_out = np.zeros((nRows, nCols), dtype=np.float32) - cdef float[:,:] _conv_out = conv_out + conv_out = np.zeros((nFrames, nRows, nCols), dtype=np.float32) + cdef float[:,:,:] _conv_out = conv_out with nogil: - % if sch=='unthreaded': - for r in range(nRows): - for c in range(nCols): - % elif sch=='threaded': - for r in prange(nRows): - for c in prange(nCols): - % else: - for r in prange(nRows,schedule="${sch.split('_')[1]}"): - for c in prange(nCols,schedule="${sch.split('_')[1]}"): - % endif - acc = 0 - for kr in range(nRows_kernel): - for kc in range(nCols_kernel): - local_row = min(max(r+(kr-center_r),0),nRows-1) - local_col = min(max(c+(kc-center_c),0),nCols-1) - acc = acc + kernel[kr,kc] * image[local_row, local_col] - _conv_out[r,c] = acc - - return conv_out + for f in range(nFrames): + % if sch=='unthreaded': + for r in range(nRows): + for c in range(nCols): + % elif sch=='threaded': + for r in prange(nRows): + for c in prange(nCols): + % else: + for r in prange(nRows,schedule="${sch.split('_')[1]}"): + for c in prange(nCols,schedule="${sch.split('_')[1]}"): + % endif + acc = 0 + for kr in range(nRows_kernel): + for kc in range(nCols_kernel): + local_row = min(max(r+(kr-center_r),0),nRows-1) + local_col = min(max(c+(kc-center_c),0),nCols-1) + acc = acc + kernel[kr,kc] * image[f,local_row, local_col] + _conv_out[f,r,c] = acc + + return np.squeeze(conv_out) % endfor @@ -99,55 +102,62 @@ class Convolution(LiquidEngine): dc = device['device'] cl_queue = cl.CommandQueue(cl_ctx) - image_out = np.zeros((image.shape[0], image.shape[1]), dtype=np.float32) + image_out = np.zeros((image.shape[0], image.shape[1], image.shape[2]), dtype=np.float32) mf = cl.mem_flags - input_image = cl.image_from_array(cl_ctx, image, mode='r') - input_kernel = cl.image_from_array(cl_ctx, kernel, mode='r') - output_opencl = cl.image_from_array(cl_ctx, image_out, mode='w') + input_image = cl.Buffer(cl_ctx, mf.READ_ONLY, image.nbytes) + cl.enqueue_copy(cl_queue, input_image, image).wait() + + input_kernel = cl.Buffer(cl_ctx, mf.READ_ONLY, kernel.nbytes) + cl.enqueue_copy(cl_queue, input_kernel, kernel).wait() + + output_opencl = cl.Buffer(cl_ctx, mf.WRITE_ONLY, image_out.nbytes) + + kernelsize = kernel.shape[0] cl_queue.finish() code = self._get_cl_code("_le_convolution.cl", device['DP']) prg = cl.Program(cl_ctx, code).build() - knl = prg.conv2d + knl = prg.conv2d_2 knl(cl_queue, - (1,image.shape[0], image.shape[1]), - self.get_work_group(device['device'],(1,image.shape[0], image.shape[1])), + (image.shape[0],image.shape[1],image.shape[2]), + None,#self.get_work_group(device['device'],(image.shape[0], image.shape[1], image.shape[2])), input_image, output_opencl, - input_kernel).wait() + input_kernel, + np.int32(kernelsize)).wait() - cl.enqueue_copy(cl_queue, image_out, output_opencl,origin=(0,0), region=(image.shape[0], image.shape[1])).wait() + cl.enqueue_copy(cl_queue, image_out, output_opencl).wait() cl_queue.finish() - return image_out + return np.squeeze(image_out) def _run_python(self, image, kernel): """ @cpu """ - return convolution2D_python(image, kernel).astype(np.float32) + return np.squeeze(convolution2D_python(image, kernel)) def _run_transonic(self, image, kernel): """ @cpu @threaded """ - return convolution2D_transonic(image, kernel).astype(np.float32) + return np.squeeze(convolution2D_transonic(image, kernel)) def _run_dask(self, image, kernel): """ @cpu @threaded """ - return convolution2D_dask(image, kernel).astype(np.float32) + return np.squeeze(convolution2D_dask(image, kernel)) def _run_cuda(self, image, kernel): """ @gpu """ - return convolution2D_cuda(image, kernel).astype(np.float32) + return np.squeeze(convolution2D_cuda(image, kernel)) def _run_numba(self, image, kernel): """ @@ -155,4 +165,4 @@ class Convolution(LiquidEngine): @threaded @numba """ - return convolution2D_numba(image, kernel).astype(np.float32) + return np.squeeze(convolution2D_numba(image, kernel)) diff --git a/src/mako_templates/nanopyx.core.transform._le_esrrf.pyx b/src/mako_templates/nanopyx.core.transform._le_esrrf.pyx index 09710307..d1c5b98e 100644 --- a/src/mako_templates/nanopyx.core.transform._le_esrrf.pyx +++ b/src/mako_templates/nanopyx.core.transform._le_esrrf.pyx @@ -14,6 +14,7 @@ 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 @@ -28,15 +29,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, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True, run_type=None): + 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, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting, run_type=run_type) + 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, radius: float = 1.5, sensitivity: float = 1, doIntensityWeighting: bool = True): + 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, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) + return super().benchmark(image, magnification=magnification, grad_magnification=grad_magnification, radius=radius, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) - def _run_opencl(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, mem_div=1): + def _run_opencl(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True, device=None, mem_div=1): """ @gpu """ @@ -50,36 +51,53 @@ class eSRRF(LiquidEngine): output_shape = (image.shape[0], int(image.shape[1]*magnification), int(image.shape[2]*magnification)) - # needs input image, first interpolation output, roberts cross output, magnified gradients and output image + # 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((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) + 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 + 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((2,2)).nbytes, dc, max_slices)) + row_kernel_cl = cl.Buffer(cl_ctx, mf.READ_ONLY, self._check_max_buffer_size(np.empty((2,2)).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*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)) + 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 - 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 + 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*2*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 @@ -88,7 +106,7 @@ class eSRRF(LiquidEngine): cr_knl(cl_queue, (n_slices, int(image.shape[1]*magnification), int(image.shape[2]*magnification)), - self.get_work_group(dc, (n_slices, image.shape[1]*magnification, image.shape[2]*magnification)), + None, input_cl, magnified_cl, np.float32(0), @@ -96,51 +114,61 @@ class eSRRF(LiquidEngine): np.float32(magnification), np.float32(magnification)).wait() - rc_knl(cl_queue, - (n_slices,), - None, + 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, - np.int32(image.shape[1]), - np.int32(image.shape[2])).wait() + row_kernel_cl, + np.int32(2)).wait() cr_knl(cl_queue, - (n_slices, int(image.shape[1]*magnification*2), int(image.shape[2]*magnification*2)), - self.get_work_group(dc, (n_slices, image.shape[1]*magnification*2, image.shape[2]*magnification*2)), + (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*2), - np.float32(magnification*2)).wait() + np.float32(magnification*grad_magnification), + np.float32(magnification*grad_magnification)).wait() cr_knl(cl_queue, - (n_slices, int(image.shape[1]*magnification*2), int(image.shape[2]*magnification*2)), - self.get_work_group(dc, (n_slices, image.shape[1]*magnification*2, image.shape[2]*magnification*2)), + (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*2), - np.float32(magnification*2)).wait() + np.float32(magnification*grad_magnification), + np.float32(magnification*grad_magnification)).wait() rgc_knl(cl_queue, - (n_slices, (image.shape[1]*magnification-magnification*2) - magnification*2, image.shape[2]*magnification-magnification*2 - magnification*2), - self.get_work_group(dc, (n_slices, (image.shape[1]*magnification-magnification*2) - magnification*2, image.shape[2]*magnification-magnification*2 - magnification*2)), + (n_slices, highest_row - lowest_row, highest_col - lowest_col), + None, col_magnified_gradients_cl, row_magnified_gradients_cl, magnified_cl, output_cl, - np.int32(image.shape[2]*magnification), - np.int32(image.shape[1]*magnification), + np.int32(output_shape[2]), # Ensure correct dimensions + np.int32(output_shape[1]), np.int32(magnification), - np.float32(2), + np.float32(grad_magnification), np.float32(radius), - np.float32(2 * (radius / 2.355) + 1), - np.float32(2 * (radius / 2.355) * (radius / 2.355)), + np.float32(2 * (radius / 2.355) + 1), # Match sigma calculation + np.float32(2 * (radius / 2.355) ** 2), np.float32(sensitivity), - np.int32(doIntensityWeighting)).wait() + 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() @@ -153,7 +181,7 @@ class eSRRF(LiquidEngine): return output_image % for sch in schedulers: - def _run_${sch}(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_${sch}(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -161,32 +189,40 @@ class eSRRF(LiquidEngine): """ runtype = "${sch}".capitalize() crsm = ShiftAndMagnify(verbose=False) - rbc = GradientRobertsCross(verbose=False) + conv = Convolution(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, 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) + 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) return radial_gradients % endfor - def _run_unthreaded(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_unthreaded(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @cython """ runtype = "Unthreaded" crsm = ShiftAndMagnify(verbose=False) - rbc = GradientRobertsCross(verbose=False) + conv = Convolution(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, 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) + 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) return radial_gradients \ No newline at end of file diff --git a/src/mako_templates/nanopyx.core.transform._le_esrrf3d.pyx b/src/mako_templates/nanopyx.core.transform._le_esrrf3d.pyx index d46f2979..78cb801c 100644 --- a/src/mako_templates/nanopyx.core.transform._le_esrrf3d.pyx +++ b/src/mako_templates/nanopyx.core.transform._le_esrrf3d.pyx @@ -1,6 +1,4 @@ -<%! -schedulers = ['threaded','threaded_guided','threaded_dynamic','threaded_static', 'unthreaded'] -%># cython: infer_types=True, wraparound=False, nonecheck=False, boundscheck=False, cdivision=True, language_level=3, profile=False, autogen_pxd=False +# cython: infer_types=True, wraparound=False, nonecheck=False, boundscheck=False, cdivision=True, language_level=3, profile=False, autogen_pxd=False import numpy as np import math @@ -57,13 +55,10 @@ class eSRRF3D(LiquidEngine): image = image.reshape((1, image.shape[0], image.shape[1], image.shape[2])) return super().benchmark(image, magnification_xy=magnification_xy, magnification_z=magnification_z, radius=radius, PSF_voxel_ratio=PSF_voxel_ratio, sensitivity=sensitivity, doIntensityWeighting=doIntensityWeighting) - % for sch in schedulers: - def _run_${sch}(self, float[:,:,:,:] image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_voxel_ratio: float = 4.0, sensitivity: float = 1, doIntensityWeighting: bool = True): + def _run_threaded(self, float[:,:,:,:] image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_voxel_ratio: float = 4.0, sensitivity: float = 1, doIntensityWeighting: bool = True): """ @cpu - % if sch!='unthreaded': @threaded - % endif @cython """ cdef float sigma = radius / 2.355 @@ -118,13 +113,7 @@ class eSRRF3D(LiquidEngine): with nogil: for sM in range(0, n_slices_mag): - % if sch=="unthreaded": - for rM in range(0, n_rows_mag): - % elif sch=="threaded": for rM in prange(0, n_rows_mag): - % else: - for rM in prange(0, n_rows_mag, schedule="${sch.split('_')[1]}"): - % endif for cM in range(0, n_cols_mag): if _doIntensityWeighting: rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) @@ -134,7 +123,277 @@ class eSRRF3D(LiquidEngine): rgc_map[f, sM, rM, cM] = rgc_val return np.asarray(rgc_map) - % endfor + def _run_threaded_guided(self, float[:,:,:,:] image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_voxel_ratio: float = 4.0, sensitivity: float = 1, doIntensityWeighting: bool = True): + """ + @cpu + @threaded + @cython + """ + cdef float sigma = radius / 2.355 + cdef float fwhm = radius + cdef float tSS = 2 * sigma * sigma + cdef float tSO = 2 * sigma + 1 + cdef float sigma_z = radius * PSF_voxel_ratio / 2.355 # Taking voxel size into account + cdef float fwhm_z = radius * PSF_voxel_ratio + cdef float tSS_z = 2 * sigma_z * sigma_z + cdef float tSO_z = 2 * sigma_z + 1 + cdef int Gx_Gy_MAGNIFICATION = 2 + cdef int _magnification_xy = magnification_xy + cdef int _magnification_z = magnification_z + cdef int _doIntensityWeighting = doIntensityWeighting + + cdef int n_frames, n_slices, n_rows, n_cols, n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum + n_frames, n_slices, n_rows, n_cols = image.shape[0], image.shape[1], image.shape[2], image.shape[3] + + cdef float[:, :, :, :] rgc_map = np.zeros((n_frames, n_slices*magnification_z, n_rows*magnification_xy, n_cols*magnification_xy), dtype=np.float32) + + cdef float[:, :, :] image_interpolated, gradients_s, gradients_r, gradients_c, gradients_s_interpolated, gradients_r_interpolated, gradients_c_interpolated, padded, img_dum + + cdef int f, n_slices_mag, n_rows_mag, n_cols_mag, sM, rM, cM, z0 + + cdef float rgc_val, zcof + + for f in range(n_frames): + image_interpolated = interpolate_3d_zlinear(image[f,:,:,:], _magnification_xy, _magnification_z) + + n_slices_mag, n_rows_mag, n_cols_mag = image_interpolated.shape[0], image_interpolated.shape[1], image_interpolated.shape[2] + + img_dum = interpolate_3d_zlinear(image[f], Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION) + n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum = img_dum.shape[0], img_dum.shape[1], img_dum.shape[2] + + gradients_c = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + gradients_r = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + gradients_s = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + with nogil: + _c_gradient_3d(&img_dum[0, 0, 0], &gradients_c[0, 0, 0], &gradients_r[0, 0, 0], &gradients_s[0, 0, 0], n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum) + + gradients_s_interpolated = interpolate_3d_zlinear(gradients_s, _magnification_xy, _magnification_z) + gradients_r_interpolated = interpolate_3d_zlinear(gradients_r, _magnification_xy, _magnification_z) + gradients_c_interpolated = interpolate_3d_zlinear(gradients_c, _magnification_xy, _magnification_z) + + if self.keep_gradients: + self._gradients_s_interpolated = gradients_s_interpolated.copy() + self._gradients_r_interpolated = gradients_r_interpolated.copy() + self._gradients_c_interpolated = gradients_c_interpolated.copy() + + if self.keep_interpolated: + self._img_interpolated = img_dum.copy() + + with nogil: + for sM in range(0, n_slices_mag): + for rM in prange(0, n_rows_mag, schedule="guided"): + for cM in range(0, n_cols_mag): + if _doIntensityWeighting: + rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) + rgc_map[f, sM, rM, cM] = rgc_val * image_interpolated[sM, rM, cM] + else: + rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) + rgc_map[f, sM, rM, cM] = rgc_val + + return np.asarray(rgc_map) + def _run_threaded_dynamic(self, float[:,:,:,:] image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_voxel_ratio: float = 4.0, sensitivity: float = 1, doIntensityWeighting: bool = True): + """ + @cpu + @threaded + @cython + """ + cdef float sigma = radius / 2.355 + cdef float fwhm = radius + cdef float tSS = 2 * sigma * sigma + cdef float tSO = 2 * sigma + 1 + cdef float sigma_z = radius * PSF_voxel_ratio / 2.355 # Taking voxel size into account + cdef float fwhm_z = radius * PSF_voxel_ratio + cdef float tSS_z = 2 * sigma_z * sigma_z + cdef float tSO_z = 2 * sigma_z + 1 + cdef int Gx_Gy_MAGNIFICATION = 2 + cdef int _magnification_xy = magnification_xy + cdef int _magnification_z = magnification_z + cdef int _doIntensityWeighting = doIntensityWeighting + + cdef int n_frames, n_slices, n_rows, n_cols, n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum + n_frames, n_slices, n_rows, n_cols = image.shape[0], image.shape[1], image.shape[2], image.shape[3] + + cdef float[:, :, :, :] rgc_map = np.zeros((n_frames, n_slices*magnification_z, n_rows*magnification_xy, n_cols*magnification_xy), dtype=np.float32) + + cdef float[:, :, :] image_interpolated, gradients_s, gradients_r, gradients_c, gradients_s_interpolated, gradients_r_interpolated, gradients_c_interpolated, padded, img_dum + + cdef int f, n_slices_mag, n_rows_mag, n_cols_mag, sM, rM, cM, z0 + + cdef float rgc_val, zcof + + for f in range(n_frames): + image_interpolated = interpolate_3d_zlinear(image[f,:,:,:], _magnification_xy, _magnification_z) + + n_slices_mag, n_rows_mag, n_cols_mag = image_interpolated.shape[0], image_interpolated.shape[1], image_interpolated.shape[2] + + img_dum = interpolate_3d_zlinear(image[f], Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION) + n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum = img_dum.shape[0], img_dum.shape[1], img_dum.shape[2] + + gradients_c = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + gradients_r = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + gradients_s = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + with nogil: + _c_gradient_3d(&img_dum[0, 0, 0], &gradients_c[0, 0, 0], &gradients_r[0, 0, 0], &gradients_s[0, 0, 0], n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum) + + gradients_s_interpolated = interpolate_3d_zlinear(gradients_s, _magnification_xy, _magnification_z) + gradients_r_interpolated = interpolate_3d_zlinear(gradients_r, _magnification_xy, _magnification_z) + gradients_c_interpolated = interpolate_3d_zlinear(gradients_c, _magnification_xy, _magnification_z) + + if self.keep_gradients: + self._gradients_s_interpolated = gradients_s_interpolated.copy() + self._gradients_r_interpolated = gradients_r_interpolated.copy() + self._gradients_c_interpolated = gradients_c_interpolated.copy() + + if self.keep_interpolated: + self._img_interpolated = img_dum.copy() + + with nogil: + for sM in range(0, n_slices_mag): + for rM in prange(0, n_rows_mag, schedule="dynamic"): + for cM in range(0, n_cols_mag): + if _doIntensityWeighting: + rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) + rgc_map[f, sM, rM, cM] = rgc_val * image_interpolated[sM, rM, cM] + else: + rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) + rgc_map[f, sM, rM, cM] = rgc_val + + return np.asarray(rgc_map) + def _run_threaded_static(self, float[:,:,:,:] image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_voxel_ratio: float = 4.0, sensitivity: float = 1, doIntensityWeighting: bool = True): + """ + @cpu + @threaded + @cython + """ + cdef float sigma = radius / 2.355 + cdef float fwhm = radius + cdef float tSS = 2 * sigma * sigma + cdef float tSO = 2 * sigma + 1 + cdef float sigma_z = radius * PSF_voxel_ratio / 2.355 # Taking voxel size into account + cdef float fwhm_z = radius * PSF_voxel_ratio + cdef float tSS_z = 2 * sigma_z * sigma_z + cdef float tSO_z = 2 * sigma_z + 1 + cdef int Gx_Gy_MAGNIFICATION = 2 + cdef int _magnification_xy = magnification_xy + cdef int _magnification_z = magnification_z + cdef int _doIntensityWeighting = doIntensityWeighting + + cdef int n_frames, n_slices, n_rows, n_cols, n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum + n_frames, n_slices, n_rows, n_cols = image.shape[0], image.shape[1], image.shape[2], image.shape[3] + + cdef float[:, :, :, :] rgc_map = np.zeros((n_frames, n_slices*magnification_z, n_rows*magnification_xy, n_cols*magnification_xy), dtype=np.float32) + + cdef float[:, :, :] image_interpolated, gradients_s, gradients_r, gradients_c, gradients_s_interpolated, gradients_r_interpolated, gradients_c_interpolated, padded, img_dum + + cdef int f, n_slices_mag, n_rows_mag, n_cols_mag, sM, rM, cM, z0 + + cdef float rgc_val, zcof + + for f in range(n_frames): + image_interpolated = interpolate_3d_zlinear(image[f,:,:,:], _magnification_xy, _magnification_z) + + n_slices_mag, n_rows_mag, n_cols_mag = image_interpolated.shape[0], image_interpolated.shape[1], image_interpolated.shape[2] + + img_dum = interpolate_3d_zlinear(image[f], Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION) + n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum = img_dum.shape[0], img_dum.shape[1], img_dum.shape[2] + + gradients_c = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + gradients_r = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + gradients_s = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + with nogil: + _c_gradient_3d(&img_dum[0, 0, 0], &gradients_c[0, 0, 0], &gradients_r[0, 0, 0], &gradients_s[0, 0, 0], n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum) + + gradients_s_interpolated = interpolate_3d_zlinear(gradients_s, _magnification_xy, _magnification_z) + gradients_r_interpolated = interpolate_3d_zlinear(gradients_r, _magnification_xy, _magnification_z) + gradients_c_interpolated = interpolate_3d_zlinear(gradients_c, _magnification_xy, _magnification_z) + + if self.keep_gradients: + self._gradients_s_interpolated = gradients_s_interpolated.copy() + self._gradients_r_interpolated = gradients_r_interpolated.copy() + self._gradients_c_interpolated = gradients_c_interpolated.copy() + + if self.keep_interpolated: + self._img_interpolated = img_dum.copy() + + with nogil: + for sM in range(0, n_slices_mag): + for rM in prange(0, n_rows_mag, schedule="static"): + for cM in range(0, n_cols_mag): + if _doIntensityWeighting: + rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) + rgc_map[f, sM, rM, cM] = rgc_val * image_interpolated[sM, rM, cM] + else: + rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) + rgc_map[f, sM, rM, cM] = rgc_val + + return np.asarray(rgc_map) + def _run_unthreaded(self, float[:,:,:,:] image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_voxel_ratio: float = 4.0, sensitivity: float = 1, doIntensityWeighting: bool = True): + """ + @cpu + @cython + """ + cdef float sigma = radius / 2.355 + cdef float fwhm = radius + cdef float tSS = 2 * sigma * sigma + cdef float tSO = 2 * sigma + 1 + cdef float sigma_z = radius * PSF_voxel_ratio / 2.355 # Taking voxel size into account + cdef float fwhm_z = radius * PSF_voxel_ratio + cdef float tSS_z = 2 * sigma_z * sigma_z + cdef float tSO_z = 2 * sigma_z + 1 + cdef int Gx_Gy_MAGNIFICATION = 2 + cdef int _magnification_xy = magnification_xy + cdef int _magnification_z = magnification_z + cdef int _doIntensityWeighting = doIntensityWeighting + + cdef int n_frames, n_slices, n_rows, n_cols, n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum + n_frames, n_slices, n_rows, n_cols = image.shape[0], image.shape[1], image.shape[2], image.shape[3] + + cdef float[:, :, :, :] rgc_map = np.zeros((n_frames, n_slices*magnification_z, n_rows*magnification_xy, n_cols*magnification_xy), dtype=np.float32) + + cdef float[:, :, :] image_interpolated, gradients_s, gradients_r, gradients_c, gradients_s_interpolated, gradients_r_interpolated, gradients_c_interpolated, padded, img_dum + + cdef int f, n_slices_mag, n_rows_mag, n_cols_mag, sM, rM, cM, z0 + + cdef float rgc_val, zcof + + for f in range(n_frames): + image_interpolated = interpolate_3d_zlinear(image[f,:,:,:], _magnification_xy, _magnification_z) + + n_slices_mag, n_rows_mag, n_cols_mag = image_interpolated.shape[0], image_interpolated.shape[1], image_interpolated.shape[2] + + img_dum = interpolate_3d_zlinear(image[f], Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION) + n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum = img_dum.shape[0], img_dum.shape[1], img_dum.shape[2] + + gradients_c = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + gradients_r = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + gradients_s = np.zeros((n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum), dtype=np.float32) + with nogil: + _c_gradient_3d(&img_dum[0, 0, 0], &gradients_c[0, 0, 0], &gradients_r[0, 0, 0], &gradients_s[0, 0, 0], n_slices_mag_dum, n_rows_mag_dum, n_cols_mag_dum) + + gradients_s_interpolated = interpolate_3d_zlinear(gradients_s, _magnification_xy, _magnification_z) + gradients_r_interpolated = interpolate_3d_zlinear(gradients_r, _magnification_xy, _magnification_z) + gradients_c_interpolated = interpolate_3d_zlinear(gradients_c, _magnification_xy, _magnification_z) + + if self.keep_gradients: + self._gradients_s_interpolated = gradients_s_interpolated.copy() + self._gradients_r_interpolated = gradients_r_interpolated.copy() + self._gradients_c_interpolated = gradients_c_interpolated.copy() + + if self.keep_interpolated: + self._img_interpolated = img_dum.copy() + + with nogil: + for sM in range(0, n_slices_mag): + for rM in range(0, n_rows_mag): + for cM in range(0, n_cols_mag): + if _doIntensityWeighting: + rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) + rgc_map[f, sM, rM, cM] = rgc_val * image_interpolated[sM, rM, cM] + else: + rgc_val = _c_calculate_rgc3D(cM, rM, sM, &gradients_c_interpolated[0,0,0], &gradients_r_interpolated[0,0,0], &gradients_s_interpolated[0,0,0], n_cols_mag, n_rows_mag, n_slices_mag, _magnification_xy, _magnification_z, PSF_voxel_ratio, Gx_Gy_MAGNIFICATION, Gx_Gy_MAGNIFICATION, fwhm, fwhm_z, tSO, tSO_z, tSS, tSS_z, sensitivity) + rgc_map[f, sM, rM, cM] = rgc_val + + return np.asarray(rgc_map) def get_gradients(self): if self._gradients_c_interpolated is None or self._gradients_r_interpolated is None or self._gradients_s_interpolated is None: @@ -143,4 +402,4 @@ class eSRRF3D(LiquidEngine): return self._gradients_c_interpolated, self._gradients_r_interpolated, self._gradients_s_interpolated def get_interpolated_image(self): - return self._img_interpolated + return self._img_interpolated \ No newline at end of file diff --git a/src/mako_templates/nanopyx.core.transform._le_radial_gradient_convergence.pyx b/src/mako_templates/nanopyx.core.transform._le_radial_gradient_convergence.pyx index 3ad05ca6..fe64b25a 100644 --- a/src/mako_templates/nanopyx.core.transform._le_radial_gradient_convergence.pyx +++ b/src/mako_templates/nanopyx.core.transform._le_radial_gradient_convergence.pyx @@ -13,7 +13,7 @@ from ...__liquid_engine__ import LiquidEngine from .__interpolation_tools__ import check_image cdef extern from "_c_sr_radial_gradient_convergence.h": - 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) nogil + 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) nogil class RadialGradientConvergence(LiquidEngine): """ @@ -27,29 +27,30 @@ class RadialGradientConvergence(LiquidEngine): verbose=verbose) - def run(self, gradient_col_interp, gradient_row_interp, image_interp, magnification: int = 5, radius: float = 1.5, sensitivity: float = 1 , doIntensityWeighting: bool = True, run_type = None): + def run(self, gradient_col_interp, gradient_row_interp, image_interp, magnification: int = 5, grad_magnification: int = 1, radius: float = 1.5, sensitivity: float = 1 , doIntensityWeighting: bool = True, offset: float = 0, xyoffset: float = 0, angle: float = 0, run_type = None): gradient_col_interp = check_image(gradient_col_interp) gradient_row_interp = check_image(gradient_row_interp) image_interp = check_image(image_interp) - return self._run(gradient_col_interp, gradient_row_interp, image_interp, magnification, radius, sensitivity, doIntensityWeighting, run_type=run_type) + return self._run(gradient_col_interp, gradient_row_interp, image_interp, magnification, grad_magnification, radius, sensitivity, doIntensityWeighting, offset, xyoffset, angle, run_type=run_type) - def benchmark(self, gradient_col_interp, gradient_row_interp, image_interp, magnification: int = 5, radius: float = 1.5, sensitivity: float = 1 , doIntensityWeighting: bool = True): + def benchmark(self, gradient_col_interp, gradient_row_interp, image_interp, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1 , doIntensityWeighting: bool = True, offset: float = 0, xyoffset: float = 0, angle: float = 0): gradient_col_interp = check_image(gradient_col_interp) gradient_row_interp = check_image(gradient_row_interp) image_interp = check_image(image_interp) - return super().benchmark(gradient_col_interp, gradient_row_interp, image_interp, magnification, radius, sensitivity, doIntensityWeighting) + return super().benchmark(gradient_col_interp, gradient_row_interp, image_interp, magnification, grad_magnification, radius, sensitivity, doIntensityWeighting, offset, xyoffset, angle) - def _run_unthreaded(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_unthreaded(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): """ @cpu @cython """ cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting @@ -63,17 +64,17 @@ class RadialGradientConvergence(LiquidEngine): cdef int f, rM, cM with nogil: for f in range(nFrames): - for rM in range(_magnification*2, rowsM - _magnification*2): - for cM in range(_magnification*2, colsM - _magnification*2): + 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) * 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, offset, xyoffset, angle) * 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) + 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) 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, radius=1.5, sensitivity=1, doIntensityWeighting=True): + 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): """ @cpu @threaded @@ -81,9 +82,10 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting @@ -98,20 +100,20 @@ class RadialGradientConvergence(LiquidEngine): with nogil: for f in range(nFrames): % if sch=='threaded': - for rM in prange(_magnification*2, rowsM - _magnification*2): + for rM in prange(margin, rowsM - margin): % else: - for rM in prange(_magnification*2, rowsM - _magnification*2, schedule="${sch.split('_')[1]}"): + for rM in prange(margin, rowsM - margin, schedule="${sch.split('_')[1]}"): % endif - for cM in range(_magnification*2, colsM - _magnification*2): + 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) * 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, offset, xyoffset, angle) * 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) + 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) return np.asarray(rgc_map,dtype=np.float32) % endfor - 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): + 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): """ @gpu """ @@ -129,12 +131,16 @@ class RadialGradientConvergence(LiquidEngine): # Parameters cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting + cdef float _offset = offset + cdef float _xyoffset = xyoffset + cdef float _angle = angle # Sizes cdef int nFrames = gradient_col_interp.shape[0] @@ -142,10 +148,10 @@ class RadialGradientConvergence(LiquidEngine): cdef int cols_interpolated = (gradient_row_interp.shape[2] / Gx_Gy_MAGNIFICATION) # Grid size of the global work space - lowest_row = _magnification*2 - highest_row = rows_interpolated - _magnification*2 - lowest_col = _magnification*2 - highest_col = cols_interpolated - _magnification*2 + lowest_row = margin # TODO discuss edges calculation + highest_row = rows_interpolated - margin + lowest_col = margin + highest_col = cols_interpolated - margin # Output rgc_map = np.zeros((nFrames, rows_interpolated, cols_interpolated), dtype=np.float32) @@ -194,7 +200,10 @@ class RadialGradientConvergence(LiquidEngine): np.float32(tSO), np.float32(tSS), np.float32(_sensitivity), - np.int32(_doIntensityWeighting)).wait() + np.int32(_doIntensityWeighting), + np.float32(_offset), + np.float32(_xyoffset), + np.float32(_angle)).wait() # Copy output cl.enqueue_copy(cl_queue, rgc_map[i:i+n_slices,:,:], rgc_map_out).wait() diff --git a/src/nanopyx/core/analysis/_le_channel_registration.pyx b/src/nanopyx/core/analysis/_le_channel_registration.pyx index 09c6e4c0..bbf79d8d 100644 --- a/src/nanopyx/core/analysis/_le_channel_registration.pyx +++ b/src/nanopyx/core/analysis/_le_channel_registration.pyx @@ -83,10 +83,9 @@ def _gaussian_filter(inpimg, _runtype, sigma): knl1 = np.exp(-0.5 / (sigma*sigma) * x1 ** 2) knl1 = knl1 / knl1.sum() knl1 = knl1.astype(np.float32) - img1 = conv.run(inpimg,knl1.reshape((len(x1), 1)),run_type=_runtype) + img1 = np.expand_dims(img1, axis=0) img2 = conv.run(img1,knl1.reshape((1, len(x1))),run_type=_runtype) - return img2 class ChannelRegistrationEstimator(LiquidEngine): diff --git a/src/nanopyx/core/transform/_le_convolution.cl b/src/nanopyx/core/transform/_le_convolution.cl index 61913d38..62380cf6 100644 --- a/src/nanopyx/core/transform/_le_convolution.cl +++ b/src/nanopyx/core/transform/_le_convolution.cl @@ -25,3 +25,36 @@ conv2d(__read_only image2d_t image_in, __write_only image2d_t image_out, __read_ write_imagef(image_out,image_coords,acc); } +__kernel void +conv2d_2(__global float *image, __global float *image_out, __global float *kernel_array, int kernel_size){ + + int frame = get_global_id(0); + int row = get_global_id(1); + int col = get_global_id(2); + + int nframes = get_global_size(0); + int nrows = get_global_size(1); + int ncols = get_global_size(2); + + int kernel_center = (kernel_size-1)/2; + + float acc = 0; + int localrow; + int localcol; + + for (int kr = 0; kr= max_slices: n_slices = max_slices @@ -86,7 +104,7 @@ class eSRRF(LiquidEngine): cr_knl(cl_queue, (n_slices, int(image.shape[1]*magnification), int(image.shape[2]*magnification)), - self.get_work_group(dc, (n_slices, image.shape[1]*magnification, image.shape[2]*magnification)), + None, input_cl, magnified_cl, np.float32(0), @@ -94,51 +112,61 @@ class eSRRF(LiquidEngine): np.float32(magnification), np.float32(magnification)).wait() - rc_knl(cl_queue, - (n_slices,), - None, + 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, - np.int32(image.shape[1]), - np.int32(image.shape[2])).wait() + row_kernel_cl, + np.int32(2)).wait() cr_knl(cl_queue, - (n_slices, int(image.shape[1]*magnification*2), int(image.shape[2]*magnification*2)), - self.get_work_group(dc, (n_slices, image.shape[1]*magnification*2, image.shape[2]*magnification*2)), + (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*2), - np.float32(magnification*2)).wait() + np.float32(magnification*grad_magnification), + np.float32(magnification*grad_magnification)).wait() cr_knl(cl_queue, - (n_slices, int(image.shape[1]*magnification*2), int(image.shape[2]*magnification*2)), - self.get_work_group(dc, (n_slices, image.shape[1]*magnification*2, image.shape[2]*magnification*2)), + (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*2), - np.float32(magnification*2)).wait() + np.float32(magnification*grad_magnification), + np.float32(magnification*grad_magnification)).wait() rgc_knl(cl_queue, - (n_slices, (image.shape[1]*magnification-magnification*2) - magnification*2, image.shape[2]*magnification-magnification*2 - magnification*2), - self.get_work_group(dc, (n_slices, (image.shape[1]*magnification-magnification*2) - magnification*2, image.shape[2]*magnification-magnification*2 - magnification*2)), + (n_slices, highest_row - lowest_row, highest_col - lowest_col), + None, col_magnified_gradients_cl, row_magnified_gradients_cl, magnified_cl, output_cl, - np.int32(image.shape[2]*magnification), - np.int32(image.shape[1]*magnification), + np.int32(output_shape[2]), # Ensure correct dimensions + np.int32(output_shape[1]), np.int32(magnification), - np.float32(2), + np.float32(grad_magnification), np.float32(radius), - np.float32(2 * (radius / 2.355) + 1), - np.float32(2 * (radius / 2.355) * (radius / 2.355)), + np.float32(2 * (radius / 2.355) + 1), # Match sigma calculation + np.float32(2 * (radius / 2.355) ** 2), np.float32(sensitivity), - np.int32(doIntensityWeighting)).wait() + 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() @@ -150,7 +178,7 @@ class eSRRF(LiquidEngine): return output_image - def _run_threaded(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_threaded(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -158,17 +186,21 @@ class eSRRF(LiquidEngine): """ runtype = "threaded".capitalize() crsm = ShiftAndMagnify(verbose=False) - rbc = GradientRobertsCross(verbose=False) + conv = Convolution(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, 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) + 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) return radial_gradients - def _run_threaded_guided(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_threaded_guided(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -176,17 +208,21 @@ class eSRRF(LiquidEngine): """ runtype = "threaded_guided".capitalize() crsm = ShiftAndMagnify(verbose=False) - rbc = GradientRobertsCross(verbose=False) + conv = Convolution(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, 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) + 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) return radial_gradients - def _run_threaded_dynamic(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_threaded_dynamic(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -194,17 +230,21 @@ class eSRRF(LiquidEngine): """ runtype = "threaded_dynamic".capitalize() crsm = ShiftAndMagnify(verbose=False) - rbc = GradientRobertsCross(verbose=False) + conv = Convolution(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, 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) + 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) return radial_gradients - def _run_threaded_static(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_threaded_static(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @threaded @@ -212,31 +252,39 @@ class eSRRF(LiquidEngine): """ runtype = "threaded_static".capitalize() crsm = ShiftAndMagnify(verbose=False) - rbc = GradientRobertsCross(verbose=False) + conv = Convolution(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, 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) + 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) return radial_gradients - def _run_unthreaded(self, image, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_unthreaded(self, image, magnification=5, grad_magnification=2, radius=1.5, sensitivity=1, doIntensityWeighting=True): """ @cpu @cython """ runtype = "Unthreaded" crsm = ShiftAndMagnify(verbose=False) - rbc = GradientRobertsCross(verbose=False) + conv = Convolution(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, 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) + 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) return radial_gradients \ No newline at end of file diff --git a/src/nanopyx/core/transform/_le_esrrf3d.pyx b/src/nanopyx/core/transform/_le_esrrf3d.pyx index 73177569..78cb801c 100644 --- a/src/nanopyx/core/transform/_le_esrrf3d.pyx +++ b/src/nanopyx/core/transform/_le_esrrf3d.pyx @@ -402,4 +402,4 @@ class eSRRF3D(LiquidEngine): return self._gradients_c_interpolated, self._gradients_r_interpolated, self._gradients_s_interpolated def get_interpolated_image(self): - return self._img_interpolated + return self._img_interpolated \ No newline at end of file diff --git a/src/nanopyx/core/transform/_le_radial_gradient_convergence.cl b/src/nanopyx/core/transform/_le_radial_gradient_convergence.cl index 258e473c..746b34ef 100644 --- a/src/nanopyx/core/transform/_le_radial_gradient_convergence.cl +++ b/src/nanopyx/core/transform/_le_radial_gradient_convergence.cl @@ -1,7 +1,9 @@ -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 _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 @@ -19,38 +21,62 @@ double _c_calculate_dk(float Gx, float Gy, float dx, float dy, float 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 _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 xc = (xM + 0.5) / magnification; - float yc = (yM + 0.5) / magnification; + 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)(Gx_Gy_MAGNIFICATION * fwhm); - int _end = (int)(Gx_Gy_MAGNIFICATION * fwhm + 1); + int _start = -(int)(2 * fwhm); //changed to have "Factor Cagança" but with properly iterating over the desired range + int _end = (int)(2 * fwhm + 1); // TODO discuss with Ricardo for (int j = _start; j < _end; j++) { - vy = (int)(Gx_Gy_MAGNIFICATION * yc) + j; - vy /= Gx_Gy_MAGNIFICATION; + vy = yc + j; - if (0 < vy && vy <= rowsM - 1) { + if (0 < vy && vy <= (rowsM/magnification) - 1) { for (int i = _start; i < _end; i++) { - vx = (int)(Gx_Gy_MAGNIFICATION * xc) + i; - vx /= Gx_Gy_MAGNIFICATION; + vx = xc + i; - if (0 < vx && vx <= colsM - 1) { + if (0 < vx && vx <= (colsM/magnification) - 1) { dx = vx - xc; dy = vy - yc; distance = sqrt(dx * dx + dy * dy); if (distance != 0 && distance <= tSO) { - Gx = imIntGx[(int)(vy * magnification * Gx_Gy_MAGNIFICATION * colsM * Gx_Gy_MAGNIFICATION) + (int)(vx * magnification * Gx_Gy_MAGNIFICATION)]; - Gy = imIntGy[(int)(vy * magnification * Gx_Gy_MAGNIFICATION * colsM * Gx_Gy_MAGNIFICATION) + (int)(vx * magnification * Gx_Gy_MAGNIFICATION)]; + Gx = imIntGx[(int)((vy+xyoffset) * magnification * Gx_Gy_MAGNIFICATION * colsM * Gx_Gy_MAGNIFICATION) + (int)((vx+xyoffset) * magnification * Gx_Gy_MAGNIFICATION)]; + Gy = imIntGy[(int)((vy+xyoffset) * magnification * Gx_Gy_MAGNIFICATION * colsM * Gx_Gy_MAGNIFICATION) + (int)((vx+xyoffset) * magnification * Gx_Gy_MAGNIFICATION)]; + + // Rotate the gradient components + float2 rotatedG = _rotate_vector(Gx, Gy, angle); + Gx = rotatedG.x; + Gy = rotatedG.y; distanceWeight = _c_calculate_dw(distance, tSS); distanceWeightSum += distanceWeight; @@ -78,7 +104,7 @@ float _c_calculate_rgc(int xM, int yM, __global float* imIntGx, __global float* } - __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) { + __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); @@ -91,13 +117,13 @@ float _c_calculate_rgc(int xM, int yM, __global float* imIntGx, __global float* // gradient image dimensions int nPixels_grad = nRows*Gx_Gy_MAGNIFICATION * nCols*Gx_Gy_MAGNIFICATION; - row = row + magnification*2; - col = col + magnification*2; + row = row + fwhm*2*magnification; + col = col + fwhm*2*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) * imInt[f * nPixels_out + row * nCols + col]; + 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); + 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/core/transform/_le_radial_gradient_convergence.pyx b/src/nanopyx/core/transform/_le_radial_gradient_convergence.pyx index 0e9d1ab1..f12e82d3 100644 --- a/src/nanopyx/core/transform/_le_radial_gradient_convergence.pyx +++ b/src/nanopyx/core/transform/_le_radial_gradient_convergence.pyx @@ -11,7 +11,7 @@ from ...__liquid_engine__ import LiquidEngine from .__interpolation_tools__ import check_image cdef extern from "_c_sr_radial_gradient_convergence.h": - 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) nogil + 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) nogil class RadialGradientConvergence(LiquidEngine): """ @@ -25,29 +25,30 @@ class RadialGradientConvergence(LiquidEngine): verbose=verbose) - def run(self, gradient_col_interp, gradient_row_interp, image_interp, magnification: int = 5, radius: float = 1.5, sensitivity: float = 1 , doIntensityWeighting: bool = True, run_type = None): + def run(self, gradient_col_interp, gradient_row_interp, image_interp, magnification: int = 5, grad_magnification: int = 1, radius: float = 1.5, sensitivity: float = 1 , doIntensityWeighting: bool = True, offset: float = 0, xyoffset: float = 0, angle: float = 0, run_type = None): gradient_col_interp = check_image(gradient_col_interp) gradient_row_interp = check_image(gradient_row_interp) image_interp = check_image(image_interp) - return self._run(gradient_col_interp, gradient_row_interp, image_interp, magnification, radius, sensitivity, doIntensityWeighting, run_type=run_type) + return self._run(gradient_col_interp, gradient_row_interp, image_interp, magnification, grad_magnification, radius, sensitivity, doIntensityWeighting, offset, xyoffset, angle, run_type=run_type) - def benchmark(self, gradient_col_interp, gradient_row_interp, image_interp, magnification: int = 5, radius: float = 1.5, sensitivity: float = 1 , doIntensityWeighting: bool = True): + def benchmark(self, gradient_col_interp, gradient_row_interp, image_interp, magnification: int = 5, grad_magnification: int = 2, radius: float = 1.5, sensitivity: float = 1 , doIntensityWeighting: bool = True, offset: float = 0, xyoffset: float = 0, angle: float = 0): gradient_col_interp = check_image(gradient_col_interp) gradient_row_interp = check_image(gradient_row_interp) image_interp = check_image(image_interp) - return super().benchmark(gradient_col_interp, gradient_row_interp, image_interp, magnification, radius, sensitivity, doIntensityWeighting) + return super().benchmark(gradient_col_interp, gradient_row_interp, image_interp, magnification, grad_magnification, radius, sensitivity, doIntensityWeighting, offset, xyoffset, angle) - def _run_unthreaded(self, float[:,:,:] gradient_col_interp, float[:,:,:] gradient_row_interp, float[:,:,:] image_interp, magnification=5, radius=1.5, sensitivity=1, doIntensityWeighting=True): + def _run_unthreaded(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): """ @cpu @cython """ cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting @@ -61,16 +62,16 @@ class RadialGradientConvergence(LiquidEngine): cdef int f, rM, cM with nogil: for f in range(nFrames): - for rM in range(_magnification*2, rowsM - _magnification*2): - for cM in range(_magnification*2, colsM - _magnification*2): + 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) * 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, offset, xyoffset, angle) * 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) + 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) 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, radius=1.5, sensitivity=1, doIntensityWeighting=True): + 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): """ @cpu @threaded @@ -78,9 +79,10 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting @@ -94,14 +96,14 @@ class RadialGradientConvergence(LiquidEngine): cdef int f, rM, cM with nogil: for f in range(nFrames): - for rM in prange(_magnification*2, rowsM - _magnification*2): - for cM in range(_magnification*2, colsM - _magnification*2): + 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) * 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, offset, xyoffset, angle) * 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) + 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) 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, radius=1.5, sensitivity=1, doIntensityWeighting=True): + 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): """ @cpu @threaded @@ -109,9 +111,10 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting @@ -125,14 +128,14 @@ class RadialGradientConvergence(LiquidEngine): cdef int f, rM, cM with nogil: for f in range(nFrames): - for rM in prange(_magnification*2, rowsM - _magnification*2, schedule="guided"): - for cM in range(_magnification*2, colsM - _magnification*2): + 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) * 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, offset, xyoffset, angle) * 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) + 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) 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, radius=1.5, sensitivity=1, doIntensityWeighting=True): + 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): """ @cpu @threaded @@ -140,9 +143,10 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting @@ -156,14 +160,14 @@ class RadialGradientConvergence(LiquidEngine): cdef int f, rM, cM with nogil: for f in range(nFrames): - for rM in prange(_magnification*2, rowsM - _magnification*2, schedule="dynamic"): - for cM in range(_magnification*2, colsM - _magnification*2): + 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) * 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, offset, xyoffset, angle) * 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) + 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) 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, radius=1.5, sensitivity=1, doIntensityWeighting=True): + 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): """ @cpu @threaded @@ -171,9 +175,10 @@ class RadialGradientConvergence(LiquidEngine): """ cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting @@ -187,16 +192,16 @@ class RadialGradientConvergence(LiquidEngine): cdef int f, rM, cM with nogil: for f in range(nFrames): - for rM in prange(_magnification*2, rowsM - _magnification*2, schedule="static"): - for cM in range(_magnification*2, colsM - _magnification*2): + 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) * 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, offset, xyoffset, angle) * 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) + 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) return np.asarray(rgc_map,dtype=np.float32) - 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): + 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): """ @gpu """ @@ -214,12 +219,16 @@ class RadialGradientConvergence(LiquidEngine): # Parameters cdef float sigma = radius / 2.355 cdef float fwhm = radius + cdef int margin = int(fwhm*2) * magnification cdef float tSS = 2 * sigma * sigma cdef float tSO = 2 * sigma + 1 - cdef float Gx_Gy_MAGNIFICATION = 2.0 + cdef float Gx_Gy_MAGNIFICATION = grad_magnification cdef int _magnification = magnification cdef float _sensitivity = sensitivity cdef int _doIntensityWeighting = doIntensityWeighting + cdef float _offset = offset + cdef float _xyoffset = xyoffset + cdef float _angle = angle # Sizes cdef int nFrames = gradient_col_interp.shape[0] @@ -227,10 +236,10 @@ class RadialGradientConvergence(LiquidEngine): cdef int cols_interpolated = (gradient_row_interp.shape[2] / Gx_Gy_MAGNIFICATION) # Grid size of the global work space - lowest_row = _magnification*2 - highest_row = rows_interpolated - _magnification*2 - lowest_col = _magnification*2 - highest_col = cols_interpolated - _magnification*2 + lowest_row = margin # TODO discuss edges calculation + highest_row = rows_interpolated - margin + lowest_col = margin + highest_col = cols_interpolated - margin # Output rgc_map = np.zeros((nFrames, rows_interpolated, cols_interpolated), dtype=np.float32) @@ -279,7 +288,10 @@ class RadialGradientConvergence(LiquidEngine): np.float32(tSO), np.float32(tSS), np.float32(_sensitivity), - np.int32(_doIntensityWeighting)).wait() + np.int32(_doIntensityWeighting), + np.float32(_offset), + np.float32(_xyoffset), + np.float32(_angle)).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/_transonic.py b/src/nanopyx/core/transform/_transonic.py index 75663b06..96212ed5 100644 --- a/src/nanopyx/core/transform/_transonic.py +++ b/src/nanopyx/core/transform/_transonic.py @@ -4,8 +4,9 @@ @jit(backend="numba") def convolution2D(image, kernel): - nRows = image.shape[0] - nCols = image.shape[1] + nFrames = image.shape[0] + nRows = image.shape[1] + nCols = image.shape[2] nRows_kernel = kernel.shape[0] nCols_kernel = kernel.shape[1] @@ -15,16 +16,17 @@ def convolution2D(image, kernel): acc = 0.0 - conv_out = np.zeros((nRows, nCols), dtype=np.float32) - - for r in range(nRows): - for c in range(nCols): - acc = 0 - for kr in range(nRows_kernel): - for kc in range(nCols_kernel): - local_row = min(max(r+(kr-center_r), 0), nRows-1) - local_col = min(max(c+(kc-center_c), 0), nCols-1) - acc = acc + kernel[kr, kc] * image[local_row, local_col] - conv_out[r, c] = acc + conv_out = np.zeros((nFrames, nRows, nCols), dtype=np.float32) + + for f in range(nFrames): + for r in range(nRows): + for c in range(nCols): + acc = 0 + for kr in range(nRows_kernel): + for kc in range(nCols_kernel): + local_row = min(max(r+(kr-center_r), 0), nRows-1) + local_col = min(max(c+(kc-center_c), 0), nCols-1) + acc = acc + kernel[kr, kc] * image[f, local_row, local_col] + conv_out[f, r, c] = acc return conv_out diff --git a/src/nanopyx/core/transform/convolution.py b/src/nanopyx/core/transform/convolution.py index 3775911c..e850dd98 100644 --- a/src/nanopyx/core/transform/convolution.py +++ b/src/nanopyx/core/transform/convolution.py @@ -32,33 +32,34 @@ def wrapper(func): print("Optional dependecy Dask_image is not installed. Implementations using it will be ignored.") -def check_array(image: np.ndarray): - """ - Check the given image and ensure it meets the required conditions. - - Parameters: - image (numpy.ndarray): The image to be checked. - - Returns: - numpy.ndarray: The checked and potentially modified image. - - Raises: - TypeError: If the image is not of type numpy.ndarray. - ValueError: If the image is not 2D or 3D. - """ - image = np.asarray(image) - if type(image) is not np.ndarray: - raise TypeError("Image must be of type np.ndarray") - if image.ndim != 2: - raise ValueError("Image must be 2D") - if image.dtype != np.float32: - image = image.astype(np.float32, copy=False) - return image +# def check_array(image: np.ndarray): +# """ +# Check the given image and ensure it meets the required conditions. + +# Parameters: +# image (numpy.ndarray): The image to be checked. + +# Returns: +# numpy.ndarray: The checked and potentially modified image. + +# Raises: +# TypeError: If the image is not of type numpy.ndarray. +# ValueError: If the image is not 2D or 3D. +# """ +# image = np.asarray(image) +# if type(image) is not np.ndarray: +# raise TypeError("Image must be of type np.ndarray") +# if image.ndim != 3: +# raise ValueError("Image must be 2D") +# if image.dtype != np.float32: +# image = image.astype(np.float32, copy=False) +# return image def convolution2D_python(image: np.ndarray, kernel: np.ndarray): - nRows = image.shape[0] - nCols = image.shape[1] + nFrames = image.shape[0] + nRows = image.shape[1] + nCols = image.shape[2] nRows_kernel = kernel.shape[0] nCols_kernel = kernel.shape[1] @@ -68,17 +69,18 @@ def convolution2D_python(image: np.ndarray, kernel: np.ndarray): acc = 0.0 - conv_out = np.zeros((nRows, nCols), dtype=np.float32) - - for r in range(nRows): - for c in range(nCols): - acc = 0 - for kr in range(nRows_kernel): - for kc in range(nCols_kernel): - local_row = min(max(r+(kr-center_r), 0), nRows-1) - local_col = min(max(c+(kc-center_c), 0), nCols-1) - acc = acc + kernel[kr, kc] * image[local_row, local_col] - conv_out[r, c] = acc + conv_out = np.zeros((nFrames, nRows, nCols), dtype=np.float32) + + for f in range(nFrames): + for r in range(nRows): + for c in range(nCols): + acc = 0 + for kr in range(nRows_kernel): + for kc in range(nCols_kernel): + local_row = min(max(r+(kr-center_r), 0), nRows-1) + local_col = min(max(c+(kc-center_c), 0), nCols-1) + acc = acc + kernel[kr, kc] * image[f,local_row, local_col] + conv_out[f,r, c] = acc return conv_out @@ -98,8 +100,9 @@ def convolution2D_transonic(image, kernel): @njit(cache=True, parallel=True) def convolution2D_numba(image, kernel): - nRows = image.shape[0] - nCols = image.shape[1] + nFrames = image.shape[0] + nRows = image.shape[1] + nCols = image.shape[2] nRows_kernel = kernel.shape[0] nCols_kernel = kernel.shape[1] @@ -109,23 +112,28 @@ def convolution2D_numba(image, kernel): acc = 0.0 - conv_out = np.zeros((nRows, nCols), dtype=np.float32) + conv_out = np.zeros((nFrames, nRows, nCols), dtype=np.float32) - for r in range(nRows): - for c in range(nCols): - acc = 0 - for kr in range(nRows_kernel): - for kc in range(nCols_kernel): - local_row = min(max(r+(kr-center_r), 0), nRows-1) - local_col = min(max(c+(kc-center_c), 0), nCols-1) - acc = acc + kernel[kr, kc] * image[local_row, local_col] - conv_out[r, c] = acc + for f in range(nFrames): + for r in range(nRows): + for c in range(nCols): + acc = 0 + for kr in range(nRows_kernel): + for kc in range(nCols_kernel): + local_row = min(max(r+(kr-center_r), 0), nRows-1) + local_col = min(max(c+(kc-center_c), 0), nCols-1) + acc = acc + kernel[kr, kc] * image[f,local_row, local_col] + conv_out[f,r, c] = acc return conv_out def convolution2D_dask(image, kernel): - return np.asarray(dask_convolve(da.from_array(image), da.from_array(kernel))) + + conv_out = np.zeros_like(image) + for i in range(image.shape[0]): + conv_out[i] = dask_convolve(da.from_array(image[i]), da.from_array(kernel)) + return conv_out def convolution2D_cuda(image, kernel): diff --git a/src/nanopyx/methods/channel_registration/estimator.py b/src/nanopyx/methods/channel_registration/estimator.py index 95d81e16..a2f698cc 100644 --- a/src/nanopyx/methods/channel_registration/estimator.py +++ b/src/nanopyx/methods/channel_registration/estimator.py @@ -3,7 +3,9 @@ from skimage.io import imsave from .corrector import ChannelRegistrationCorrector -from ...core.analysis._le_channel_registration import ChannelRegistrationEstimator as leChannelRegistrationEstimator +from ...core.analysis._le_channel_registration import ( + ChannelRegistrationEstimator as leChannelRegistrationEstimator, +) # this class assumes that the image is a numpy array with shape = [n_channels, height, width] @@ -31,10 +33,11 @@ class ChannelRegistrationEstimator(object): It calculates translation masks for channel alignment and provides options for saving results to files. """ - def __init__(self) -> None: + def __init__(self, verbose=True) -> None: """ Initialize a ChannelRegistrationEstimator instance. """ + self.verbose = verbose self.translation_masks = None def apply_elastic_transform(self, img_stack): @@ -56,8 +59,9 @@ def apply_elastic_transform(self, img_stack): for aligning the channels in the provided image stack based on the stored translation masks. """ corrector = ChannelRegistrationCorrector() - return corrector.align_channels(img_stack, translation_masks=self.translation_masks) - + return corrector.align_channels( + img_stack, translation_masks=self.translation_masks + ) def save_translation_mask(self, path=None): """ @@ -77,7 +81,12 @@ def save_translation_mask(self, path=None): it prompts the user to input a file path and appends "_translation_masks.tif" to it as the default file name. """ if path is None: - path = input("Please provide a filepath to save the translation masks") + "_translation_masks.tif" + path = ( + input( + "Please provide a filepath to save the translation masks" + ) + + "_translation_masks.tif" + ) imsave(path + "_translation_masks.tif", self.translation_masks) @@ -123,11 +132,19 @@ def estimate( """ if ref_channel > img_stack.shape[0]: - print("Reference channel number cannot be bigger than number of channels!") + print( + "Reference channel number cannot be bigger than number of channels!" + ) return None - estimator = leChannelRegistrationEstimator() - self.translation_masks = estimator.run(np.asarray(img_stack, dtype=np.float32), ref_channel, max_shift, blocks_per_axis, min_similarity) + estimator = leChannelRegistrationEstimator(verbose=self.verbose) + self.translation_masks = estimator.run( + np.asarray(img_stack, dtype=np.float32), + ref_channel, + max_shift, + blocks_per_axis, + min_similarity, + ) if save_translation_masks: self.save_translation_mask(path=translation_mask_save_path) diff --git a/src/nanopyx/methods/drift_alignment/estimator.py b/src/nanopyx/methods/drift_alignment/estimator.py index 5945bfda..78fb396a 100644 --- a/src/nanopyx/methods/drift_alignment/estimator.py +++ b/src/nanopyx/methods/drift_alignment/estimator.py @@ -8,7 +8,9 @@ from ...core.utils.timeit import timeit from ...core.analysis.ccm import calculate_ccm from ...core.analysis.rcc import rcc -from ...core.analysis._le_drift_calculator import DriftEstimator as leDriftEstimator +from ...core.analysis._le_drift_calculator import ( + DriftEstimator as leDriftEstimator, +) class DriftEstimator(object): @@ -60,7 +62,7 @@ class DriftEstimator(object): It provides methods for estimating drift using cross-correlation and applying drift correction to an image stack. """ - def __init__(self): + def __init__(self, verbose=True): """ Initialize the `DriftEstimator` object. @@ -73,6 +75,7 @@ def __init__(self): Example: estimator = DriftEstimator() """ + self.verbose = verbose self.estimator_table = DriftEstimatorTable() self.cross_correlation_map = None self.drift_xy = None @@ -109,20 +112,25 @@ def estimate(self, image_array, **kwargs): # x0, y0, x1, y1 correspond to the exact coordinates of the roi to be used or full image dims and should be a tuple if ( - self.estimator_table.params["use_roi"] and self.estimator_table.params["roi"] is not None + self.estimator_table.params["use_roi"] + and self.estimator_table.params["roi"] is not None ): # crops image to roi - print(self.estimator_table.params["use_roi"], self.estimator_table.params["roi"]) + print( + self.estimator_table.params["use_roi"], + self.estimator_table.params["roi"], + ) x0, y0, x1, y1 = tuple(self.estimator_table.params["roi"]) image_arr = image_array[:, y0 : y1 + 1, x0 : x1 + 1] else: image_arr = image_array - estimator = leDriftEstimator() + estimator = leDriftEstimator(verbose=self.verbose) self.estimator_table.drift_table = estimator.run( np.asarray(image_arr, dtype=np.float32), time_averaging=self.estimator_table.params["time_averaging"], max_drift=self.estimator_table.params["max_expected_drift"], - ref_option=self.estimator_table.params["ref_option"]) + ref_option=self.estimator_table.params["ref_option"], + ) if self.estimator_table.params["apply"]: drift_corrector = DriftCorrector() diff --git a/src/nanopyx/methods/esrrf_3d/run.py b/src/nanopyx/methods/esrrf_3d/run.py index ab6c63d6..4ac718b7 100644 --- a/src/nanopyx/methods/esrrf_3d/run.py +++ b/src/nanopyx/methods/esrrf_3d/run.py @@ -1,7 +1,12 @@ from ...core.transform._le_esrrf3d import eSRRF3D -from ...core.transform.sr_temporal_correlations import calculate_eSRRF3d_temporal_correlations +from ...core.transform.sr_temporal_correlations import ( + calculate_eSRRF3d_temporal_correlations, +) -def run_esrrf3d(img,correlation="AVG", framewindow=5, rollingoverlap=2, **kwargs): + +def run_esrrf3d( + img, correlation="AVG", framewindow=5, rollingoverlap=2, **kwargs +): """ Calculate the eSRRF3D temporal correlations for the given image. @@ -24,4 +29,9 @@ def run_esrrf3d(img,correlation="AVG", framewindow=5, rollingoverlap=2, **kwargs The calculated eSRRF3D temporal correlations. """ esrrf_calculator = eSRRF3D() - return calculate_eSRRF3d_temporal_correlations(esrrf_calculator.run(img, **kwargs), correlation=correlation, framewindow=framewindow, rollingoverlap=rollingoverlap) \ No newline at end of file + return calculate_eSRRF3d_temporal_correlations( + esrrf_calculator.run(img, **kwargs), + correlation=correlation, + framewindow=framewindow, + rollingoverlap=rollingoverlap, + ) diff --git a/tests/test_conv2d.py b/tests/test_conv2d.py index 023969dd..145284d3 100644 --- a/tests/test_conv2d.py +++ b/tests/test_conv2d.py @@ -5,14 +5,22 @@ def test_conv_small(): - small = np.random.random((100, 100)).astype(np.float32) + small = np.random.random((10,100, 100)).astype(np.float32) kernel = np.ones((23, 23)).astype(np.float32) conv = Convolution() conv.benchmark(small, kernel) def test_convolution2D_numba(): - img = np.random.random((100, 100)).astype(np.float32) + img = np.random.random((10,100, 100)).astype(np.float32) kernel = np.random.random((3, 3)).astype(np.float32) convolution2D_numba(img, kernel) + + +def test_conv_small_2d(): + small = np.random.random((100, 100)).astype(np.float32) + kernel = np.ones((23, 23)).astype(np.float32) + conv = Convolution() + conv.benchmark(small, kernel) +