diff --git a/docs/source/_static/figures/center_find/find_center_entropy_1253_smoothing3.jpg b/docs/source/_static/figures/center_find/find_center_entropy_1253_smoothing3.jpg new file mode 100755 index 00000000..f3505e32 Binary files /dev/null and b/docs/source/_static/figures/center_find/find_center_entropy_1253_smoothing3.jpg differ diff --git a/docs/source/_static/figures/center_find/find_center_entropy_1263_no_smoothing.jpg b/docs/source/_static/figures/center_find/find_center_entropy_1263_no_smoothing.jpg new file mode 100755 index 00000000..caaca4f5 Binary files /dev/null and b/docs/source/_static/figures/center_find/find_center_entropy_1263_no_smoothing.jpg differ diff --git a/docs/source/reference/api.rst b/docs/source/reference/api.rst index 5f5d29f0..8ea94152 100644 --- a/docs/source/reference/api.rst +++ b/docs/source/reference/api.rst @@ -4,7 +4,7 @@ API reference ========================== -This section contains the API reference and usage information for HttomolibGPU. +This section contains the API reference and usage information for HTTomolibGPU. HTTomolibGPU Modules --------------------- diff --git a/docs/source/reference/methods.rst b/docs/source/reference/methods.rst index ec31ad7d..b4008f86 100644 --- a/docs/source/reference/methods.rst +++ b/docs/source/reference/methods.rst @@ -11,6 +11,7 @@ Here we present a list of methods of the HTTomolibGPU library with more detailed methods_list/correction_methods methods_list/stripe_removal_methods methods_list/phase_contrast_methods + methods_list/cor_finder_methods methods_list/reconstruction_methods methods_list/denoising_methods methods_list/rescale_methods diff --git a/docs/source/reference/methods_list/cor_finder/find_center_metric_recon.rst b/docs/source/reference/methods_list/cor_finder/find_center_metric_recon.rst new file mode 100644 index 00000000..972921ed --- /dev/null +++ b/docs/source/reference/methods_list/cor_finder/find_center_metric_recon.rst @@ -0,0 +1,77 @@ +.. _find_center_metric_recon_doc: + +Find center by metric +^^^^^^^^^^^^^^^^^^^^^ + +**Description** + +The optimal centre of rotation is found by minimising a metric of the reconstructed image: + +.. math:: + + c^{*} = \arg\min_{c\in\mathcal{C}} M(c), + +where :math:`\mathcal{C}` is the set of candidate centres and :math:`M` is an image-quality metric. + +The choice of metric is important because different metrics respond differently to noise, background, contrast, and reconstruction artefacts. It is therefore useful to compare several metrics rather than relying on a single measure. + + +**Where and how to use it:** + +Use it when the centre of rotation is unknown and needs to be found automatically. + +This method can be applied to any data, however the main benefit of this method +is when the data is incomplete, i.e., limited angle/missing wedge data. + +.. note:: For parallel-beam tomography, projections separated by approximately 180 degrees contain complementary information. With a correctly determined COR, the projections should exhibit the expected geometrical relationship. Other methods, such as, :mod:`httomolibgpu.recon.rotation.find_center_vo` and :mod:`httomolibgpu.recon.rotation.find_center_pc`, rely on that assumption. So it is advisable to use other methods first, if the data is complete. + +**What are the adjustable parameters:** + +There are quite a lot of parameters that can be adjusted to improve the performance of this method. +The most important ones are: + +* :code:`metric_type` the choice of metrics to be computed on the reconstructed image. Choose from: + + * :code:`'entropy'` - Shannon entropy measures the information content or complexity of an image intensity distribution. For a discrete intensity distribution with probabilities :math:`p_i`, the entropy is :math:`H = -\sum_i p_i \log_2(p_i)`. A higher entropy indicates a broader or more complex intensity distribution. For an image, the intensity distribution can be estimated using a histogram This metric is sensitive to noise and background, but less sensitive to contrast and reconstruction artefacts. Normally a second best choice after :code:`'tv'` metric. + + * :code:`'tv'` - the total variation (the sum of the magnitude of the image gradient) of the reconstructed image. :math:`|\nabla I| = \sqrt{\left(\frac{\partial I}{\partial x}\right)^2 + \left(\frac{\partial I}{\partial y}\right)^2}`. The :code:`'tv'` metric is the sum of gradient's magnitude: :math:`M_\mathrm{grad} = \sum{(|\nabla I|)}`. Default choice and of the most robust metrics. + + * :code:`'sharpness'` - the variance of of the Laplacian. The metric is defined as :math:`M_\mathrm{Lap} = \operatorname{Var}(\nabla^2 I)`, where :math:`I` is the reconstructed image. Very sensitive to high-frequency noise. + + +.. note:: We found that smoothing the image before computing any of the metrics with :code:`'gaussian_filter_sigma'` can improve the performance of the correct CoR estimation significantly. Usually the range between 1 and 3 is a good choice for the :code:`'gaussian_filter_sigma'` parameter (see images below). + +* :code:`reconstruction_method` is set to :code:`'LPRec3d_tomobar'` by default, however the iterative methods such as :code:`'SIRT3d_tomobar'`, :code:`'CGLS3d_tomobar'` are also available. In some cases (especially limited angle data) the iterative methods can provide better results than the default direct methods. Set :code:`'gaussian_filter_sigma' = None'` if iterative methods are enabled. + +* :code:`cor_initialisation_value` is the initial guess of the centre of rotation. The search for the optimal CoR will be performed in the range :code:`[cor_initialisation_value - cor_search_range, cor_initialisation_value + cor_search_range]`. The default value is set to :code:`None`, which means that the initialisation value will be set to the middle of the detector. + +* :code:`squared_mask_x_y_size` is a very important parameter as it defines the squared mask that is applied to the reconstructed image BEFORE computing the metric. The mask can be placed in any position on the reconstructed image. Very frequently, the optimal location of the mask leads to optimal performance of the estimation. However, placing the mask on the ring/streak artifacts or at the edge of the reconstruction leads to incorrect results. + + The mask is defined as a :code:`list` with 3 values :code:`[X, Y, size]`. The positive offsets :code:`X, Y` place the mask with respect to the left top corner of the reconstructed image and the third parameter is the size of the cropped image in percents with respect to the whole reconstructed image. + + Example: :code:`squared_mask_x_y_size = [10, 20, 50]` will apply a mask that is 10 pixels away from the left top corner in the horizontal direction and 20 pixels away in the vertical direction, the mask will be 50% in size of the reconstructed image size. + + The default value when :code:`squared_mask_x_y_size = None` is the mask in the center of the reconstructed image (not always the best position as many ring artifacts are present) with a size of 50% of the reconstructed image size. + +* :code:`save_recon_tiff` - this feature is useful for debugging. When set to a valid path, the reconstructed image together with the mask after smoothing are saved in to a specified folder. This allows for visual inspection of the reconstruction and mask placement. + +**Practical example:** + +In this example we demonstrate how the choice of smoothing can affect the performance of the CoR estimation. + +.. list-table:: + + + * - .. figure:: ../../../_static/figures/center_find/find_center_entropy_1263_no_smoothing.jpg + + :code:`metric_type = entropy`, :code:`squared_mask_x_y_size = None` (center), :code:`gaussian_filter_sigma = None`. Incorrect CoR found as the masked image is not smoothed. + + - .. figure:: ../../../_static/figures/center_find/find_center_entropy_1253_smoothing3.jpg + + :code:`metric_type = entropy`, :code:`squared_mask_x_y_size = None` (center), :code:`gaussian_filter_sigma = 3`. Correct CoR found as the masked image is smoothed. + + + + + + diff --git a/docs/source/reference/methods_list/cor_finder_methods.rst b/docs/source/reference/methods_list/cor_finder_methods.rst new file mode 100644 index 00000000..45af416c --- /dev/null +++ b/docs/source/reference/methods_list/cor_finder_methods.rst @@ -0,0 +1,14 @@ +.. _cor_finder_module: + +Centre of Rotation +****************** + +Methods from :mod:`httomolibgpu.recon.rotation` module are used for automatic determination of the centre of rotation (COR). The centre of rotation (COR) is a critical parameter in tomographic reconstruction. An incorrect COR introduces characteristic reconstruction artefacts, including blurring, double edges, and duplicated structures. See more about CoR `here `_ + + +.. toctree:: + :maxdepth: 2 + + cor_finder/find_center_metric_recon + + diff --git a/httomolibgpu/__init__.py b/httomolibgpu/__init__.py index 36aa0a91..202ef190 100644 --- a/httomolibgpu/__init__.py +++ b/httomolibgpu/__init__.py @@ -26,4 +26,9 @@ ADMM3d_tomobar, ) -from httomolibgpu.recon.rotation import find_center_vo, find_center_360, find_center_pc +from httomolibgpu.recon.rotation import ( + find_center_vo, + find_center_360, + find_center_pc, + find_center_metric_recon, +) diff --git a/httomolibgpu/misc/blend.py b/httomolibgpu/misc/blend.py index 4ab84d1f..bdfbf471 100644 --- a/httomolibgpu/misc/blend.py +++ b/httomolibgpu/misc/blend.py @@ -18,7 +18,11 @@ # Created By : Tomography Team at DLS # Created Date: 5 August 2026 # --------------------------------------------------------------------------- -"""Module for data type morphing functions""" +"""Module for data type morphing functions: + +* :mod:`httomolibgpu.misc.blend.seam_blend_stitched_data` + +""" import numpy as np from httomolibgpu import cupywrapper diff --git a/httomolibgpu/misc/corr.py b/httomolibgpu/misc/corr.py index f303e79c..85a84d23 100644 --- a/httomolibgpu/misc/corr.py +++ b/httomolibgpu/misc/corr.py @@ -18,7 +18,13 @@ # Created By : Tomography Team at DLS # Created Date: 21/October/2022 # --------------------------------------------------------------------------- -"""Module for data correction. For more detailed information see :ref:`data_correction_module`.""" +"""Module for data correction. For more detailed information see :ref:`data_correction_module`. + +* :mod:`httomolibgpu.misc.corr.median_filter` + +* :mod:`httomolibgpu.misc.corr.remove_outlier` + +""" from httomolibgpu import cupywrapper diff --git a/httomolibgpu/misc/denoise.py b/httomolibgpu/misc/denoise.py index 8e289367..dd306893 100644 --- a/httomolibgpu/misc/denoise.py +++ b/httomolibgpu/misc/denoise.py @@ -18,7 +18,13 @@ # Created By : Tomography Team at DLS # Created Date: 18/December/2024 # --------------------------------------------------------------------------- -"""Module for data denoising. For more detailed information see :ref:`data_denoising_module`.""" +"""Module for data denoising. For more detailed information see :ref:`data_denoising_module`: + +* :mod:`httomolibgpu.misc.denoise.total_variation_ROF` + +* :mod:`httomolibgpu.misc.denoise.total_variation_PD` + +""" from httomolibgpu import cupywrapper diff --git a/httomolibgpu/misc/morph.py b/httomolibgpu/misc/morph.py index e58d32d6..f9ac00cd 100644 --- a/httomolibgpu/misc/morph.py +++ b/httomolibgpu/misc/morph.py @@ -18,7 +18,15 @@ # Created By : Tomography Team at DLS # Created Date: 23 March 2023 # --------------------------------------------------------------------------- -"""Module for data type morphing functions""" +"""Module for data type morphing functions: + +* :mod:`httomolibgpu.misc.morph.sino_360_to_180` + +* :mod:`httomolibgpu.misc.morph.data_resampler` + +* :mod:`httomolibgpu.misc.morph.average_projection_frames` + +""" import numpy as np from httomolibgpu import cupywrapper diff --git a/httomolibgpu/misc/rescale.py b/httomolibgpu/misc/rescale.py index 845e31bb..2602d381 100644 --- a/httomolibgpu/misc/rescale.py +++ b/httomolibgpu/misc/rescale.py @@ -18,7 +18,11 @@ # Created By : Tomography Team at DLS # Created Date: 1 March 2024 # --------------------------------------------------------------------------- -"""Module for data rescaling. For more detailed information see :ref:`data_rescale_module`.""" +"""Module for data rescaling. For more detailed information see :ref:`data_rescale_module` : + +* :mod:`httomolibgpu.misc.rescale.rescale_to_int` + +""" import numpy as np from httomolibgpu import cupywrapper diff --git a/httomolibgpu/prep/alignment.py b/httomolibgpu/prep/alignment.py index 8f282bba..6d825cb2 100644 --- a/httomolibgpu/prep/alignment.py +++ b/httomolibgpu/prep/alignment.py @@ -18,7 +18,11 @@ # Created By : Tomography Team at DLS # Created Date: 01 November 2022 # --------------------------------------------------------------------------- -"""Modules for data correction""" +"""Modules for data correction: + +* :mod:`httomolibgpu.prep.alignment.distortion_correction_proj_discorpy` + +""" from httomolibgpu import cupywrapper diff --git a/httomolibgpu/prep/normalize.py b/httomolibgpu/prep/normalize.py index a6224946..7c7df152 100644 --- a/httomolibgpu/prep/normalize.py +++ b/httomolibgpu/prep/normalize.py @@ -18,7 +18,13 @@ # Created By : Tomography Team at DLS # Created Date: 01 November 2022 # --------------------------------------------------------------------------- -"""Modules for raw projection data normalization""" +"""Modules for flat/dark-field correction of projection data and normalisation: + +* :mod:`httomolibgpu.prep.normalize.dark_flat_field_correction` + +* :mod:`httomolibgpu.prep.normalize.minus_log` + +""" from httomolibgpu import cupywrapper import numpy as np diff --git a/httomolibgpu/prep/phase.py b/httomolibgpu/prep/phase.py index 087aea9d..962cdc99 100644 --- a/httomolibgpu/prep/phase.py +++ b/httomolibgpu/prep/phase.py @@ -18,7 +18,13 @@ # Created By : Tomography Team at DLS # Created Date: 01 November 2022 # --------------------------------------------------------------------------- -"""Modules for phase retrieval and phase-contrast enhancement. For more detailed information, see :ref:`phase_contrast_module`.""" +"""Modules for phase retrieval and phase-contrast enhancement. For more detailed information, see :ref:`phase_contrast_module`. + +* :mod:`httomolibgpu.prep.phase.paganin_filter` + +* :mod:`httomolibgpu.prep.phase.paganin_filter_savu_legacy` + +""" import numpy as np from httomolibgpu import cupywrapper diff --git a/httomolibgpu/prep/stripe.py b/httomolibgpu/prep/stripe.py index 719b4b97..2e250c6d 100644 --- a/httomolibgpu/prep/stripe.py +++ b/httomolibgpu/prep/stripe.py @@ -18,7 +18,19 @@ # Created By : Tomography Team at DLS # Created Date: 01 November 2022 # --------------------------------------------------------------------------- -"""Module for stripes removal""" +"""Module for stripes removal (ring artifacts suppression), see more in :ref:`stripes_removal_module`. + +* :mod:`httomolibgpu.prep.stripe.remove_stripe_based_sorting` + +* :mod:`httomolibgpu.prep.stripe.remove_stripe_fw` + +* :mod:`httomolibgpu.prep.stripe.remove_stripe_ti` + +* :mod:`httomolibgpu.prep.stripe.remove_all_stripe` + +* :mod:`httomolibgpu.prep.stripe.raven_filter` + +""" import numpy as np import pywt diff --git a/httomolibgpu/recon/algorithm.py b/httomolibgpu/recon/algorithm.py index 48959e85..70181219 100644 --- a/httomolibgpu/recon/algorithm.py +++ b/httomolibgpu/recon/algorithm.py @@ -18,7 +18,25 @@ # Created By : Tomography Team at DLS # Changes relative to ToMoBAR 2024.01 version # --------------------------------------------------------------------------- -"""Module for tomographic reconstruction. For more detailed information, see :ref:`image_reconstruction_module`""" +"""Module for tomographic reconstruction. For more detailed information, see :ref:`image_reconstruction_module` + +* :mod:`httomolibgpu.recon.algorithm.FBP2d_astra` + +* :mod:`httomolibgpu.recon.algorithm.FBP3d_tomobar` + +* :mod:`httomolibgpu.recon.algorithm.LPRec3d_tomobar` + +* :mod:`httomolibgpu.recon.algorithm.SIRT3d_tomobar` + +* :mod:`httomolibgpu.recon.algorithm.CGLS3d_tomobar` + +* :mod:`httomolibgpu.recon.algorithm.FISTA3d_tomobar` + +* :mod:`httomolibgpu.recon.algorithm.ADMM3d_tomobar` + +* :mod:`httomolibgpu.recon.algorithm.OSEM3d_tomobar` + +""" import numpy as np from httomolibgpu import cupywrapper diff --git a/httomolibgpu/recon/rotation.py b/httomolibgpu/recon/rotation.py index 6cab12c0..806c724a 100644 --- a/httomolibgpu/recon/rotation.py +++ b/httomolibgpu/recon/rotation.py @@ -18,19 +18,40 @@ # Created By : Tomography Team at DLS # Created Date: 01 November 2022 # --------------------------------------------------------------------------- -"""Modules for finding the axis of rotation for 180 and 360 degrees scans""" +"""Modules for finding the axis of rotation for 180 or 360 degrees scans. + +* :mod:`httomolibgpu.recon.rotation.find_center_vo` + +* :mod:`httomolibgpu.recon.rotation.find_center_360` + +* :mod:`httomolibgpu.recon.rotation.find_center_pc` + +* :mod:`httomolibgpu.recon.rotation.find_center_metric_recon` +""" import numpy as np +import os +import pathlib +from pathlib import Path from numpy.polynomial import Polynomial from httomolibgpu import cupywrapper +from PIL import Image +from scipy.ndimage import laplace cp = cupywrapper.cp cupy_run = cupywrapper.cupy_run + from unittest.mock import Mock if cupy_run: from httomolibgpu.cuda_kernels import load_cuda_module + from httomolibgpu.recon.algorithm import ( + FBP3d_tomobar, + LPRec3d_tomobar, + SIRT3d_tomobar, + CGLS3d_tomobar, + ) from cupyx.scipy.ndimage import shift, gaussian_filter from ._phase_cross_correlation import phase_cross_correlation from cupyx.scipy.fftpack import get_fft_plan @@ -60,6 +81,7 @@ "find_center_vo", "find_center_360", "find_center_pc", + "find_center_metric_recon", ] @@ -863,3 +885,287 @@ def find_center_pc( ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +## %%%%%%%%%%%%%%%%%%%find_center_metric_recon%%%%%%%%%%%%%%%%%%%%% +def find_center_metric_recon( + data: cp.ndarray, + angles: np.ndarray, + metric_type: Literal["entropy", "tv", "sharpness"] = "tv", + range: int = 10, + step: Union[float, int] = 0.5, + reconstruction_method: Literal[ + "LPRec3d_tomobar", "FBP3d_tomobar", "SIRT3d_tomobar", "CGLS3d_tomobar" + ] = "LPRec3d_tomobar", + ind: Optional[int] = None, + recon_iterations: Optional[int] = None, + cor_initialisation_value: Optional[float] = None, + squared_mask_x_y_size: Optional[list] = None, + gaussian_filter_sigma: Optional[float] = None, + save_recon_tiff: Optional[os.PathLike] = None, +) -> np.float32: + """ + Find the rotation axis location using different metrics that are applied to the result of the reconstruction. This method iteratively + assesses the quality of the reconstruction while changing the CoR value within the provided range. This approach for the CoR finding is suitable + for limited angle (missing wedge) data, when the symmetry and consistency of the sinogram is lost. + See more about the method and its parameters in :ref:`find_center_metric_recon_doc`. + + + Parameters + ---------- + data : cp.ndarray + 3D [angles, 1, detX] tomographic data as a CuPy array. + angles : np.ndarray + An array of angles given in radians. + metric_type : str, + Type of image quality metric to use on the reconstructed image. Available metrics are :code:`'entropy'`, :code:`'tv'`, :code:`'sharpness'`. + range : int + CoR search range/radius. The search will be performed in the range: :code:`[-range:cor_initialisation_value:range]`. + step : float + Step for CoR value. + reconstruction_method : str, + Type of the reconstruction method to be used. Choose from: :code:`"LPRec3d_tomobar"`, :code:`"FBP3d_tomobar"`, :code:`"SIRT3d_tomobar"`, :code:`"CGLS3d_tomobar"`. Default :code:`'LPRec3d_tomobar'`. + ind : int, optional + Index of the slice to be used for estimate the CoR. If 'None' is given, the zero slice will be used. + recon_iterations: int, optional + Set only for iterative methods: :code:`'SIRT3d_tomobar'`, :code:`'CGLS3d_tomobar'`. + cor_initialisation_value : float, optional + The initial approximation for the centre of rotation. If the value is None, use the horizontal centre of the projection/sinogram image. + squared_mask_x_y_size: list, optional + Apply a square mask to the reconstructed image. Selection of the mask is crucial for this algorithm to work successfully. + The mask is defined as a :code:`list` with 3 values :code:`[X, Y, size]`. The positive offsets :code:`X, Y` place the mask with respect to the left top corner of the reconstructed image and the third parameter is the size of the cropped image in percents with respect to the whole reconstructed image. + Example: :code:`squared_mask_x_y_size = [10, 20, 50]` will apply a mask that is 10 pixels away from the left top corner in the horizontal direction and 20 pixels away in the vertical direction, the mask will be 50% in size of the reconstructed image size. + The default value when :code:`squared_mask_x_y_size = None` is the mask in the center of the reconstructed image (not always the best position as many ring artifacts are present) with a size of 50% of the reconstructed image size. + gaussian_filter_sigma: float, optional + Enable gaussian filtering of the reconstructed image, highly recommended for noisy data. Good range of values 1.0-4.0. + save_recon_tiff: path, optional + Path to output directory for the saved reconstruction image when :code:`CoR = cor_initialisation_value` and also the masked image. Useful for debugging. + + Returns + ------- + float32 + Rotation axis location with a subpixel precision. + """ + ### Data and parameters checks ### + if ind is None: + data = data[:, 0:1, :] + else: + data = data[:, ind : ind + 1, :] + methods_name = "find_center_metric_recon" + __check_if_data_correct_type( + data, accepted_type=["float32"], methods_name=methods_name + ) + __check_variable_type( + metric_type, [str], "metric_type", ["entropy", "tv", "sharpness"], methods_name + ) + __check_variable_type(range, [int], "range", [], methods_name) + __check_variable_type(step, [int, float], "step", [], methods_name) + __check_variable_type( + reconstruction_method, + [str], + "reconstruction_method", + ["LPRec3d_tomobar", "FBP3d_tomobar", "SIRT3d_tomobar", "CGLS3d_tomobar"], + methods_name, + ) + __check_variable_type( + recon_iterations, [int, type(None)], "recon_iterations", [], methods_name + ) + if recon_iterations is None and reconstruction_method in [ + "SIRT3d_tomobar", + "CGLS3d_tomobar", + ]: + recon_iterations_n = 15 + __check_variable_type( + gaussian_filter_sigma, + [float, type(None)], + "gaussian_filter_sigma", + [], + methods_name, + ) + __check_variable_type( + squared_mask_x_y_size, + [list, type(None)], + "squared_mask_x_y_size", + [], + methods_name, + ) + __check_variable_type( + save_recon_tiff, + [str, pathlib.PosixPath, type(None)], + "save_recon_tiff", + [], + methods_name, + ) + ################################################################### + detectorX_size = data.shape[2] + if cor_initialisation_value is None: + center = detectorX_size // 2 + else: + center = cor_initialisation_value + + if squared_mask_x_y_size is None: + squared_mask_x_y_size = [ + detectorX_size // 2, + detectorX_size // 2, + 50, + ] # default mask is the center of the reconstructed image with a size of 50% of the reconstructed image size + + centres = np.arange(center - range, center + range + step, step) + + values = [] + + for i, center in enumerate(centres): + # ---------- RECONSTRUCTION ---------------- + if reconstruction_method == "LPRec3d_tomobar": + recon = LPRec3d_tomobar( + data, angles, center=center, detector_pad=True, recon_mask_radius=2.0 + ) + elif reconstruction_method == "FBP3d_tomobar": + recon = FBP3d_tomobar( + data, angles, center=center, detector_pad=True, recon_mask_radius=2.0 + ) + elif reconstruction_method == "SIRT3d_tomobar": + recon = SIRT3d_tomobar( + data, + angles, + center=center, + detector_pad=True, + recon_mask_radius=2.0, + iterations=recon_iterations_n, + ) + elif reconstruction_method == "CGLS3d_tomobar": + recon = CGLS3d_tomobar( + data, + angles, + center=center, + detector_pad=True, + recon_mask_radius=2.0, + iterations=recon_iterations_n, + ) + recon = recon[:, 0, :] + # ----CROPPING with the position defined by the squared mask----- + x_start, x_end, y_start, y_end = square_mask_bounds( + recon, + centre_x=squared_mask_x_y_size[0], + centre_y=squared_mask_x_y_size[1], + size_pct=squared_mask_x_y_size[2], + ) + image = recon[y_start:y_end, x_start:x_end] + + cp.nan_to_num(image, copy=False, nan=0.0, posinf=0, neginf=0) + + # ---------- SMOOTHING ----------- + image = ( + gaussian_filter(image, sigma=gaussian_filter_sigma, mode="reflect") + if gaussian_filter_sigma is not None + else image + ) + + image = cp.asnumpy(image) + + # SAVING reconstruction and masked image as TIFF + if save_recon_tiff is not None: + save_dir = Path(save_recon_tiff) + save_dir.mkdir(parents=True, exist_ok=True) + + save_path_recon = save_dir / f"reconstruction_cor_{center:.3f}.tiff" + save_path_masked = save_dir / f"masked_recon_{center:.3f}.tiff" + + Image.fromarray(to_uint16(cp.asnumpy(recon))).save(save_path_recon) + Image.fromarray(to_uint16(image)).save(save_path_masked) + + # ------- METRIC CALCULATION -------- + + if metric_type == "entropy": + metric_val = image_entropy(image) + elif metric_type == "tv": + metric_val = gradient_energy(image) + elif metric_type == "sharpness": + metric_val = sharpness_metric(image) + + values.append(metric_val) + + print( + f"{i+1}/{len(centres)} " f"center={center:.3f}, " f"value={metric_val:.6f}" + ) + + values_arr = np.asarray(values) + + best_index = np.argmin(values_arr) + + best_center = centres[best_index] + value_function_min = values_arr[best_index] + + return best_center + ################################### + + +def image_entropy(image: np.ndarray) -> float: + """Compute the Shannon entropy of a reconstructed (cropped) 2D image.""" + + bins = ( + np.shape(image)[0] // 10 + ) # choosing the number of bins based on the image size + x = image.flatten() + # Histogram -> probability distribution + counts, _ = np.histogram(x, bins=bins) + p = counts[counts > 0] / counts.sum() + + # Shannon entropy + return float(-np.sum(p * np.log2(p))) + + +def gradient_energy(image: np.ndarray) -> float: + """Compute the gradient energy of a reconstructed (cropped) 2D image.""" + + # gradients + gx = np.zeros_like(image) + gy = np.zeros_like(image) + + gx[:-1, :] = np.diff(image, axis=0) + gy[:, :-1] = np.diff(image, axis=1) + + # gradient magnitude + return float(np.sum(np.sqrt(gx**2 + gy**2))) + + +def sharpness_metric(image: np.ndarray) -> float: + """Variance of the Laplacian; higher values indicate sharper images.""" + + lap = laplace(image) + + return float(np.var(lap)) + + +def square_mask_bounds( + image: cp.ndarray, + centre_x: float, + centre_y: float, + size_pct: float, +) -> tuple[int, int, int, int]: + height, width = image.shape[:2] + + # Square side length as percentage of image size + side = int(round(min(height, width) * size_pct / 100)) + + half_side = side // 2 + + x_start = max(0, int(round(centre_x)) - half_side) + x_end = min(width, int(round(centre_x)) + half_side) + + y_start = max(0, int(round(centre_y)) - half_side) + y_end = min(height, int(round(centre_y)) + half_side) + + return x_start, x_end, y_start, y_end + + +# Rescale each image independently to uint16 +def to_uint16(arr: np.ndarray) -> np.ndarray: + arr_min = np.nanmin(arr) + arr_max = np.nanmax(arr) + + if arr_max == arr_min: + return np.zeros_like(arr, dtype=np.uint16) + + arr = (arr - arr_min) / (arr_max - arr_min) + return np.round(arr * 65535).astype(np.uint16) diff --git a/tests/test_recon/test_rotation.py b/tests/test_recon/test_rotation.py index 273977a4..bc1e7719 100644 --- a/tests/test_recon/test_rotation.py +++ b/tests/test_recon/test_rotation.py @@ -9,6 +9,7 @@ from httomolibgpu.recon.rotation import ( _calculate_chunks, find_center_360, + find_center_metric_recon, find_center_vo, find_center_pc, ) @@ -28,6 +29,39 @@ def test_find_center_vo(data, flats, darks): assert_allclose(cor, 79.5) +@pytest.mark.parametrize( + "metric", + [ + "tv", + "entropy", + "sharpness", + ], +) +def test_find_center_metric_recon(data, host_angles, flats, darks, metric): + data_normalize = dark_flat_field_correction(cp.copy(data), flats, darks) + + # --- testing the center of rotation on tomo_standard ---# + cor = find_center_metric_recon( + data_normalize, + angles=host_angles, + metric_type=metric, + range=10, + step=0.5, + reconstruction_method="LPRec3d_tomobar", + ind=64, + gaussian_filter_sigma=2.0, + squared_mask_x_y_size=[20, 20, 20], + ) + + data = None + if metric == "tv": + assert_allclose(cor, 88.0) + elif metric == "entropy": + assert_allclose(cor, 81.5) + elif metric == "sharpness": + assert_allclose(cor, 84.5) + + def test_find_center_vo_ones(ensure_clean_memory): mat = cp.ones(shape=(103, 450, 230), dtype=cp.float32) cor = find_center_vo(mat)