Merge pull request #81 from DiamondLightSource/paganin_changes

Paganin filter related changes

Author: dkazanc

httomo_backends/methods_database/packages/backends/httomolibgpu/httomolibgpu.yaml

@@ -78,7 +78,7 @@ prep:
         multiplier: None
         method: module
   phase:
-    paganin_filter_tomopy:
+    paganin_filter:
       pattern: projection
       output_dims_change: False
       implementation: gpu_cupy

httomo_backends/methods_database/packages/backends/httomolibgpu/supporting_funcs/prep/phase.py

@@ -20,18 +20,17 @@
 # ---------------------------------------------------------------------------
 """Modules for memory estimation for phase retrieval and phase-contrast enhancement"""
 
-import math
 from typing import Tuple
 import numpy as np
 
 from httomo_backends.cufft import CufftType, cufft_estimate_2d
 
 __all__ = [
-    "_calc_memory_bytes_paganin_filter_tomopy",
+    "_calc_memory_bytes_paganin_filter",
 ]
 
 
-def _calc_memory_bytes_paganin_filter_tomopy(
+def _calc_memory_bytes_paganin_filter(
     non_slice_dims_shape: Tuple[int, int],
     dtype: np.dtype,
     **kwargs,

httomo_backends/methods_database/packages/backends/httomolibgpu/supporting_funcs/recon/algorithm.py

@@ -189,7 +189,7 @@ def _calc_memory_bytes_LPRec3d_tomobar(
     if detector_pad is True:
         detector_pad = __estimate_detectorHoriz_padding(non_slice_dims_shape[1])
     elif detector_pad is False:
-        detector_pad = 0        
+        detector_pad = 0
 
     min_mem_usage_filter = False
     if "min_mem_usage_filter" in kwargs:
@@ -208,9 +208,9 @@ def _calc_memory_bytes_LPRec3d_tomobar(
     _CENTER_SIZE_MIN = 192  # must be divisible by 8
 
     n = DetectorsLengthH
-    if(power_of_2_cropping):
+    if power_of_2_cropping:
         n_pow2 = 2 ** math.ceil(math.log2(n))
-        if( 0.9 < n / n_pow2 ):
+        if 0.9 < n / n_pow2:
             n = n_pow2
 
     odd_horiz = False
@@ -309,10 +309,7 @@ def _calc_memory_bytes_LPRec3d_tomobar(
         DetectorsLengthH_prepad * DetectorsLengthH_prepad * np.float32().itemsize
     )
     ifft2_plan_slice_size = (
-        cufft_estimate_2d(
-            nx=2 * n, ny=2 * n, fft_type=CufftType.CUFFT_C2C
-        )
-        / 2
+        cufft_estimate_2d(nx=2 * n, ny=2 * n, fft_type=CufftType.CUFFT_C2C) / 2
     )
     circular_mask_size = np.prod(output_dims) / 2 * np.int64().itemsize * 4
     after_recon_swapaxis_slice = recon_output_size
@@ -351,38 +348,74 @@ def add_to_memory_counters(amount, per_slice: bool):
     add_to_memory_counters(tmp_p_input_slice, True)
     if min_mem_usage_filter:
         add_to_memory_counters(rfft_plan_slice_size / 4 / projection_chunk_count, False)
-        add_to_memory_counters(irfft_plan_slice_size / 4 / projection_chunk_count, False)
+        add_to_memory_counters(
+            irfft_plan_slice_size / 4 / projection_chunk_count, False
+        )
         add_to_memory_counters(padded_tmp_p_input_slice / projection_chunk_count, False)
 
         add_to_memory_counters(rfft_result_size / projection_chunk_count, False)
-        add_to_memory_counters(filtered_rfft_result_size / projection_chunk_count, False)
+        add_to_memory_counters(
+            filtered_rfft_result_size / projection_chunk_count, False
+        )
         add_to_memory_counters(-rfft_result_size / projection_chunk_count, False)
-        add_to_memory_counters(-padded_tmp_p_input_slice / projection_chunk_count, False)
+        add_to_memory_counters(
+            -padded_tmp_p_input_slice / projection_chunk_count, False
+        )
 
-        add_to_memory_counters(irfft_scratch_memory_size / projection_chunk_count, False)
-        add_to_memory_counters(-irfft_scratch_memory_size / projection_chunk_count, False)
+        add_to_memory_counters(
+            irfft_scratch_memory_size / projection_chunk_count, False
+        )
+        add_to_memory_counters(
+            -irfft_scratch_memory_size / projection_chunk_count, False
+        )
         add_to_memory_counters(irfft_result_size / projection_chunk_count, False)
-        add_to_memory_counters(-filtered_rfft_result_size / projection_chunk_count, False)
+        add_to_memory_counters(
+            -filtered_rfft_result_size / projection_chunk_count, False
+        )
 
         add_to_memory_counters(-irfft_result_size / projection_chunk_count, False)
     else:
-        add_to_memory_counters(rfft_plan_slice_size / chunk_count / projection_chunk_count * 2, True)
-        add_to_memory_counters(irfft_plan_slice_size / chunk_count / projection_chunk_count * 2, True)
+        add_to_memory_counters(
+            rfft_plan_slice_size / chunk_count / projection_chunk_count * 2, True
+        )
+        add_to_memory_counters(
+            irfft_plan_slice_size / chunk_count / projection_chunk_count * 2, True
+        )
         # add_to_memory_counters(irfft_scratch_memory_size / chunk_count / projection_chunk_count, True)
         for _ in range(0, chunk_count):
-            add_to_memory_counters(padded_tmp_p_input_slice / chunk_count / projection_chunk_count, True)
-
-            add_to_memory_counters(rfft_result_size / chunk_count / projection_chunk_count, True)
-            add_to_memory_counters(filtered_rfft_result_size / chunk_count / projection_chunk_count, True)
-            add_to_memory_counters(-rfft_result_size / chunk_count / projection_chunk_count, True)
-            add_to_memory_counters(-padded_tmp_p_input_slice / chunk_count / projection_chunk_count, True)
-
-            add_to_memory_counters(irfft_scratch_memory_size / chunk_count / projection_chunk_count, True)
-            add_to_memory_counters(-irfft_scratch_memory_size / chunk_count / projection_chunk_count, True)
-            add_to_memory_counters(irfft_result_size / chunk_count / projection_chunk_count, True)
-            add_to_memory_counters(-filtered_rfft_result_size / chunk_count / projection_chunk_count, True)
-
-            add_to_memory_counters(-irfft_result_size / chunk_count / projection_chunk_count, True)
+            add_to_memory_counters(
+                padded_tmp_p_input_slice / chunk_count / projection_chunk_count, True
+            )
+
+            add_to_memory_counters(
+                rfft_result_size / chunk_count / projection_chunk_count, True
+            )
+            add_to_memory_counters(
+                filtered_rfft_result_size / chunk_count / projection_chunk_count, True
+            )
+            add_to_memory_counters(
+                -rfft_result_size / chunk_count / projection_chunk_count, True
+            )
+            add_to_memory_counters(
+                -padded_tmp_p_input_slice / chunk_count / projection_chunk_count, True
+            )
+
+            add_to_memory_counters(
+                irfft_scratch_memory_size / chunk_count / projection_chunk_count, True
+            )
+            add_to_memory_counters(
+                -irfft_scratch_memory_size / chunk_count / projection_chunk_count, True
+            )
+            add_to_memory_counters(
+                irfft_result_size / chunk_count / projection_chunk_count, True
+            )
+            add_to_memory_counters(
+                -filtered_rfft_result_size / chunk_count / projection_chunk_count, True
+            )
+
+            add_to_memory_counters(
+                -irfft_result_size / chunk_count / projection_chunk_count, True
+            )
 
     add_to_memory_counters(-padded_in_slice_size, True)
     add_to_memory_counters(-filter_size, False)
@@ -433,7 +466,7 @@ def _calc_memory_bytes_SIRT3d_tomobar(
     if detector_pad is True:
         detector_pad = __estimate_detectorHoriz_padding(non_slice_dims_shape[1])
     elif detector_pad is False:
-        detector_pad = 0        
+        detector_pad = 0
 
     anglesnum = non_slice_dims_shape[0]
     DetectorsLengthH_padded = non_slice_dims_shape[1] + 2 * detector_pad
@@ -456,7 +489,7 @@ def _calc_memory_bytes_SIRT3d_tomobar(
     Res_times_R = Res
     C_times_res = out_data_size
 
-    astra_projection = (in_data_size + out_data_size)
+    astra_projection = in_data_size + out_data_size
 
     tot_memory_bytes = int(
         recon_data_size_original
@@ -483,7 +516,7 @@ def _calc_memory_bytes_CGLS3d_tomobar(
     if detector_pad is True:
         detector_pad = __estimate_detectorHoriz_padding(non_slice_dims_shape[1])
     elif detector_pad is False:
-        detector_pad = 0        
+        detector_pad = 0
 
     anglesnum = non_slice_dims_shape[0]
     DetectorsLengthH_padded = non_slice_dims_shape[1] + 2 * detector_pad
@@ -535,7 +568,7 @@ def _calc_memory_bytes_FISTA3d_tomobar(
     if detector_pad is True:
         detector_pad = __estimate_detectorHoriz_padding(non_slice_dims_shape[1])
     elif detector_pad is False:
-        detector_pad = 0        
+        detector_pad = 0
 
     anglesnum = non_slice_dims_shape[0]
     DetectorsLengthH_padded = non_slice_dims_shape[1] + 2 * detector_pad

httomo_backends/pipelines_full/deg360_paganin_FBP3d_tomobar_directive.yaml

@@ -8,7 +8,7 @@
   module_path: httomolibgpu.misc.morph
 - method: remove_stripe_based_sorting
   module_path: httomolibgpu.prep.stripe
-- method: paganin_filter_tomopy
+- method: paganin_filter
   module_path: httomolibgpu.prep.phase
 - method: FBP3d_tomobar
   module_path: httomolibgpu.recon.algorithm

httomo_backends/pipelines_full/sweep_paganin_FBP3d_tomobar_directive.yaml

@@ -4,9 +4,9 @@
   module_path: httomolibgpu.recon.rotation  
 - method: normalize
   module_path: httomolibgpu.prep.normalize
-- method: paganin_filter_tomopy
+- method: paganin_filter
   module_path: httomolibgpu.prep.phase
-  sweep_parameter: alpha
-  sweep_values: [0.01, 0.001, 0.0001]
+  sweep_parameter: ratio_delta_beta
+  sweep_values: [10, 150, 350]
 - method: FBP3d_tomobar
-  module_path: httomolibgpu.recon.algorithm
+  module_path: httomolibgpu.recon.algorithm

httomo_backends/scripts/yaml_pipelines_generator.py

@@ -138,7 +138,7 @@ def yaml_pipelines_generator(
 
             pipeline_full.yaml_set_start_comment(
                 "This pipeline should be supported by the latest developments of HTTomo. Use module load httomo/latest module at Diamond."
-            )            
+            )
 
             if "loaders" in module_name:
                 # should be the first method in the list
@@ -246,8 +246,20 @@ def yaml_pipelines_generator(
                 )
                 pipeline_full += yaml_template_method
                 pipeline_full[i]["parameters"].yaml_add_eol_comment(
-                    key="alpha",
-                    comment="Controls the balance between the strength of the filter and the amount of noise reduction. Smaller values lead to less noise and more blur.",
+                    key="ratio_delta_beta",
+                    comment="The ratio of delta/beta for filter strength control. Larger values lead to more smoothing.",
+                )
+                pipeline_full[i]["parameters"].yaml_add_eol_comment(
+                    key="pixel_size",
+                    comment="Detector pixel size (resolution) in MICRON units.",
+                )
+                pipeline_full[i]["parameters"].yaml_add_eol_comment(
+                    key="distance",
+                    comment="Propagation distance of the wavefront from sample to detector in METRE units.",
+                )
+                pipeline_full[i]["parameters"].yaml_add_eol_comment(
+                    key="energy",
+                    comment="Beam energy in keV.",
                 )
             elif "stripe" in module_name:
                 pipeline_full.yaml_set_comment_before_after_key(
@@ -273,7 +285,7 @@ def yaml_pipelines_generator(
                 )
                 pipeline_full[i]["parameters"].yaml_add_eol_comment(
                     key="recon_mask_radius",
-                    comment="Zero pixels outside the mask-circle radius.",
+                    comment="Zero pixels outside the mask-circle radius. Make radius equal to 2.0 to remove the mask effect.",
                 )
                 pipeline_full[i]["parameters"].yaml_add_eol_comment(
                     key="neglog",

tests/test_httomolibgpu.py

@@ -15,7 +15,7 @@
 
 from httomolibgpu.misc.morph import data_resampler, sino_360_to_180
 from httomolibgpu.prep.normalize import normalize
-from httomolibgpu.prep.phase import paganin_filter_tomopy
+from httomolibgpu.prep.phase import paganin_filter
 from httomolibgpu.prep.alignment import distortion_correction_proj_discorpy
 from httomolibgpu.prep.stripe import (
     remove_stripe_based_sorting,
@@ -228,19 +228,19 @@ def test_denoiser_PD_TV_memoryhook(ensure_clean_memory, slices):
 @pytest.mark.parametrize("slices", [64, 128])
 @pytest.mark.parametrize("dim_x", [81, 260, 320])
 @pytest.mark.parametrize("dim_y", [340, 135, 96])
-def test_paganin_filter_tomopy_memoryhook(slices, dim_x, dim_y, ensure_clean_memory):
+def test_paganin_filter_memoryhook(slices, dim_x, dim_y, ensure_clean_memory):
     data = cp.random.random_sample((slices, dim_x, dim_y), dtype=np.float32)
     hook = MaxMemoryHook()
     with hook:
-        data_filtered = paganin_filter_tomopy(cp.copy(data)).get()
+        data_filtered = paganin_filter(cp.copy(data)).get()
 
     # make sure estimator function is within range (80% min, 100% max)
     max_mem = (
         hook.max_mem
     )  # the amount of memory in bytes needed for the method according to memoryhook
 
     # now we estimate how much of the total memory required for this data
-    (estimated_memory_bytes, subtract_bytes) = _calc_memory_bytes_paganin_filter_tomopy(
+    (estimated_memory_bytes, subtract_bytes) = _calc_memory_bytes_paganin_filter(
         (dim_x, dim_y), dtype=np.float32()
     )
     estimated_memory_mb = round(slices * estimated_memory_bytes / (1024**2), 2)
@@ -574,7 +574,9 @@ def test_recon_LPRec3d_tomobar_0_pi_memoryhook(
 
 @pytest.mark.full
 @pytest.mark.cupy
-@pytest.mark.parametrize("min_mem_usage_filter_ifft2", [(False, False), (True, False), (True, True)])
+@pytest.mark.parametrize(
+    "min_mem_usage_filter_ifft2", [(False, False), (True, False), (True, True)]
+)
 @pytest.mark.parametrize("power_of_2_cropping", [False, True])
 @pytest.mark.parametrize("padding_detx", [0, 10, 50, 100, 800])
 @pytest.mark.parametrize("projections", [1500, 1801, 2560, 3601])
@@ -605,7 +607,9 @@ def test_recon_LPRec3d_tomobar_0_pi_memoryhook_full(
 
 @pytest.mark.full
 @pytest.mark.cupy
-@pytest.mark.parametrize("min_mem_usage_filter_ifft2", [(False, False), (True, False), (True, True)])
+@pytest.mark.parametrize(
+    "min_mem_usage_filter_ifft2", [(False, False), (True, False), (True, True)]
+)
 @pytest.mark.parametrize("power_of_2_cropping", [False, True])
 @pytest.mark.parametrize("padding_detx", [0, 10, 50, 100, 800])
 @pytest.mark.parametrize("projections", [1500, 1801, 2560, 3601])
@@ -742,10 +746,13 @@ def test_recon_SIRT3d_tomobar_memoryhook(slices, recon_size_it, ensure_clean_mem
     assert estimated_memory_mb >= max_mem_mb
     assert percents_relative_maxmem <= 100
 
+
 @pytest.mark.cupy
 @pytest.mark.parametrize("slices", [3])
 @pytest.mark.parametrize("recon_size_it", [2560])
-def test_recon_SIRT3d_tomobar_autopad_memoryhook(slices, recon_size_it, ensure_clean_memory):
+def test_recon_SIRT3d_tomobar_autopad_memoryhook(
+    slices, recon_size_it, ensure_clean_memory
+):
     angles_tot = 901
     det_size = 2560
     data = cp.random.random_sample((angles_tot, slices, det_size), dtype=np.float32)
@@ -758,8 +765,8 @@ def test_recon_SIRT3d_tomobar_autopad_memoryhook(slices, recon_size_it, ensure_c
         recon_data = SIRT3d_tomobar(
             cp.copy(data),
             np.linspace(0.0 * np.pi / 180.0, 180.0 * np.pi / 180.0, data.shape[0]),
-            center = 1200,
-            detector_pad = True,
+            center=1200,
+            detector_pad=True,
             recon_size=recon_size_it,
             iterations=2,
             nonnegativity=True,
@@ -803,8 +810,8 @@ def test_recon_CGLS3d_tomobar_memoryhook(slices, recon_size_it, ensure_clean_mem
         recon_data = CGLS3d_tomobar(
             cp.copy(data),
             np.linspace(0.0 * np.pi / 180.0, 180.0 * np.pi / 180.0, data.shape[0]),
-            center = 1200,
-            detector_pad = False,
+            center=1200,
+            detector_pad=False,
             recon_size=recon_size_it,
             iterations=2,
             nonnegativity=True,
@@ -835,7 +842,9 @@ def test_recon_CGLS3d_tomobar_memoryhook(slices, recon_size_it, ensure_clean_mem
 @pytest.mark.cupy
 @pytest.mark.parametrize("slices", [3])
 @pytest.mark.parametrize("recon_size_it", [2560])
-def test_recon_CGLS3d_tomobar_autopad_memoryhook(slices, recon_size_it, ensure_clean_memory):
+def test_recon_CGLS3d_tomobar_autopad_memoryhook(
+    slices, recon_size_it, ensure_clean_memory
+):
     angles_tot = 901
     det_size = 2560
     data = cp.random.random_sample((angles_tot, slices, det_size), dtype=np.float32)
@@ -848,8 +857,8 @@ def test_recon_CGLS3d_tomobar_autopad_memoryhook(slices, recon_size_it, ensure_c
         recon_data = CGLS3d_tomobar(
             cp.copy(data),
             np.linspace(0.0 * np.pi / 180.0, 180.0 * np.pi / 180.0, data.shape[0]),
-            center = 1200,
-            detector_pad = True,
+            center=1200,
+            detector_pad=True,
             recon_size=recon_size_it,
             iterations=2,
             nonnegativity=True,
@@ -876,6 +885,7 @@ def test_recon_CGLS3d_tomobar_autopad_memoryhook(slices, recon_size_it, ensure_c
     assert estimated_memory_mb >= max_mem_mb
     assert percents_relative_maxmem <= 85
 
+
 @pytest.mark.cupy
 @pytest.mark.parametrize("slices", [3, 5])
 @pytest.mark.parametrize("recon_size_it", [2560])