Update FluxNormaliser for Flex geometry (#2347)

* Update FluxNormaliser for Flex geometry

* simplify tests

* projection_index as an alias for now deprecated angles

* Update changelog

Author: gfardell

CHANGELOG.md

@@ -3,6 +3,7 @@
     - LSQR algorithm added to the CIL algorithm class (#1975)
     - Add `VolumeShrinker` tool to reduce the size of the reconstruction volume from an `AcquisitionData` (#2221)
     - LaminographyGeometryCorrector tool added to processors (#2259)
+    - `FluxNormaliser` can be used on `Cone3D_Flex` data (#2347)
   - Bug fixes:
     - `CentreOfRotationCorrector.image_sharpness` data is now correctly smoothed to reduce aliasing artefacts and improve robustness. (#2202)
     - `PaganinProcessor` now correctly applies scaling with magnification for cone-beam geometry (#2225)

Wrappers/Python/cil/processors/FluxNormaliser.py

@@ -116,10 +116,7 @@ def check_input(self, dataset):
         if not (type(dataset), AcquisitionData):
             raise TypeError("Expected AcquistionData, found {}"
                             .format(type(dataset)))
-        
-        if dataset.geometry.geom_type & AcquisitionType.CONE_FLEX:
-            raise NotImplementedError("FluxNormaliser does not yet support CONE3D_FLEX data")
-        
+
         image_axes = 0
         if 'vertical' in dataset.dimension_labels:
             self.v_axis = dataset.get_dimension_axis('vertical')
@@ -248,7 +245,7 @@ def _calculate_target(self):
             raise TypeError("Target must be string or a number, found {}"
                             .format(type(self.target)))
 
-    def preview_configuration(self, angle=None, channel=None, log=False):
+    def preview_configuration(self, projection_index=None, channel=None, log=False, **kwargs):
         '''
         Preview the FluxNormalisation processor configuration for roi mode.
         Plots the region of interest on the image and the mean, maximum and 
@@ -258,8 +255,8 @@ def preview_configuration(self, angle=None, channel=None, log=False):
         
         Parameters:
         -----------
-        angle: float, optional
-            Index of the angle to plot, default=None displays the data with the 
+        projection_index: int, optional
+            Index of the projection to plot, default=None displays the data with the 
             minimum and maximum pixel values in the roi. For 2D data, the roi is 
             plotted on the sinogram.
 
@@ -270,12 +267,26 @@ def preview_configuration(self, angle=None, channel=None, log=False):
         log: bool, default=False
             If True, plot the image with a log scale, default is False
 
+        **kwargs:
+            angle: int, optional
+                Deprecated alias for `projection_index`. Use `projection_index` instead.
+
         Returns:
         --------
         matplotlib.figure.Figure
             The figure object created to plot the configuration
         '''
         import matplotlib.pyplot as plt
+        
+        if 'angle' in kwargs:
+            if projection_index is not None:
+                raise TypeError("Both projection_index and angle were specified; angle is deprecated, use projection_index instead")
+            projection_index = kwargs.pop('angle')
+            warnings.warn(
+                "The 'angle' keyword argument is deprecated and will be removed in a future version; use 'projection_index' instead.",
+                DeprecationWarning, stacklevel=2)
+        if kwargs:
+            raise TypeError(f"preview_configuration() got unexpected keyword arguments {list(kwargs)}")
 
         self._calculate_flux()
 
@@ -307,48 +318,49 @@ def preview_configuration(self, angle=None, channel=None, log=False):
 
             plt.figure(figsize=(8,8))
             if data.geometry.dimension == '3D':
-                if angle is None:
-                    if 'angle' in data.dimension_labels:
+                if projection_index is None:
+                    if 'angle' in data.dimension_labels or 'projection' in data.dimension_labels:
                         self._plot_slice_roi(angle_index=numpy.argmin(min), channel_index=channel, log=log, ax=221)
                         self._plot_slice_roi(angle_index=numpy.argmax(max), channel_index=channel, log=log, ax=222)
                     else:
                         self._plot_slice_roi(log=log, channel_index=channel, ax=211)
                 else:
-                    if 'angle' in data.dimension_labels:
-                        self._plot_slice_roi(angle_index=angle, channel_index=channel, log=log, ax=211)
+                    if 'angle' in data.dimension_labels or 'projection' in data.dimension_labels:
+                        self._plot_slice_roi(angle_index=projection_index, channel_index=channel, log=log, ax=211)
                     else:
                         self._plot_slice_roi(log=log, channel_index=channel, ax=211)
 
             # if data is 2D plot roi on all angles
             elif data.geometry.dimension == '2D':
-                if angle is None:
+                if projection_index is None:
                     self._plot_slice_roi(channel_index=channel, log=log, ax=211)
                 else:
-                    raise ValueError("Cannot plot ROI for a single angle on 2D data, please specify angle=None to plot ROI on the sinogram")
+                    raise ValueError("Cannot plot ROI for a single projection on 2D data, please specify projection_index=None to plot ROI on the sinogram")
 
             plt.subplot(212)
             if data.geometry.num_projections==1:
                 plt.plot(0, flux_array, '.r', label='Mean')
                 plt.plot(0, min,'.k', label='Minimum')
                 plt.plot(0, max,'.k', label='Maximum')
             else:
-                indices = range(data.get_dimension_size('angle'))
+                indices = range(data.geometry.num_projections)
                 plt.plot(indices, flux_array, 'r', label='Mean')
                 plt.plot(indices, min,'--k', label='Minimum')
                 plt.plot(indices, max,'--k', label='Maximum')
 
             plt.legend()
-            plt.xlabel('angle index')
+            plt.xlabel('projection index')
             plt.ylabel('Intensity in roi')
             plt.grid()
 
             ax1 = plt.gca()
-            ax2 = ax1.twiny()
-            valid_ticks = [int(tick) for tick in ax1.get_xticks() if 0 <= tick < data.geometry.num_projections]
-            ax2.set_xticks(valid_ticks)
-            ax2.set_xbound(ax1.get_xbound())
-            ax2.set_xticklabels([data.geometry.angles[tick] for tick in valid_ticks])
-            ax2.set_xlabel('angle')
+            if 'angle' in data.geometry.dimension_labels:
+                ax2 = ax1.twiny()
+                valid_ticks = [int(tick) for tick in ax1.get_xticks() if 0 <= tick < data.geometry.num_projections]
+                ax2.set_xticks(valid_ticks)
+                ax2.set_xbound(ax1.get_xbound())
+                ax2.set_xticklabels([data.geometry.angles[tick] for tick in valid_ticks])
+                ax2.set_xlabel('angle')
 
             plt.tight_layout()
 
@@ -363,7 +375,7 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11
         Parameters:
         -----------
         angle_index: int, optional
-            Index of the angle to plot
+            Index of the projection to plot
         channel_index: int, optional
             Index of the channel to plot
         log: bool, optional
@@ -376,14 +388,16 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11
         data = self.get_input()
         if angle_index is not None and 'angle' in data.dimension_labels:
             data_slice = data.get_slice(angle=angle_index)
+        elif angle_index is not None and 'projection' in data.dimension_labels:
+            data_slice = data.get_slice(projection=angle_index)
         else:
             data_slice = data
 
         if 'channel' in data.dimension_labels:
             data_slice = data_slice.get_slice(channel=channel_index)
 
         if len(data_slice.shape) != 2:
-            raise ValueError("Data shape not compatible with preview_configuration(), data must have at least two of 'horizontal', 'vertical' and 'angle'")
+            raise ValueError("Data shape not compatible with preview_configuration(), data must have at least two of 'horizontal', 'vertical' and 'angle'/'projection'")
 
         # if horizontal and vertical are not specified in the roi, get the
         # min and max extent from the full size of the dimension
@@ -409,14 +423,14 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11
         v = data_slice.dimension_labels[0]
 
         # get the box to plot from the roi
-        if h == 'angle':
+        if h == 'angle' or h == 'projection':
             h_min = min_angle
             h_max = max_angle
         else:
             h_min = self.roi[h][0]
             h_max = self.roi[h][1]
 
-        if v == 'angle':
+        if v == 'angle' or v == 'projection':
             v_min = min_angle
             v_max = max_angle
         else:
@@ -431,8 +445,12 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11
         ax1.plot([h_max, h_max],[v_min, v_max],'--r')
 
         title = 'ROI'
-        if angle_index is not None:
+        if angle_index is not None and 'angle' in  data_slice.dimension_labels:
             title += ' angle = ' + str(data.geometry.angles[angle_index])
+
+        if angle_index is not None and 'projection' in  data_slice.dimension_labels:
+            title += ' projection = ' + str(angle_index)
+
         if channel_index is not None:
             title += ' channel = ' + str(channel_index)
         ax1.set_title(title)

Wrappers/Python/test/test_DataProcessor.py

@@ -3588,15 +3588,15 @@ def setUp(self):
             .set_panel([3,3])
         arr = numpy.array([[[1,2,3],[1,2,3],[1,2,3]],
                         [[4,5,6],[4,5,6],[4,5,6]], 
-                        [[7,8,9],[7,8,9],[7,8,9]]])
+                        [[7,8,9],[7,8,9],[7,8,9]]],dtype=numpy.float32)
         self.data_simple = AcquisitionData(arr, geometry=ag)
 
         source_position_set=[[0,-100000,0]]*3
         detector_position_set=[[0,0,0]]*3
         detector_direction_x_set=[[1, 0, 0]]*3
         detector_direction_y_set=[[0, 0, 1]]*3
-        cone_flex_ag = AcquisitionGeometry.create_Cone3D_Flex(source_position_set, detector_position_set, detector_direction_x_set, detector_direction_y_set).set_panel([3,3])
-        self.cone_flex  = AcquisitionData(arr, geometry=cone_flex_ag)
+        cone_flex_ag = AcquisitionGeometry.create_Cone3D_Flex(source_position_set, detector_position_set, detector_direction_x_set, detector_direction_y_set).set_panel([128,128])
+        self.cone_flex  = cone_flex_ag.allocate(1)
 
     def error_message(self,processor, test_parameter):
             return "Failed with processor " + str(processor) + " on test parameter " + test_parameter
@@ -3623,11 +3623,6 @@ def test_check_input(self):
         with self.assertRaises(ValueError):
             processor.check_input(self.data_cone)
 
-        # check there's a not implemented error if cone flex geom is used:
-        processor = FluxNormaliser(flux=[1,2,3])
-        with self.assertRaises(NotImplementedError):
-            processor.check_input(self.cone_flex)
-
     def test_calculate_flux(self):
         # check there is an error if flux array size is not equal to the number of angles in data
         processor = FluxNormaliser(flux = [1,2,3])
@@ -3814,7 +3809,8 @@ def test_preview_configuration(self, mock_show):
 
         # Test no error with preview_configuration with different data shapes
         for data in [self.data_cone, self.data_parallel, self.data_multichannel, 
-                     self.data_slice, self.data_reorder, self.data_single_angle]:
+                     self.data_slice, self.data_reorder, self.data_single_angle,
+                     self.cone_flex]:
             mock_show.reset_mock()
 
             roi = {'horizontal':(25,40)}
@@ -3826,15 +3822,15 @@ def test_preview_configuration(self, mock_show):
 
             # for 3D, check no error specifying a single angle to plot
             if data.geometry.dimension == '3D':
-                processor.preview_configuration(angle=1)
+                processor.preview_configuration(projection_index=1)
             # if 2D, attempt to plot single angle should cause error
             else:
                 with self.assertRaises(ValueError):
-                    processor.preview_configuration(angle=1)
+                    processor.preview_configuration(projection_index=1)
 
             # if data is multichannel, check no error specifying a single channel to plot
             if 'channel' in data.dimension_labels:
-                processor.preview_configuration(angle=1, channel=1)
+                processor.preview_configuration(projection_index=1, channel=1)
                 processor.preview_configuration(channel=1)
             # if single channel, check specifying channel causes an error
             else:
@@ -3844,6 +3840,25 @@ def test_preview_configuration(self, mock_show):
         # Re-enable logging
         logging.disable(logging.NOTSET)
 
+    @unittest.skipIf(not has_matplotlib, "matplotlib not installed")
+    @patch('matplotlib.pyplot.show')
+    def test_preview_configuration_deprecated_angle(self, mock_show):
+        logging.disable(logging.CRITICAL)
+
+        roi = {'horizontal':(25,40)}
+        processor = FluxNormaliser(roi=roi)
+        processor.set_input(self.data_cone)
+
+        with self.assertWarns(DeprecationWarning):
+            processor.preview_configuration(angle=1)
+
+        processor.preview_configuration(projection_index=1)
+
+        with self.assertRaises(TypeError):
+            processor.preview_configuration(projection_index=1, angle=1)
+
+        logging.disable(logging.NOTSET)
+
     def test_FluxNormaliser(self, accelerated=False):
 
         # Suppress backround range warning
@@ -3855,6 +3870,12 @@ def test_FluxNormaliser(self, accelerated=False):
         data_norm = processor.get_output()
         numpy.testing.assert_allclose(data_norm.array, self.data_cone.array)
 
+        #Test flux with no target on ConeFlex
+        processor = FluxNormaliser(flux=1, accelerated=accelerated)
+        processor.set_input(self.cone_flex)
+        data_norm = processor.get_output()
+        numpy.testing.assert_allclose(data_norm.array, self.cone_flex.array)
+
         #Test flux with target
         processor = FluxNormaliser(flux=10, target=5.0, accelerated=accelerated)
         processor.set_input(self.data_cone)