Make Paganin processor work with one angle (#1920)

Make Paganin processor work with one angle Signed-off-by: Hannah Robarts <77114597+hrobarts@users.noreply.github.com>
TomographicImaging · Dec 13, 2024 · cdb2bc9 · cdb2bc9
1 parent 554a7fd
commit cdb2bc9
Show file tree

Hide file tree

Showing 3 changed files with 60 additions and 97 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -3,6 +3,7 @@
   - Fix bug with 'median' and 'mean' methods in Masker averaging over the wrong axes.
   - `SPDHG` `gamma` parameter is now applied correctly so that the product of the dual and primal step sizes remains constant as `gamma` varies (#1644)
   - Allow MaskGenerator to be run on DataContainers (#2001)
+  - Make Paganin Processor work with AcquistionData with one angle (#1920)
   - Fix bug passing `kwargs` to PDHG (#2010)
 - Enhancements:
   - Removed multiple exits from numba implementation of KullbackLeibler divergence (#1901)
@@ -13,6 +14,7 @@
 - Changes that break backwards compatibility:
     - Deprecated `norms` and `prob` in the `SPDHG` algorithm to be set in the `BlockOperator` and `Sampler` respectively (#1644)
     - The `run` method in the cil algorithm class will no longer run if a number of iterations is not passed (#1940)
+    - Paganin processor now requires the CIL data order (#1920)
 
 
 * 24.2.0
@@ -48,7 +50,6 @@
     - Armijo step size rule now by default initialises the search for a step size from the previously calculated step size (#1934)
   - Changes that break backwards compatibility:
     - CGLS will no longer automatically stop iterations once a default tolerance is reached. The option to pass `tolerance` will be deprecated to be replaced by `optimisation.utilities.callbacks` (#1892)
-
 
 * 24.1.0
   - New Features:

diff --git a/Wrappers/Python/cil/processors/PaganinProcessor.py b/Wrappers/Python/cil/processors/PaganinProcessor.py
@@ -204,101 +204,68 @@ def check_input(self, data):
 
         if data.dtype!=np.float32:
             raise TypeError('Processor only support dtype=float32')
+
+        cil_order = tuple(AcquisitionDimension.get_order_for_engine('cil',data.geometry))
+        if data.dimension_labels != cil_order:
+            raise ValueError("This processor requires CIL data order, consider using `data.reorder('cil')`")
 
         return True
 
     def process(self, out=None):
 
-        data  = self.get_input()
-        cil_order = tuple(AcquisitionDimension.get_order_for_engine('cil',data.geometry))
-        if data.dimension_labels != cil_order:
-            log.warning(msg="This processor will work most efficiently using\
-                        \nCIL data order, consider using `data.reorder('cil')`")
-
+        data  = self.get_input()            
         # set the geometry parameters to use from data.geometry unless the
         # geometry is overridden with an override_geometry
         self._set_geometry(data.geometry, self.override_geometry)
 
         if out is None:
             out = data.geometry.allocate(None)
 
-        # make slice indices to get the projection
-        slice_proj = [slice(None)]*len(data.shape)
-        angle_axis = data.get_dimension_axis('angle')
-        slice_proj[angle_axis] = 0
-
-        if data.geometry.channels>1:
-            channel_axis = data.get_dimension_axis('channel')
-            slice_proj[channel_axis] = 0
-        else:
-            channel_axis = None
-
-        data_proj = data.as_array()[tuple(slice_proj)]
-
-        # create an empty axis if the data is 2D
-        if len(data.shape) == 2:
-            data.array = np.expand_dims(data.array, len(data.shape))
-            slice_proj.append(slice(None))
-            data_proj = data.as_array()[tuple(slice_proj)]
-
-        elif len(data_proj.shape) == 2:
-            pass
-        else:
-            raise(ValueError('Data must be 2D or 3D per channel'))
-
-        if len(out.shape) == 2:
-            out.array = np.expand_dims(out.array, len(out.shape))
+        target_shape = (
+            data.get_dimension_size('channel')  if 'channel' in data.dimension_labels else 1,
+            data.get_dimension_size('angle') if 'angle' in data.dimension_labels else 1,
+            data.get_dimension_size('vertical') if 'vertical' in data.dimension_labels else 1, 
+            data.get_dimension_size('horizontal') if 'horizontal' in data.dimension_labels else 1)
+
+        data.array = np.reshape(data.array, target_shape)
+        out.array = np.reshape(out.array, target_shape)
 
         # create a filter based on the shape of the data
-        filter_shape = np.shape(data_proj)
-        self.filter_Nx = filter_shape[0]+self.pad*2
-        self.filter_Ny = filter_shape[1]+self.pad*2
+        proj_shape = data.array.shape[-2:]
+        self.filter_Nx = proj_shape[0]+self.pad*2
+        self.filter_Ny = proj_shape[1]+self.pad*2
         self._create_filter(self.override_filter)
 
+        # allocate padded buffer
+        padded_buffer = np.zeros((self.filter_Nx, self.filter_Ny), dtype=data.dtype)
+        buffer_slice = (slice(self.pad, self.pad + proj_shape[0]), 
+                    slice(self.pad, self.pad + proj_shape[1]))
+
         # pre-calculate the scaling factor
         scaling_factor = -(1/self.mu)
-
-        # allocate padded buffer
-        padded_buffer = np.zeros(tuple(x+self.pad*2 for x in data_proj.shape), dtype=data.dtype)
+        mag2 = self.magnification ** 2
 
-        # make slice indices to unpad the data
-        if self.pad>0:
-            slice_pad = tuple([slice(self.pad,-self.pad)]
-                                *len(padded_buffer.shape))
-        else:
-            slice_pad = tuple([slice(None)]*len(padded_buffer.shape))
         # loop over the channels
-        mag2 = self.magnification**2
         for j in range(data.geometry.channels):
-            if channel_axis is not None:
-                slice_proj[channel_axis] = j
-            # loop over the projections
-            for i in tqdm(range(len(data.geometry.angles))):
-
-                slice_proj[angle_axis] = i
-                padded_buffer.fill(0)
-                padded_buffer[slice_pad] = data.array[(tuple(slice_proj))]
+            # loop over the angles
+            for i in tqdm(range(target_shape[1])):
+                padded_buffer[:] = 0
+                padded_buffer[buffer_slice] = data.array[j, i, :, :]
 
                 if self.full_retrieval==True:
                     # apply the filter in fourier space, apply log and scale
-                    # by magnification
                     padded_buffer*=mag2
                     fI = fft2(padded_buffer)
                     iffI = ifft2(fI*self.filter).real
                     np.log(iffI, out=padded_buffer)
-                    # apply scaling factor
                     np.multiply(scaling_factor, padded_buffer, out=padded_buffer)
                 else:
                     # apply the filter in fourier space
                     fI = fft2(padded_buffer)
                     padded_buffer[:] = ifft2(fI*self.filter).real
-
-                if data.geometry.channels>1:
-                    out.fill(padded_buffer[slice_pad], angle = i, 
-                             channel=j)
-                else:
-                    out.fill(padded_buffer[slice_pad], angle = i)
-
+
+                out.array[j, i, :, :] = padded_buffer[buffer_slice]
+
         data.array = np.squeeze(data.array)
         out.array = np.squeeze(out.array)
         return out

diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py
@@ -18,6 +18,9 @@
 
 import unittest
 import numpy
+import sys
+import os
+
 from cil.framework import DataContainer, ImageGeometry, ImageData, VectorGeometry, AcquisitionData, AcquisitionGeometry
 from cil.utilities import dataexample
 from cil.utilities import quality_measures
@@ -2904,6 +2907,11 @@ def test_PaganinProcessor_check_input(self):
             with self.assertRaises(TypeError):
                 processor.set_input(dc)
 
+            # check with different data order
+            data.reorder('astra')
+            with self.assertRaises(ValueError):
+                processor.set_input(data)
+
 
     def test_PaganinProcessor_set_geometry(self):
         processor = PaganinProcessor()
@@ -3059,9 +3067,11 @@ def test_PaganinProcessor_create_filter(self):
     def test_PaganinProcessor(self):
 
         wavelength = (constants.h*constants.speed_of_light)/(40000*constants.electron_volt)
-        mu = 4.0*numpy.pi*1e-2/(wavelength)        
+        mu = 4.0*numpy.pi*1e-2/(wavelength)
 
         data_array = [self.data_cone, self.data_parallel, self.data_multichannel]
+        sys.stdout = open(os.devnull, "w")
+        sys.stderr = open(os.devnull, "w")
         for data in data_array:
             data.geometry.config.units = 'm'
             data_abs = -(1/mu)*numpy.log(data)
@@ -3099,37 +3109,12 @@ def test_PaganinProcessor(self):
             numpy.testing.assert_allclose(thickness.as_array(), thickness_inline.as_array())
             filtered_image_inline = PaganinProcessor(full_retrieval=False, pad=10)(data, override_geometry={'propagation_distance':1})
             numpy.testing.assert_allclose(filtered_image.as_array(), filtered_image_inline.as_array())
-
-            # check with different data order
-            data.reorder('astra')
-            data_abs = -(1/mu)*numpy.log(data)
-            processor = PaganinProcessor(full_retrieval=True, pad=10, return_units='m')
-            processor.set_input(data)
-            with self.assertLogs(level='WARN') as log:
-                thickness = processor.get_output(override_geometry={'propagation_distance':1})
-            self.assertLessEqual(quality_measures.mse(thickness, data_abs), 1e-5)
-            processor = PaganinProcessor(full_retrieval=False, pad=10)
-            processor.set_input(data)
-            with self.assertLogs(level='WARN') as log:
-                filtered_image = processor.get_output(override_geometry={'propagation_distance':1})
-            self.assertLessEqual(quality_measures.mse(filtered_image, data), 1e-5)
-
-            # check with different channel data order
-            if data.geometry.channels>1:
-                data.reorder(('vertical','channel','horizontal','angle'))
-                data_abs = -(1/mu)*numpy.log(data)
-                processor = PaganinProcessor(full_retrieval=True, pad=10, return_units='m')
-                processor.set_input(data)
-                with self.assertLogs(level='WARN') as log:
-                    thickness = processor.get_output(override_geometry={'propagation_distance':1})
-                self.assertLessEqual(quality_measures.mse(thickness, data_abs), 1e-5)
-                processor = PaganinProcessor(full_retrieval=False, pad=10)
-                processor.set_input(data)
-                with self.assertLogs(level='WARN') as log:
-                    filtered_image = processor.get_output(override_geometry={'propagation_distance':1})
-                self.assertLessEqual(quality_measures.mse(filtered_image, data), 1e-5)
+
+        sys.stdout = sys.__stdout__
+        sys.stderr = sys.__stderr__
 
     def test_PaganinProcessor_2D(self):
+
         self.data_parallel.geometry.config.units = 'm'
         data_slice = self.data_parallel.get_slice(vertical=10)
         wavelength = (constants.h*constants.speed_of_light)/(40000*constants.electron_volt)
@@ -3138,21 +3123,31 @@ def test_PaganinProcessor_2D(self):
 
         processor = PaganinProcessor(pad=10)
         processor.set_input(data_slice)
+        sys.stdout = open(os.devnull, "w")
+        sys.stderr = open(os.devnull, "w")
         output = processor.get_output(override_geometry={'propagation_distance':1})
+        sys.stdout = sys.__stdout__
+        sys.stderr = sys.__stderr__
         self.assertLessEqual(quality_measures.mse(output, thickness), 0.05)
 
-        # check with different data order
-        data_slice.reorder(('horizontal','angle'))
+    def test_PaganinProcessor_1angle(self):
+        data = self.data_cone.get_slice(angle=1)
+        data.geometry.config.units = 'm'
         wavelength = (constants.h*constants.speed_of_light)/(40000*constants.electron_volt)
         mu = 4.0*numpy.pi*1e-2/(wavelength) 
-        thickness = -(1/mu)*numpy.log(data_slice)
+        thickness = -(1/mu)*numpy.log(data)
 
         processor = PaganinProcessor(pad=10)
-        processor.set_input(data_slice)
+        processor.set_input(data)
+        sys.stdout = open(os.devnull, "w")
+        sys.stderr = open(os.devnull, "w")
         output = processor.get_output(override_geometry={'propagation_distance':1})
+        sys.stdout = sys.__stdout__
+        sys.stderr = sys.__stderr__
         self.assertLessEqual(quality_measures.mse(output, thickness), 0.05)
 
 
+
 if __name__ == "__main__":
 
     d = TestDataProcessor()