diff --git a/setup.py b/setup.py index fade518b..442af95e 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ ] EXTRA_LING_ARGS = [] -VERSION = "1.2.1" # sets version number for whole package +VERSION = "1.2.2" # sets version number for whole package def run_command(command: str) -> str: diff --git a/src/include/_c_sr_radial_gradient_convergence.c b/src/include/_c_sr_radial_gradient_convergence.c index 7162dd46..e7b9cee9 100644 --- a/src/include/_c_sr_radial_gradient_convergence.c +++ b/src/include/_c_sr_radial_gradient_convergence.c @@ -140,7 +140,7 @@ float _c_get_bound_value(float* im, int slices, int rows, int cols, int s, int r float _c_calculate_rgc3D(int xM, int yM, int sliceM, float* imIntGx, float* imIntGy, float* imIntGz, int colsM, int rowsM, int slicesM, int magnification_xy, int magnification_z, float voxel_ratio, float fwhm, float fwhm_z, float tSO, float tSO_z, float tSS, float tSS_z, float sensitivity) { - float vx, vy, vz, Gx, Gy, Gz, dx, dy, dz, dz_real, distance, distance_xy, distance_z, distanceWeight, distanceWeight_xy, distanceWeight_z, GdotR, Dk; + float vx, vy, vz, Gx, Gy, Gz, dx, dy, dz, distance, distance_xy, distance_z, distanceWeight, distanceWeight_xy, distanceWeight_z, GdotR, Dk; float xc = (float)(xM) / magnification_xy; float yc = (float)(yM) / magnification_xy; @@ -171,11 +171,9 @@ float _c_calculate_rgc3D(int xM, int yM, int sliceM, float* imIntGx, float* imIn if (0 < vz && vz <= (slicesM/magnification_z) - 1) { dx = vx - xc; dy = vy - yc; - dz = vz - zc; - dz_real = dz * voxel_ratio; - distance = sqrt(dx * dx + dy * dy + dz_real * dz_real); + distance_z = vz - zc; + distance = sqrt(dx * dx + dy * dy + distance_z * distance_z); distance_xy = sqrt(dx * dx + dy * dy); - distance_z = dz_real; if (distance != 0 && distance_xy <= tSO && distance_z <= tSO_z) { int linear_index = (int)(vz * magnification_z) * rowsM * colsM + @@ -189,10 +187,10 @@ float _c_calculate_rgc3D(int xM, int yM, int sliceM, float* imIntGx, float* imIn distanceWeight = _c_calculate_dw3D(distance, distance_xy, distance_z, tSS, tSS_z); distanceWeightSum += distanceWeight; - GdotR = Gx*dx + Gy*dy + Gz*dz_real; + GdotR = Gx*dx + Gy*dy + Gz*distance_z; if (GdotR < 0) { - Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, dz_real, distance); + Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, distance_z, distance); RGC += (Dk * distanceWeight); } } diff --git a/src/liquid_benchmarks/_le_esrrf3d/eSRRF3D.yml b/src/liquid_benchmarks/_le_esrrf3d/eSRRF3D.yml index cff089c5..6ed306c8 100644 --- a/src/liquid_benchmarks/_le_esrrf3d/eSRRF3D.yml +++ b/src/liquid_benchmarks/_le_esrrf3d/eSRRF3D.yml @@ -1,48 +1,53 @@ opencl: - ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'': - 2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'': - 1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})' - : - 12363612.0 - - 0.534501292015193 - - 0.4527143339801114 - - 0.4526337919814978 + ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'': + 5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'': + 1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})' + : - 50399999.99999999 + - 6.163972874986939 + - 6.691725915996358 + - 5.600578582962044 + ? '([''shape(1000, 30, 29, 30)''], {''magnification_xy'': 2, ''magnification_z'': + 2, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'': + 1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})' + : - 657719999.9999999 + - 51.01166029192973 threaded: - ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'': - 2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'': - 1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})' - : - 12363612.0 - - 1.9790987920132466 - - 1.9240770420001354 - - 1.895047333004186 + ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'': + 5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'': + 1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})' + : - 50399999.99999999 + - 19.27102233399637 + - 18.66210575005971 + - 18.99253195791971 threaded_dynamic: - ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'': - 2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'': - 1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})' - : - 12363612.0 - - 1.7660593329928815 - - 1.719053625012748 - - 1.7240022089972626 + ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'': + 5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'': + 1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})' + : - 50399999.99999999 + - 17.991310499957763 + - 18.140572750009596 + - 18.62613470805809 threaded_guided: - ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'': - 2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'': - 1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})' - : - 12363612.0 - - 1.781112333002966 - - 1.8617710830003489 - - 2.0126889999955893 + ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'': + 5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'': + 1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})' + : - 50399999.99999999 + - 18.974687915993854 + - 18.85600141703617 + - 18.480262708035298 threaded_static: - ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'': - 2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'': - 1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})' - : - 12363612.0 - - 1.8515684169833548 - - 1.8730584580043796 - - 1.8657334169838578 + ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'': + 5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'': + 1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})' + : - 50399999.99999999 + - 19.723863166058436 + - 19.532742790994234 + - 19.2554613329703 unthreaded: - ? '([''shape(1, 101, 101, 101)''], {''magnification_xy'': 2, ''magnification_z'': - 2, ''radius'': 1.5, ''radius_z'': 0.5, ''voxel_ratio'': 4.0, ''sensitivity'': - 1.0, ''mode'': ''average'', ''doIntensityWeighting'': False})' - : - 12363612.0 - - 4.863078458001837 - - 4.9081600419885945 - - 4.920950291998452 + ? '([''shape(10, 20, 40, 40)''], {''magnification_xy'': 5, ''magnification_z'': + 5, ''radius'': 1.5, ''radius_z'': 1.0499999999999998, ''voxel_ratio'': 4.0, ''sensitivity'': + 1.0, ''mode'': ''average'', ''doIntensityWeighting'': True})' + : - 50399999.99999999 + - 35.281799208023585 + - 35.78980629099533 + - 35.422467958997004 diff --git a/src/mako_templates/nanopyx.core.transform._le_esrrf3d.pyx b/src/mako_templates/nanopyx.core.transform._le_esrrf3d.pyx index d35aa728..d5aacec9 100644 --- a/src/mako_templates/nanopyx.core.transform._le_esrrf3d.pyx +++ b/src/mako_templates/nanopyx.core.transform._le_esrrf3d.pyx @@ -34,25 +34,28 @@ class eSRRF3D(LiquidEngine): self._designation = "eSRRF_3D" super().__init__(clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose) - def run(self, image, magnification_xy: int = 2, magnification_z: int = 2, radius: float = 1.5, radius_z: float = 1.5, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True, run_type=None): + def run(self, image, magnification_xy: int = 2, magnification_z: int = 2, radius: float = 1.5, PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True, run_type=None): # TODO: complete and check _run inputs, need to complete variables? + radius_z = radius * PSF_ratio / voxel_ratio if len(image.shape) == 3: image = image.reshape((1, image.shape[0], image.shape[1], image.shape[2])) if len(image.shape) != 4: print("Warning:image must either be 3D or 4D. If 3D, it will be reshaped to 4D.") return None - if radius * 2 > (image.shape[2]) / 2 or radius * 2 > (image.shape[3] / 2): - print("Warning: Radius is too big for the image. Half the radius must be smaller than both half the number of columns and half number of rows of the image.") + if radius > (image.shape[2]) / 2 or radius * 2 > (image.shape[3] / 2): + print("Warning: Radius is too big for the image. The radius must be smaller than half the number of columns and half number of rows of the image.") return None - if radius_z * 2 > image.shape[1] / 2: - print("Warning: Radius_z is too big for the image. Half the radius_z must be smaller than half of number of Z planes.") + if radius_z > image.shape[1] / 2: + print("Warning: Radius_z is too big for the image. The radius_z must be smaller than half of number of Z planes. Radius_z is automatically calculated from radius * PSF_ratio / voxel_ratio.") return None if image.dtype != np.float32: image = image.astype(np.float32) return self._run(image, magnification_xy=magnification_xy, magnification_z=magnification_z, radius=radius, radius_z=radius_z, voxel_ratio=voxel_ratio, sensitivity=sensitivity, mode=mode, doIntensityWeighting=doIntensityWeighting, run_type=run_type) - def benchmark(self, image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, radius_z: float = 1.5, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True): + def benchmark(self, image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True): + + radius_z = radius * PSF_ratio / voxel_ratio if image.dtype != np.float32: image = image.astype(np.float32) if len(image.shape) == 4: @@ -70,7 +73,6 @@ class eSRRF3D(LiquidEngine): % endif @cython """ - time_start = time.time() # calculate all constants cdef float sigma = radius / 2.355 diff --git a/src/nanopyx/core/transform/_le_esrrf3d.pyx b/src/nanopyx/core/transform/_le_esrrf3d.pyx index 7af43854..b4a8e5df 100644 --- a/src/nanopyx/core/transform/_le_esrrf3d.pyx +++ b/src/nanopyx/core/transform/_le_esrrf3d.pyx @@ -32,25 +32,28 @@ class eSRRF3D(LiquidEngine): self._designation = "eSRRF_3D" super().__init__(clear_benchmarks=clear_benchmarks, testing=testing, verbose=verbose) - def run(self, image, magnification_xy: int = 2, magnification_z: int = 2, radius: float = 1.5, radius_z: float = 1.5, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True, run_type=None): + def run(self, image, magnification_xy: int = 2, magnification_z: int = 2, radius: float = 1.5, PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True, run_type=None): # TODO: complete and check _run inputs, need to complete variables? + radius_z = radius * PSF_ratio / voxel_ratio if len(image.shape) == 3: image = image.reshape((1, image.shape[0], image.shape[1], image.shape[2])) if len(image.shape) != 4: print("Warning:image must either be 3D or 4D. If 3D, it will be reshaped to 4D.") return None - if radius * 2 > (image.shape[2]) / 2 or radius * 2 > (image.shape[3] / 2): - print("Warning: Radius is too big for the image. Half the radius must be smaller than both half the number of columns and half number of rows of the image.") + if radius > (image.shape[2]) / 2 or radius * 2 > (image.shape[3] / 2): + print("Warning: Radius is too big for the image. The radius must be smaller than half the number of columns and half number of rows of the image.") return None - if radius_z * 2 > image.shape[1] / 2: - print("Warning: Radius_z is too big for the image. Half the radius_z must be smaller than half of number of Z planes.") + if radius_z > image.shape[1] / 2: + print("Warning: Radius_z is too big for the image. The radius_z must be smaller than half of number of Z planes. Radius_z is automatically calculated from radius * PSF_ratio / voxel_ratio.") return None if image.dtype != np.float32: image = image.astype(np.float32) return self._run(image, magnification_xy=magnification_xy, magnification_z=magnification_z, radius=radius, radius_z=radius_z, voxel_ratio=voxel_ratio, sensitivity=sensitivity, mode=mode, doIntensityWeighting=doIntensityWeighting, run_type=run_type) - def benchmark(self, image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, radius_z: float = 1.5, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True): + def benchmark(self, image, magnification_xy: int = 5, magnification_z: int = 5, radius: float = 1.5, PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, mode: str = "average", doIntensityWeighting: bool = True): + + radius_z = radius * PSF_ratio / voxel_ratio if image.dtype != np.float32: image = image.astype(np.float32) if len(image.shape) == 4: @@ -65,7 +68,6 @@ class eSRRF3D(LiquidEngine): @threaded @cython """ - time_start = time.time() # calculate all constants cdef float sigma = radius / 2.355 @@ -142,7 +144,6 @@ class eSRRF3D(LiquidEngine): @threaded @cython """ - time_start = time.time() # calculate all constants cdef float sigma = radius / 2.355 @@ -219,7 +220,6 @@ class eSRRF3D(LiquidEngine): @threaded @cython """ - time_start = time.time() # calculate all constants cdef float sigma = radius / 2.355 @@ -296,7 +296,6 @@ class eSRRF3D(LiquidEngine): @threaded @cython """ - time_start = time.time() # calculate all constants cdef float sigma = radius / 2.355 @@ -372,7 +371,6 @@ class eSRRF3D(LiquidEngine): @cpu @cython """ - time_start = time.time() # calculate all constants cdef float sigma = radius / 2.355 diff --git a/src/nanopyx/core/transform/_le_esrrf3d_.cl b/src/nanopyx/core/transform/_le_esrrf3d_.cl index f5d65dfb..4ae6b83d 100644 --- a/src/nanopyx/core/transform/_le_esrrf3d_.cl +++ b/src/nanopyx/core/transform/_le_esrrf3d_.cl @@ -239,11 +239,9 @@ float _c_calculate_rgc3D(int xM, int yM, int sliceM, __global float* imIntGx, __ if (0 < vz && vz <= (slicesM/magnification_z) - 1) { dx = vx - xc; dy = vy - yc; - dz = vz - zc; - dz_real = dz * voxel_ratio; - distance = sqrt(dx * dx + dy * dy + dz_real * dz_real); + distance_z = vz - zc; + distance = sqrt(dx * dx + dy * dy + distance_z * distance_z); distance_xy = sqrt(dx * dx + dy * dy); - distance_z = dz_real; if (distance != 0 && distance_xy <= tSO && distance_z <= tSO_z) { int linear_index = (int)(vz * magnification_z) * rowsM * colsM + @@ -264,10 +262,10 @@ float _c_calculate_rgc3D(int xM, int yM, int sliceM, __global float* imIntGx, __ distanceWeightSum = distanceWeightSum + distanceWeight; // distanceWeightSum_xy += distanceWeight_xy; // distanceWeightSum_z += distanceWeight_z; - GdotR = Gx*dx + Gy*dy + Gz*dz_real; + GdotR = Gx*dx + Gy*dy + Gz*distance_z; if (GdotR < 0) { - Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, dz_real, distance); + Dk = _c_calculate_dk3D(Gx, Gy, Gz, dx, dy, distance_z, distance); RGC = RGC + (Dk * distanceWeight); } } diff --git a/src/nanopyx/methods/esrrf_3d/eSRRF3D_workflow.py b/src/nanopyx/methods/esrrf_3d/eSRRF3D_workflow.py index c524e3cc..58f128fe 100644 --- a/src/nanopyx/methods/esrrf_3d/eSRRF3D_workflow.py +++ b/src/nanopyx/methods/esrrf_3d/eSRRF3D_workflow.py @@ -10,9 +10,10 @@ def eSRRF3D( magnification_xy=2, magnification_z=2, radius: float = 1.5, - radius_z: float = 0.5, + PSF_ratio: float = 2.8, voxel_ratio: float = 4.0, sensitivity: float = 1, + frames_per_timepoint: int = 0, mode: str = "average", doIntensityWeighting: bool = True, macro_pixel_correction: bool = True, @@ -26,8 +27,9 @@ def eSRRF3D( magnification_xy (int, optional): Magnification factor in XY plane (default is 2). magnification_z (int, optional): Magnification factor in Z plane (default is 2). radius (float, optional): Radius parameter for eSRRF3D analysis (default is 1.5). - radius_z (float, optional): Radius parameter in Z direction for eSRRF3D analysis (default is 0.5). + PSF_ratio (float, optional): Ratio of PSF shape as Z/XY for eSRRF3D analysis (default is 2.8). voxel_ratio (float, optional): Ratio of voxel size in XY to Z direction (default is 4.0). + frames_per_timepoint (int, optional): Number of frames per timepoint (default is 0, which means all frames are used). sensitivity (float, optional): Sensitivity parameter for eSRRF3D analysis (default is 1). mode (str, optional): Time projection mode (default is "average"). doIntensityWeighting (bool, optional): Enable intensity weighting (default is True). @@ -43,33 +45,63 @@ def eSRRF3D( Note: - eSRRF3D (enhanced Super-Resolution Radial Fluctuations 3D) is a method for super-resolution localization microscopy in three dimensions. - This function sets up a workflow to perform eSRRF3D analysis on the input image. - + See Also: - eSRRF3D_ST: The eSRRF3D step that performs the actual analysis. - - Workflow: The class used to define and run analysis workflows. + - Workflow: The class used to define and run analysis workflows. """ - _eSRRF3D = Workflow( + if frames_per_timepoint == 0: + frames_per_timepoint = img.shape[0] + elif frames_per_timepoint > img.shape[0]: + frames_per_timepoint = img.shape[0] + + number_of_timepoints = img.shape[0] // frames_per_timepoint + if img.shape[0] % frames_per_timepoint != 0: + number_of_timepoints += 1 + + output_array = np.zeros( ( - eSRRF3D_ST(verbose=False), - (img,), - { - "magnification_xy": magnification_xy, - "magnification_z": magnification_z, - "radius": radius, - "radius_z": radius_z, - "voxel_ratio": voxel_ratio, - "sensitivity": sensitivity, - "mode": mode, - "doIntensityWeighting": doIntensityWeighting, - }, - ) + number_of_timepoints, + img.shape[1] * magnification_z, + img.shape[2] * magnification_xy, + img.shape[3] * magnification_xy, + ), + dtype=np.float32, ) - if macro_pixel_correction: - return macro_pixel_corrector( - _eSRRF3D.calculate(_force_run_type=_force_run_type)[0], - magnification=magnification_xy, + + for i in range(number_of_timepoints): + _eSRRF3D = Workflow( + ( + eSRRF3D_ST(verbose=False), + ( + img[ + frames_per_timepoint + * i : frames_per_timepoint + * (i + 1) + ], + ), + { + "magnification_xy": magnification_xy, + "magnification_z": magnification_z, + "radius": radius, + "PSF_ratio": PSF_ratio, + "voxel_ratio": voxel_ratio, + "sensitivity": sensitivity, + "mode": mode, + "doIntensityWeighting": doIntensityWeighting, + }, + ) ) - else: - return _eSRRF3D.calculate(_force_run_type=_force_run_type)[0] + if macro_pixel_correction: + output_array[i] = macro_pixel_corrector( + _eSRRF3D.calculate(_force_run_type=_force_run_type)[0], + magnification=magnification_xy, + ) + else: + output_array[i] = _eSRRF3D.calculate( + _force_run_type=_force_run_type + )[0] + + return output_array.astype(np.float32) diff --git a/tests/test_esrrf3d.py b/tests/test_esrrf3d.py index 69d46546..05ec48ae 100644 --- a/tests/test_esrrf3d.py +++ b/tests/test_esrrf3d.py @@ -10,7 +10,7 @@ def test_esrrf3d_workflow(downloader): ) small_dataset = dataset[:10, :20, :20] - dataset = dataset[np.newaxis, ...] + small_dataset = small_dataset[np.newaxis, ...] eSRRF3D_w(small_dataset) @@ -96,4 +96,4 @@ def test_esrrf3d_radiusz_too_high(downloader): e3d = eSRRF3D_ST(testing=True, clear_benchmarks=True) - assert e3d.run(small_dataset, radius_z=15) is None + assert e3d.run(small_dataset, PSF_ratio=15) is None