Tigre algorithm wrapping

Wrappers/Python/cil/plugins/tigre/Algorithms.py

@@ -0,0 +1,219 @@
+#  Copyright 2025 United Kingdom Research and Innovation
+#  Copyright 2025 The University of Manchester
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+#
+# Authors:
+# CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
+
+from cil.recon import Reconstructor
+try:
+    import tigre.algorithms as algs
+    from cil.plugins.tigre import CIL2TIGREGeometry
+except ImportError:
+    raise ImportError("TIGRE is not installed. Please install it to use this module.")
+from cil.framework import ImageData
+import logging
+import numpy as np
+import warnings
+from cil.framework.labels import AcquisitionDimension
+
+log = logging.getLogger(__name__)
+
+try:
+    from tigre.utilities.gpu import GpuIds
+    has_gpu_sel = True
+except ModuleNotFoundError:
+    has_gpu_sel = False
+
+import weakref
+
+class tigre_algo_wrapper(Reconstructor):
+
+    def __init__(self, algorithm_name=None,  initial=None, image_geometry=None,  data=None, number_iterations=0,  **kwargs):
+        """
+        A wrapper for TIGRE algorithms, allowing the use of CIL geometries and data.
+
+        Parameters
+        ----------
+        algorithm_name : str
+            algorithm_Name of the TIGRE algorithm to use (e.g., 'ART', 'SART', 'SIRT', 'OSSART').
+        initial : ImageData, optional
+            Initial guess for the reconstruction. If None, a zero-initialized image is used.
+        image_geometry : ImageGeometry
+            The geometry of the image to be reconstructed.
+        data : AcquisitionData
+            The measured projection data.
+        number_iterations : int, default=0
+            Number of iterations for the reconstruction algorithm.
+        **kwargs : dict
+            Additional keyword arguments passed to the TIGRE reconstruction algorithm.
+
+        Returns
+        -------
+        ImageData
+            The reconstructed image.
+        quality : float
+            Quality measures computed by the algorithm, if applicable. See the tigre algorithm documentation for details.
+
+        Raises
+        ------
+        ValueError
+            If  `data` is None.
+
+        Notes
+        -----
+        This class is designed to facilitate the use of TIGRE algorithms within the CIL framework,
+        allowing for the use of CIL's `ImageGeometry` and `AcquisitionData` classes. It handles the conversion
+        of CIL geometries to TIGRE geometries and prepares the data for the specified algorithm.
+        The `algorithm_name` parameter should match one of the available TIGRE algorithms for example: 'art', 'sirt', 'sart', 'ossart', 'cgls', 'lsmr', 'hybrid_lsqr', 'ista', 'fista', 'sart_tv', 'ossart_tv'.
+
+        Note
+        ----
+        We are aware that running the TIGRE algorithms: ISTA, FISTA, SART_TV, OSSART_TV using 2D data can lead to incorrect restults in the TV denoising step, particularly when using more than one GPU. https://github.com/CERN/TIGRE/issues/681
+        You can change the gpuids by passing the `gpuids` keyword argument, for example:
+        ```python
+        from tigre.utilities.gpu import GpuIds
+        gpuids = GpuIds()
+        gpuids.devices = [0]  # Specify the GPU device IDs you want to use
+        algo = tigre_algo_wrapper(algorithm_name='fista', initial=initial_image, image_geometry=image_geom, data=acquisition_data, number_iterations=10, gpuids=gpuids)
+        ```
+
+
+        Example
+        -------
+        >>> from cil.plugins.tigre import tigre_algo_wrapping 
+        >>> algo = tigre_algo_wrapper(algorithm_name='SART', initial=initial_image, image_geometry=image_geom, data=acquisition_data, number_iterations=10)
+        >>> reconstructed_image, quality = algo.run()
+
+        """
+
+
+        if data is None:
+            raise ValueError("`data` is required")
+        if image_geometry is None and initial is None:
+            image_geometry = data.geometry.get_ImageGeometry()
+        elif image_geometry is None and initial is not None:
+            image_geometry = initial.geometry
+
+        log.info("%s setting up tigre geometry", self.__class__.__name__)
+
+        if initial is None:
+            initial = image_geometry.allocate(0)
+
+        self.tigre_initial = initial.copy().as_array()
+        ig = image_geometry
+        ag = data.geometry
+        self.tigre_geom, self.tigre_angles = CIL2TIGREGeometry.getTIGREGeometry(
+            ig, ag)
+        self.tigre_projections = data.as_array()
+
+        if self.tigre_projections.ndim == 2:
+            if any( a==algorithm_name for a in ['ista', 'fista', 'sart_tv', 'ossart_tv']):
+                warnings.warn(
+                    "We are aware that the TIGRE algorithms: ISTA, FISTA, SART_TV, OSSART_TV using 2D data can lead to incorrect results in the TV denoising step, particularly when using more than one GPU.", UserWarning, stacklevel=2)
+
+
+        if data.dimension_labels[0] != AcquisitionDimension.ANGLE:
+            self.tigre_projections = np.expand_dims(self.tigre_projections, axis=0)
+
+        if self.tigre_geom.is2D:
+            self.tigre_projections = np.expand_dims(self.tigre_projections, axis=1)
+            self.tigre_initial = np.expand_dims(self.tigre_initial, axis=0)
+
+        self.tigre_algo = getattr(algs, algorithm_name)
+        self.number_iterations = number_iterations
+        self.kwargs = kwargs
+        if has_gpu_sel:
+            self.gpuids = self.kwargs.pop('gpuids', None)
+            if self.gpuids is None: 
+                self.gpuids =  GpuIds()
+            log.info("Using GPU ids:", self.gpuids)
+
+        self._input = None  
+
+        super(tigre_algo_wrapper, self).__init__(data, image_geometry=ig, backend='tigre')
+
+        log.info("%s configured", self.__class__.__name__)
+
+    def set_input(self, input):
+        """
+        When called by the parent class during initialisation, sets the input data to run the reconstructor on. The geometry of the dataset must be compatible with the reconstructor.
+        When called after initialisation, raises NotImplementedError as changing the input is not currently supported.
+        Parameters
+        ----------
+        input : AcquisitionData
+            A dataset with a compatible geometry
+        """
+        if self._input is None:
+            if input.geometry != self.acquisition_geometry:
+                raise ValueError ("Input not compatible with configured reconstructor. Initialise a new reconstructor with this geometry")
+            else:
+                self._input = weakref.ref(input)
+
+        else:
+            raise NotImplementedError("Setting the input after initialisation is not currently supported.")
+
+    def run(self, out=None):
+        """
+        Run the specified TIGRE algorithm with the provided parameters.
+
+        Parameters
+        ----------
+        out : ImageData, optional
+            Output image data to store the result. If None, a new ImageData object is created.
+
+        Returns
+        -------
+        out : ImageData
+            The reconstructed image data.
+        quality : float
+            Quality measures computed by the tigre algorithm, if applicable.
+
+        """
+
+        log.info("%s passing to the tigre algorithm", self.__class__.__name__)
+        if has_gpu_sel:
+            result = self.tigre_algo(
+                proj=self.tigre_projections,
+                geo=self.tigre_geom,
+                angles=self.tigre_angles,
+                init=self.tigre_initial,
+                niter=self.number_iterations,
+                gpuids=self.gpuids,
+                **self.kwargs
+            )
+        else:
+            result = self.tigre_algo(
+                proj=self.tigre_projections,
+                geo=self.tigre_geom,
+                angles=self.tigre_angles,
+                init=self.tigre_initial,
+                niter=self.number_iterations,
+                **self.kwargs
+            )
+        if isinstance(result, tuple):
+            img = result[0]
+            quality = result[1]
+        else:
+            img = result
+            quality = None
+
+        if out is None:
+            out = self._image_geometry.allocate(0)
+
+        out.fill(img)
+
+        log.info("%s completed", self.__class__.__name__)
+
+        return out, quality

Wrappers/Python/cil/plugins/tigre/__init__.py

@@ -19,3 +19,4 @@
 from .Geometry import CIL2TIGREGeometry
 from .ProjectionOperator import ProjectionOperator
 from .FBP import FBP
+from .Algorithms import tigre_algo_wrapper

Wrappers/Python/test/test_PluginsTigre_Algorithms.py

@@ -0,0 +1,147 @@
+#  Copyright 2025 United Kingdom Research and Innovation
+#  Copyright 2025 The University of Manchester
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+#
+# Authors:
+# CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt
+import unittest
+from utils import has_tigre, has_nvidia, initialise_tests
+
+from cil.utilities.dataexample import (
+    SIMULATED_PARALLEL_BEAM_DATA,
+    SIMULATED_CONE_BEAM_DATA,
+    SIMULATED_SPHERE_VOLUME
+)
+from cil.processors import TransmissionAbsorptionConverter, Slicer
+from cil.framework import ImageData
+
+import numpy as np
+from unittest_parametrize import parametrize
+from unittest_parametrize import ParametrizedTestCase
+
+initialise_tests()
+
+import warnings
+
+from testclass import CCPiTestClass
+
+if has_tigre:
+    from tigre.utilities.gpu import GpuIds
+    from cil.plugins.tigre import ProjectionOperator, tigre_algo_wrapper
+
+class TestTigreReconstructionAlgorithms(ParametrizedTestCase,  unittest.TestCase):
+
+
+    def get_geometry_data(self, geometry_type):
+        if geometry_type == "parallel_2d":
+            data = SIMULATED_PARALLEL_BEAM_DATA.get().get_slice(vertical='centre')
+            gt = SIMULATED_SPHERE_VOLUME.get().get_slice(vertical='centre')
+        elif geometry_type == "parallel_3d":
+            data = SIMULATED_PARALLEL_BEAM_DATA.get()
+            gt = SIMULATED_SPHERE_VOLUME.get()
+        elif geometry_type == "cone_2d":
+            gt = SIMULATED_SPHERE_VOLUME.get().get_slice(vertical='centre')
+            data = SIMULATED_CONE_BEAM_DATA.get().get_slice(vertical='centre')
+        elif geometry_type == "cone_3d":
+            gt = SIMULATED_SPHERE_VOLUME.get()
+            data = SIMULATED_CONE_BEAM_DATA.get()
+        else:
+            raise ValueError(f"Unknown geometry type: {geometry_type}")
+
+        absorption = TransmissionAbsorptionConverter()(data)
+        ig = gt.geometry
+        return ig, absorption, gt
+
+
+
+
+    def run_algorithm(self, algorithm_name, geometry_type, expect_warning=False, **kwargs):
+        ig, absorption, gt = self.get_geometry_data(geometry_type)
+
+        if expect_warning:
+            with warnings.catch_warnings(record=True) as w:
+                warnings.simplefilter("always")
+                algo = tigre_algo_wrapper(
+                    algorithm_name=algorithm_name,
+                    initial=None,
+                    image_geometry=ig,
+                    data=absorption,
+                    number_iterations=2,
+                    **kwargs
+                )
+                img, qual = algo.run()
+                warning_msgs = [str(warn.message) for warn in w]
+                self.assertTrue(
+                    any("incorrect results in the TV denoising step" in msg for msg in warning_msgs),
+                    f"Expected warning not raised for {algorithm_name} with {geometry_type}"
+                )
+        else:
+            algo = tigre_algo_wrapper(
+                algorithm_name=algorithm_name,
+                initial=None,
+                image_geometry=ig,
+                data=absorption,
+                number_iterations=2,
+                **kwargs
+            )
+            img, qual = algo.run()
+
+        self.assertIsInstance(img, ImageData)
+        self.assertEqual(img.shape, ig.shape)
+        if qual is not None:
+            self.assertTrue(isinstance(qual, (float, int, np.ndarray)))
+
+
+    @parametrize(
+        ("algorithm_name", "kwargs", "expect_warning", "geometry_type"),
+        [
+            ("sart", {}, False, "parallel_2d"),
+            ("sirt", {}, False, "parallel_3d"),
+            ("ossart", {}, False, "cone_2d"),
+            ("lsmr", {}, False, "cone_3d"),
+            ("cgls", {}, False, "parallel_2d"),
+            ('hybrid_lsqr', {}, False, "cone_2d"),
+             ("ista", {
+                "hyper": lambda self, ig, ag: ProjectionOperator(ig, ag).norm()**2,
+                "Quameasopts": ['RMSE'],
+                "tvlambda": 0.01
+            }, True, "cone_2d"),
+            ("fista", {
+                "hyper": lambda self, ig, ag: 2 * ProjectionOperator(ig, ag).norm()**2,
+                "Quameasopts": ['RMSE'],
+                "tvlambda": 0.001
+            }, True, "cone_2d"),
+            ("sart_tv", {"tvlambda": 50}, True, "parallel_2d"),
+            ("ossart_tv", {"tvlambda": 0.005}, True, "parallel_2d"),
+        ]
+    )
+    @unittest.skipUnless(has_tigre, "Requires TIGRE")
+    @unittest.skipUnless(has_nvidia, "Requires NVIDIA GPU for TIGRE")
+    def test_tigre_algorithms_with_geometries(self, algorithm_name, kwargs, expect_warning, geometry_type):
+        ig, absorption, _ = self.get_geometry_data(geometry_type)
+
+
+        gpuids = GpuIds()
+        if expect_warning:
+            gpuids.devices = [0]
+        kwargs['gpuids'] = gpuids
+
+
+        resolved_kwargs = {
+            k: v(self, ig, absorption.geometry) if callable(v) else v
+            for k, v in kwargs.items()
+        }
+        self.run_algorithm(algorithm_name, geometry_type, expect_warning=expect_warning, **resolved_kwargs)
+
+

docs/source/plugins.rst

@@ -83,7 +83,18 @@ Projection Operator
 .. autoclass:: cil.plugins.tigre.ProjectionOperator
    :members:
 
+Algorithms
+----------------
+This class is designed to facilitate the use of TIGRE algorithms within the CIL framework,
+allowing for the use of CIL's `ImageGeometry` and `AcquisitionData` classes. It handles the conversion
+of CIL geometries to TIGRE geometries and prepares the data for the specified algorithm.
+The `name` parameter should match one of the available TIGRE algorithms for example: 'art', 'sirt', 'sart', 'ossart', 'cgls', 'lsmr', 'hybrid_lsqr', 'ista', 'fista', 'sart_tv', 'ossart_tv'.
 
+Note that we provide this wrapper for the convenience of our users who use Tigre or want a specific algorithm and have tested the CIL side of converting the geomtries and passing to Tigre. However, we do not provide support for the algorithms themselves or guarantee that they will work as expected. Please refer to the TIGRE documentation for details on the algorithms and their parameters.
+
+.. autoclass:: cil.plugins.tigre.Algorithms
+   :members:
+   :inherited-members:
 
 ASTRA
 =====
@@ -110,3 +121,5 @@ Projection Operator
    :members:
 
 :ref:`Return Home <mastertoc>`
+
+