diff --git a/CHANGELOG.md b/CHANGELOG.md index a226ea1b5..41ed4e4e7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ - New Features: - Added method to plot filter in `GenericFilteredBackProjection` (#1667) - Added wavelet operator, wrapping PyWavelets operator as a CIL operator (#1618) + - Added PaganinProcessor processor, to perform phase retrieval from phase contrast images (#1737) - Added L1Sparsity function, allowing calculations of `|Ax-b|_1` and its proximal, in the case of orthogonal operators, `A` (#1618) - Options in algorithms GD, ISTA and FISTA to pass a `cil.optimisation.utilities.StepSizeRule` or a `cil.optimisation.utilities.Preconditioner`(#1768) - an implementation of the Armijo Rule as a child class of `cil.optimisation.utilities.StepSizeRule` (#1768) diff --git a/Wrappers/Python/cil/framework/framework.py b/Wrappers/Python/cil/framework/framework.py index da02038b2..0ea32c538 100644 --- a/Wrappers/Python/cil/framework/framework.py +++ b/Wrappers/Python/cil/framework/framework.py @@ -967,7 +967,22 @@ def get_centre_slice(self): return self def calculate_magnification(self): - return [None, None, 1.0] + '''Method to calculate magnification and distance from the sample to + the detector using the detector positions and the rotation axis. + For parallel beam geometry magnification = 1 + + Returns + ------- + list + A list containing the [0] distance from the source to the rotate + axis, [1] distance from the rotate axis to the detector, + [2] magnification of the system + + ''' + ab = (self.rotation_axis.position - self.detector.position) + dist_center_detector = float(numpy.sqrt(ab.dot(ab))) + + return [None, dist_center_detector, 1.0] class Parallel3D(SystemConfiguration): r'''This class creates the SystemConfiguration of a parallel beam 3D tomographic system @@ -1116,7 +1131,22 @@ def __eq__(self, other): return False def calculate_magnification(self): - return [None, None, 1.0] + '''Method to calculate magnification and distance from the sample to + the detector using the detector positions and the rotation axis. + For parallel beam geometry magnification = 1 + + Returns + ------- + list + A list containing the [0] distance from the source to the rotate + axis, [1] distance from the rotate axis to the detector, + [2] magnification of the system + + ''' + ab = (self.rotation_axis.position - self.detector.position) + dist_center_detector = float(numpy.sqrt(ab.dot(ab))) + + return [None, dist_center_detector, 1.0] def get_centre_slice(self): """Returns the 2D system configuration corresponding to the centre slice diff --git a/Wrappers/Python/cil/processors/PaganinProcessor.py b/Wrappers/Python/cil/processors/PaganinProcessor.py new file mode 100644 index 000000000..33ad87841 --- /dev/null +++ b/Wrappers/Python/cil/processors/PaganinProcessor.py @@ -0,0 +1,574 @@ +# Copyright 2024 United Kingdom Research and Innovation +# Copyright 2024 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.framework import Processor, AcquisitionData, DataOrder + +import numpy as np +from scipy.fft import fft2 +from scipy.fft import ifft2 +from scipy.fft import ifftshift +from scipy import constants +from tqdm import tqdm +import logging + +log = logging.getLogger(__name__) + +class PaganinProcessor(Processor): + + r""" + Processor to retrieve quantitative information from phase contrast images + using the Paganin phase retrieval algorithm described in [1] + + Parameters + ---------- + delta: float (optional) + Real part of the deviation of the material refractive index from 1, + where refractive index :math:`n = (1 - \delta) + i \beta` energy- + dependent refractive index information for x-ray wavelengths can be + found at [2], default is 1 + + beta: float (optional) + Complex part of the material refractive index, where refractive index + :math:`n = (1 - \delta) + i \beta` energy-dependent refractive index + information for x-ray wavelengths can be found at [2], default is 1e-2 + + energy: float (optional) + Energy of the incident photon, default is 40000 + + energy_units: string (optional) + Energy units, default is 'eV' + + full_retrieval : bool, optional + If True, perform the full phase retrieval and return the thickness. If + False, return a filtered image, default is True + + filter_type: string (optional) + The form of the Paganin filter to use, either 'paganin_method' + (default) or 'generalised_paganin_method' as described in [3] + + pad: int (optional) + Number of pixels to pad the image in Fourier space to reduce aliasing, + default is 0 + + return_units: string (optional) + The distance units to return the sample thickness in, must be one of + 'm', 'cm', 'mm' or 'um'. Only applies if full_retrieval=True (default + is'cm') + + Returns + ------- + AcquisitionData + AcquisitionData corrected for phase effects, retrieved sample thickness + or (if :code:`full_retrieval=False`) filtered data + + Example + ------- + >>> processor = PaganinProcessor(delta=5, beta=0.05, energy=18000) + >>> processor.set_input(data) + >>> thickness = processor.get_output() + + Example + ------- + >>> processor = PaganinProcessor(delta=1,beta=10e2, full_retrieval=False) + >>> processor.set_input(data) + >>> filtered_image = processor.get_output() + + Example + ------- + >>> processor = PaganinProcessor() + >>> processor.set_input(data) + >>> thickness = processor.get_output(override_filter={'alpha':10}) + >>> phase_retrieved_image = thickness*processor.mu + + Notes + ----- + This processor will work most efficiently using the cil data order with + `data.reorder('cil')` + + Notes + ----- + This processor uses the phase retrieval algorithm described by Paganin et + al. [1] to retrieve the sample thickness + + .. math:: T(x,y) = - \frac{1}{\mu}\ln\left (\mathcal{F}^{-1}\left + (\frac{\mathcal{F}\left ( M^2I_{norm}(x, y,z = \Delta) \right )}{1 + + \alpha\left ( k_x^2 + k_y^2 \right )} \right )\right ), + + where + + - :math:`T`, is the sample thickness, + - :math:`\mu = \frac{4\pi\beta}{\lambda}` is the material linear + attenuation coefficient where :math:`\beta` is the complex part of the + material refractive index and :math:`\lambda=\frac{hc}{E}` is the probe + wavelength, + - :math:`M` is the magnification at the detector, + - :math:`I_{norm}` is the input image which is expected to be the + normalised transmission data, + - :math:`\Delta` is the propagation distance, + - :math:`\alpha = \frac{\Delta\delta}{\mu}` is a parameter determining + the strength of the filter to be applied in Fourier space where + :math:`\delta` is the real part of the deviation of the material + refractive index from 1 + - :math:`k_x, k_y = \left ( \frac{2\pi p}{N_xW}, \frac{2\pi q}{N_yW} + \right )` where :math:`p` and :math:`q` are co-ordinates in a Fourier + mesh in the range :math:`-N_x/2` to :math:`N_x/2` for an image with + size :math:`N_x, N_y` and pixel size :math:`W`. + + A generalised form of the Paganin phase retrieval method can be called + using :code:`filter_type='generalised_paganin_method'`, which uses the + form of the algorithm described in [2] + + .. math:: T(x,y) = -\frac{1}{\mu}\ln\left (\mathcal{F}^{-1}\left (\frac{ + \mathcal{F}\left ( M^2I_{norm}(x, y,z = \Delta) \right )}{1 - \frac{2 + \alpha}{W^2}\left ( \cos(Wk_x) + \cos(Wk_y) -2 \right )} \right ) + \right ) + + The phase retrieval is valid under the following assumptions + + - used with paraxial propagation-induced phase contrast images which + can be assumed to be single-material locally + - using intensity data which has been flat field corrected + - and under the assumption that the Fresnel number + :math:`F_N = W^2/(\lambda\Delta) >> 1` + + To apply a filter to images using the Paganin method, call + :code:`full_retrieval=False`. In this case the pre-scaling and conversion + to absorption is not applied so the requirement to supply flat field + corrected intensity data is relaxed, + + .. math:: I_{filt} = \mathcal{F}^{-1}\left (\frac{\mathcal{F}\left ( + I(x, y,z = \Delta) \right )} + {1 - \alpha\left ( k_x^2 + k_y^2 \right )} \right ) + + References + --------- + - [1] https://doi.org/10.1046/j.1365-2818.2002.01010.x + - [2] https://henke.lbl.gov/optical_constants/getdb2.html + - [3] https://iopscience.iop.org/article/10.1088/2040-8986/abbab9 + With thanks to colleagues at DTU for help with the initial implementation + of the phase retrieval algorithm + + """ + + def __init__(self, delta=1, beta=1e-2, energy=40000, + energy_units='eV', full_retrieval=True, + filter_type='paganin_method', pad=0, + return_units='cm'): + + kwargs = { + 'energy' : energy, + 'wavelength' : self._energy_to_wavelength(energy, energy_units, + return_units), + 'delta': delta, + 'beta': beta, + '_delta_user' : delta, + '_beta_user' : beta, + 'filter_Nx' : None, + 'filter_Ny' : None, + 'filter_type' : filter_type, + 'mu' : None, + 'alpha' : None, + 'pixel_size' : None, + 'propagation_distance' : None, + 'magnification' : None, + 'filter' : None, + 'full_retrieval' : full_retrieval, + 'pad' : pad, + 'override_geometry' : None, + 'override_filter' : None, + 'return_units' : return_units + } + + super(PaganinProcessor, self).__init__(**kwargs) + + def check_input(self, data): + if not isinstance(data, (AcquisitionData)): + raise TypeError('Processor only supports AcquisitionData') + + return True + + def process(self, out=None): + + data = self.get_input() + cil_order = tuple(DataOrder.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')`") + + # 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_proj.shape) == 1: + 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')) + + # 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 + self._create_filter(self.override_filter) + + # 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)) + + # 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 + 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(out.geometry.angles))): + + slice_proj[angle_axis] = i + padded_buffer[slice_pad] = data.array[(tuple(slice_proj))] + + if self.full_retrieval==True: + # apply the filter in fourier space, apply log and scale + # by magnification + fI = fft2(self.magnification**2*padded_buffer) + iffI = ifft2(fI*self.filter) + # apply scaling factor + padded_buffer = scaling_factor*np.log(iffI) + else: + # apply the filter in fourier space + fI = fft2(padded_buffer) + padded_buffer = ifft2(fI*self.filter) + if data.geometry.channels>1: + out.fill(np.squeeze(padded_buffer[slice_pad]), angle = i, + channel=j) + else: + out.fill(np.squeeze(padded_buffer[slice_pad]), angle = i) + data.array = np.squeeze(data.array) + return out + + def set_input(self, dataset): + """ + Set the input data to the processor + + Parameters + ---------- + dataset : AcquisitionData + The input AcquisitionData + """ + return super().set_input(dataset) + + def get_output(self, out=None, override_geometry=None, + override_filter=None): + r''' + Function to get output from the PaganinProcessor + + Parameters + ---------- + out : DataContainer, optional + Fills the referenced DataContainer with the processed data + + override_geometry: dict, optional + Geometry parameters to use in the phase retrieval if you want to + over-ride values found in `data.geometry`. Specify parameters as a + dictionary :code:`{'parameter':value}` where parameter is + :code:`'magnification', 'propagation_distance'` or + :code:`'pixel_size'` and value is the new value to use. Specify + distance parameters in the same units as :code:`return_units` + (default is cm). + + override_filter: dict, optional + Over-ride the filter parameters to use in the phase retrieval. + Specify parameters as :code:`{'parameter':value}` where parameter + is :code:`'delta', 'beta'` or :code:`'alpha'` and value is the new + value to use. + + Returns + ------- + AcquisitionData + AcquisitionData corrected for phase effects, retrieved sample + thickness or (if :code:`full_retrieval=False`) filtered data + + Example + ------- + >>> processor = PaganinProcessor(delta=5, beta=0.05, energy=18000) + >>> processor.set_input(data) + >>> thickness = processor.get_output() + + Example + ------- + >>> processor = PaganinProcessor(delta=1,beta=10e2, + full_retrieval=False) + >>> processor.set_input(data) + >>> filtered_image = processor.get_output() + + Example + ------- + >>> processor = PaganinProcessor() + >>> processor.set_input(data) + >>> thickness = processor.get_output(override_filter={'alpha':10}) + >>> phase_retrieved_image = thickness*processor.mu + + Notes + ----- + If :code:`'alpha'` is specified in override_filter the new value will + be used and delta will be ignored but beta will still be used to + calculate :math:`\mu = \frac{4\pi\beta}{\lambda}` which is used for + scaling the thickness, therefore it is only recommended to specify + alpha when also using :code:`get_output(full_retrieval=False)`, or + re-scaling the result by :math:`\mu` e.g. + :code:`thickness*processor.mu` If :code:`alpha` is not specified, + it will be calculated :math:`\frac{\Delta\delta\lambda}{4\pi\beta}` + + ''' + self.override_geometry = override_geometry + self.override_filter = override_filter + + return super().get_output(out) + + def __call__(self, x, out=None, override_geometry=None, + override_filter=None): + self.set_input(x) + + if out is None: + out = self.get_output(override_geometry=override_geometry, + override_filter=override_filter) + else: + self.get_output(out=out, override_geometry=override_geometry, + override_filter=override_filter) + + return out + + def _set_geometry(self, geometry, override_geometry=None): + ''' + Function to set the geometry parameters for the processor. Values are + from the data geometry unless the geometry is overridden with an + override_geometry dictionary. + ''' + + parameters = ['magnification', 'propagation_distance', 'pixel_size'] + # specify parameter names as defined in geometry + geometry_parameters = ['magnification', 'dist_center_detector', + ('pixel_size_h', 'pixel_size_v')] + # specify if parameter requires unit conversion + convert_units = [False, True, True] + + if override_geometry is None: + override_geometry = {} + + # get and check parameters from over-ride geometry dictionary + for parameter in override_geometry.keys(): + if parameter not in parameters: + raise ValueError('Parameter {} not recognised, expected one of\ + {}.'.format(parameter, parameters)) + elif (override_geometry[parameter] is None) \ + | (override_geometry[parameter] == 0): + raise ValueError("Parameter {} cannot be {}, please update \ + data.geometry.{} or over-ride with \ + processor.get_output(override_geometry= \ + {{ '{}' : value }} )"\ + .format(parameter, str(getattr(self, parameter)), + geometry_parameters[i], parameter)) + else: + self.__setattr__(parameter, override_geometry[parameter]) + + + # get and check parameters from geometry if they are not in the + # over-ride geometry dictionary + for i, parameter in enumerate(parameters): + if parameter not in override_geometry: + if type(geometry_parameters[i])==tuple: + param1 = getattr(geometry, geometry_parameters[i][0]) + param2 = getattr(geometry, geometry_parameters[i][1]) + if (param1 - param2) / (param1 + param2) >= 1e-5: + raise ValueError("Parameter {} is not homogeneous up \ + to 1e-5: got {} and {}, please update\ + geometry using data.geometry.{} and \ + data.geometry.{} or over-ride with \ + processor.get_output(\ + override_geometry={{ '{}' : value }})" + .format(parameter, str(param1), + str(param2), + geometry_parameters[i][0], + geometry_parameters[i][1], + parameter)) + else: + param1 = getattr(geometry, geometry_parameters[i]) + + if (param1 is None) | (param1 == 0): + raise ValueError("Parameter {} cannot be {}, please update\ + data.geometry.{} or over-ride with \ + processor.get_output(override_geometry\ + ={{ '{}' : value }} )" + .format(parameter, str(param1), + str(geometry_parameters[i]), + parameter)) + else: + if convert_units[i]: + param1 = self._convert_units(param1, 'distance', + geometry.config.units, + self.return_units) + self.__setattr__(parameter, param1) + + + def _create_filter(self, override_filter=None): + ''' + Function to create the Paganin filter, either using the paganin [1] or + generalised paganin [2] method + The filter is created on a mesh in Fourier space kx, ky + [1] https://doi.org/10.1046/j.1365-2818.2002.01010.x + [2] https://iopscience.iop.org/article/10.1088/2040-8986/abbab9 + ''' + if override_filter is None: + override_filter = {} + + # update any parameter which has been over-ridden with override_filter + if ('alpha' in override_filter) & ('delta' in override_filter): + log.warning(msg="Because you specified alpha, it will not be \ + calculated and therefore delta will be ignored") + + if ('delta' in override_filter): + self.delta = override_filter['delta'] + else: + self.delta = self._delta_user + + if ('beta' in override_filter): + self.beta = override_filter['beta'] + else: + self.beta = self._beta_user + + self._calculate_mu() + + if ('alpha' in override_filter): + self.alpha = override_filter['alpha'] + else: + self._calculate_alpha() + + # create the Fourier mesh + kx,ky = np.meshgrid( + np.arange(-self.filter_Nx/2, self.filter_Nx/2, 1, dtype=np.float64) + * (2*np.pi)/(self.filter_Nx*self.pixel_size), + np.arange(-self.filter_Ny/2, self.filter_Ny/2, 1, dtype=np.float64) + * (2*np.pi)/(self.filter_Ny*self.pixel_size), + sparse=False, + indexing='ij' + ) + + # create the filter using either paganin or generalised paganin method + if self.filter_type == 'paganin_method': + self.filter = ifftshift(1/(1. + self.alpha*(kx**2 + ky**2))) + elif self.filter_type == 'generalised_paganin_method': + self.filter = ifftshift(1/(1. - (2*self.alpha/self.pixel_size**2) + *(np.cos(self.pixel_size*kx) + + np.cos(self.pixel_size*ky) -2))) + else: + raise ValueError("filter_type not recognised: got {0} expected one\ + of 'paganin_method' or \ + 'generalised_paganin_method'" + .format(self.filter_type)) + + def _calculate_mu(self): + ''' + Function to calculate the linear attenutation coefficient mu + ''' + self.mu = 4.0*np.pi*self.beta/self.wavelength + + def _calculate_alpha(self): + ''' + Function to calculate alpha, a constant defining the Paganin filter + strength + ''' + self.alpha = self.propagation_distance*self.delta/self.mu + + def _energy_to_wavelength(self, energy, energy_units, return_units): + ''' + Function to convert photon energy in eV to wavelength in return_units + + Parameters + ---------- + energy: float + Photon energy + + energy_units + Energy units + + return_units + Distance units in which to return the wavelength + + Returns + ------- + float + Photon wavelength in return_units + ''' + top = self._convert_units(constants.h*constants.speed_of_light, + 'distance', 'm', return_units) + bottom = self._convert_units(energy, 'energy', energy_units, 'J') + + return top/bottom + + def _convert_units(self, value, unit_type, input_unit, output_unit): + unit_types = ['distance','energy','angle'] + + if unit_type == unit_types[0]: + unit_list = ['m','cm','mm','um'] + unit_multipliers = [1.0, 1e-2, 1e-3, 1e-6] + elif unit_type == unit_types[1]: + unit_list = ['meV', 'eV', 'keV', 'MeV', 'J'] + unit_multipliers = [1e-3, 1, 1e3, 1e6, 1/constants.eV] + elif unit_type == unit_types[2]: + unit_list = ['deg', 'rad'] + unit_multipliers = [1, np.rad2deg(1)] + else: + raise ValueError("Unit type '{}' not recognised, must be one of {}" + .format(unit_type, unit_types)) + + for x in [input_unit, output_unit]: + if x not in unit_list: + raise ValueError("Unit '{}' not recognised, must be one of {}.\ + \nGeometry units can be updated using geometry.config.units" + .format(x, unit_list)) + + return value*unit_multipliers[unit_list.index(input_unit)]\ + /unit_multipliers[unit_list.index(output_unit)] \ No newline at end of file diff --git a/Wrappers/Python/cil/processors/__init__.py b/Wrappers/Python/cil/processors/__init__.py index baae130b2..15ba249e7 100644 --- a/Wrappers/Python/cil/processors/__init__.py +++ b/Wrappers/Python/cil/processors/__init__.py @@ -26,3 +26,4 @@ from .TransmissionAbsorptionConverter import TransmissionAbsorptionConverter from .Masker import Masker from .Padder import Padder +from .PaganinProcessor import PaganinProcessor \ No newline at end of file diff --git a/Wrappers/Python/test/test_AcquisitionGeometry.py b/Wrappers/Python/test/test_AcquisitionGeometry.py index c4fd16208..30a8853dc 100644 --- a/Wrappers/Python/test/test_AcquisitionGeometry.py +++ b/Wrappers/Python/test/test_AcquisitionGeometry.py @@ -704,7 +704,8 @@ def test_get_centre_slice(self): def test_calculate_magnification(self): AG = AcquisitionGeometry.create_Parallel2D() out = AG.config.system.calculate_magnification() - self.assertEqual(out, [None, None, 1]) + detector_position = np.array(AG.config.system.detector.position) + self.assertEqual(out, [None, float(np.sqrt(detector_position.dot(detector_position))), 1]) def test_calculate_centre_of_rotation(self): AG = AcquisitionGeometry.create_Parallel2D() @@ -858,7 +859,8 @@ def test_get_centre_slice(self): def test_calculate_magnification(self): AG = AcquisitionGeometry.create_Parallel3D() out = AG.config.system.calculate_magnification() - self.assertEqual(out, [None, None, 1]) + detector_position = np.array(AG.config.system.detector.position) + self.assertEqual(out, [None, float(np.sqrt(detector_position.dot(detector_position))), 1]) def test_calculate_centre_of_rotation(self): diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index a7516d8a6..a821bba10 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -22,6 +22,7 @@ from cil.framework import ImageGeometry, VectorGeometry, AcquisitionGeometry from cil.framework import ImageData, AcquisitionData from cil.utilities import dataexample +from cil.utilities import quality_measures from cil.framework import AX, CastDataContainer, PixelByPixelDataProcessor from cil.recon import FBP @@ -29,9 +30,12 @@ from cil.processors import CentreOfRotationCorrector from cil.processors.CofR_xcorrelation import CofR_xcorrelation from cil.processors import TransmissionAbsorptionConverter, AbsorptionTransmissionConverter -from cil.processors import Slicer, Binner, MaskGenerator, Masker, Padder +from cil.processors import Slicer, Binner, MaskGenerator, Masker, Padder, PaganinProcessor import gc +from scipy import constants +from scipy.fft import ifftshift + from utils import has_astra, has_tigre, has_nvidia, has_tomophantom, initialise_tests, has_ipp initialise_tests() @@ -2756,11 +2760,319 @@ def Masker_check(self, mask, data, data_init): data_test = data.copy().as_array() data_test[2,3] = (data_test[1,3] + data_test[3,3]) / 2 data_test[4,5] = (data_test[3,5] + data_test[5,5]) / 2 + + numpy.testing.assert_allclose(res.as_array(), data_test, rtol=1E-6) - numpy.testing.assert_allclose(res.as_array(), data_test, rtol=1E-6) + +class TestPaganinProcessor(unittest.TestCase): + + def setUp(self): + self.data_parallel = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get() + self.data_cone = dataexample.SIMULATED_CONE_BEAM_DATA.get() + ag = AcquisitionGeometry.create_Parallel3D()\ + .set_angles(numpy.linspace(0,360,360,endpoint=False))\ + .set_panel([128,128],0.1)\ + .set_channels(4) + + self.data_multichannel = ag.allocate('random') + + def error_message(self,processor, test_parameter): + return "Failed with processor " + str(processor) + " on test parameter " + test_parameter + + def test_PaganinProcessor_init(self): + # test default values are initialised + processor = PaganinProcessor() + test_parameter = ['energy', 'wavelength', 'delta', 'beta', 'full_retrieval', + 'filter_type', 'pad', 'return_units'] + test_value = [40000, 1e2*(constants.h*constants.speed_of_light)/(40000*constants.electron_volt), + 1, 1e-2, True, 'paganin_method', 0, 'cm'] + + for i in numpy.arange(len(test_value)): + self.assertEqual(getattr(processor,test_parameter[i]), test_value[i], msg=self.error_message(processor, test_parameter[i])) + + # test non-default values are initialised + processor = PaganinProcessor(1, 2, 3, 'keV', False, 'string', 19, 'mm') + test_value = [3, 1e3*(constants.h*constants.speed_of_light)/(3000*constants.electron_volt), 1, 2, False, 'string', 19, 'mm'] + for i in numpy.arange(len(test_value)): + self.assertEqual(getattr(processor,test_parameter[i]), test_value[i], msg=self.error_message(processor, test_parameter[i])) + + with self.assertRaises(ValueError): + processor = PaganinProcessor(return_units='string') + + def test_PaganinProcessor_energy_to_wavelength(self): + processor = PaganinProcessor() + wavelength = processor._energy_to_wavelength(10, 'meV', 'mm') + self.assertAlmostEqual(wavelength, 0.12398419) + def test_PaganinProcessor_check_input(self): + processor = PaganinProcessor() + for data in [self.data_cone, self.data_parallel, self.data_multichannel]: + processor.set_input(data) + data2 = processor.get_input() + numpy.testing.assert_allclose(data2.as_array(), data.as_array()) + # check there is an error when the wrong data type is input + with self.assertRaises(TypeError): + processor.set_input(data.geometry) + + with self.assertRaises(TypeError): + processor.set_input(data.as_array()) + + dc = DataContainer(data.as_array()) + with self.assertRaises(TypeError): + processor.set_input(dc) + + + def test_PaganinProcessor_set_geometry(self): + processor = PaganinProcessor() + data = self.data_cone + # check there is an error when the data geometry does not have units + processor.set_input(data) + with self.assertRaises(ValueError): + processor._set_geometry(data.geometry, None) + + # check there is no error when the geometry unit is provided + data.geometry.config.units = 'um' + processor._set_geometry(data.geometry, None) + multiplier = 1e-4 # convert um to return units cm + + # check the processor finds the correct geometry values, scaled by the units + self.assertAlmostEqual(processor.propagation_distance, data.geometry.dist_center_detector*multiplier, + msg=self.error_message(processor, 'propagation_distance')) + self.assertEqual(processor.magnification, data.geometry.magnification, + msg=self.error_message(processor, 'magnification')) + self.assertAlmostEqual(processor.pixel_size, data.geometry.pixel_size_h*multiplier, + msg=self.error_message(processor, 'pixel_size')) + + # check there is an error when the data geometry does not have propagation distance, and it is not provided in override geometry + processor.set_input(self.data_parallel) + with self.assertRaises(ValueError): + processor._set_geometry(self.data_parallel.geometry, None) + + # check override_geometry + for data in [self.data_parallel, self.data_cone, self.data_multichannel]: + processor.set_input(data) + processor._set_geometry(self.data_cone.geometry, override_geometry={'propagation_distance':1,'magnification':2, 'pixel_size':3}) + + self.assertEqual(processor.propagation_distance, 1, + msg=self.error_message(processor, 'propagation_distance')) + self.assertEqual(processor.magnification, 2, + msg=self.error_message(processor, 'magnification')) + self.assertEqual(processor.pixel_size, 3, + msg=self.error_message(processor, 'pixel_size')) + + # check the processor goes back to values from geometry if the geometry over-ride is not passed + processor.set_input(self.data_cone) + processor._set_geometry(self.data_cone.geometry) + self.assertAlmostEqual(processor.propagation_distance, self.data_cone.geometry.dist_center_detector*multiplier, + msg=self.error_message(processor, 'propagation_distance')) + self.assertEqual(processor.magnification, self.data_cone.geometry.magnification, + msg=self.error_message(processor, 'magnification')) + self.assertAlmostEqual(processor.pixel_size, self.data_cone.geometry.pixel_size_h*multiplier, + msg=self.error_message(processor, 'pixel_size')) + + processor.set_input(self.data_parallel) + with self.assertRaises(ValueError): + processor._set_geometry(self.data_parallel.geometry) + + # check there is an error when the pixel_size_h and pixel_size_v are different + self.data_parallel.geometry.pixel_size_h = 9 + self.data_parallel.geometry.pixel_size_h = 10 + with self.assertRaises(ValueError): + processor._set_geometry(self.data_parallel.geometry, override_geometry={'propagation_distance':1}) + + def test_PaganinProcessor_create_filter(self): + image = self.data_cone.get_slice(angle=0).as_array() + Nx, Ny = image.shape + + delta = 1 + beta = 2 + energy = 3 + processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, return_units='m') + + # check alpha and mu are calculated correctly + wavelength = (constants.h*constants.speed_of_light)/(energy*constants.electron_volt) + mu = 4.0*numpy.pi*beta/(wavelength) + alpha = 60000*delta/mu + + self.data_cone.geometry.config.units='m' + processor.set_input(self.data_cone) + processor._set_geometry(self.data_cone.geometry) + processor.filter_Nx = Nx + processor.filter_Ny = Ny + processor._create_filter() + + self.assertEqual(processor.alpha, alpha, msg=self.error_message(processor, 'alpha')) + self.assertEqual(processor.mu, mu, msg=self.error_message(processor, 'mu')) + + kx,ky = numpy.meshgrid( + numpy.arange(-Nx/2, Nx/2, 1, dtype=numpy.float64) * (2*numpy.pi)/(Nx*self.data_cone.geometry.pixel_size_h), + numpy.arange(-Ny/2, Ny/2, 1, dtype=numpy.float64) * (2*numpy.pi)/(Nx*self.data_cone.geometry.pixel_size_h), + sparse=False, + indexing='ij' + ) + + # check default filter is created with paganin_method + filter = ifftshift(1/(1. + alpha*(kx**2 + ky**2))) + numpy.testing.assert_allclose(processor.filter, filter) + + # check generalised_paganin_method + processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, filter_type='generalised_paganin_method', return_units='m') + processor.set_input(self.data_cone) + processor._set_geometry(self.data_cone.geometry) + processor.filter_Nx = Nx + processor.filter_Ny = Ny + processor._create_filter() + filter = ifftshift(1/(1. - (2*alpha/self.data_cone.geometry.pixel_size_h**2)*(numpy.cos(self.data_cone.geometry.pixel_size_h*kx) + numpy.cos(self.data_cone.geometry.pixel_size_h*ky) -2))) + numpy.testing.assert_allclose(processor.filter, filter) + + # check unknown method raises error + processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, filter_type='unknown_method', return_units='m') + processor.set_input(self.data_cone) + processor._set_geometry(self.data_cone.geometry) + processor.filter_Nx = Nx + processor.filter_Ny = Ny + with self.assertRaises(ValueError): + processor._create_filter() + + # check parameter override + processor = PaganinProcessor(delta=delta, beta=beta, energy=energy, return_units='m') + processor.set_input(self.data_cone) + processor._set_geometry(self.data_cone.geometry) + delta = 100 + beta=200 + processor.filter_Nx = Nx + processor.filter_Ny = Ny + processor._create_filter(override_filter={'delta':delta, 'beta':beta}) + + # check alpha and mu are calculated correctly + wavelength = (constants.h*constants.speed_of_light)/(energy*constants.electron_volt) + mu = 4.0*numpy.pi*beta/(wavelength) + alpha = 60000*delta/mu + self.assertEqual(processor.delta, delta, msg=self.error_message(processor, 'delta')) + self.assertEqual(processor.beta, beta, msg=self.error_message(processor, 'beta')) + self.assertEqual(processor.alpha, alpha, msg=self.error_message(processor, 'alpha')) + self.assertEqual(processor.mu, mu, msg=self.error_message(processor, 'mu')) + filter = ifftshift(1/(1. + alpha*(kx**2 + ky**2))) + numpy.testing.assert_allclose(processor.filter, filter) + + # test specifying alpha, delta and beta + delta = 12 + beta = 13 + alpha = 14 + processor.filter_Nx = Nx + processor.filter_Ny = Ny + with self.assertLogs(level='WARN') as log: + processor._create_filter(override_filter = {'delta':delta, 'beta':beta, 'alpha':alpha}) + wavelength = (constants.h*constants.speed_of_light)/(energy*constants.electron_volt) + mu = 4.0*numpy.pi*beta/(wavelength) + + self.assertEqual(processor.delta, delta, msg=self.error_message(processor, 'delta')) + self.assertEqual(processor.beta, beta, msg=self.error_message(processor, 'beta')) + self.assertEqual(processor.alpha, alpha, msg=self.error_message(processor, 'alpha')) + self.assertEqual(processor.mu, mu, msg=self.error_message(processor, 'mu')) + filter = ifftshift(1/(1. + alpha*(kx**2 + ky**2))) + numpy.testing.assert_allclose(processor.filter, filter) + + def test_PaganinProcessor(self): + + wavelength = (constants.h*constants.speed_of_light)/(40000*constants.electron_volt) + mu = 4.0*numpy.pi*1e-2/(wavelength) + + data_array = [self.data_cone, self.data_parallel, self.data_multichannel] + for data in data_array: + data.geometry.config.units = 'm' + data_abs = -(1/mu)*numpy.log(data) + processor = PaganinProcessor(full_retrieval=True) + processor.set_input(data) + thickness = processor.get_output(override_geometry={'propagation_distance':1}) + self.assertLessEqual(quality_measures.mse(thickness, data_abs), 1e-5) + processor = PaganinProcessor(full_retrieval=False) + processor.set_input(data) + filtered_image = processor.get_output(override_geometry={'propagation_distance':1}) + self.assertLessEqual(quality_measures.mse(filtered_image, data), 1e-5) + + # test with GPM + processor = PaganinProcessor(full_retrieval=True, filter_type='generalised_paganin_method') + processor.set_input(data) + thickness = processor.get_output(override_geometry={'propagation_distance':1}) + self.assertLessEqual(quality_measures.mse(thickness, data_abs), 1e-5) + processor = PaganinProcessor(full_retrieval=False, filter_type='generalised_paganin_method') + processor.set_input(data) + filtered_image = processor.get_output(override_geometry={'propagation_distance':1}) + self.assertLessEqual(quality_measures.mse(filtered_image, data), 1e-5) + + # test with padding + processor = PaganinProcessor(full_retrieval=True, pad=10) + processor.set_input(data) + 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) + filtered_image = processor.get_output(override_geometry={'propagation_distance':1}) + self.assertLessEqual(quality_measures.mse(filtered_image, data), 1e-5) + + # test in-line + thickness_inline = PaganinProcessor(full_retrieval=True, pad=10)(data, override_geometry={'propagation_distance':1}) + 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) + 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) + 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) + + 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) + mu = 4.0*numpy.pi*1e-2/(wavelength) + thickness = -(1/mu)*numpy.log(data_slice) + + processor = PaganinProcessor(pad=10) + processor.set_input(data_slice) + output = processor.get_output(override_geometry={'propagation_distance':1}) + self.assertLessEqual(quality_measures.mse(output, thickness), 0.05) + + # check with different data order + data_slice.reorder(('horizontal','angle')) + 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) + + processor = PaganinProcessor(pad=10) + processor.set_input(data_slice) + output = processor.get_output(override_geometry={'propagation_distance':1}) + self.assertLessEqual(quality_measures.mse(output, thickness), 0.05) + + # 'horizontal, vertical, angles + if __name__ == "__main__": d = TestDataProcessor() diff --git a/docs/source/processors.rst b/docs/source/processors.rst index 052589a59..1381a065f 100644 --- a/docs/source/processors.rst +++ b/docs/source/processors.rst @@ -140,4 +140,12 @@ Ring Remover :inherited-members: set_input, get_output +Paganin Processor +----------------- + +.. autoclass:: cil.processors.PaganinProcessor + :exclude-members: check_input, get_input + :members: + :inherited-members: + :ref:`Return Home `