diff --git a/CHANGELOG.md b/CHANGELOG.md index 63e12e6484..075d199617 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,7 @@ * XX.X.X - New features: - LSQR algorithm added to the CIL algorithm class (#1975) + - LaminographyGeometryCorrector tool added to processors (#2259) - Bug fixes: - `CentreOfRotationCorrector.image_sharpness` data is now correctly smoothed to reduce aliasing artefacts and improve robustness. (#2202) - `PaganinProcessor` now correctly applies scaling with magnification for cone-beam geometry (#2225) diff --git a/Wrappers/Python/cil/processors/LaminographyGeometryCorrector.py b/Wrappers/Python/cil/processors/LaminographyGeometryCorrector.py new file mode 100644 index 0000000000..54a455367d --- /dev/null +++ b/Wrappers/Python/cil/processors/LaminographyGeometryCorrector.py @@ -0,0 +1,411 @@ +# Copyright 2026 United Kingdom Research and Innovation +# Copyright 2026 The University of Manchester +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# Authors: +# CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt +import numpy as np +from scipy.spatial.transform import Rotation as R +from scipy.optimize import minimize +from scipy.ndimage import gaussian_filter, sobel +import matplotlib.pyplot as plt +import numpy as np +import importlib + +from cil.framework import Processor +from cil.processors import Binner, Slicer +from cil.framework import AcquisitionData +from cil.framework.labels import AcquisitionType, AcquisitionDimension + +import logging +log = logging.getLogger(__name__) + + + +class LaminographyGeometryCorrector(Processor): + """ + LaminographyGeometryCorrector processor to fit the geometry of a + parallel beam laminography dataset to find tilt and center-of-rotation. + + Parameters + ---------- + parameter_bounds : list of tuple of float, optional + Bounds for the parameters [(tilt_min_deg, tilt_max_deg), (CoR_min_pix, CoR_max_pix)]. + Defaults to [(-10, 10), (-20, 20)]. + + parameter_tolerance : tuple of float, optional + Convergence tolerance for optimisation of parameters, (tilt_tol_deg, CoR_tol_pix). + Defaults to (0.01, 0.01). + + coarse_binning : int, optional + Initial binning factor applied to the input dataset for coarse optimisation. + If None, a value based on dataset size is used. + + final_binning : int, optional + Final binning factor applied for fine optimisation. + If None, no binning is applied in the fine optimisation step. + + angle_subsampling : float, optional + Subsampling factor for the angle dimension during optimisation. + If None, automatically determined based on input dataset. + + image_geometry : ImageGeometry, optional + Pass a reduced volume ImageGeometry to be used for fitting. + If None, the full dataset is used. + + backend : {'astra'}, optional + The backend to use for the reconstruction. Currently only 'astra' + is supported. + + + Example + ------- + + >>> processor = LaminographyGeometryCorrector(parameter_bounds=(tilt_bounds, CoR_bounds), + parameter_tolerance=(tilt_tol, CoR_tol)) + >>> processor.set_input(data) + >>> data_corrected = processor.get_output() + + """ + + _supported_backends = ['astra'] + + def __init__(self, parameter_bounds=[(-10, 10),(-20, 20)], parameter_tolerance=(0.01, 0.01), + coarse_binning=None, final_binning = None, angle_subsampling = None, image_geometry = None, backend='astra'): + + + FBP, ProjectionOperator = self._configure_FBP(backend) + kwargs = { + 'initial_parameters' : None, + 'parameter_bounds' : parameter_bounds, + 'parameter_tolerance' : parameter_tolerance, + 'image_geometry' : image_geometry, + 'coarse_binning' : coarse_binning, + 'final_binning' : final_binning, + 'angle_subsampling' : angle_subsampling, + 'backend' : backend, + 'FBP' : FBP, + 'ProjectionOperator' : ProjectionOperator, + 'evaluations' : [] + } + super(LaminographyGeometryCorrector, self).__init__(**kwargs) + + def check_input(self, dataset): + if not isinstance(dataset, (AcquisitionData)): + raise TypeError('Processor only supports AcquisitionData') + + if dataset.geometry.geom_type & AcquisitionType.CONE_FLEX: + raise NotImplementedError("Processor not implemented for CONE_FLEX geometry.") + + if dataset.geometry.geom_type & AcquisitionType.CONE: + raise NotImplementedError("LaminographyGeometryCorrector does not yet support CONE data") + + if not AcquisitionDimension.check_order_for_engine('astra', dataset.geometry): + raise ValueError("LaminographyGeometryCorrector must be used with astra data order, try `data.reorder('astra')`") + + if not dataset.geometry.dimension & AcquisitionType.DIM3: + raise ValueError("LaminographyGeometryCorrector must be used 3D data") + + return True + + def _calculate_initial_parameters(self): + """ + Get the current tilt and centre of rotation from the geometry + """ + # get initial parameters from geometry + dataset = self.get_input() + U = dataset.geometry.config.system.rotation_axis.direction + V = dataset.geometry.config.system.detector.direction_y + c = np.cross(U, V) + d = np.dot(U, V) + c_norm = np.linalg.norm(c) + tilt_deg = np.rad2deg(np.arctan2(c_norm, d)) + CoR_pix = dataset.geometry.get_centre_of_rotation('pixels')['offset'][0] + self.initial_parameters = (tilt_deg, CoR_pix) + + def _update_geometry(self, ag, tilt_deg, cor_pix, + tilt_direction_vector = np.array([1, 0, 0]), + original_rotation_axis=np.array([0, 0, 1])): + """ + Update the rotation matrix direction and centre of rotation from a tilt + in degrees and centre of rotation offset in pixels + """ + tilt_rad = np.deg2rad(tilt_deg) + rotation_matrix = R.from_rotvec(tilt_rad * tilt_direction_vector) + tilted_rotation_axis = rotation_matrix.apply(original_rotation_axis) + + ag.set_centre_of_rotation(offset=cor_pix, distance_units='pixels') + ag.config.system.rotation_axis.direction = tilted_rotation_axis + + return ag + + def _projection_reprojection(self, data, recon_buffer, ig, ag, ag_downsampled, data_downsampled, residual_buffer, tilt_deg, cor_pix): + """ + Reconstruct the data then re-project and calculate the residual. Then + filter the residual and calculate the L2Norm loss. + + Parameters + ---------- + data: AcquisitionData + The full size dataset + recon_buffer: ImageData + Pre-allocated buffer for the reconstruction volume + ig: ImageGeometry + Reconstruction volume geometry + ag: AcquisitionGeometry + Copy of full data geometry + ag_downsampled: AcquisitionGeometry + Geometry of the downsampled dataset used for reprojection + data_downsampled: AcquisitionData + Downsampled dataset used for reprojection + residual_buffer: AcquisitionData + Pre-allocated buffer, size of the downsampled data + tilt_deg: float + Latest tilt angle in degrees + cor_pix: float + Latest centre of rotation offset in pixels + + """ + + # update the geometry with latest values and get reconstruction + ag = self._update_geometry(ag, tilt_deg, cor_pix) + FBP = self.FBP(ig, ag) + FBP.set_input(data) + FBP.get_output(out=recon_buffer) + recon_buffer.apply_circular_mask(0.9) + + # update the downsampled data geometry and get forward projection + ag_downsampled = self._update_geometry(ag_downsampled, tilt_deg, cor_pix) + A = self.ProjectionOperator(ig, ag_downsampled) + A.direct(recon_buffer, out=residual_buffer) + # subtract the downsampled reference data + residual_buffer.subtract(data_downsampled, out=residual_buffer) + + # apply Gaussian and Sobel filter - note the axes are hard coded here for astra, this would need to be updated for tigre + residual_buffer.subtract(gaussian_filter(residual_buffer.as_array(), sigma=3.0, axes=(0,2)), out=residual_buffer) + np.sqrt((sobel(residual_buffer.array, axis=0))**2 + (sobel(residual_buffer.array, axis=2))**2, out=residual_buffer.array) + + loss = float(np.sum(residual_buffer**2)) + + return loss + + def _minimise_geometry(self, data, binning, p0, bounds): + """ + Setup and run the scipy Powell minimize method + + Parameters + ---------- + data: AcquistionData + Full dataset + binning: int + Current detector binning value + p0: list or tuple + Initial start parameters for search (tilt_deg, cor_pix) + bounds: list of tuple of floats + Bounds for the parameters [(tilt_min_deg, tilt_max_deg), (CoR_min_pix, CoR_max_pix)] + """ + + current_run_evaluations = [] + xtol = self.parameter_tolerance + + # scale the start values and bounds by the binning + p0_binned = (p0[0], p0[1]/binning) + bounds_binned = (bounds[0], (bounds[1][0]/binning, bounds[1][1]/binning)) + + # scale the start values and bounds so xtol can be 1 + p0_scaled = np.array([p0_binned[0] / xtol[0], + p0_binned[1] / xtol[1]], dtype=float) + + bounds_scaled = [(bounds_binned[0][0] / xtol[0], bounds_binned[0][1] / xtol[0]), + (bounds_binned[1][0] / xtol[1], bounds_binned[1][1] / xtol[1])] + + direc = np.diag(np.asarray(xtol) / np.min(xtol)) + + # get y_ref: a subset of the real data to compare with the reprojections + target = max(np.ceil(data.get_dimension_size('angle') / 10), 36) + divider = np.floor(data.get_dimension_size('angle') / target) + data_downsampled = Slicer(roi={'angle':(None, None, divider)})(data) + + # also get a matching reference geometry + ag = data.geometry.copy() + ag_downsampled = Slicer(roi={'angle':(None, None, divider)})(ag) + + if self.image_geometry is None: + ig = ag.get_ImageGeometry() + else: + ig = Binner(roi={'horizontal_x':(None, None,binning), 'horizontal_y':(None, None,binning), 'vertical':(None, None,binning)})(self.image_geometry) + + # pre-allocate reconstruction volume and residual array + recon = ig.allocate(0) + residual = ag_downsampled.allocate() + + def loss_function_wrapper(p): + """ + Function wrapper for the loss function, self._projection_reprojection + to be called by scipy.minimize. Rescales tilt and cor by xtol so a + single tolerance can be used for each parameter. + """ + tilt = p[0] * xtol[0] + cor = p[1] * xtol[1] + loss = self._projection_reprojection(data, recon, ig, ag, ag_downsampled, data_downsampled, residual, tilt, cor) + + current_run_evaluations.append((tilt, cor * binning, loss)) + if log.isEnabledFor(logging.DEBUG): + print(f"tilt: {tilt:.3f}, CoR: {cor*binning:.3f}, loss: {loss:.3e}") + + return loss + + # call minimize + res_scaled = minimize(loss_function_wrapper, p0_scaled, + method='Powell', + bounds=bounds_scaled, + options={'maxiter': 5, 'disp': False, 'xtol': 1.0, 'direc': direc}) + + # re-scale the results + res_real = res_scaled + res_real.x = np.array([res_scaled.x[0] * xtol[0], + res_scaled.x[1] * xtol[1] * binning]) + + # save information about the minimisation + self.evaluations.append({ + "p0": p0, + "bounds": bounds, + "binning": binning, + "xtol": xtol, + "result": res_real, + "evaluations": current_run_evaluations + }) + + return res_real + + def process(self, out=None): + + data = self.get_input() + self._calculate_initial_parameters() + + # apply coarse binning to the data + if self.coarse_binning is None: + # if no coarse binning provided, get a default binning based on the size of the panel + self.coarse_binning = min(int(np.ceil(data.geometry.config.panel.num_pixels[0] / 256)),5) + binning = self.coarse_binning + roi = { + 'horizontal': (None, None, binning), + 'vertical': (None, None, binning) + } + data_binned = Binner(roi)(data) + + # sub-sample the angles + if self.angle_subsampling is None: + # if no sub-sampling value is provided, get a default subsampling based on the Nyquist criteria + self.angle_subsampling = np.ceil(data.get_dimension_size('angle')/(data.get_dimension_size('horizontal')*(np.pi/2))) + roi={'angle':(None, None, self.angle_subsampling*binning)} + data_binned = Slicer(roi)(data_binned) + + # run coarse minimisation + coarse_tolerance = (self.parameter_tolerance[0], self.parameter_tolerance[1]) + res = self._minimise_geometry(data_binned, binning=binning, + p0=self.initial_parameters, + bounds=self.parameter_bounds) + + tilt_min = res.x[0] + cor_min = res.x[1] + if log.isEnabledFor(logging.DEBUG): + print(f"Coarse scan optimised tilt = {tilt_min:.3f}, CoR = {cor_min:.3f}") + + # apply final binning + if self.final_binning is None: + binning = 1 + else: + binning = self.final_binning + roi = { + 'horizontal': (None, None, binning), + 'vertical': (None, None, binning), + 'angle': (None, None, self.angle_subsampling) + } + + data_binned = Binner(roi)(data) + + # calculate new search ranges based on coarse minimisation results + search_factor = 2 # multiplier on parameter_tolerance + min_search_range_tilt = 1.0 + min_search_range_cor = 1.0 + + half_width_tilt = max(search_factor * coarse_tolerance[0], min_search_range_tilt/2) + fine_bounds_tilt = (max((tilt_min - half_width_tilt), self.parameter_bounds[0][0]), # whichever is higher out of the calculated fine lower bound and coarse bound + min((tilt_min + half_width_tilt), self.parameter_bounds[0][1])) # whichever is lower out of the calculated fine upper bound and coarse bound + + half_width_cor = max(search_factor * coarse_tolerance[1], min_search_range_cor/2) + fine_bounds_cor = (max((cor_min - half_width_cor), self.parameter_bounds[1][0]), # whichever is higher out of the calculated fine lower bound and coarse bound + min((cor_min + half_width_cor), self.parameter_bounds[1][1])) # whichever is lower out of the calculated fine upper bound and coarse bound + + # run fine minimisation + res = self._minimise_geometry(data_binned, binning=binning, + p0=(tilt_min, cor_min), + bounds=[fine_bounds_tilt, fine_bounds_cor]) + tilt_min = res.x[0] + cor_min = res.x[1] + print(f"Fine scan optimised tilt = {tilt_min:.3f}, CoR ={cor_min:.3f}") + + if log.isEnabledFor(logging.DEBUG): + self.plot_evaluations() + + new_geometry = data.geometry.copy() + self._update_geometry(new_geometry, tilt_min, cor_min) + + if out is None: + return AcquisitionData(array=data.as_array(), deep_copy=True, geometry=new_geometry) + else: + out.geometry = new_geometry + return out + + def plot_evaluations(self): + """ + Plot results from the minimisation. Plots the loss function value as a + function of tilt and centre of rotation offset position. + """ + num_evals = len(self.evaluations) + if num_evals > 0: + fig, axs = plt.subplots(nrows=1, ncols=num_evals, figsize=(14, 6)) + + for i in np.arange(num_evals): + eval = self.evaluations[i] + tilts = [t[0] for t in eval['evaluations']] + cors = [t[1] for t in eval['evaluations']] + losses = [t[2] for t in eval['evaluations']] + + ax = axs[i] + scatter = ax.scatter(tilts, cors, c=losses, cmap='viridis_r', s=100, edgecolors='k') + fig.colorbar(scatter, label='Loss value', ax=ax) + ax.set_xlabel('Tilt') + ax.set_ylabel('Cor') + ax.set_title('bounds = ({:.2f}:{:.2f}), ({:.2f}:{:.2f}), binning = {}, xtol = ({}, {}) \n result = ({:.3f}, {:.3f})' + .format(*eval['bounds'][0], *eval['bounds'][1], eval['binning'], *eval['xtol'], eval['result'].x[0], eval['result'].x[1])) + ax.grid() + plt.tight_layout() + else: + raise ValueError("No evaluation available to plot. Run processor.process() first.") + + def _configure_FBP(self, backend='astra'): + """ + Configures the recon and projection operator for the right engine. + """ + if backend not in self._supported_backends: + raise ValueError("Backend unsupported. Supported backends: {}".format(self._supported_backends)) + + module = importlib.import_module(f'cil.plugins.{backend}') + + return module.FBP, module.ProjectionOperator + + + diff --git a/Wrappers/Python/cil/processors/__init__.py b/Wrappers/Python/cil/processors/__init__.py index c86bb50cb7..6a81d00cf3 100644 --- a/Wrappers/Python/cil/processors/__init__.py +++ b/Wrappers/Python/cil/processors/__init__.py @@ -27,4 +27,5 @@ from .Masker import Masker from .Padder import Padder from .PaganinProcessor import PaganinProcessor -from .FluxNormaliser import FluxNormaliser \ No newline at end of file +from .FluxNormaliser import FluxNormaliser +from .LaminographyGeometryCorrector import LaminographyGeometryCorrector \ No newline at end of file diff --git a/Wrappers/Python/test/test_DataProcessor.py b/Wrappers/Python/test/test_DataProcessor.py index 32914873c2..c6dfe3d497 100644 --- a/Wrappers/Python/test/test_DataProcessor.py +++ b/Wrappers/Python/test/test_DataProcessor.py @@ -36,7 +36,7 @@ from cil.processors.CofR_xcorrelation import CofR_xcorrelation from cil.processors.CofR_image_sharpness import CofR_image_sharpness from cil.processors import TransmissionAbsorptionConverter, AbsorptionTransmissionConverter -from cil.processors import Slicer, Binner, MaskGenerator, Masker, Padder, PaganinProcessor, FluxNormaliser, Normaliser +from cil.processors import Slicer, Binner, MaskGenerator, Masker, Padder, PaganinProcessor, FluxNormaliser, Normaliser, LaminographyGeometryCorrector import gc from utils import has_numba @@ -45,6 +45,7 @@ from scipy import constants from scipy.fft import ifftshift +from scipy.spatial.transform import Rotation as R from utils import has_astra, has_tigre, has_nvidia, has_tomophantom, initialise_tests, has_ipp, has_matplotlib @@ -3261,6 +3262,313 @@ def test_PaganinProcessor_1angle(self): sys.stderr = sys.__stderr__ self.assertLessEqual(quality_measures.mse(output, thickness), 0.05) +class TestLaminographyGeometryCorrector(unittest.TestCase): + + def setUp(self): + self.data_parallel = dataexample.SIMULATED_PARALLEL_BEAM_DATA.get() + self.data_parallel.reorder('astra') + + # Additional test data + self.data_cone = dataexample.SIMULATED_CONE_BEAM_DATA.get() + self.data_cone.reorder('astra') + + # Create cone flex geometry data + source_position_set = [[0, -100000, 0]] * 3 + detector_position_set = [[0, 0, 0]] * 3 + detector_direction_x_set = [[1, 0, 0]] * 3 + detector_direction_y_set = [[0, 0, 1]] * 3 + cone_flex_ag = AcquisitionGeometry.create_Cone3D_Flex( + source_position_set, detector_position_set, detector_direction_x_set, detector_direction_y_set + ).set_panel([32, 32]) + arr = numpy.random.rand(3, 32, 32).astype(numpy.float32) + self.data_cone_flex = AcquisitionData(arr, geometry=cone_flex_ag) + + def error_message(self, processor, test_parameter): + return "Failed with processor " + str(processor) + " on test parameter " + test_parameter + + @unittest.skipUnless(has_astra, "Astra not installed") + def test_LaminographyGeometryCorrector_init(self): + # test default values are initialised + processor = LaminographyGeometryCorrector() + test_parameter = ['parameter_bounds', 'parameter_tolerance', + 'coarse_binning', 'final_binning', 'angle_subsampling', 'image_geometry', 'evaluations'] + test_value = [[(-10, 10), (-20, 20)], (0.01, 0.01), + None, 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 = LaminographyGeometryCorrector( + parameter_bounds=[(20, 35), (-5, 15)], + parameter_tolerance=(0.1, 0.1), + coarse_binning=2, + final_binning=1, + angle_subsampling=2, + image_geometry=None) + test_value = [[(20, 35), (-5, 15)], (0.1, 0.1), + 2, 1, 2, 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])) + + @unittest.skipUnless(has_astra, "Astra not installed") + def test_LaminographyGeometryCorrector_check_input(self): + processor = LaminographyGeometryCorrector() + + # test with parallel beam data - should work + processor.set_input(self.data_parallel) + data2 = processor.get_input() + numpy.testing.assert_allclose(data2.as_array(), self.data_parallel.as_array()) + + # check there is an error when the wrong data type is input + with self.assertRaises(TypeError): + processor.set_input(self.data_parallel.geometry) + + with self.assertRaises(TypeError): + processor.set_input(self.data_parallel.as_array()) + + dc = DataContainer(self.data_parallel.as_array()) + with self.assertRaises(TypeError): + processor.set_input(dc) + + # check with different data order - should raise error + data_reorder = self.data_parallel.copy() + data_reorder.reorder('astra') + data_reorder.reorder(['angle','horizontal','vertical']) + with self.assertRaises(ValueError): + processor.set_input(data_reorder) + + # Test that cone beam data raises NotImplementedError + processor = LaminographyGeometryCorrector() + with self.assertRaises(NotImplementedError): + processor.set_input(self.data_cone) + processor.check_input(self.data_cone) + + # Test that cone flex geometry raises NotImplementedError + processor = LaminographyGeometryCorrector() + with self.assertRaises(NotImplementedError): + processor.set_input(self.data_cone_flex) + + @unittest.skipUnless(has_astra, "Astra not installed") + def test_LaminographyGeometryCorrector_update_geometry(self): + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + ag.set_angles(numpy.linspace(0, numpy.pi, 10)) + ag.set_panel([512, 512], pixel_size=(1.0, 1.0)) + + tilt_deg = 35.0 + cor_pix = 5.0 + + ag_updated = processor._update_geometry(ag, tilt_deg, cor_pix) + + # Verify CoR was updated + self.assertAlmostEqual(ag_updated.config.system.rotation_axis.position[0], cor_pix) + + # Verify rotation axis direction was modified (should be tilted) + original_axis = numpy.array([0, 0, 1]) + tilted_axis = ag_updated.config.system.rotation_axis.direction + assert not numpy.allclose(tilted_axis, original_axis) + + # Test tilt around x-axis (default) + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + ag.set_angles(numpy.linspace(0, numpy.pi, 10)) + ag.set_panel([512, 512], pixel_size=(1.0, 1.0)) + + tilt_deg = 30.0 + cor_pix = 0.0 + tilt_x = numpy.array([1, 0, 0]) + + ag_updated_x = processor._update_geometry(ag, tilt_deg, cor_pix, tilt_direction_vector=tilt_x) + + updated_axis = ag_updated_x.config.system.rotation_axis.direction + # When tilting around x-axis, z-component should change, x-component should stay ~1 + self.assertAlmostEqual(updated_axis[0], 0, places=5) + + # Test tilt around y-axis + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + ag.set_angles(numpy.linspace(0, numpy.pi, 10)) + ag.set_panel([512, 512], pixel_size=(1.0, 1.0)) + + tilt_deg = 30.0 + cor_pix = 0.0 + tilt_y = numpy.array([0, 1, 0]) + + ag_updated_y = processor._update_geometry(ag, tilt_deg, cor_pix, tilt_direction_vector=tilt_y) + + updated_axis = ag_updated_y.config.system.rotation_axis.direction + # When tilting around y-axis, x-component should change + original_axis = numpy.array([0, 0, 1]) + self.assertFalse(numpy.allclose(updated_axis, original_axis)) + + # Test tilt around z-axis + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + ag.set_angles(numpy.linspace(0, numpy.pi, 10)) + ag.set_panel([512, 512], pixel_size=(1.0, 1.0)) + + tilt_deg = 45.0 + cor_pix = 0.0 + tilt_z = numpy.array([0, 0, 1]) + + ag_updated_z = processor._update_geometry(ag, tilt_deg, cor_pix, tilt_direction_vector=tilt_z) + + # Tilting around the rotation axis itself should not change the direction + original_axis = numpy.array([0, 0, 1]) + updated_axis = ag_updated_z.config.system.rotation_axis.direction + numpy.testing.assert_allclose(updated_axis, original_axis, atol=1e-5) + + # Test with zero tilt + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + ag.set_angles(numpy.linspace(0, numpy.pi, 10)) + ag.set_panel([512, 512], pixel_size=(1.0, 1.0)) + + tilt_deg = 0.0 + cor_pix = 5.0 + + ag_updated = processor._update_geometry(ag, tilt_deg, cor_pix) + + # With zero tilt, rotation axis should remain unchanged + original_axis = numpy.array([0, 0, 1]) + updated_axis = ag_updated.config.system.rotation_axis.direction + numpy.testing.assert_allclose(updated_axis, original_axis, atol=1e-5) + # CoR should still be updated + self.assertAlmostEqual(ag_updated.config.system.rotation_axis.position[0], cor_pix) + + # Test that panel size is preserved after geometry update + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + ag.set_angles(numpy.linspace(0, numpy.pi, 10)) + panel_size = (512, 512) + ag.set_panel(panel_size, pixel_size=(1.0, 1.0)) + + original_panel = ag.config.panel.num_pixels + + ag_updated = processor._update_geometry(ag, 35.0, 5.0) + updated_panel = ag_updated.config.panel.num_pixels + + # Panel size should be unchanged + self.assertEqual(original_panel[0], updated_panel[0]) + self.assertEqual(original_panel[1], updated_panel[1]) + + # Test that number of angles is preserved + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + num_angles = 15 + ag.set_angles(numpy.linspace(0, numpy.pi, num_angles)) + ag.set_panel([128, 128]) + + original_num_angles = len(ag.angles) + + ag_updated = processor._update_geometry(ag, 30.0, 3.0) + updated_num_angles = len(ag_updated.angles) + + # Number of angles should be preserved + self.assertEqual(original_num_angles, updated_num_angles) + + # Test that rotation axis direction remains normalized (unit vector) + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + ag.set_angles(numpy.linspace(0, numpy.pi, 10)) + ag.set_panel([256, 256]) + + tilt_deg = 45.0 + cor_pix = 10.0 + + ag_updated = processor._update_geometry(ag, tilt_deg, cor_pix) + + # Check that rotation axis is a unit vector + axis = ag_updated.config.system.rotation_axis.direction + axis_norm = numpy.linalg.norm(axis) + self.assertAlmostEqual(axis_norm, 1.0, places=5) + + # Test with extreme tilt angles + processor = LaminographyGeometryCorrector() + + ag = AcquisitionGeometry.create_Parallel3D() + ag.set_angles(numpy.linspace(0, numpy.pi, 10)) + ag.set_panel([128, 128]) + + original_panel = ag.config.panel.num_pixels + + # Test with 90 degree tilt + ag_updated_90 = processor._update_geometry(ag, 90.0, 0.0) + self.assertEqual(ag_updated_90.config.panel.num_pixels[0], original_panel[0]) + + # Test with 180 degree tilt + ag_updated_180 = processor._update_geometry(ag, 180.0, 0.0) + self.assertEqual(ag_updated_180.config.panel.num_pixels[0], original_panel[0]) + + @unittest.skipUnless(has_astra and has_nvidia, "ASTRA GPU not installed") + def test_LaminographyGeometryCorrector_process(self): + # get volume + vol = dataexample.SIMULATED_SPHERE_VOLUME.get() + vol.reorder('astra') + ig = vol.geometry + + # forward project with tilt and cor offset + tilt = 37.7 # degrees + cor_offset = -2.4 # pixels + tilt_direction = numpy.array([1, 0, 0]) + untilted_rotation_axis = numpy.array([0, 0, 1]) + + rotation_matrix = R.from_rotvec(numpy.deg2rad(tilt) * tilt_direction) + tilted_rotation_axis = rotation_matrix.apply(untilted_rotation_axis) + angles = numpy.arange(0,360,5) + ag = AcquisitionGeometry.create_Parallel3D(rotation_axis_direction=tilted_rotation_axis)\ + .set_angles(angles)\ + .set_panel(vol.shape[0:2], (vol.geometry.voxel_size_x, vol.geometry.voxel_size_y)) + ag.set_centre_of_rotation(cor_offset, distance_units='pixels') + ag.set_labels(['vertical', 'angle','horizontal']) + + A = AstraProjectionOperator(ig, ag) + proj = A.direct(vol) + + # update geometry with guess values + tilt_guess = 35 # degrees + cor_guess = 0 + rotation_matrix_guess = R.from_rotvec(numpy.deg2rad(tilt_guess) * tilt_direction) + tilted_rotation_axis_guess = rotation_matrix_guess.apply(untilted_rotation_axis) + proj.geometry.config.system.rotation_axis.direction = tilted_rotation_axis_guess + proj.geometry.set_centre_of_rotation(cor_guess, distance_units='pixels') + + # run LaminographyGeometryCorrector + tolerance = 0.5 # choose high tolerance for quick convergence + processor = LaminographyGeometryCorrector(parameter_bounds=[(30,40),(-5,5)], + parameter_tolerance=(tolerance, tolerance)) + processor.set_input(proj) + + sys.stdout = open(os.devnull, "w") + sys.stderr = open(os.devnull, "w") + corrected_proj = processor.get_output() + sys.stdout = sys.__stdout__ + sys.stderr = sys.__stderr__ + + # check the fitted values are within the tolerance of the true values + U = corrected_proj.geometry.config.system.rotation_axis.direction + V = corrected_proj.geometry.config.system.detector.direction_y + c = numpy.cross(U, V) + d = numpy.dot(U, V) + c_norm = numpy.linalg.norm(c) + tilt_fit = numpy.rad2deg(numpy.arctan2(c_norm, d)) + cor_fit = corrected_proj.geometry.get_centre_of_rotation(distance_units='pixels')['offset'][0] + + self.assertLessEqual(abs(tilt_fit - tilt), tolerance) + self.assertLessEqual(abs(cor_fit - cor_offset), tolerance) + class TestFluxNormaliser(unittest.TestCase): def setUp(self): diff --git a/docs/source/images/lamino_tilt1.png b/docs/source/images/lamino_tilt1.png new file mode 100644 index 0000000000..cc39cc1b4d Binary files /dev/null and b/docs/source/images/lamino_tilt1.png differ diff --git a/docs/source/images/lamino_tilt2.png b/docs/source/images/lamino_tilt2.png new file mode 100644 index 0000000000..06179d90fd Binary files /dev/null and b/docs/source/images/lamino_tilt2.png differ diff --git a/docs/source/processors.rst b/docs/source/processors.rst index 12b26f9c6b..6825bfc645 100644 --- a/docs/source/processors.rst +++ b/docs/source/processors.rst @@ -156,4 +156,37 @@ Paganin Processor :members: :inherited-members: set_input, get_output +Laminography Geometry Corrector +----------------------------- + +Laminography is a technique similar to tomography in which the rotation axis +is tilted relative to the beam. This can be setup as below, by tilting the +sample and rotation axis by an angle :math:`\phi`: + +.. figure:: images/lamino_tilt1.png + :align: center + :alt: laminography geometry showing tilted rotation axis + :figclass: align-center + + Laminography geometry from tilted rotation axis + +or equivalently by offsetting the source and detector: + +.. figure:: images/lamino_tilt2.png + :align: center + :alt: laminography geometry showing tilted rotation axis + :figclass: align-center + + Laminography geometry from offset source and detector + +This processor fits the tilt angle :math:`\phi` as well as the centre of +rotation offset (defined above) for parallel beam laminography data. + + +.. autoclass:: cil.processors.LaminographyGeometryCorrector + :exclude-members: check_input, get_input + :members: + :inherited-members: set_input, get_output + + :ref:`Return Home `