From c477096e1382872d7e97243ea31a6ac494c6eba9 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Thu, 25 Jul 2024 11:04:49 +0000 Subject: [PATCH 01/25] Add flux normaliser processor, init tests --- .../Python/cil/processors/FluxNormaliser.py | 79 +++++++++++++++++++ Wrappers/Python/cil/processors/__init__.py | 3 +- Wrappers/Python/test/test_DataProcessor.py | 34 +++++++- 3 files changed, 113 insertions(+), 3 deletions(-) create mode 100644 Wrappers/Python/cil/processors/FluxNormaliser.py diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py new file mode 100644 index 0000000000..c791ed832c --- /dev/null +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -0,0 +1,79 @@ +# 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, DataContainer, AcquisitionData,\ + AcquisitionGeometry, ImageGeometry, ImageData +import numpy + +class FluxNormaliser(Processor): + '''Flux normalisation based on float or region of interest + + This processor reads in a AcquisitionData and normalises it based on + a float or array of float values, or a region of interest. + + Input: AcquisitionData + Parameter: 2D projection with flat field (or stack) + 2D projection with dark field (or stack) + Output: AcquisitionDataSet + ''' + + def __init__(self, flux=None, roi=None, tolerance=1e-5): + kwargs = { + 'flux' : flux, + 'roi' : roi, + # very small number. Used when there is a division by zero + 'tolerance' : tolerance + } + super(FluxNormaliser, self).__init__(**kwargs) + + def check_input(self, dataset): + flux_size = (numpy.shape(self.flux)) + if len(flux_size) > 0: + data_size = numpy.shape(dataset.geometry.angles) + if data_size != flux_size: + raise ValueError("Flux must be a scalar or array with length \ + \n = data.geometry.angles, found {} and {}" + .format(flux_size, data_size)) + + return True + + def process(self, out=None): + + data = self.get_input() + + if out is None: + out = data.copy() + + flux_size = (numpy.shape(self.flux)) + + proj_axis = data.get_dimension_axis('angle') + slice_proj = [slice(None)]*len(data.shape) + slice_proj[proj_axis] = 0 + + f = self.flux + for i in range(numpy.shape(data)[proj_axis]): + if len(flux_size) > 0: + f = self.flux[i] + + slice_proj[proj_axis] = i + with numpy.errstate(divide='ignore', invalid='ignore'): + out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]/f + + out.array[ ~ numpy.isfinite( out.array )] = self.tolerance + + return out diff --git a/Wrappers/Python/cil/processors/__init__.py b/Wrappers/Python/cil/processors/__init__.py index 15ba249e74..c86bb50cb7 100644 --- a/Wrappers/Python/cil/processors/__init__.py +++ b/Wrappers/Python/cil/processors/__init__.py @@ -26,4 +26,5 @@ from .TransmissionAbsorptionConverter import TransmissionAbsorptionConverter from .Masker import Masker from .Padder import Padder -from .PaganinProcessor import PaganinProcessor \ No newline at end of file +from .PaganinProcessor import PaganinProcessor +from .FluxNormaliser import FluxNormaliser \ No newline at end of file diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index a821bba105..6a063596ca 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -30,7 +30,7 @@ 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, PaganinProcessor +from cil.processors import Slicer, Binner, MaskGenerator, Masker, Padder, PaganinProcessor, FluxNormaliser import gc from scipy import constants @@ -3072,7 +3072,37 @@ def test_PaganinProcessor_2D(self): self.assertLessEqual(quality_measures.mse(output, thickness), 0.05) # 'horizontal, vertical, angles - + +class TestFluxNormaliser(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_init(self): + # test default values are initialised + processor = FluxNormaliser() + test_parameter = ['flux','roi','tolerance'] + test_value = [None, None, 1e-5] + + 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 = FluxNormaliser(1,2,3) + test_value = [1, 2, 3] + 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])) + if __name__ == "__main__": d = TestDataProcessor() From 45cb0a5c9c34cda589778484c16a97d5bc898cb3 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Fri, 26 Jul 2024 09:16:35 +0000 Subject: [PATCH 02/25] Add roi check input logic --- Wrappers/Python/test/test_DataProcessor.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 6a063596ca..52c2ad7071 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3103,6 +3103,12 @@ def test_init(self): 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])) + def test_FluxNormaliser(self): + processor = FluxNormaliser(flux=10) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, self.data_cone.array/10) + if __name__ == "__main__": d = TestDataProcessor() From d9cc16749d9db5d05816d7f4fe359dd5330fbcaa Mon Sep 17 00:00:00 2001 From: hrobarts Date: Wed, 31 Jul 2024 14:53:15 +0000 Subject: [PATCH 03/25] Add roi functionality --- .../Python/cil/processors/FluxNormaliser.py | 75 ++++++++++++++++++- 1 file changed, 73 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index c791ed832c..9a49327eca 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -18,7 +18,11 @@ from cil.framework import Processor, DataContainer, AcquisitionData,\ AcquisitionGeometry, ImageGeometry, ImageData +from cil.utilities.display import show2D import numpy +import logging + +log = logging.getLogger(__name__) class FluxNormaliser(Processor): '''Flux normalisation based on float or region of interest @@ -35,13 +39,44 @@ class FluxNormaliser(Processor): def __init__(self, flux=None, roi=None, tolerance=1e-5): kwargs = { 'flux' : flux, - 'roi' : roi, + 'roi' : roi, + 'roi_slice' : None, # very small number. Used when there is a division by zero 'tolerance' : tolerance } super(FluxNormaliser, self).__init__(**kwargs) def check_input(self, dataset): + + if not issubclass(type(dataset), DataContainer): + raise TypeError("Expected DataContainer or subclass, found {}" + .format(type(dataset))) + if self.roi is not None: + if self.flux is not None: + raise ValueError("Please specify either flux or roi, not both") + else: + if isinstance(self.roi,dict): + slc = [slice(None)]*len(data.shape) + axes=[] + for r in self.roi: + if (r == 'horizontal' or r == 'vertical') and (r in data.dimension_labels): + if not all(isinstance(i, int) for i in self.roi[r]): + raise ValueError("roi values must be int, found {}" + .format(str(type(self.roi[r])))) + else: + ax = data.get_dimension_axis(r) + slc[ax] = slice(self.roi[r][0],self.roi[r][1]) + axes.append(ax) + else: + raise ValueError("roi must be 'horizontal' or 'vertical', and in data.dimension_labels, found {}" + .format(str(r))) + + self.flux = np.mean(data.array[tuple(slc)],axis=tuple(axes)) + + else: + raise TypeError("roi must be a dictionary, found {}" + .format(str(type(self.roi)))) + flux_size = (numpy.shape(self.flux)) if len(flux_size) > 0: data_size = numpy.shape(dataset.geometry.angles) @@ -52,8 +87,44 @@ def check_input(self, dataset): return True - def process(self, out=None): + def show_roi(self, angle='range'): + if angle=='range': + plot_slice_roi_angle(angle=0, ax=231) + plot_slice_roi_angle(angle=len(data.geometry.angles)//2, ax=232) + plot_slice_roi_angle(angle=len(data.geometry.angles)-1, ax=233) + + plt.subplot(212) + plt.plot(self.flux) + plt.xlabel('Angle index') + plt.ylabel('Mean intensity in roi') + + else: + data = self.get_input() + plt.imshow(data.array[1,:,:], cmap='gray') + plt.plot([0,120],[0,120],'--k') + + def plot_slice_roi_angle(angle, ax): + data_slice = data.get_slice(angle=angle) + plt.subplot(ax) + plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower') + h = data_slice.dimension_labels[0] + v = data_slice.dimension_labels[1] + plt.plot([self.roi[h][0],self.roi[h][0]], [self.roi[v][0],self.roi[v][1]],'--r') + plt.plot([self.roi[h][1],self.roi[h][1]], [self.roi[v][0],self.roi[v][1]],'--r') + + plt.plot([self.roi[h][0],self.roi[h][1]], [self.roi[v][0],self.roi[v][0]],'--r') + plt.plot([self.roi[h][0],self.roi[h][1]], [self.roi[v][1],self.roi[v][1]],'--r') + + plt.xlabel(h) + plt.ylabel(v) + plt.title('Angle = ' + str(angle)) + + + + + def process(self, out=None): + data = self.get_input() if out is None: From fd53e4d71111e383e6dc2b405e0289c5d5600ef9 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Mon, 5 Aug 2024 09:29:20 +0000 Subject: [PATCH 04/25] Add show_roi and deal with no roi specified --- .../Python/cil/processors/FluxNormaliser.py | 98 +++++++++++++------ Wrappers/Python/test/test_DataProcessor.py | 49 ++++++++++ 2 files changed, 115 insertions(+), 32 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 9a49327eca..813e262235 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -21,6 +21,8 @@ from cil.utilities.display import show2D import numpy import logging +import matplotlib.pyplot as plt +from scipy import stats log = logging.getLogger(__name__) @@ -51,31 +53,58 @@ def check_input(self, dataset): if not issubclass(type(dataset), DataContainer): raise TypeError("Expected DataContainer or subclass, found {}" .format(type(dataset))) + if self.roi is not None: if self.flux is not None: raise ValueError("Please specify either flux or roi, not both") + else: if isinstance(self.roi,dict): - slc = [slice(None)]*len(data.shape) + if not all (r in dataset.dimension_labels for r in self.roi): + raise ValueError("roi must be 'horizontal' or 'vertical', found '{}'" + .format(str(self.roi))) + + slc = [slice(None)]*len(dataset.shape) axes=[] - for r in self.roi: - if (r == 'horizontal' or r == 'vertical') and (r in data.dimension_labels): - if not all(isinstance(i, int) for i in self.roi[r]): - raise ValueError("roi values must be int, found {}" - .format(str(type(self.roi[r])))) + roi_list = list(dataset.dimension_labels) + + # loop through all dimensions in the dataset apart from angle + roi_list.remove('angle') + for r in roi_list: + # check the dimension is in the user specified roi + if r in self.roi: + # only allow horizontal and vertical roi + if (r == 'horizontal' or r == 'vertical'): + # check indices are ints + if not all(isinstance(i, int) for i in self.roi[r]): + raise TypeError("roi values must be int, found {} and {}" + .format(str(type(self.roi[r][0])), str(type(self.roi[r][1])))) + # check indices are in range + elif (self.roi[r][0] >= self.roi[r][1]) or (self.roi[r][0] < 0) or self.roi[r][1] > dataset.get_dimension_size(r): + raise ValueError("roi values must be start > stop and between 0 and {}, found start={} and stop={} for direction '{}'" + .format(str(dataset.get_dimension_size(r)), str(self.roi[r][0]), str(self.roi[r][1]), r )) + # create slice + else: + ax = dataset.get_dimension_axis(r) + slc[ax] = slice(self.roi[r][0], self.roi[r][1]) + axes.append(ax) else: - ax = data.get_dimension_axis(r) - slc[ax] = slice(self.roi[r][0],self.roi[r][1]) - axes.append(ax) + raise ValueError("roi must be 'horizontal' or 'vertical', found '{}'" + .format(str(r))) + # if the dimension is not in roi, average across the whole dimension else: - raise ValueError("roi must be 'horizontal' or 'vertical', and in data.dimension_labels, found {}" - .format(str(r))) + ax = dataset.get_dimension_axis(r) + axes.append(ax) + self.roi.update({r:(0,dataset.get_dimension_size(r))}) - self.flux = np.mean(data.array[tuple(slc)],axis=tuple(axes)) + self.flux = numpy.mean(dataset.array[tuple(slc)], axis=tuple(axes)) else: raise TypeError("roi must be a dictionary, found {}" .format(str(type(self.roi)))) + else: + if self.flux is None: + raise ValueError("Please specify flux or roi") flux_size = (numpy.shape(self.flux)) if len(flux_size) > 0: @@ -87,41 +116,46 @@ def check_input(self, dataset): return True - def show_roi(self, angle='range'): - if angle=='range': - plot_slice_roi_angle(angle=0, ax=231) - plot_slice_roi_angle(angle=len(data.geometry.angles)//2, ax=232) - plot_slice_roi_angle(angle=len(data.geometry.angles)-1, ax=233) + def show_roi(self, angle_index='range'): + if angle_index=='range': + data = self.get_input() + max_zscore = numpy.argmax(numpy.abs(stats.zscore(self.flux))) + self._plot_slice_roi_angle(angle_index=0, ax=231) + self._plot_slice_roi_angle(angle_index=max_zscore, ax=232) + self._plot_slice_roi_angle(angle_index=len(data.geometry.angles)-1, ax=233) plt.subplot(212) plt.plot(self.flux) + plt.plot(max_zscore, self.flux[max_zscore],'rx') + + plt.legend(['Mean intensity in roi','Max Z-score']) plt.xlabel('Angle index') plt.ylabel('Mean intensity in roi') + plt.tight_layout() + + elif isinstance(angle_index, int): + self._plot_slice_roi_angle(angle_index=angle_index) else: - data = self.get_input() - plt.imshow(data.array[1,:,:], cmap='gray') - plt.plot([0,120],[0,120],'--k') + raise TypeError("angle must be an int or 'range'") - def plot_slice_roi_angle(angle, ax): - data_slice = data.get_slice(angle=angle) + def _plot_slice_roi_angle(self, angle_index, ax=111): + data = self.get_input() + data_slice = data.get_slice(angle=angle_index) plt.subplot(ax) plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower') h = data_slice.dimension_labels[0] v = data_slice.dimension_labels[1] - plt.plot([self.roi[h][0],self.roi[h][0]], [self.roi[v][0],self.roi[v][1]],'--r') - plt.plot([self.roi[h][1],self.roi[h][1]], [self.roi[v][0],self.roi[v][1]],'--r') + plt.plot([self.roi[v][0],self.roi[v][0]], [self.roi[h][0],self.roi[h][1]],'--r') + plt.plot([self.roi[v][1],self.roi[v][1]], [self.roi[h][0],self.roi[h][1]],'--r') - plt.plot([self.roi[h][0],self.roi[h][1]], [self.roi[v][0],self.roi[v][0]],'--r') - plt.plot([self.roi[h][0],self.roi[h][1]], [self.roi[v][1],self.roi[v][1]],'--r') - - plt.xlabel(h) - plt.ylabel(v) - plt.title('Angle = ' + str(angle)) - - + plt.plot([self.roi[v][0],self.roi[v][1]], [self.roi[h][0],self.roi[h][0]],'--r') + plt.plot([self.roi[v][0],self.roi[v][1]], [self.roi[h][1],self.roi[h][1]],'--r') + plt.xlabel(v) + plt.ylabel(h) + plt.title('Angle index = ' + str(angle_index)) def process(self, out=None): diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 52c2ad7071..84c4967c52 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3103,11 +3103,60 @@ def test_init(self): 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])) + def test_check_input(self): + + # check there is an error if no flux or roi is specified + processor = FluxNormaliser() + with self.assertRaises(ValueError): + processor.check_input(self.data_cone) + + # check there is an error if no flux array size is not equal to the number of angles in data + processor = FluxNormaliser(flux = [1,2,3]) + with self.assertRaises(ValueError): + processor.check_input(self.data_cone) + + # check there is an error if roi is not specified as a dictionary + processor = FluxNormaliser(roi='string') + with self.assertRaises(TypeError): + processor.check_input(self.data_cone) + + # check there is an error if roi is specified with float values + processor = FluxNormaliser(roi={'horizontal':(1.5, 6.5)}) + with self.assertRaises(TypeError): + processor.check_input(self.data_cone) + + # check there is an error if roi stop is greater than start + processor = FluxNormaliser(roi={'horizontal':(10, 5)}) + with self.assertRaises(ValueError): + processor.check_input(self.data_cone) + + # check there is an error if roi stop is greater than the size of the axis + processor = FluxNormaliser(roi={'horizontal':(0, self.data_cone.get_dimension_size('horizontal')+1)}) + with self.assertRaises(ValueError): + processor.check_input(self.data_cone) + def test_FluxNormaliser(self): processor = FluxNormaliser(flux=10) processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, self.data_cone.array/10) + + processor = FluxNormaliser(flux=10*numpy.ones(self.data_cone.get_dimension_size('angle'))) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, self.data_cone.array/10) + + roi = {'vertical':(0,10), 'horizontal':(0,10)} + processor = FluxNormaliser(roi=roi) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) + + roi = {'vertical':(0,2)} + processor = FluxNormaliser(roi=roi) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) if __name__ == "__main__": From 36cfa783c43c90e4390025ff7519b51dcec15711 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Wed, 14 Aug 2024 16:31:30 +0000 Subject: [PATCH 05/25] Preview min and max, log plot, plot angle --- .../Python/cil/processors/FluxNormaliser.py | 100 ++++++++++++------ 1 file changed, 70 insertions(+), 30 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 813e262235..409e081359 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -27,15 +27,28 @@ log = logging.getLogger(__name__) class FluxNormaliser(Processor): - '''Flux normalisation based on float or region of interest + ''' + Flux normalisation based on float or region of interest This processor reads in a AcquisitionData and normalises it based on a float or array of float values, or a region of interest. - Input: AcquisitionData - Parameter: 2D projection with flat field (or stack) - 2D projection with dark field (or stack) - Output: AcquisitionDataSet + Parameters: + ----------- + flux: float or list of floats (optional) + Divide the image by the flux value. If flux is a list it must have length + equal to the number of angles in the dataset. + + roi: dict (optional) + Dictionary describing the region of interest containing the background + in the image. The image is divided by the mean value in the roi. + + tolerance: float (optional) + Small number to set to when there is a division by zero, default is 1e-5. + + Returns: + -------- + Output: AcquisitionData normalised by flux or mean intensity in roi ''' def __init__(self, flux=None, roi=None, tolerance=1e-5): @@ -43,7 +56,7 @@ def __init__(self, flux=None, roi=None, tolerance=1e-5): 'flux' : flux, 'roi' : roi, 'roi_slice' : None, - # very small number. Used when there is a division by zero + 'roi_axes' : None, 'tolerance' : tolerance } super(FluxNormaliser, self).__init__(**kwargs) @@ -98,7 +111,9 @@ def check_input(self, dataset): self.roi.update({r:(0,dataset.get_dimension_size(r))}) self.flux = numpy.mean(dataset.array[tuple(slc)], axis=tuple(axes)) - + self.roi_slice = slc + self.roi_axes = axes + else: raise TypeError("roi must be a dictionary, found {}" .format(str(type(self.roi)))) @@ -116,34 +131,59 @@ def check_input(self, dataset): return True - def show_roi(self, angle_index='range'): - if angle_index=='range': - data = self.get_input() - max_zscore = numpy.argmax(numpy.abs(stats.zscore(self.flux))) - self._plot_slice_roi_angle(angle_index=0, ax=231) - self._plot_slice_roi_angle(angle_index=max_zscore, ax=232) - self._plot_slice_roi_angle(angle_index=len(data.geometry.angles)-1, ax=233) - - plt.subplot(212) - plt.plot(self.flux) - plt.plot(max_zscore, self.flux[max_zscore],'rx') + def preview_configuration(self, angle='min_and_max', log=False): + ''' + Preview the FluxNormalisation processor configuration for roi mode. + Plots the region of interest on the image and the mean, maximum and + minimum intensity in the roi. the angles with minimum and maximum intensity in the roi, - plt.legend(['Mean intensity in roi','Max Z-score']) - plt.xlabel('Angle index') - plt.ylabel('Mean intensity in roi') - plt.tight_layout() - - elif isinstance(angle_index, int): - self._plot_slice_roi_angle(angle_index=angle_index) - + Parameters: + ----------- + angle: float or str (optional) + Angle to plot, default='min_and_max' displays the data with the minimum + and maximum pixel values in the roi, otherwise the angle to display + can be specified as a float (closest angle is displayed). + + log: bool (optional) + If True, plot the image with a log scale, default is False + ''' + if self.roi_slice is None: + raise ValueError('Preview available with roi, run `processor= FluxNormaliser(roi=roi)` then `set_input(data)`') else: - raise TypeError("angle must be an int or 'range'") + data = self.get_input() + min = numpy.min(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) + max = numpy.max(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) + + if angle=='min_and_max': + self._plot_slice_roi_angle(angle_index=numpy.argmin(min), log=log, ax=221) + self._plot_slice_roi_angle(angle_index=numpy.argmax(max), log=log, ax=222) + + plt.subplot(212) + + else: + plt.figure(figsize=(3,6)) + angle_index = numpy.argmin(numpy.abs(angle-data.geometry.angles)) + self._plot_slice_roi_angle(angle_index=angle_index, log=log, ax=211) - def _plot_slice_roi_angle(self, angle_index, ax=111): + plt.subplot(212) + + plt.plot(data.geometry.angles, self.flux, 'r', label='Mean') + plt.plot(data.geometry.angles, min,'--k', label='Minimum') + plt.plot(data.geometry.angles, max,'--k', label='Maximum') + plt.legend() + plt.xlabel('Angle') + plt.ylabel('Intensity in roi') + plt.grid() + plt.tight_layout() + + def _plot_slice_roi_angle(self, angle_index, log=False, ax=111): data = self.get_input() data_slice = data.get_slice(angle=angle_index) plt.subplot(ax) - plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower') + if log: + plt.imshow(numpy.log(data_slice.array), cmap='gray',aspect='equal', origin='lower') + else: + plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower') h = data_slice.dimension_labels[0] v = data_slice.dimension_labels[1] @@ -155,7 +195,7 @@ def _plot_slice_roi_angle(self, angle_index, ax=111): plt.xlabel(v) plt.ylabel(h) - plt.title('Angle index = ' + str(angle_index)) + plt.title('Angle = ' + str(data.geometry.angles[angle_index])) def process(self, out=None): From 07964e9bd1c0e81bfe20d4ef4b163fb1f6c56ad7 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Mon, 19 Aug 2024 13:01:40 +0000 Subject: [PATCH 06/25] Update preview_configuration() to deal with different data shapes --- .../Python/cil/processors/FluxNormaliser.py | 178 ++++++++++++------ Wrappers/Python/test/test_DataProcessor.py | 9 + 2 files changed, 133 insertions(+), 54 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 409e081359..4386699c5d 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -74,41 +74,45 @@ def check_input(self, dataset): else: if isinstance(self.roi,dict): if not all (r in dataset.dimension_labels for r in self.roi): - raise ValueError("roi must be 'horizontal' or 'vertical', found '{}'" + raise ValueError("roi labels must be in the dataset dimension_labels, found {}" .format(str(self.roi))) slc = [slice(None)]*len(dataset.shape) axes=[] - roi_list = list(dataset.dimension_labels) + + for r in self.roi: + # only allow roi to be specified in horizontal and vertical + if (r != 'horizontal' and r != 'vertical'): + raise ValueError("roi must be 'horizontal' or 'vertical', found '{}'" + .format(str(r))) + dimension_label_list = list(dataset.dimension_labels) # loop through all dimensions in the dataset apart from angle - roi_list.remove('angle') - for r in roi_list: + if 'angle' in dimension_label_list: + dimension_label_list.remove('angle') + + for d in dimension_label_list: # check the dimension is in the user specified roi - if r in self.roi: - # only allow horizontal and vertical roi - if (r == 'horizontal' or r == 'vertical'): - # check indices are ints - if not all(isinstance(i, int) for i in self.roi[r]): - raise TypeError("roi values must be int, found {} and {}" - .format(str(type(self.roi[r][0])), str(type(self.roi[r][1])))) - # check indices are in range - elif (self.roi[r][0] >= self.roi[r][1]) or (self.roi[r][0] < 0) or self.roi[r][1] > dataset.get_dimension_size(r): - raise ValueError("roi values must be start > stop and between 0 and {}, found start={} and stop={} for direction '{}'" - .format(str(dataset.get_dimension_size(r)), str(self.roi[r][0]), str(self.roi[r][1]), r )) - # create slice - else: - ax = dataset.get_dimension_axis(r) - slc[ax] = slice(self.roi[r][0], self.roi[r][1]) - axes.append(ax) + if d in self.roi: + # check indices are ints + if not all(isinstance(i, int) for i in self.roi[d]): + raise TypeError("roi values must be int, found {} and {}" + .format(str(type(self.roi[d][0])), str(type(self.roi[d][1])))) + # check indices are in range + elif (self.roi[d][0] >= self.roi[d][1]) or (self.roi[d][0] < 0) or self.roi[d][1] > dataset.get_dimension_size(d): + raise ValueError("roi values must be start > stop and between 0 and {}, found start={} and stop={} for direction '{}'" + .format(str(dataset.get_dimension_size(d)), str(self.roi[d][0]), str(self.roi[d][1]), d )) + # create slice else: - raise ValueError("roi must be 'horizontal' or 'vertical', found '{}'" - .format(str(r))) + ax = dataset.get_dimension_axis(d) + slc[ax] = slice(self.roi[d][0], self.roi[d][1]) + axes.append(ax) + # if the dimension is not in roi, average across the whole dimension else: - ax = dataset.get_dimension_axis(r) + ax = dataset.get_dimension_axis(d) axes.append(ax) - self.roi.update({r:(0,dataset.get_dimension_size(r))}) + self.roi.update({d:(0,dataset.get_dimension_size(d))}) self.flux = numpy.mean(dataset.array[tuple(slc)], axis=tuple(axes)) self.roi_slice = slc @@ -131,7 +135,7 @@ def check_input(self, dataset): return True - def preview_configuration(self, angle='min_and_max', log=False): + def preview_configuration(self, angle=None, channel=None, log=False): ''' Preview the FluxNormalisation processor configuration for roi mode. Plots the region of interest on the image and the mean, maximum and @@ -140,9 +144,14 @@ def preview_configuration(self, angle='min_and_max', log=False): Parameters: ----------- angle: float or str (optional) - Angle to plot, default='min_and_max' displays the data with the minimum + Angle to plot, default=None displays the data with the minimum and maximum pixel values in the roi, otherwise the angle to display - can be specified as a float (closest angle is displayed). + can be specified as a float and the closest angle will be displayed. + For 2D data, the roi is plotted on the sinogram. + + channel: int (optional) + The channel to plot, default=None displays the central channel if + the data has channels log: bool (optional) If True, plot the image with a log scale, default is False @@ -151,51 +160,112 @@ def preview_configuration(self, angle='min_and_max', log=False): raise ValueError('Preview available with roi, run `processor= FluxNormaliser(roi=roi)` then `set_input(data)`') else: data = self.get_input() + min = numpy.min(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) max = numpy.max(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) - if angle=='min_and_max': - self._plot_slice_roi_angle(angle_index=numpy.argmin(min), log=log, ax=221) - self._plot_slice_roi_angle(angle_index=numpy.argmax(max), log=log, ax=222) + if data.geometry.dimension == '3D': + if angle is None: + if 'angle' 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: + angle_index = numpy.argmin(numpy.abs(angle-data.geometry.angles)) + self._plot_slice_roi(angle_index=angle_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: + 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") - plt.subplot(212) - + plt.subplot(212) + if len(data.geometry.angles)==1: + plt.plot(data.geometry.angles, self.flux, '.r', label='Mean') + plt.plot(data.geometry.angles, min,'.k', label='Minimum') + plt.plot(data.geometry.angles, max,'.k', label='Maximum') else: - plt.figure(figsize=(3,6)) - angle_index = numpy.argmin(numpy.abs(angle-data.geometry.angles)) - self._plot_slice_roi_angle(angle_index=angle_index, log=log, ax=211) - - plt.subplot(212) - - plt.plot(data.geometry.angles, self.flux, 'r', label='Mean') - plt.plot(data.geometry.angles, min,'--k', label='Minimum') - plt.plot(data.geometry.angles, max,'--k', label='Maximum') + plt.plot(data.geometry.angles, self.flux, 'r', label='Mean') + plt.plot(data.geometry.angles, min,'--k', label='Minimum') + plt.plot(data.geometry.angles, max,'--k', label='Maximum') plt.legend() plt.xlabel('Angle') plt.ylabel('Intensity in roi') plt.grid() plt.tight_layout() - def _plot_slice_roi_angle(self, angle_index, log=False, ax=111): + def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=111): + data = self.get_input() - data_slice = data.get_slice(angle=angle_index) + if angle_index is not None and 'angle' in data.dimension_labels: + data_slice = data.get_slice(angle=angle_index) + else: + data_slice = data + + if 'channel' in data.dimension_labels: + if channel_index is None: + channel_index = int(numpy.round(data_slice.get_dimension_size('channel')/2)) + data_slice = data_slice.get_slice(channel=channel_index) + else: + if channel_index is not None: + raise ValueError("Channel not found") + + 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'") + + extent = [0, data_slice.shape[1], 0, data_slice.shape[0]] + if 'angle' in data_slice.dimension_labels: + min_angle = data_slice.geometry.angles[0] + max_angle = data_slice.geometry.angles[-1] + for i, d in enumerate(data_slice.dimension_labels): + if d !='angle': + extent[i*2]=min_angle + extent[i*2+1]=max_angle + plt.subplot(ax) if log: - plt.imshow(numpy.log(data_slice.array), cmap='gray',aspect='equal', origin='lower') + plt.imshow(numpy.log(data_slice.array), cmap='gray',aspect='equal', origin='lower', extent=extent) else: - plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower') + plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower', extent=extent) - h = data_slice.dimension_labels[0] - v = data_slice.dimension_labels[1] - plt.plot([self.roi[v][0],self.roi[v][0]], [self.roi[h][0],self.roi[h][1]],'--r') - plt.plot([self.roi[v][1],self.roi[v][1]], [self.roi[h][0],self.roi[h][1]],'--r') + h = data_slice.dimension_labels[1] + v = data_slice.dimension_labels[0] - plt.plot([self.roi[v][0],self.roi[v][1]], [self.roi[h][0],self.roi[h][0]],'--r') - plt.plot([self.roi[v][0],self.roi[v][1]], [self.roi[h][1],self.roi[h][1]],'--r') + if h == 'angle': + h_min = min_angle + h_max = max_angle + else: + h_min = self.roi[h][0] + h_max = self.roi[h][1] + + if v == 'angle': + v_min = min_angle + v_max = max_angle + else: + v_min = self.roi[v][0] + v_max = self.roi[v][1] + + plt.plot([h_min, h_max],[v_min, v_min],'--r') + plt.plot([h_min, h_max],[v_max, v_max],'--r') + + plt.plot([h_min, h_min],[v_min, v_max],'--r') + plt.plot([h_max, h_max],[v_min, v_max],'--r') + + title = 'ROI' + if angle_index is not None: + title += ' angle = ' + str(data.geometry.angles[angle_index]) + if channel_index is not None: + title += ' channel = ' + str(channel_index) + plt.title(title) - plt.xlabel(v) - plt.ylabel(h) - plt.title('Angle = ' + str(data.geometry.angles[angle_index])) + plt.xlabel(h) + plt.ylabel(v) def process(self, out=None): diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 84c4967c52..5bb1add3e1 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3135,6 +3135,15 @@ def test_check_input(self): with self.assertRaises(ValueError): processor.check_input(self.data_cone) + def test_preview_configuration(self): + # Test error in preview configuration if there is no roi + processor = FluxNormaliser(flux=10) + processor.set_input(self.data_cone) + with self.assertRaises(ValueError): + processor.preview_configuration() + + + def test_FluxNormaliser(self): processor = FluxNormaliser(flux=10) processor.set_input(self.data_cone) From e724a3a20740b39aa825025af50e4671f5d68c0b Mon Sep 17 00:00:00 2001 From: hrobarts Date: Mon, 19 Aug 2024 13:43:55 +0000 Subject: [PATCH 07/25] Make process work with one angle --- .../Python/cil/processors/FluxNormaliser.py | 48 +++++++++++++------ 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 4386699c5d..44a397f570 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -32,17 +32,28 @@ class FluxNormaliser(Processor): This processor reads in a AcquisitionData and normalises it based on a float or array of float values, or a region of interest. + + The normalised image :math:`I_{norm}` is calculated from the original image + :math:`I` by + I_{norm} = I*(norm_value/flux) + Parameters: ----------- flux: float or list of floats (optional) - Divide the image by the flux value. If flux is a list it must have length - equal to the number of angles in the dataset. + The value to divide the image by. If flux is a list it must have length + equal to the number of angles in the dataset. If flux=None, calculate + flux from the roi, default is None. roi: dict (optional) Dictionary describing the region of interest containing the background - in the image. The image is divided by the mean value in the roi. + in the image. The image is divided by the mean value in the roi. If None, + specify flux direcrtly, default is None. + norm_value: float (optional) + The value to multiply the image by. If None, the mean flux value will + be used, default is None. + tolerance: float (optional) Small number to set to when there is a division by zero, default is 1e-5. @@ -51,12 +62,13 @@ class FluxNormaliser(Processor): Output: AcquisitionData normalised by flux or mean intensity in roi ''' - def __init__(self, flux=None, roi=None, tolerance=1e-5): + def __init__(self, flux=None, roi=None, norm_value=None, tolerance=1e-5): kwargs = { 'flux' : flux, 'roi' : roi, 'roi_slice' : None, 'roi_axes' : None, + 'norm_value' : norm_value, 'tolerance' : tolerance } super(FluxNormaliser, self).__init__(**kwargs) @@ -133,6 +145,9 @@ def check_input(self, dataset): \n = data.geometry.angles, found {} and {}" .format(flux_size, data_size)) + if self.norm_value is None: + self.norm_value = numpy.mean(self.flux) + return True def preview_configuration(self, angle=None, channel=None, log=False): @@ -275,20 +290,23 @@ def process(self, out=None): out = data.copy() flux_size = (numpy.shape(self.flux)) + f = self.flux - proj_axis = data.get_dimension_axis('angle') - slice_proj = [slice(None)]*len(data.shape) - slice_proj[proj_axis] = 0 + if 'angle' in data.dimension_labels: + proj_axis = data.get_dimension_axis('angle') + slice_proj = [slice(None)]*len(data.shape) + slice_proj[proj_axis] = 0 - f = self.flux - for i in range(numpy.shape(data)[proj_axis]): - if len(flux_size) > 0: - f = self.flux[i] + for i in range(len(data.geometry.angles)): + if len(flux_size) > 0: + f = self.flux[i] + slice_proj[proj_axis] = i + with numpy.errstate(divide='ignore', invalid='ignore'): + out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]*self.norm_value/f + else: + with numpy.errstate(divide='ignore', invalid='ignore'): + out.array = data.array*self.norm_value/f - slice_proj[proj_axis] = i - with numpy.errstate(divide='ignore', invalid='ignore'): - out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]/f - out.array[ ~ numpy.isfinite( out.array )] = self.tolerance return out From b4f516c56385e4ae4dea5efc14fb78a6c93e45b2 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Mon, 19 Aug 2024 16:32:09 +0000 Subject: [PATCH 08/25] Update tests --- Wrappers/Python/test/test_DataProcessor.py | 83 +++++++++++++++++++--- 1 file changed, 74 insertions(+), 9 deletions(-) diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 5bb1add3e1..6160f308ac 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3084,6 +3084,10 @@ def setUp(self): .set_channels(4) self.data_multichannel = ag.allocate('random') + self.data_slice = self.data_parallel.get_slice(vertical=1) + self.data_reorder = self.data_cone + self.data_reorder.reorder(['horizontal','angle','vertical']) + self.data_single_angle = self.data_cone.get_slice(angle=1) def error_message(self,processor, test_parameter): return "Failed with processor " + str(processor) + " on test parameter " + test_parameter @@ -3091,15 +3095,15 @@ def error_message(self,processor, test_parameter): def test_init(self): # test default values are initialised processor = FluxNormaliser() - test_parameter = ['flux','roi','tolerance'] - test_value = [None, None, 1e-5] + test_parameter = ['flux','roi','norm_value', 'tolerance'] + test_value = [None, None, None, 1e-5] 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 = FluxNormaliser(1,2,3) - test_value = [1, 2, 3] + processor = FluxNormaliser(1,2,3,4) + test_value = [1, 2, 3, 4] 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])) @@ -3142,30 +3146,91 @@ def test_preview_configuration(self): with self.assertRaises(ValueError): processor.preview_configuration() - + # Test error in preview configuration if set_input not called + roi = {'horizontal':(25,40)} + processor = FluxNormaliser(roi=roi) + with self.assertRaises(ValueError): + processor.preview_configuration() + + # 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]: + roi = {'horizontal':(25,40)} + processor = FluxNormaliser(roi=roi) + processor.set_input(data) + processor.preview_configuration() + + # for 3D, check no error specifying a single angle to plot + if data.geometry.dimension == '3D': + processor.preview_configuration(angle=1) + # if 2D, attempt to plot single angle should cause error + else: + with self.assertRaises(ValueError): + processor.preview_configuration(angle=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(channel=1) + # if single channel, check specifying channel causes an error + else: + with self.assertRaises(ValueError): + processor.preview_configuration(channel=1) + def test_FluxNormaliser(self): + #Test flux with no norm_value processor = FluxNormaliser(flux=10) processor.set_input(self.data_cone) data_norm = processor.get_output() - numpy.testing.assert_allclose(data_norm.array, self.data_cone.array/10) + numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) + + #Test flux with norm_value + processor = FluxNormaliser(flux=10, norm_value=5) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) + #Test flux array with no norm_value processor = FluxNormaliser(flux=10*numpy.ones(self.data_cone.get_dimension_size('angle'))) processor.set_input(self.data_cone) data_norm = processor.get_output() - numpy.testing.assert_allclose(data_norm.array, self.data_cone.array/10) + numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) + + #Test flux array with norm_value + processor = FluxNormaliser(flux=10*numpy.ones(self.data_cone.get_dimension_size('angle')), norm_value=5) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) + #Test roi with no norm_value roi = {'vertical':(0,10), 'horizontal':(0,10)} processor = FluxNormaliser(roi=roi) processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) + #Test roi with norm_value + roi = {'vertical':(0,10), 'horizontal':(0,10)} + processor = FluxNormaliser(roi=roi, norm_value=5) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) + + # Test roi with just one dimension roi = {'vertical':(0,2)} - processor = FluxNormaliser(roi=roi) + processor = FluxNormaliser(roi=roi, norm_value=5) processor.set_input(self.data_cone) data_norm = processor.get_output() - numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) + numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) + + # for data in [self.data_cone, self.data_parallel, self.data_multichannel, + # self.data_slice, self.data_reorder, self.data_single_angle]: + # roi = {'horizontal':(25,40)} + # processor = FluxNormaliser(roi=roi) + # processor.set_input(data) + # data_norm = processor.get_output() + # numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) if __name__ == "__main__": From c4b39cce3eb057379589eaff43c7cd75a218fb19 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 20 Aug 2024 09:57:18 +0000 Subject: [PATCH 09/25] FluxNormaliser unit tests --- Wrappers/Python/test/test_DataProcessor.py | 49 ++++++++++++++++++---- 1 file changed, 41 insertions(+), 8 deletions(-) diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 6160f308ac..4003e2e1e0 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3223,14 +3223,47 @@ def test_FluxNormaliser(self): processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) - - # for data in [self.data_cone, self.data_parallel, self.data_multichannel, - # self.data_slice, self.data_reorder, self.data_single_angle]: - # roi = {'horizontal':(25,40)} - # processor = FluxNormaliser(roi=roi) - # processor.set_input(data) - # data_norm = processor.get_output() - # numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) + def test_FluxNormaliser2(self): + # test roi with different data shapes + for data in [self.data_cone, self.data_parallel, self.data_multichannel, + self.data_slice, self.data_reorder]: + roi = {'horizontal':(25,40)} + processor = FluxNormaliser(roi=roi, norm_value=5) + processor.set_input(data) + data_norm = processor.get_output() + + ax = data.get_dimension_axis('horizontal') + slc = [slice(None)]*len(data.shape) + slc[ax] = slice(25,40) + axes=[ax] + if 'vertical' in data.dimension_labels: + axes.append(data.get_dimension_axis('vertical')) + if 'channel' in data.dimension_labels: + axes.append(data.get_dimension_axis('channel')) + flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) + slice_proj = [slice(None)]*len(data.shape) + proj_axis = data.get_dimension_axis('angle') + data_norm_test = data.copy() + for i in range(len(data.geometry.angles)): + f = flux[i] + slice_proj[proj_axis] = i + with numpy.errstate(divide='ignore', invalid='ignore'): + data_norm_test.array[tuple(slice_proj)] = 5/f*data.array[tuple(slice_proj)] + numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6, + err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) + + data = self.data_single_angle + processor = FluxNormaliser(roi=roi, norm_value=5) + processor.set_input(data) + data_norm = processor.get_output() + ax = data.get_dimension_axis('horizontal') + slc = [slice(None)]*len(data.shape) + slc[ax] = slice(25,40) + axes=[ax,data.get_dimension_axis('vertical')] + flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) + + numpy.testing.assert_allclose(data_norm.array, 5/flux*data.array, atol=1e-6, + err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) if __name__ == "__main__": From ece454f114e7534a0c10752b1525cde4f15c0867 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 20 Aug 2024 10:15:21 +0000 Subject: [PATCH 10/25] Docs update --- Wrappers/Python/cil/processors/FluxNormaliser.py | 6 ++++-- docs/source/processors.rst | 10 +++++++++- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 44a397f570..4a44af03d2 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -27,7 +27,7 @@ log = logging.getLogger(__name__) class FluxNormaliser(Processor): - ''' + r''' Flux normalisation based on float or region of interest This processor reads in a AcquisitionData and normalises it based on @@ -35,8 +35,10 @@ class FluxNormaliser(Processor): The normalised image :math:`I_{norm}` is calculated from the original image :math:`I` by - I_{norm} = I*(norm_value/flux) + .. math:: I_{norm} = I\frac{n}{F} + + where :math:`F` is the flux and :math:`n` is the norm_value Parameters: ----------- diff --git a/docs/source/processors.rst b/docs/source/processors.rst index 1381a065fc..bd80c074e0 100644 --- a/docs/source/processors.rst +++ b/docs/source/processors.rst @@ -146,6 +146,14 @@ Paganin Processor .. autoclass:: cil.processors.PaganinProcessor :exclude-members: check_input, get_input :members: - :inherited-members: + :inherited-members: set_input, get_output + +Flux Normaliser +----------------- + +.. autoclass:: cil.processors.FluxNormaliser + :exclude-members: check_input, get_input + :members: + :inherited-members: set_input, get_output :ref:`Return Home ` From 19f494aa269091ea35a12395474a4c71bcacdb8d Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 20 Aug 2024 13:11:02 +0000 Subject: [PATCH 11/25] Remove redundant dependency --- Wrappers/Python/cil/processors/FluxNormaliser.py | 1 - 1 file changed, 1 deletion(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 4a44af03d2..19ed31bd17 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -22,7 +22,6 @@ import numpy import logging import matplotlib.pyplot as plt -from scipy import stats log = logging.getLogger(__name__) From 608ffee0d5f447ad1a5becd663b55e4f80898549 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 20 Aug 2024 13:45:26 +0000 Subject: [PATCH 12/25] Remove plotting in tests --- .../Python/cil/processors/FluxNormaliser.py | 2 +- Wrappers/Python/test/test_DataProcessor.py | 74 +++++++++---------- 2 files changed, 38 insertions(+), 38 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 19ed31bd17..f9b4ea726e 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -226,7 +226,7 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11 if 'channel' in data.dimension_labels: if channel_index is None: - channel_index = int(numpy.round(data_slice.get_dimension_size('channel')/2)) + channel_index = int(data_slice.get_dimension_size('channel')/2) data_slice = data_slice.get_slice(channel=channel_index) else: if channel_index is not None: diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index c3a677bf50..ba5e76d761 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3142,43 +3142,43 @@ def test_check_input(self): with self.assertRaises(ValueError): processor.check_input(self.data_cone) - def test_preview_configuration(self): - # Test error in preview configuration if there is no roi - processor = FluxNormaliser(flux=10) - processor.set_input(self.data_cone) - with self.assertRaises(ValueError): - processor.preview_configuration() - - # Test error in preview configuration if set_input not called - roi = {'horizontal':(25,40)} - processor = FluxNormaliser(roi=roi) - with self.assertRaises(ValueError): - processor.preview_configuration() - - # 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]: - roi = {'horizontal':(25,40)} - processor = FluxNormaliser(roi=roi) - processor.set_input(data) - processor.preview_configuration() - - # for 3D, check no error specifying a single angle to plot - if data.geometry.dimension == '3D': - processor.preview_configuration(angle=1) - # if 2D, attempt to plot single angle should cause error - else: - with self.assertRaises(ValueError): - processor.preview_configuration(angle=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(channel=1) - # if single channel, check specifying channel causes an error - else: - with self.assertRaises(ValueError): - processor.preview_configuration(channel=1) + # def test_preview_configuration(self): + # # Test error in preview configuration if there is no roi + # processor = FluxNormaliser(flux=10) + # processor.set_input(self.data_cone) + # with self.assertRaises(ValueError): + # processor.preview_configuration() + + # # Test error in preview configuration if set_input not called + # roi = {'horizontal':(25,40)} + # processor = FluxNormaliser(roi=roi) + # with self.assertRaises(ValueError): + # processor.preview_configuration() + + # # 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]: + # roi = {'horizontal':(25,40)} + # processor = FluxNormaliser(roi=roi) + # processor.set_input(data) + # processor.preview_configuration() + + # # for 3D, check no error specifying a single angle to plot + # if data.geometry.dimension == '3D': + # processor.preview_configuration(angle=1) + # # if 2D, attempt to plot single angle should cause error + # else: + # with self.assertRaises(ValueError): + # processor.preview_configuration(angle=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(channel=1) + # # if single channel, check specifying channel causes an error + # else: + # with self.assertRaises(ValueError): + # processor.preview_configuration(channel=1) def test_FluxNormaliser(self): From 2985e6214bc74c4293fc6b74e65efbaafbf22a56 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 20 Aug 2024 14:09:45 +0000 Subject: [PATCH 13/25] Mock plot --- Wrappers/Python/test/test_DataProcessor.py | 77 +++++++++++----------- 1 file changed, 40 insertions(+), 37 deletions(-) diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index ba5e76d761..846f825026 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -18,6 +18,8 @@ import unittest import numpy +from unittest.mock import patch + from cil.framework import DataContainer from cil.framework import ImageGeometry, VectorGeometry, AcquisitionGeometry from cil.framework import ImageData, AcquisitionData @@ -3142,43 +3144,44 @@ def test_check_input(self): with self.assertRaises(ValueError): processor.check_input(self.data_cone) - # def test_preview_configuration(self): - # # Test error in preview configuration if there is no roi - # processor = FluxNormaliser(flux=10) - # processor.set_input(self.data_cone) - # with self.assertRaises(ValueError): - # processor.preview_configuration() - - # # Test error in preview configuration if set_input not called - # roi = {'horizontal':(25,40)} - # processor = FluxNormaliser(roi=roi) - # with self.assertRaises(ValueError): - # processor.preview_configuration() - - # # 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]: - # roi = {'horizontal':(25,40)} - # processor = FluxNormaliser(roi=roi) - # processor.set_input(data) - # processor.preview_configuration() - - # # for 3D, check no error specifying a single angle to plot - # if data.geometry.dimension == '3D': - # processor.preview_configuration(angle=1) - # # if 2D, attempt to plot single angle should cause error - # else: - # with self.assertRaises(ValueError): - # processor.preview_configuration(angle=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(channel=1) - # # if single channel, check specifying channel causes an error - # else: - # with self.assertRaises(ValueError): - # processor.preview_configuration(channel=1) + @patch("matplotlib.pyplot") + def test_preview_configuration(self, mock_plot): + # Test error in preview configuration if there is no roi + processor = FluxNormaliser(flux=10) + processor.set_input(self.data_cone) + with self.assertRaises(ValueError): + processor.preview_configuration() + + # Test error in preview configuration if set_input not called + roi = {'horizontal':(25,40)} + processor = FluxNormaliser(roi=roi) + with self.assertRaises(ValueError): + processor.preview_configuration() + + # 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]: + roi = {'horizontal':(25,40)} + processor = FluxNormaliser(roi=roi) + processor.set_input(data) + processor.preview_configuration() + + # for 3D, check no error specifying a single angle to plot + if data.geometry.dimension == '3D': + processor.preview_configuration(angle=1) + # if 2D, attempt to plot single angle should cause error + else: + with self.assertRaises(ValueError): + processor.preview_configuration(angle=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(channel=1) + # if single channel, check specifying channel causes an error + else: + with self.assertRaises(ValueError): + processor.preview_configuration(channel=1) def test_FluxNormaliser(self): From cf4b56f197fc32b6ad1310be63c22c9184d52f82 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 20 Aug 2024 14:20:46 +0000 Subject: [PATCH 14/25] Changelog and test format update --- CHANGELOG.md | 1 + Wrappers/Python/test/test_DataProcessor.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index eee67f4cc8..f3cd4a5450 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ - Enhancements: - Use ravel instead of flat in KullbackLeibler numba backend (#1874) - Upgrade Python wrapper (#1873, #1875) + - Add FluxNormaliser processor (#1878) - Bug fixes: - `ImageData` removes dimensions of size 1 from the input array. This fixes an issue where single slice reconstructions from 3D data would fail due to shape mismatches (#1885) - Make Binner accept accelerated=False (#1887) diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 846f825026..4f37fd336c 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3229,7 +3229,7 @@ def test_FluxNormaliser(self): processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) - def test_FluxNormaliser2(self): + # test roi with different data shapes for data in [self.data_cone, self.data_parallel, self.data_multichannel, self.data_slice, self.data_reorder]: From 3a0f7bb8341191339849dae11fa46aea5fa7f5d6 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 20 Aug 2024 14:30:56 +0000 Subject: [PATCH 15/25] Update mock --- Wrappers/Python/cil/processors/FluxNormaliser.py | 5 +++-- Wrappers/Python/test/test_DataProcessor.py | 3 +-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index f9b4ea726e..da8970aff1 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -200,7 +200,7 @@ def preview_configuration(self, angle=None, channel=None, log=False): 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") - + plt.figure() plt.subplot(212) if len(data.geometry.angles)==1: plt.plot(data.geometry.angles, self.flux, '.r', label='Mean') @@ -215,7 +215,8 @@ def preview_configuration(self, angle=None, channel=None, log=False): plt.ylabel('Intensity in roi') plt.grid() plt.tight_layout() - + plt.show() + def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=111): data = self.get_input() diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 4f37fd336c..a9b6091f61 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3144,7 +3144,7 @@ def test_check_input(self): with self.assertRaises(ValueError): processor.check_input(self.data_cone) - @patch("matplotlib.pyplot") + @patch("matplotlib.pyplot.figure") def test_preview_configuration(self, mock_plot): # Test error in preview configuration if there is no roi processor = FluxNormaliser(flux=10) @@ -3183,7 +3183,6 @@ def test_preview_configuration(self, mock_plot): with self.assertRaises(ValueError): processor.preview_configuration(channel=1) - def test_FluxNormaliser(self): #Test flux with no norm_value processor = FluxNormaliser(flux=10) From 787b6bdc396980405138bb6b45fed286075502fe Mon Sep 17 00:00:00 2001 From: hrobarts Date: Wed, 18 Sep 2024 15:04:24 +0000 Subject: [PATCH 16/25] Add colorbar to plots and warning based on data range --- .../Python/cil/processors/FluxNormaliser.py | 20 +++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index da8970aff1..583c6f121a 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -128,6 +128,16 @@ def check_input(self, dataset): self.roi.update({d:(0,dataset.get_dimension_size(d))}) self.flux = numpy.mean(dataset.array[tuple(slc)], axis=tuple(axes)) + + dataset_range = dataset.max(axis=tuple(axes)) - dataset.min(axis=tuple(axes)) + + if (numpy.mean(self.flux) > dataset.mean()): + if numpy.mean(self.flux/dataset_range.array) < 0.9: + log.warning('Warning: mean value in selected roi is more than 10 percent of data range - may not represent the background') + else: + if numpy.mean(self.flux/dataset_range.array) > 0.1: + log.warning('Warning: mean value in selected roi is more than 10 percent of data range - may not represent the background') + self.roi_slice = slc self.roi_axes = axes @@ -179,7 +189,7 @@ def preview_configuration(self, angle=None, channel=None, log=False): min = numpy.min(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) max = numpy.max(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) - + plt.figure() if data.geometry.dimension == '3D': if angle is None: if 'angle' in data.dimension_labels: @@ -200,7 +210,7 @@ def preview_configuration(self, angle=None, channel=None, log=False): 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") - plt.figure() + plt.subplot(212) if len(data.geometry.angles)==1: plt.plot(data.geometry.angles, self.flux, '.r', label='Mean') @@ -247,9 +257,11 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11 plt.subplot(ax) if log: - plt.imshow(numpy.log(data_slice.array), cmap='gray',aspect='equal', origin='lower', extent=extent) + im = plt.imshow(numpy.log(data_slice.array), cmap='gray',aspect='equal', origin='lower', extent=extent) + plt.gcf().colorbar(im, ax=plt.gca()) else: - plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower', extent=extent) + im = plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower', extent=extent) + plt.gcf().colorbar(im, ax=plt.gca()) h = data_slice.dimension_labels[1] v = data_slice.dimension_labels[0] From 21e773cddc014644a2246f883347561963ec2faf Mon Sep 17 00:00:00 2001 From: hrobarts Date: Wed, 18 Sep 2024 16:05:42 +0000 Subject: [PATCH 17/25] Review changes --- CHANGELOG.md | 2 +- .../Python/cil/processors/FluxNormaliser.py | 76 +++++++++++++------ Wrappers/Python/test/test_DataProcessor.py | 10 +-- docs/source/processors.rst | 16 ++-- 4 files changed, 67 insertions(+), 37 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2ec4bd3ed5..d4ad09745e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ - Added callback `optimisation.utilities.callbacks.EarlyStoppingObjectiveValue` which stops iterations if an algorithm objective changes less than a provided threshold (#1892) - Added callback `optimisation.utilities.callbacks.CGLSEarlyStopping` which replicates the automatic behaviour of CGLS in CIL versions <=24. (#1892) - Added `labels` module with `ImageDimension`, `AcquisitionDimension`, `AcquisitionType`, `AngleUnit`, `FillType` (#1692) + - Add FluxNormaliser processor (#1878) - Enhancements: - Use ravel instead of flat in KullbackLeibler numba backend (#1874) - Upgrade Python wrapper (#1873, #1875) @@ -16,7 +17,6 @@ - Internal refactor: Replaced string-based label checks with enum-based checks for improved type safety and consistency (#1692) - Internal refactor: Separate framework into multiple files (#1692) - Allow the SIRT algorithm to take `initial=None` (#1906) - - Add FluxNormaliser processor (#1878) - Testing: - New unit tests for operators and functions to check for in place errors and the behaviour of `out` (#1805) - Updates in SPDHG vs PDHG unit test to reduce test time and adjustments to parameters (#1898) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 583c6f121a..927be20f33 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -16,8 +16,7 @@ # Authors: # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt -from cil.framework import Processor, DataContainer, AcquisitionData,\ - AcquisitionGeometry, ImageGeometry, ImageData +from cil.framework import Processor, AcquisitionData from cil.utilities.display import show2D import numpy import logging @@ -29,7 +28,7 @@ class FluxNormaliser(Processor): r''' Flux normalisation based on float or region of interest - This processor reads in a AcquisitionData and normalises it based on + This processor reads in an AcquisitionData and normalises it based on a float or array of float values, or a region of interest. The normalised image :math:`I_{norm}` is calculated from the original image @@ -48,36 +47,64 @@ class FluxNormaliser(Processor): roi: dict (optional) Dictionary describing the region of interest containing the background - in the image. The image is divided by the mean value in the roi. If None, - specify flux direcrtly, default is None. + in the image. The roi is specified as `{‘axis_name1’:(start,stop), + ‘axis_name2’:(start,stop)}`, where the key is the axis name to calculate + the flux from. If the dataset has an axis which is not specified in the + roi dictionary, all data from that axis will be used in the calculation. + The image is divided by the mean value in the roi. If None, specify flux + directly, default is None. norm_value: float (optional) - The value to multiply the image by. If None, the mean flux value will - be used, default is None. - - tolerance: float (optional) - Small number to set to when there is a division by zero, default is 1e-5. + The value to multiply the image by. If None, the mean flux value will be used, + either the flux provided or the mean value in the roi across all + projections, default is None. Returns: -------- Output: AcquisitionData normalised by flux or mean intensity in roi + + Example + ------- + >>> from cil.processors import FluxNormaliser + >>> processor = FluxNormaliser(flux=10) + >>> processor.set_input(data) + >>> data_norm = processor.get_output() + + Example + ------- + >>> from cil.processors import FluxNormaliser + >>> processor = FluxNormaliser(flux=np.arange(1,2,(2-1)/(data.get_dimension_size('angle')))) + >>> processor.set_input(data) + >>> data_norm = processor.get_output() + + Example + ------- + >>> from cil.processors import FluxNormaliser + >>> processor = FluxNormaliser(flux=10) + >>> processor.set_input(data) + >>> data_norm = processor.get_output() + + Note + ---- + The roi indices provided are start inclusive, stop exclusive. + All elements along a dimension will be included if the axis does not appear + in the roi dictionary ''' - def __init__(self, flux=None, roi=None, norm_value=None, tolerance=1e-5): + def __init__(self, flux=None, roi=None, norm_value=None): kwargs = { 'flux' : flux, 'roi' : roi, 'roi_slice' : None, 'roi_axes' : None, - 'norm_value' : norm_value, - 'tolerance' : tolerance + 'norm_value' : norm_value } super(FluxNormaliser, self).__init__(**kwargs) def check_input(self, dataset): - if not issubclass(type(dataset), DataContainer): - raise TypeError("Expected DataContainer or subclass, found {}" + if not (type(dataset), AcquisitionData): + raise TypeError("Expected AcquistionData, found {}" .format(type(dataset))) if self.roi is not None: @@ -85,7 +112,7 @@ def check_input(self, dataset): raise ValueError("Please specify either flux or roi, not both") else: - if isinstance(self.roi,dict): + if isinstance(self.roi, dict): if not all (r in dataset.dimension_labels for r in self.roi): raise ValueError("roi labels must be in the dataset dimension_labels, found {}" .format(str(self.roi))) @@ -155,6 +182,13 @@ def check_input(self, dataset): raise ValueError("Flux must be a scalar or array with length \ \n = data.geometry.angles, found {} and {}" .format(flux_size, data_size)) + if numpy.any(self.flux==0): + raise ValueError('Flux value can\'t be 0, provide a different flux\ + or region of interest with non-zero values') + else: + if self.flux==0: + raise ValueError('Flux value can\'t be 0, provide a different flux\ + or region of interest with non-zero values') if self.norm_value is None: self.norm_value = numpy.mean(self.flux) @@ -165,7 +199,7 @@ def preview_configuration(self, angle=None, channel=None, log=False): ''' Preview the FluxNormalisation processor configuration for roi mode. Plots the region of interest on the image and the mean, maximum and - minimum intensity in the roi. the angles with minimum and maximum intensity in the roi, + minimum intensity in the roi. Parameters: ----------- @@ -315,12 +349,8 @@ def process(self, out=None): if len(flux_size) > 0: f = self.flux[i] slice_proj[proj_axis] = i - with numpy.errstate(divide='ignore', invalid='ignore'): - out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]*self.norm_value/f + out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]*self.norm_value/f else: - with numpy.errstate(divide='ignore', invalid='ignore'): - out.array = data.array*self.norm_value/f - - out.array[ ~ numpy.isfinite( out.array )] = self.tolerance + out.array = data.array*self.norm_value/f return out diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 0fdbabb60d..630e6d7180 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3100,7 +3100,7 @@ def setUp(self): self.data_multichannel = ag.allocate('random') self.data_slice = self.data_parallel.get_slice(vertical=1) - self.data_reorder = self.data_cone + self.data_reorder = self.data_cone.copy() self.data_reorder.reorder(['horizontal','angle','vertical']) self.data_single_angle = self.data_cone.get_slice(angle=1) @@ -3129,7 +3129,7 @@ def test_check_input(self): with self.assertRaises(ValueError): processor.check_input(self.data_cone) - # check there is an error if no flux array size is not equal to the number of angles in data + # 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]) with self.assertRaises(ValueError): processor.check_input(self.data_cone) @@ -3207,13 +3207,13 @@ def test_FluxNormaliser(self): numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) #Test flux array with no norm_value - processor = FluxNormaliser(flux=10*numpy.ones(self.data_cone.get_dimension_size('angle'))) + processor = FluxNormaliser(flux=numpy.arange(1,2,(2-1)/(data.get_dimension_size('angle')))) processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) #Test flux array with norm_value - processor = FluxNormaliser(flux=10*numpy.ones(self.data_cone.get_dimension_size('angle')), norm_value=5) + processor = FluxNormaliser(flux=numpy.arange(1,2,(2-1)/(data.get_dimension_size('angle'))), norm_value=5) processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) @@ -3239,7 +3239,7 @@ def test_FluxNormaliser(self): data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) - # test roi with different data shapes + # test roi with different data shapes and different flux values per projection for data in [self.data_cone, self.data_parallel, self.data_multichannel, self.data_slice, self.data_reorder]: roi = {'horizontal':(25,40)} diff --git a/docs/source/processors.rst b/docs/source/processors.rst index bd80c074e0..12b26f9c6b 100644 --- a/docs/source/processors.rst +++ b/docs/source/processors.rst @@ -112,6 +112,14 @@ Data Normaliser :members: :inherited-members: set_input, get_output +Flux Normaliser +----------------- + +.. autoclass:: cil.processors.FluxNormaliser + :exclude-members: check_input, get_input + :members: + :inherited-members: set_input, get_output + Transmission to Absorption Converter ------------------------------------- @@ -148,12 +156,4 @@ Paganin Processor :members: :inherited-members: set_input, get_output -Flux Normaliser ------------------ - -.. autoclass:: cil.processors.FluxNormaliser - :exclude-members: check_input, get_input - :members: - :inherited-members: set_input, get_output - :ref:`Return Home ` From e8423bbf88c4c3322b0543e3c403ba9ae493acb5 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Mon, 23 Sep 2024 08:39:19 +0000 Subject: [PATCH 18/25] Fix tests after removing tolerance --- Wrappers/Python/cil/processors/FluxNormaliser.py | 6 +++--- Wrappers/Python/test/test_DataProcessor.py | 12 ++++++------ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 927be20f33..8a91cb41a8 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -156,13 +156,13 @@ def check_input(self, dataset): self.flux = numpy.mean(dataset.array[tuple(slc)], axis=tuple(axes)) - dataset_range = dataset.max(axis=tuple(axes)) - dataset.min(axis=tuple(axes)) + dataset_range = numpy.max(dataset.array, axis=tuple(axes)) - numpy.min(dataset.array, axis=tuple(axes)) if (numpy.mean(self.flux) > dataset.mean()): - if numpy.mean(self.flux/dataset_range.array) < 0.9: + if numpy.mean(self.flux/dataset_range) < 0.9: log.warning('Warning: mean value in selected roi is more than 10 percent of data range - may not represent the background') else: - if numpy.mean(self.flux/dataset_range.array) > 0.1: + if numpy.mean(self.flux/dataset_range) > 0.1: log.warning('Warning: mean value in selected roi is more than 10 percent of data range - may not represent the background') self.roi_slice = slc diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 630e6d7180..3e23a11526 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3110,15 +3110,15 @@ def error_message(self,processor, test_parameter): def test_init(self): # test default values are initialised processor = FluxNormaliser() - test_parameter = ['flux','roi','norm_value', 'tolerance'] - test_value = [None, None, None, 1e-5] + test_parameter = ['flux','roi','norm_value'] + test_value = [None, None, None] 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 = FluxNormaliser(1,2,3,4) - test_value = [1, 2, 3, 4] + processor = FluxNormaliser(1,2,3) + test_value = [1, 2, 3] 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])) @@ -3207,13 +3207,13 @@ def test_FluxNormaliser(self): numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) #Test flux array with no norm_value - processor = FluxNormaliser(flux=numpy.arange(1,2,(2-1)/(data.get_dimension_size('angle')))) + processor = FluxNormaliser(flux=numpy.arange(1,2,(2-1)/(self.data_cone.get_dimension_size('angle')))) processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) #Test flux array with norm_value - processor = FluxNormaliser(flux=numpy.arange(1,2,(2-1)/(data.get_dimension_size('angle'))), norm_value=5) + processor = FluxNormaliser(flux=numpy.arange(1,2,(2-1)/(self.data_cone.get_dimension_size('angle'))), norm_value=5) processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) From 1d7e15793cbc9ce93f5bca9a580e254e84f29427 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 24 Sep 2024 09:27:12 +0000 Subject: [PATCH 19/25] Update tests --- Wrappers/Python/test/test_DataProcessor.py | 143 +++++++++++---------- 1 file changed, 77 insertions(+), 66 deletions(-) diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 3e23a11526..33afc7e330 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3207,78 +3207,89 @@ def test_FluxNormaliser(self): numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) #Test flux array with no norm_value - processor = FluxNormaliser(flux=numpy.arange(1,2,(2-1)/(self.data_cone.get_dimension_size('angle')))) + flux = numpy.arange(1,2,(2-1)/(self.data_cone.get_dimension_size('angle'))) + processor = FluxNormaliser(flux=flux) processor.set_input(self.data_cone) data_norm = processor.get_output() - numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) - - #Test flux array with norm_value - processor = FluxNormaliser(flux=numpy.arange(1,2,(2-1)/(self.data_cone.get_dimension_size('angle'))), norm_value=5) - processor.set_input(self.data_cone) - data_norm = processor.get_output() - numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) - - #Test roi with no norm_value - roi = {'vertical':(0,10), 'horizontal':(0,10)} - processor = FluxNormaliser(roi=roi) - processor.set_input(self.data_cone) - data_norm = processor.get_output() - numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) - - #Test roi with norm_value - roi = {'vertical':(0,10), 'horizontal':(0,10)} - processor = FluxNormaliser(roi=roi, norm_value=5) + data_norm_test = self.data_cone.copy() + for a in range(data_norm_test.get_dimension_size('angle')): + data_norm_test.array[a,:,:] /= flux[a] + data_norm_test.array[a,:,:]*= numpy.mean(flux) + numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6) + + # #Test flux array with norm_value + flux = numpy.arange(1,2,(2-1)/(self.data_cone.get_dimension_size('angle'))) + norm_value = 5 + processor = FluxNormaliser(flux=flux, norm_value=norm_value) processor.set_input(self.data_cone) data_norm = processor.get_output() - numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) - - # Test roi with just one dimension - roi = {'vertical':(0,2)} - processor = FluxNormaliser(roi=roi, norm_value=5) - processor.set_input(self.data_cone) - data_norm = processor.get_output() - numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) + data_norm_test = self.data_cone.copy() + for a in range(data_norm_test.get_dimension_size('angle')): + data_norm_test.array[a,:,:] /= flux[a] + data_norm_test.array[a,:,:]*= norm_value + numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6) + + # #Test roi with no norm_value + # roi = {'vertical':(0,10), 'horizontal':(0,10)} + # processor = FluxNormaliser(roi=roi) + # processor.set_input(self.data_cone) + # data_norm = processor.get_output() + # numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) + + # #Test roi with norm_value + # roi = {'vertical':(0,10), 'horizontal':(0,10)} + # processor = FluxNormaliser(roi=roi, norm_value=5) + # processor.set_input(self.data_cone) + # data_norm = processor.get_output() + # numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) + + # # Test roi with just one dimension + # roi = {'vertical':(0,2)} + # processor = FluxNormaliser(roi=roi, norm_value=5) + # processor.set_input(self.data_cone) + # data_norm = processor.get_output() + # numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) # test roi with different data shapes and different flux values per projection - for data in [self.data_cone, self.data_parallel, self.data_multichannel, - self.data_slice, self.data_reorder]: - roi = {'horizontal':(25,40)} - processor = FluxNormaliser(roi=roi, norm_value=5) - processor.set_input(data) - data_norm = processor.get_output() - - ax = data.get_dimension_axis('horizontal') - slc = [slice(None)]*len(data.shape) - slc[ax] = slice(25,40) - axes=[ax] - if 'vertical' in data.dimension_labels: - axes.append(data.get_dimension_axis('vertical')) - if 'channel' in data.dimension_labels: - axes.append(data.get_dimension_axis('channel')) - flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) - slice_proj = [slice(None)]*len(data.shape) - proj_axis = data.get_dimension_axis('angle') - data_norm_test = data.copy() - for i in range(len(data.geometry.angles)): - f = flux[i] - slice_proj[proj_axis] = i - with numpy.errstate(divide='ignore', invalid='ignore'): - data_norm_test.array[tuple(slice_proj)] = 5/f*data.array[tuple(slice_proj)] - numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6, - err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) - - data = self.data_single_angle - processor = FluxNormaliser(roi=roi, norm_value=5) - processor.set_input(data) - data_norm = processor.get_output() - ax = data.get_dimension_axis('horizontal') - slc = [slice(None)]*len(data.shape) - slc[ax] = slice(25,40) - axes=[ax,data.get_dimension_axis('vertical')] - flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) - - numpy.testing.assert_allclose(data_norm.array, 5/flux*data.array, atol=1e-6, - err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) + # for data in [self.data_cone, self.data_parallel, self.data_multichannel, + # self.data_slice, self.data_reorder]: + # roi = {'horizontal':(25,40)} + # processor = FluxNormaliser(roi=roi, norm_value=5) + # processor.set_input(data) + # data_norm = processor.get_output() + + # ax = data.get_dimension_axis('horizontal') + # slc = [slice(None)]*len(data.shape) + # slc[ax] = slice(25,40) + # axes=[ax] + # if 'vertical' in data.dimension_labels: + # axes.append(data.get_dimension_axis('vertical')) + # if 'channel' in data.dimension_labels: + # axes.append(data.get_dimension_axis('channel')) + # flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) + # slice_proj = [slice(None)]*len(data.shape) + # proj_axis = data.get_dimension_axis('angle') + # data_norm_test = data.copy() + # for i in range(len(data.geometry.angles)): + # f = flux[i] + # slice_proj[proj_axis] = i + # with numpy.errstate(divide='ignore', invalid='ignore'): + # data_norm_test.array[tuple(slice_proj)] = 5/f*data.array[tuple(slice_proj)] + # numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6, + # err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) + + # data = self.data_single_angle + # processor = FluxNormaliser(roi=roi, norm_value=5) + # processor.set_input(data) + # data_norm = processor.get_output() + # ax = data.get_dimension_axis('horizontal') + # slc = [slice(None)]*len(data.shape) + # slc[ax] = slice(25,40) + # axes=[ax,data.get_dimension_axis('vertical')] + # flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) + + # numpy.testing.assert_allclose(data_norm.array, 5/flux*data.array, atol=1e-6, + # err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) if __name__ == "__main__": From ce2fe22d6f5bea2e06eb32d9123faafb335f411f Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 1 Oct 2024 13:59:51 +0000 Subject: [PATCH 20/25] Add _calculate_flux function --- .../Python/cil/processors/FluxNormaliser.py | 183 +++++++++--------- 1 file changed, 93 insertions(+), 90 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 8a91cb41a8..14f6bcfe2c 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -28,40 +28,32 @@ class FluxNormaliser(Processor): r''' Flux normalisation based on float or region of interest - This processor reads in an AcquisitionData and normalises it based on - a float or array of float values, or a region of interest. - - The normalised image :math:`I_{norm}` is calculated from the original image - :math:`I` by - - .. math:: I_{norm} = I\frac{n}{F} - - where :math:`F` is the flux and :math:`n` is the norm_value + This processor reads in an AcquisitionData and normalises it by flux from + a float or array of float values, or the mean flux in a region of interest. Parameters: ----------- flux: float or list of floats (optional) The value to divide the image by. If flux is a list it must have length - equal to the number of angles in the dataset. If flux=None, calculate + equal to the number of projections in the dataset. If flux=None, calculate flux from the roi, default is None. roi: dict (optional) Dictionary describing the region of interest containing the background - in the image. The roi is specified as `{‘axis_name1’:(start,stop), - ‘axis_name2’:(start,stop)}`, where the key is the axis name to calculate - the flux from. If the dataset has an axis which is not specified in the - roi dictionary, all data from that axis will be used in the calculation. - The image is divided by the mean value in the roi. If None, specify flux - directly, default is None. - - norm_value: float (optional) - The value to multiply the image by. If None, the mean flux value will be used, - either the flux provided or the mean value in the roi across all - projections, default is None. + in the image. The roi is specified as `{'axis_name1':(start,stop), + 'axis_name2':(start,stop)}`, where the key is the axis name 'vertical' + and/ or 'horizontal'. If an axis is not specified in the roi dictionary, + the full range will be used, default is None. + + target: string or float + The value to scale the normalised data with. If float, the data is scaled + to the float value. If string 'mean' the data is scaled to the mean value + of the input flux or flux in the roi, if 'first' the data is scaled to + the first input flux value or the flux in the roi of the first projection. Returns: -------- - Output: AcquisitionData normalised by flux or mean intensity in roi + Output: AcquisitionData normalised by flux Example ------- @@ -80,7 +72,7 @@ class FluxNormaliser(Processor): Example ------- >>> from cil.processors import FluxNormaliser - >>> processor = FluxNormaliser(flux=10) + >>> processor = FluxNormaliser(roi=(roi={'horizontal':(5, 15)) >>> processor.set_input(data) >>> data_norm = processor.get_output() @@ -91,13 +83,19 @@ class FluxNormaliser(Processor): in the roi dictionary ''' - def __init__(self, flux=None, roi=None, norm_value=None): + def __init__(self, flux=None, roi=None, target=None): + + if roi is not None and flux is not None: + raise ValueError("Please specify either flux or roi, not both") + if roi is None and flux is None: + raise ValueError("Please specify either flux or roi, found None") + kwargs = { 'flux' : flux, 'roi' : roi, 'roi_slice' : None, 'roi_axes' : None, - 'norm_value' : norm_value + 'target' : target } super(FluxNormaliser, self).__init__(**kwargs) @@ -106,81 +104,80 @@ def check_input(self, dataset): if not (type(dataset), AcquisitionData): raise TypeError("Expected AcquistionData, found {}" .format(type(dataset))) + + return True - if self.roi is not None: - if self.flux is not None: - raise ValueError("Please specify either flux or roi, not both") - else: - if isinstance(self.roi, dict): - if not all (r in dataset.dimension_labels for r in self.roi): - raise ValueError("roi labels must be in the dataset dimension_labels, found {}" - .format(str(self.roi))) - - slc = [slice(None)]*len(dataset.shape) - axes=[] - - for r in self.roi: - # only allow roi to be specified in horizontal and vertical - if (r != 'horizontal' and r != 'vertical'): - raise ValueError("roi must be 'horizontal' or 'vertical', found '{}'" - .format(str(r))) - - dimension_label_list = list(dataset.dimension_labels) - # loop through all dimensions in the dataset apart from angle - if 'angle' in dimension_label_list: - dimension_label_list.remove('angle') - - for d in dimension_label_list: - # check the dimension is in the user specified roi - if d in self.roi: - # check indices are ints - if not all(isinstance(i, int) for i in self.roi[d]): - raise TypeError("roi values must be int, found {} and {}" - .format(str(type(self.roi[d][0])), str(type(self.roi[d][1])))) - # check indices are in range - elif (self.roi[d][0] >= self.roi[d][1]) or (self.roi[d][0] < 0) or self.roi[d][1] > dataset.get_dimension_size(d): - raise ValueError("roi values must be start > stop and between 0 and {}, found start={} and stop={} for direction '{}'" - .format(str(dataset.get_dimension_size(d)), str(self.roi[d][0]), str(self.roi[d][1]), d )) - # create slice - else: - ax = dataset.get_dimension_axis(d) - slc[ax] = slice(self.roi[d][0], self.roi[d][1]) - axes.append(ax) - - # if the dimension is not in roi, average across the whole dimension + def _calculate_flux(self): + dataset = self.get_input() + # convert flux to float32 + + if self.flux is None: + + if isinstance(self.roi, dict): + if not all (r in dataset.dimension_labels for r in self.roi): + raise ValueError("roi labels must be in the dataset dimension_labels, found {}" + .format(str(self.roi))) + + slc = [slice(None)]*len(dataset.shape) + axes=[] + + for r in self.roi: + # only allow roi to be specified in horizontal and vertical + if (r != 'horizontal' and r != 'vertical'): + raise ValueError("roi must be 'horizontal' or 'vertical', found '{}'" + .format(str(r))) + + dimension_label_list = list(self.roi.keys()) + + for d in dimension_label_list: + # check the dimension is in the user specified roi + if d in self.roi: + # check indices are ints + if not all(isinstance(i, int) for i in self.roi[d]): + raise TypeError("roi values must be int, found {} and {}" + .format(str(type(self.roi[d][0])), str(type(self.roi[d][1])))) + # check indices are in range + elif (self.roi[d][0] >= self.roi[d][1]) or (self.roi[d][0] < 0) or self.roi[d][1] > dataset.get_dimension_size(d): + raise ValueError("roi values must be start > stop and between 0 and {}, found start={} and stop={} for direction '{}'" + .format(str(dataset.get_dimension_size(d)), str(self.roi[d][0]), str(self.roi[d][1]), d )) + # create slice else: ax = dataset.get_dimension_axis(d) + slc[ax] = slice(self.roi[d][0], self.roi[d][1]) axes.append(ax) - self.roi.update({d:(0,dataset.get_dimension_size(d))}) - - self.flux = numpy.mean(dataset.array[tuple(slc)], axis=tuple(axes)) - - dataset_range = numpy.max(dataset.array, axis=tuple(axes)) - numpy.min(dataset.array, axis=tuple(axes)) - - if (numpy.mean(self.flux) > dataset.mean()): - if numpy.mean(self.flux/dataset_range) < 0.9: - log.warning('Warning: mean value in selected roi is more than 10 percent of data range - may not represent the background') + + # if the dimension is not in roi, average across the whole dimension else: - if numpy.mean(self.flux/dataset_range) > 0.1: - log.warning('Warning: mean value in selected roi is more than 10 percent of data range - may not represent the background') + ax = dataset.get_dimension_axis(d) + axes.append(ax) + self.roi.update({d:(0,dataset.get_dimension_size(d))}) - self.roi_slice = slc - self.roi_axes = axes - + self.flux = numpy.mean(dataset.array[tuple(slc)], axis=tuple(axes)) + + dataset_range = numpy.max(dataset.array, axis=tuple(axes)) - numpy.min(dataset.array, axis=tuple(axes)) + + if (numpy.mean(self.flux) > dataset.mean()): + if numpy.mean(self.flux/dataset_range) < 0.9: + log.warning('Warning: mean value in selected roi is more than 10 percent of data range - may not represent the background') else: - raise TypeError("roi must be a dictionary, found {}" - .format(str(type(self.roi)))) - else: - if self.flux is None: - raise ValueError("Please specify flux or roi") + if numpy.mean(self.flux/dataset_range) > 0.1: + log.warning('Warning: mean value in selected roi is more than 10 percent of data range - may not represent the background') + + self.roi_slice = slc + self.roi_axes = axes + + else: + raise TypeError("roi must be a dictionary, found {}" + .format(str(type(self.roi)))) flux_size = (numpy.shape(self.flux)) if len(flux_size) > 0: - data_size = numpy.shape(dataset.geometry.angles) + # check this + data_size = numpy.shape(dataset.geometry.angles)*numpy.shape(dataset.geometry.channels) if data_size != flux_size: raise ValueError("Flux must be a scalar or array with length \ - \n = data.geometry.angles, found {} and {}" + \n = number of projections, found {} and {}" .format(flux_size, data_size)) if numpy.any(self.flux==0): raise ValueError('Flux value can\'t be 0, provide a different flux\ @@ -190,10 +187,15 @@ def check_input(self, dataset): raise ValueError('Flux value can\'t be 0, provide a different flux\ or region of interest with non-zero values') + self.flux = numpy.float32(self.flux) + if self.norm_value is None: self.norm_value = numpy.mean(self.flux) - return True + # def _calculate_target(): + # if isinstance(self.target, float): + + def preview_configuration(self, angle=None, channel=None, log=False): ''' @@ -219,6 +221,7 @@ def preview_configuration(self, angle=None, channel=None, log=False): if self.roi_slice is None: raise ValueError('Preview available with roi, run `processor= FluxNormaliser(roi=roi)` then `set_input(data)`') else: + self._calculate_flux() data = self.get_input() min = numpy.min(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) @@ -255,7 +258,7 @@ def preview_configuration(self, angle=None, channel=None, log=False): plt.plot(data.geometry.angles, min,'--k', label='Minimum') plt.plot(data.geometry.angles, max,'--k', label='Maximum') plt.legend() - plt.xlabel('Angle') + plt.xlabel('angle') plt.ylabel('Intensity in roi') plt.grid() plt.tight_layout() @@ -331,7 +334,7 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11 plt.ylabel(v) def process(self, out=None): - + self._calculate_flux() data = self.get_input() if out is None: From db715e6af2faa7f1b8e806650f664d7f43b83158 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Thu, 10 Oct 2024 10:48:05 +0000 Subject: [PATCH 21/25] Change target calculation --- .../Python/cil/processors/FluxNormaliser.py | 32 +++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 14f6bcfe2c..c45b1593e7 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -95,7 +95,8 @@ def __init__(self, flux=None, roi=None, target=None): 'roi' : roi, 'roi_slice' : None, 'roi_axes' : None, - 'target' : target + 'target' : target, + 'target_value' : None, } super(FluxNormaliser, self).__init__(**kwargs) @@ -174,7 +175,7 @@ def _calculate_flux(self): flux_size = (numpy.shape(self.flux)) if len(flux_size) > 0: # check this - data_size = numpy.shape(dataset.geometry.angles)*numpy.shape(dataset.geometry.channels) + data_size = numpy.shape(dataset.geometry.angles) # make this also account for channels if data_size != flux_size: raise ValueError("Flux must be a scalar or array with length \ \n = number of projections, found {} and {}" @@ -189,14 +190,26 @@ def _calculate_flux(self): self.flux = numpy.float32(self.flux) - if self.norm_value is None: - self.norm_value = numpy.mean(self.flux) + def _calculate_target(self): + if isinstance(self.target, float): + self.target_value = self.target + elif isinstance(self.target, str): + if self.target == 'first': + if len(numpy.shape(self.flux)) > 0 : + self.target_value = self.flux[0] + else: + self.target_value = self.flux + elif self.target == 'mean': + self.target_value = numpy.mean(self.flux) + else: + raise ValueError("Target string not recognised, found {}, expected 'first' or 'mean'" + .format(self.target)) + else: + raise TypeError("Target must be string or float, found {}" + .format(type(self.target))) - # def _calculate_target(): - # if isinstance(self.target, float): - def preview_configuration(self, angle=None, channel=None, log=False): ''' Preview the FluxNormalisation processor configuration for roi mode. @@ -335,6 +348,7 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11 def process(self, out=None): self._calculate_flux() + self._calculate_target() data = self.get_input() if out is None: @@ -352,8 +366,8 @@ def process(self, out=None): if len(flux_size) > 0: f = self.flux[i] slice_proj[proj_axis] = i - out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]*self.norm_value/f + out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]*self.target_value/f else: - out.array = data.array*self.norm_value/f + out.array = data.array*self.target_value/f return out From 69e4466aa839f823fe9625120b2249add548c275 Mon Sep 17 00:00:00 2001 From: hrobarts Date: Tue, 15 Oct 2024 18:20:41 +0000 Subject: [PATCH 22/25] Require h and v data order at the end --- .../Python/cil/processors/FluxNormaliser.py | 188 +++++++++++------- Wrappers/Python/test/test_DataProcessor.py | 43 ++-- 2 files changed, 144 insertions(+), 87 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index c45b1593e7..769f530e4e 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -48,8 +48,9 @@ class FluxNormaliser(Processor): target: string or float The value to scale the normalised data with. If float, the data is scaled to the float value. If string 'mean' the data is scaled to the mean value - of the input flux or flux in the roi, if 'first' the data is scaled to + of the input flux or flux across all rois, if 'first' the data is scaled to the first input flux value or the flux in the roi of the first projection. + Default is 'mean' Returns: -------- @@ -83,12 +84,7 @@ class FluxNormaliser(Processor): in the roi dictionary ''' - def __init__(self, flux=None, roi=None, target=None): - - if roi is not None and flux is not None: - raise ValueError("Please specify either flux or roi, not both") - if roi is None and flux is None: - raise ValueError("Please specify either flux or roi, found None") + def __init__(self, flux=None, roi=None, target='mean'): kwargs = { 'flux' : flux, @@ -97,22 +93,54 @@ def __init__(self, flux=None, roi=None, target=None): 'roi_axes' : None, 'target' : target, 'target_value' : None, + 'v_size' : 1, + 'v_axis' : None, + 'h_size' : 1, + 'h_axis' : None } super(FluxNormaliser, self).__init__(**kwargs) def check_input(self, dataset): + + if self.roi is not None and self.flux is not None: + raise ValueError("Please specify either flux or roi, not both") + if self.roi is None and self.flux is None: + raise ValueError("Please specify either flux or roi, found None") if not (type(dataset), AcquisitionData): raise TypeError("Expected AcquistionData, found {}" .format(type(dataset))) + image_axes = 0 + if 'vertical' in dataset.dimension_labels: + v_axis = dataset.get_dimension_axis('vertical') + self.v_size = dataset.get_dimension_size('vertical') + image_axes += 1 + + if 'horizontal' in dataset.dimension_labels: + self.h_axis = dataset.get_dimension_axis('horizontal') + self.h_size = dataset.get_dimension_size('horizontal') + image_axes += 1 + + if (( self.h_axis is not None) and (self. h_axis < (len(dataset.shape)-image_axes))) or \ + ((self.v_axis is not None) and self.v_axis < (len(dataset.shape)-image_axes)): + raise ValueError('Projections must be the last two axes of the dataset') + return True def _calculate_flux(self): - dataset = self.get_input() - # convert flux to float32 + ''' + Function to calculate flux from a region of interest in the data. If the + flux is already provided as an array, convert the array to float 32 and + check the size matches the number of projections + ''' + dataset = self.get_input() + if dataset is None: + raise ValueError('Data not found, please run `set_input(data)`') + + # Calculate the flux from the roi in the data if self.flux is None: if isinstance(self.roi, dict): @@ -128,11 +156,8 @@ def _calculate_flux(self): if (r != 'horizontal' and r != 'vertical'): raise ValueError("roi must be 'horizontal' or 'vertical', found '{}'" .format(str(r))) - - dimension_label_list = list(self.roi.keys()) - - for d in dimension_label_list: - # check the dimension is in the user specified roi + + for d in ['horizontal', 'vertical']: if d in self.roi: # check indices are ints if not all(isinstance(i, int) for i in self.roi[d]): @@ -147,15 +172,16 @@ def _calculate_flux(self): ax = dataset.get_dimension_axis(d) slc[ax] = slice(self.roi[d][0], self.roi[d][1]) axes.append(ax) - - # if the dimension is not in roi, average across the whole dimension + # if a projection dimension isn't in the roi, use the whole axis else: - ax = dataset.get_dimension_axis(d) - axes.append(ax) - self.roi.update({d:(0,dataset.get_dimension_size(d))}) + if d in dataset.dimension_labels: + ax = dataset.get_dimension_axis(d) + axes.append(ax) + self.roi.update({d:(0,dataset.get_dimension_size(d))}) self.flux = numpy.mean(dataset.array[tuple(slc)], axis=tuple(axes)) + # Warn if the flux is more than 10% of the dataset range dataset_range = numpy.max(dataset.array, axis=tuple(axes)) - numpy.min(dataset.array, axis=tuple(axes)) if (numpy.mean(self.flux) > dataset.mean()): @@ -171,36 +197,41 @@ def _calculate_flux(self): else: raise TypeError("roi must be a dictionary, found {}" .format(str(type(self.roi)))) + + - flux_size = (numpy.shape(self.flux)) - if len(flux_size) > 0: - # check this - data_size = numpy.shape(dataset.geometry.angles) # make this also account for channels - if data_size != flux_size: + # convert flux array to float32 + self.flux = numpy.array(self.flux, dtype=numpy.float32, ndmin=1) + + # check flux array is the right size + flux_size_flat = len(self.flux.ravel()) + if flux_size_flat > 1: + data_size_flat = len(dataset.geometry.angles)*dataset.geometry.channels + if data_size_flat != flux_size_flat: raise ValueError("Flux must be a scalar or array with length \ \n = number of projections, found {} and {}" - .format(flux_size, data_size)) - if numpy.any(self.flux==0): - raise ValueError('Flux value can\'t be 0, provide a different flux\ - or region of interest with non-zero values') - else: - if self.flux==0: - raise ValueError('Flux value can\'t be 0, provide a different flux\ - or region of interest with non-zero values') + .format(flux_size_flat, data_size_flat)) - self.flux = numpy.float32(self.flux) + # check if flux array contains 0s + if 0 in self.flux: + raise ValueError('Flux value can\'t be 0, provide a different flux\ + or region of interest with non-zero values') + def _calculate_target(self): + ''' + Calculate the target value for the normalisation + ''' if isinstance(self.target, float): self.target_value = self.target elif isinstance(self.target, str): if self.target == 'first': if len(numpy.shape(self.flux)) > 0 : - self.target_value = self.flux[0] + self.target_value = self.flux.flat[0] else: self.target_value = self.flux elif self.target == 'mean': - self.target_value = numpy.mean(self.flux) + self.target_value = numpy.mean(self.flux.ravel()) else: raise ValueError("Target string not recognised, found {}, expected 'first' or 'mean'" .format(self.target)) @@ -231,15 +262,29 @@ def preview_configuration(self, angle=None, channel=None, log=False): log: bool (optional) If True, plot the image with a log scale, default is False ''' + self._calculate_flux() if self.roi_slice is None: raise ValueError('Preview available with roi, run `processor= FluxNormaliser(roi=roi)` then `set_input(data)`') else: - self._calculate_flux() - data = self.get_input() + data = self.get_input() + min = numpy.min(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) max = numpy.max(data.array[tuple(self.roi_slice)], axis=tuple(self.roi_axes)) - plt.figure() + if 'channel' in data.dimension_labels: + if channel is None: + channel = int(data.get_dimension_size('channel')/2) + channel_axis = data.get_dimension_axis('channel') + flux_array = self.flux.take(indices=channel, axis=channel_axis) + min = min.take(indices=channel, axis=channel_axis) + max = max.take(indices=channel, axis=channel_axis) + else: + if channel is not None: + raise ValueError("Channel not found") + else: + flux_array = self.flux + + plt.figure(figsize=(8,8)) if data.geometry.dimension == '3D': if angle is None: if 'angle' in data.dimension_labels: @@ -263,13 +308,14 @@ def preview_configuration(self, angle=None, channel=None, log=False): plt.subplot(212) if len(data.geometry.angles)==1: - plt.plot(data.geometry.angles, self.flux, '.r', label='Mean') + plt.plot(data.geometry.angles, flux_array, '.r', label='Mean') plt.plot(data.geometry.angles, min,'.k', label='Minimum') plt.plot(data.geometry.angles, max,'.k', label='Maximum') else: - plt.plot(data.geometry.angles, self.flux, 'r', label='Mean') + plt.plot(data.geometry.angles, flux_array, 'r', label='Mean') plt.plot(data.geometry.angles, min,'--k', label='Minimum') plt.plot(data.geometry.angles, max,'--k', label='Maximum') + plt.legend() plt.xlabel('angle') plt.ylabel('Intensity in roi') @@ -286,12 +332,7 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11 data_slice = data if 'channel' in data.dimension_labels: - if channel_index is None: - channel_index = int(data_slice.get_dimension_size('channel')/2) data_slice = data_slice.get_slice(channel=channel_index) - else: - if channel_index is not None: - raise ValueError("Channel not found") 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'") @@ -305,13 +346,13 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11 extent[i*2]=min_angle extent[i*2+1]=max_angle - plt.subplot(ax) + ax1 = plt.subplot(ax) if log: - im = plt.imshow(numpy.log(data_slice.array), cmap='gray',aspect='equal', origin='lower', extent=extent) - plt.gcf().colorbar(im, ax=plt.gca()) + im = ax1.imshow(numpy.log(data_slice.array), cmap='gray',aspect='equal', origin='lower', extent=extent) + plt.gcf().colorbar(im, ax=ax1) else: - im = plt.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower', extent=extent) - plt.gcf().colorbar(im, ax=plt.gca()) + im = ax1.imshow(data_slice.array, cmap='gray',aspect='equal', origin='lower', extent=extent) + plt.gcf().colorbar(im, ax=ax1) h = data_slice.dimension_labels[1] v = data_slice.dimension_labels[0] @@ -330,21 +371,21 @@ def _plot_slice_roi(self, angle_index=None, channel_index=None, log=False, ax=11 v_min = self.roi[v][0] v_max = self.roi[v][1] - plt.plot([h_min, h_max],[v_min, v_min],'--r') - plt.plot([h_min, h_max],[v_max, v_max],'--r') + ax1.plot([h_min, h_max],[v_min, v_min],'--r') + ax1.plot([h_min, h_max],[v_max, v_max],'--r') - plt.plot([h_min, h_min],[v_min, v_max],'--r') - plt.plot([h_max, h_max],[v_min, v_max],'--r') + ax1.plot([h_min, h_min],[v_min, v_max],'--r') + ax1.plot([h_max, h_max],[v_min, v_max],'--r') title = 'ROI' if angle_index is not None: title += ' angle = ' + str(data.geometry.angles[angle_index]) if channel_index is not None: title += ' channel = ' + str(channel_index) - plt.title(title) + ax1.set_title(title) - plt.xlabel(h) - plt.ylabel(v) + ax1.set_xlabel(h) + ax1.set_ylabel(v) def process(self, out=None): self._calculate_flux() @@ -354,20 +395,29 @@ def process(self, out=None): if out is None: out = data.copy() - flux_size = (numpy.shape(self.flux)) + proj_size = self.v_size*self.h_size + num_proj = int(data.array.size / proj_size) + f = self.flux + for i in range(num_proj): + arr_proj = data.array.flat[i*proj_size:(i+1)*proj_size] + if len(self.flux.flat) > 1: + f = self.flux.flat[i] + arr_proj *= self.target_value/f + out.array.flat[i*proj_size:(i+1)*proj_size] = arr_proj + - if 'angle' in data.dimension_labels: - proj_axis = data.get_dimension_axis('angle') - slice_proj = [slice(None)]*len(data.shape) - slice_proj[proj_axis] = 0 + # if 'angle' in data.dimension_labels: + # proj_axis = data.get_dimension_axis('angle') + # slice_proj = [slice(None)]*len(data.shape) + # slice_proj[proj_axis] = 0 - for i in range(len(data.geometry.angles)): - if len(flux_size) > 0: - f = self.flux[i] - slice_proj[proj_axis] = i - out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]*self.target_value/f - else: - out.array = data.array*self.target_value/f + # for i in range(len(data.geometry.angles)): + # if len(flux_size) > 0: + # f = self.flux[i] + # slice_proj[proj_axis] = i + # out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]*self.target_value/f + # else: + # out.array = data.array*self.target_value/f return out diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 33afc7e330..288b3a4825 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3101,7 +3101,7 @@ def setUp(self): self.data_multichannel = ag.allocate('random') self.data_slice = self.data_parallel.get_slice(vertical=1) self.data_reorder = self.data_cone.copy() - self.data_reorder.reorder(['horizontal','angle','vertical']) + self.data_reorder.reorder(['angle','horizontal','vertical']) self.data_single_angle = self.data_cone.get_slice(angle=1) def error_message(self,processor, test_parameter): @@ -3110,8 +3110,8 @@ def error_message(self,processor, test_parameter): def test_init(self): # test default values are initialised processor = FluxNormaliser() - test_parameter = ['flux','roi','norm_value'] - test_value = [None, None, None] + test_parameter = ['flux','roi','target'] + test_value = [None, None, 'mean'] 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])) @@ -3129,30 +3129,36 @@ def test_check_input(self): with self.assertRaises(ValueError): processor.check_input(self.data_cone) + 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]) + processor.set_input(self.data_cone) with self.assertRaises(ValueError): - processor.check_input(self.data_cone) + processor._calculate_flux() # check there is an error if roi is not specified as a dictionary processor = FluxNormaliser(roi='string') + processor.set_input(self.data_cone) with self.assertRaises(TypeError): - processor.check_input(self.data_cone) + processor._calculate_flux() # check there is an error if roi is specified with float values processor = FluxNormaliser(roi={'horizontal':(1.5, 6.5)}) + processor.set_input(self.data_cone) with self.assertRaises(TypeError): - processor.check_input(self.data_cone) + processor._calculate_flux() # check there is an error if roi stop is greater than start processor = FluxNormaliser(roi={'horizontal':(10, 5)}) + processor.set_input(self.data_cone) with self.assertRaises(ValueError): - processor.check_input(self.data_cone) + processor._calculate_flux() # check there is an error if roi stop is greater than the size of the axis processor = FluxNormaliser(roi={'horizontal':(0, self.data_cone.get_dimension_size('horizontal')+1)}) + processor.set_input(self.data_cone) with self.assertRaises(ValueError): - processor.check_input(self.data_cone) + processor._calculate_flux() @patch("matplotlib.pyplot.figure") def test_preview_configuration(self, mock_plot): @@ -3165,7 +3171,8 @@ def test_preview_configuration(self, mock_plot): # Test error in preview configuration if set_input not called roi = {'horizontal':(25,40)} processor = FluxNormaliser(roi=roi) - with self.assertRaises(ValueError): + with self.assertRaises(TypeError): + processor.preview_configuration() # Test no error with preview_configuration with different data shapes @@ -3194,19 +3201,19 @@ def test_preview_configuration(self, mock_plot): processor.preview_configuration(channel=1) def test_FluxNormaliser(self): - #Test flux with no norm_value - processor = FluxNormaliser(flux=10) + #Test flux with no target + processor = FluxNormaliser(flux=1) processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) - #Test flux with norm_value - processor = FluxNormaliser(flux=10, norm_value=5) + #Test flux with target + processor = FluxNormaliser(flux=10, target=5.0) processor.set_input(self.data_cone) data_norm = processor.get_output() numpy.testing.assert_allclose(data_norm.array, 0.5*self.data_cone.array) - #Test flux array with no norm_value + #Test flux array with no target flux = numpy.arange(1,2,(2-1)/(self.data_cone.get_dimension_size('angle'))) processor = FluxNormaliser(flux=flux) processor.set_input(self.data_cone) @@ -3214,13 +3221,13 @@ def test_FluxNormaliser(self): data_norm_test = self.data_cone.copy() for a in range(data_norm_test.get_dimension_size('angle')): data_norm_test.array[a,:,:] /= flux[a] - data_norm_test.array[a,:,:]*= numpy.mean(flux) + data_norm_test.array[a,:,:]*= numpy.mean(flux.ravel()) numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6) - # #Test flux array with norm_value + # #Test flux array with target flux = numpy.arange(1,2,(2-1)/(self.data_cone.get_dimension_size('angle'))) - norm_value = 5 - processor = FluxNormaliser(flux=flux, norm_value=norm_value) + norm_value = 5.0 + processor = FluxNormaliser(flux=flux, target=norm_value) processor.set_input(self.data_cone) data_norm = processor.get_output() data_norm_test = self.data_cone.copy() From da5f9775d4e4647e0eb44541b98c6b0aa365958d Mon Sep 17 00:00:00 2001 From: hrobarts Date: Wed, 16 Oct 2024 12:31:20 +0000 Subject: [PATCH 23/25] New tests for calculate_target --- .../Python/cil/processors/FluxNormaliser.py | 30 +--- Wrappers/Python/test/test_DataProcessor.py | 165 ++++++++++++------ 2 files changed, 118 insertions(+), 77 deletions(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 769f530e4e..3eb3c7ecf0 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -128,7 +128,6 @@ def check_input(self, dataset): return True - def _calculate_flux(self): ''' Function to calculate flux from a region of interest in the data. If the @@ -216,13 +215,16 @@ def _calculate_flux(self): if 0 in self.flux: raise ValueError('Flux value can\'t be 0, provide a different flux\ or region of interest with non-zero values') - - + def _calculate_target(self): ''' Calculate the target value for the normalisation ''' - if isinstance(self.target, float): + + if self.flux is None: + raise ValueError('Flux not found') + + if isinstance(self.target, (int,float)): self.target_value = self.target elif isinstance(self.target, str): if self.target == 'first': @@ -236,11 +238,9 @@ def _calculate_target(self): raise ValueError("Target string not recognised, found {}, expected 'first' or 'mean'" .format(self.target)) else: - raise TypeError("Target must be string or float, found {}" + raise TypeError("Target must be string or a number, found {}" .format(type(self.target))) - - def preview_configuration(self, angle=None, channel=None, log=False): ''' Preview the FluxNormalisation processor configuration for roi mode. @@ -406,18 +406,4 @@ def process(self, out=None): arr_proj *= self.target_value/f out.array.flat[i*proj_size:(i+1)*proj_size] = arr_proj - - # if 'angle' in data.dimension_labels: - # proj_axis = data.get_dimension_axis('angle') - # slice_proj = [slice(None)]*len(data.shape) - # slice_proj[proj_axis] = 0 - - # for i in range(len(data.geometry.angles)): - # if len(flux_size) > 0: - # f = self.flux[i] - # slice_proj[proj_axis] = i - # out.array[tuple(slice_proj)] = data.array[tuple(slice_proj)]*self.target_value/f - # else: - # out.array = data.array*self.target_value/f - - return out + return out \ No newline at end of file diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 288b3a4825..c75640ac86 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3160,6 +3160,59 @@ def test_calculate_flux(self): with self.assertRaises(ValueError): processor._calculate_flux() + def test_calculate_target(self): + + # check target calculated with default method 'mean' + processor = FluxNormaliser(flux=1) + processor.set_input(self.data_cone) + processor._calculate_flux() + processor._calculate_target() + self.assertAlmostEqual(processor.target_value, 1) + + processor = FluxNormaliser(flux=numpy.linspace(1,3,len(self.data_cone.geometry.angles))) + processor.set_input(self.data_cone) + processor._calculate_flux() + processor._calculate_target() + self.assertAlmostEqual(processor.target_value, 2) + + # check target calculated with method 'first' + processor = FluxNormaliser(flux=1, target='first') + processor.set_input(self.data_cone) + processor._calculate_flux() + processor._calculate_target() + self.assertAlmostEqual(processor.target_value, 1) + + processor = FluxNormaliser(flux=numpy.linspace(1,3,len(self.data_cone.geometry.angles)), + target='first') + processor.set_input(self.data_cone) + processor._calculate_flux() + processor._calculate_target() + self.assertAlmostEqual(processor.target_value, 1) + + # check target calculated with float + processor = FluxNormaliser(flux=1, + target=55.0) + processor.set_input(self.data_cone) + processor._calculate_flux() + processor._calculate_target() + self.assertAlmostEqual(processor.target_value, 55.0) + + # check error if target is an unrecognised string + processor = FluxNormaliser(flux=1, + target='string') + processor.set_input(self.data_cone) + processor._calculate_flux() + with self.assertRaises(ValueError): + processor._calculate_target() + + # check error if target is not a string or floar + processor = FluxNormaliser(flux=1, + target={'string': 10}) + processor.set_input(self.data_cone) + processor._calculate_flux() + with self.assertRaises(TypeError): + processor._calculate_target() + @patch("matplotlib.pyplot.figure") def test_preview_configuration(self, mock_plot): # Test error in preview configuration if there is no roi @@ -3236,67 +3289,69 @@ def test_FluxNormaliser(self): data_norm_test.array[a,:,:]*= norm_value numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6) - # #Test roi with no norm_value - # roi = {'vertical':(0,10), 'horizontal':(0,10)} - # processor = FluxNormaliser(roi=roi) - # processor.set_input(self.data_cone) - # data_norm = processor.get_output() - # numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) + # #Test roi with no target + roi = {'vertical':(0,10), 'horizontal':(0,10)} + processor = FluxNormaliser(roi=roi) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, self.data_cone.array) # #Test roi with norm_value - # roi = {'vertical':(0,10), 'horizontal':(0,10)} - # processor = FluxNormaliser(roi=roi, norm_value=5) - # processor.set_input(self.data_cone) - # data_norm = processor.get_output() - # numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) + roi = {'vertical':(0,10), 'horizontal':(0,10)} + processor = FluxNormaliser(roi=roi, target=5.0) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) # # Test roi with just one dimension - # roi = {'vertical':(0,2)} - # processor = FluxNormaliser(roi=roi, norm_value=5) - # processor.set_input(self.data_cone) - # data_norm = processor.get_output() - # numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) + roi = {'vertical':(0,2)} + processor = FluxNormaliser(roi=roi, target=5) + processor.set_input(self.data_cone) + data_norm = processor.get_output() + numpy.testing.assert_allclose(data_norm.array, 5*self.data_cone.array) # test roi with different data shapes and different flux values per projection - # for data in [self.data_cone, self.data_parallel, self.data_multichannel, - # self.data_slice, self.data_reorder]: - # roi = {'horizontal':(25,40)} - # processor = FluxNormaliser(roi=roi, norm_value=5) - # processor.set_input(data) - # data_norm = processor.get_output() - - # ax = data.get_dimension_axis('horizontal') - # slc = [slice(None)]*len(data.shape) - # slc[ax] = slice(25,40) - # axes=[ax] - # if 'vertical' in data.dimension_labels: - # axes.append(data.get_dimension_axis('vertical')) - # if 'channel' in data.dimension_labels: - # axes.append(data.get_dimension_axis('channel')) - # flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) - # slice_proj = [slice(None)]*len(data.shape) - # proj_axis = data.get_dimension_axis('angle') - # data_norm_test = data.copy() - # for i in range(len(data.geometry.angles)): - # f = flux[i] - # slice_proj[proj_axis] = i - # with numpy.errstate(divide='ignore', invalid='ignore'): - # data_norm_test.array[tuple(slice_proj)] = 5/f*data.array[tuple(slice_proj)] - # numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6, - # err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) - - # data = self.data_single_angle - # processor = FluxNormaliser(roi=roi, norm_value=5) - # processor.set_input(data) - # data_norm = processor.get_output() - # ax = data.get_dimension_axis('horizontal') - # slc = [slice(None)]*len(data.shape) - # slc[ax] = slice(25,40) - # axes=[ax,data.get_dimension_axis('vertical')] - # flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) - - # numpy.testing.assert_allclose(data_norm.array, 5/flux*data.array, atol=1e-6, - # err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) + for data in [ self.data_cone, self.data_parallel, self.data_multichannel, + self.data_slice, self.data_reorder]: + roi = {'horizontal':(25,40)} + processor = FluxNormaliser(roi=roi, target=5) + processor.set_input(data) + data_norm = processor.get_output() + + ax = data.get_dimension_axis('horizontal') + slc = [slice(None)]*len(data.shape) + slc[ax] = slice(25,40) + axes=[ax] + if 'vertical' in data.dimension_labels: + axes.append(data.get_dimension_axis('vertical')) + flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) + slice_proj = [slice(None)]*len(data.shape) + proj_axis = data.get_dimension_axis('angle') + data_norm_test = data.copy() + h_size = data.get_dimension_size('horizontal') + if 'vertical' in data.dimension_labels: + v_size = data.get_dimension_size('vertical') + else: + v_size = 1 + proj_size = h_size*v_size + for i in range(len(data.geometry.angles)*data.geometry.channels): + data_norm_test.array.flat[i*proj_size:(i+1)*proj_size] /=flux.flat[i] + data_norm_test.array.flat[i*proj_size:(i+1)*proj_size] *=5 + numpy.testing.assert_allclose(data_norm.array, data_norm_test.array, atol=1e-6, + err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) + + data = self.data_single_angle + processor = FluxNormaliser(roi=roi, target=5) + processor.set_input(data) + data_norm = processor.get_output() + ax = data.get_dimension_axis('horizontal') + slc = [slice(None)]*len(data.shape) + slc[ax] = slice(25,40) + axes=[ax,data.get_dimension_axis('vertical')] + flux = numpy.mean(data.array[tuple(slc)], axis=tuple(axes)) + + numpy.testing.assert_allclose(data_norm.array, 5/flux*data.array, atol=1e-6, + err_msg='Flux Normaliser roi test failed with data shape: ' + str(data.shape) + ' and configuration:\n' + str(data.geometry.config.system)) if __name__ == "__main__": From 5780dac5983be6ad8120de2bf4b61f074e9f67bb Mon Sep 17 00:00:00 2001 From: hrobarts Date: Wed, 16 Oct 2024 14:00:26 +0000 Subject: [PATCH 24/25] Remove unused import --- Wrappers/Python/cil/processors/FluxNormaliser.py | 1 - 1 file changed, 1 deletion(-) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 3eb3c7ecf0..00fffc3f44 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -17,7 +17,6 @@ # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt from cil.framework import Processor, AcquisitionData -from cil.utilities.display import show2D import numpy import logging import matplotlib.pyplot as plt From eb9db897a4fe8b84546f9625cbb67e1798860a9e Mon Sep 17 00:00:00 2001 From: hrobarts Date: Fri, 22 Nov 2024 14:00:39 +0000 Subject: [PATCH 25/25] Review updates --- CHANGELOG.md | 2 +- .../Python/cil/processors/FluxNormaliser.py | 22 +++++++++---------- Wrappers/Python/test/test_DataProcessor.py | 21 ++++++++++++++++++ 3 files changed, 33 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0624c39793..42f77d246b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ - Fix bug with 'median' and 'mean' methods in Masker averaging over the wrong axes. - Enhancements: - Removed multiple exits from numba implementation of KullbackLeibler divergence (#1901) + - Add FluxNormaliser processor (#1878) - Dependencies: - Added scikit-image to CIL-Demos conda install command as needed for new Callbacks notebook. @@ -16,7 +17,6 @@ - Added callback `optimisation.utilities.callbacks.EarlyStoppingObjectiveValue` which stops iterations if an algorithm objective changes less than a provided threshold (#1892) - Added callback `optimisation.utilities.callbacks.CGLSEarlyStopping` which replicates the automatic behaviour of CGLS in CIL versions <=24. (#1892) - Added `labels` module with `ImageDimension`, `AcquisitionDimension`, `AcquisitionType`, `AngleUnit`, `FillType` (#1692) - - Add FluxNormaliser processor (#1878) - Enhancements: - Use ravel instead of flat in KullbackLeibler numba backend (#1874) - Upgrade Python wrapper (#1873, #1875) diff --git a/Wrappers/Python/cil/processors/FluxNormaliser.py b/Wrappers/Python/cil/processors/FluxNormaliser.py index 00fffc3f44..5b9477fa0e 100644 --- a/Wrappers/Python/cil/processors/FluxNormaliser.py +++ b/Wrappers/Python/cil/processors/FluxNormaliser.py @@ -32,22 +32,22 @@ class FluxNormaliser(Processor): Parameters: ----------- - flux: float or list of floats (optional) + flux: float or list of floats, optional The value to divide the image by. If flux is a list it must have length equal to the number of projections in the dataset. If flux=None, calculate - flux from the roi, default is None. + flux from the roi. - roi: dict (optional) + roi: dict, optional Dictionary describing the region of interest containing the background in the image. The roi is specified as `{'axis_name1':(start,stop), 'axis_name2':(start,stop)}`, where the key is the axis name 'vertical' and/ or 'horizontal'. If an axis is not specified in the roi dictionary, - the full range will be used, default is None. + the full range will be used. - target: string or float + target: {'mean', 'first'} or float, default='mean' The value to scale the normalised data with. If float, the data is scaled to the float value. If string 'mean' the data is scaled to the mean value - of the input flux or flux across all rois, if 'first' the data is scaled to + of the input flux or flux across all rois. If 'first' the data is scaled to the first input flux value or the flux in the roi of the first projection. Default is 'mean' @@ -72,7 +72,7 @@ class FluxNormaliser(Processor): Example ------- >>> from cil.processors import FluxNormaliser - >>> processor = FluxNormaliser(roi=(roi={'horizontal':(5, 15)) + >>> processor = FluxNormaliser(roi={'horizontal':(5, 15)}) >>> processor.set_input(data) >>> data_norm = processor.get_output() @@ -121,7 +121,7 @@ def check_input(self, dataset): self.h_size = dataset.get_dimension_size('horizontal') image_axes += 1 - if (( self.h_axis is not None) and (self. h_axis < (len(dataset.shape)-image_axes))) or \ + if (( self.h_axis is not None) and (self.h_axis < (len(dataset.shape)-image_axes))) or \ ((self.v_axis is not None) and self.v_axis < (len(dataset.shape)-image_axes)): raise ValueError('Projections must be the last two axes of the dataset') @@ -248,17 +248,17 @@ def preview_configuration(self, angle=None, channel=None, log=False): Parameters: ----------- - angle: float or str (optional) + angle: float, optional Angle to plot, default=None displays the data with the minimum and maximum pixel values in the roi, otherwise the angle to display can be specified as a float and the closest angle will be displayed. For 2D data, the roi is plotted on the sinogram. - channel: int (optional) + channel: int, optional The channel to plot, default=None displays the central channel if the data has channels - log: bool (optional) + log: bool, default=False If True, plot the image with a log scale, default is False ''' self._calculate_flux() diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index c51142cb25..9de8a21649 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -3205,6 +3205,27 @@ def test_calculate_flux(self): with self.assertRaises(ValueError): processor._calculate_flux() + # check error raised with 0 flux + processor = FluxNormaliser(flux = 0) + processor.set_input(self.data_cone) + with self.assertRaises(ValueError): + processor._calculate_flux() + + processor = FluxNormaliser(flux = 0.0) + processor.set_input(self.data_cone) + with self.assertRaises(ValueError): + processor._calculate_flux() + + processor = FluxNormaliser(flux=numpy.zeros(len(self.data_cone.geometry.angles))) + processor.set_input(self.data_cone) + with self.assertRaises(ValueError): + processor._calculate_flux() + + processor = FluxNormaliser(flux=numpy.zeros(len(self.data_cone.geometry.angles), dtype=numpy.uint16)) + processor.set_input(self.data_cone) + with self.assertRaises(ValueError): + processor._calculate_flux() + def test_calculate_target(self): # check target calculated with default method 'mean'