diff --git a/CHANGELOG.md b/CHANGELOG.md index 23bbdf595..aaca7e488 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,9 +1,12 @@ * XX.X.X - New features: - LSQR algorithm added to the CIL algorithm class (#1975) + - Added `get_centre_slice` method to `ImageData`, `ImageGeometry`, and `AcquisitionGeometry`. (#2235) + - Allows passing vertical='centre' to `ImageGeometry.get_slice` (#2235) - 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) + - Sets `center_x`, `center_y`, `center_z` appropriately for `ImageGeometry` returned by `get_slice` (#2235) * 25.0.0 - New features: diff --git a/Wrappers/Python/cil/framework/acquisition_data.py b/Wrappers/Python/cil/framework/acquisition_data.py index 308f0d546..406054361 100644 --- a/Wrappers/Python/cil/framework/acquisition_data.py +++ b/Wrappers/Python/cil/framework/acquisition_data.py @@ -230,6 +230,12 @@ def get_slice(self, *args, **kwargs): kwargs['force'] = args[4] return self._get_slice(**kwargs) + + def get_centre_slice(self): + ''' + Returns a new AcquisitionData of the centre slice in the vertical direction. + ''' + return self.get_slice(vertical='centre') def reorder(self, order): ''' diff --git a/Wrappers/Python/cil/framework/image_data.py b/Wrappers/Python/cil/framework/image_data.py index 67cfbff5b..13878022f 100644 --- a/Wrappers/Python/cil/framework/image_data.py +++ b/Wrappers/Python/cil/framework/image_data.py @@ -127,7 +127,7 @@ def __eq__(self, other): def get_slice(self,channel=None, vertical=None, horizontal_x=None, horizontal_y=None, force=False): ''' - Returns a new ImageData of a single slice of in the requested direction. + Returns a new ImageData of a single slice in the requested direction. ''' try: geometry_new = self.geometry.get_slice(channel=channel, vertical=vertical, horizontal_x=horizontal_x, horizontal_y=horizontal_y) @@ -156,6 +156,13 @@ def get_slice(self,channel=None, vertical=None, horizontal_x=None, horizontal_y= return out else: return ImageData(out.array, deep_copy=False, geometry=geometry_new) + + def get_centre_slice(self): + ''' + Returns a new ImageData of the centre slice in the vertical direction. + ''' + return self.get_slice(vertical='centre') + def apply_circular_mask(self, radius=0.99, in_place=True): diff --git a/Wrappers/Python/cil/framework/image_geometry.py b/Wrappers/Python/cil/framework/image_geometry.py index fc7f7771e..3e38a2314 100644 --- a/Wrappers/Python/cil/framework/image_geometry.py +++ b/Wrappers/Python/cil/framework/image_geometry.py @@ -174,27 +174,65 @@ def __init__(self, def get_slice(self,channel=None, vertical=None, horizontal_x=None, horizontal_y=None): ''' - Returns a new ImageGeometry of a single slice of in the requested direction. + Returns a new ImageGeometry of a single slice in the requested direction. + + Parameters + ---------- + channel : int or 'centre', optional + The channel index to slice. Default is None (no slicing). + vertical : int or 'centre', optional + The vertical index to slice. Default is None (no slicing). + horizontal_x : int, optional + The horizontal x index to slice. Default is None (no slicing). + horizontal_y : int, optional + The horizontal y index to slice. Default is None (no slicing). + Returns + ------- + geometry_new : ImageGeometry + A new ImageGeometry object representing the sliced geometry. + + Note + ---- + Slicing on vertical with 'centre' will return the central slice in that dimension. + Slicing on channels returns a geometry with a single channel, however the channel label is not + typically stored in the geometry. ''' + geometry_new = self.copy() if channel is not None: geometry_new.channels = 1 - try: geometry_new.channel_labels = [self.channel_labels[channel]] except: geometry_new.channel_labels = None if vertical is not None: - geometry_new.voxel_num_z = 0 + geometry_new.voxel_num_z = 1 + if vertical != 'centre': + if vertical == 0: + warnings.warn("Slicing vertical at index 0 results in a geometry \ + offset along the vertical axis. If you do not require an offset ImageGeometry, set vertical='centre", + UserWarning) + voxel_offset = (self.voxel_num_z)/2 - (vertical+0.5) + geometry_new.center_z -= voxel_offset * geometry_new.voxel_size_z if horizontal_y is not None: - geometry_new.voxel_num_y = 0 + geometry_new.voxel_num_y = 1 + voxel_offset = (self.voxel_num_y)/2 - (horizontal_y +0.5) + geometry_new.center_y -= voxel_offset * geometry_new.voxel_size_y if horizontal_x is not None: - geometry_new.voxel_num_x = 0 + geometry_new.voxel_num_x = 1 + voxel_offset = (self.voxel_num_x)/2 - (horizontal_x+0.5) + geometry_new.center_x -= voxel_offset * geometry_new.voxel_size_x return geometry_new + + def get_centre_slice(self): + ''' + Returns a new ImageGeometry of the centre slice in the vertical direction. + ''' + return self.get_slice(vertical='centre') def get_order_by_label(self, dimension_labels, default_dimension_labels): order = [] diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py index cb8bcf3d3..65cec9651 100644 --- a/Wrappers/Python/test/test_DataContainer.py +++ b/Wrappers/Python/test/test_DataContainer.py @@ -860,13 +860,20 @@ def test_ImageDataSubset(self): ss1 = vol.get_slice(horizontal_x = 0) self.assertListEqual(['channel', 'horizontal_y'], list(ss1.geometry.dimension_labels)) + ss2 = vol.get_slice(channel=0) + self.assertListEqual([ImageDimension["HORIZONTAL_Y"], ImageDimension["HORIZONTAL_X"]], list(ss2.geometry.dimension_labels)) + vg = ImageGeometry(3,4,5,channels=2) self.assertListEqual([ImageDimension["CHANNEL"], ImageDimension["VERTICAL"], ImageDimension["HORIZONTAL_Y"], ImageDimension["HORIZONTAL_X"]], list(vg.dimension_labels)) - ss2 = vg.allocate() - ss3 = vol.get_slice(channel=0) - self.assertListEqual([ImageDimension["HORIZONTAL_Y"], ImageDimension["HORIZONTAL_X"]], list(ss3.geometry.dimension_labels)) + vol2 = vg.allocate() + + ss3 = vol2.get_slice(vertical=0) + self.assertListEqual([ImageDimension["CHANNEL"], ImageDimension["HORIZONTAL_Y"], ImageDimension["HORIZONTAL_X"]], + list(ss3.geometry.dimension_labels)) + + def test_DataContainerSubset(self): dc = DataContainer(np.ones((2,3,4,5))) diff --git a/Wrappers/Python/test/test_ImageGeometry.py b/Wrappers/Python/test/test_ImageGeometry.py new file mode 100644 index 000000000..d31058379 --- /dev/null +++ b/Wrappers/Python/test/test_ImageGeometry.py @@ -0,0 +1,109 @@ +# Copyright 2025 United Kingdom Research and Innovation +# Copyright 2025 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 unittest +from utils import initialise_tests +from cil.framework import ImageGeometry +from cil.framework.labels import ImageDimension + +initialise_tests() + + +class TestImageGeometry(unittest.TestCase): + def setUp(self): + self.ig = ImageGeometry(2,3,4,channels=5) + + # get_slice tests --------------------------------------------------------------------------------------------------------------- + + def test_get_slice_vertical(self): + non_default_dimension_labels = [ImageDimension["HORIZONTAL_X"], ImageDimension["CHANNEL"], ImageDimension["HORIZONTAL_Y"], + ImageDimension["VERTICAL"]] + ig = self.ig.copy() + ig.set_labels(non_default_dimension_labels) + ig.voxel_size_z = 5.5 + sub = ig.get_slice(vertical = 1) + self.assertTrue( sub.shape == (2,5,3)) + self.assertEqual(sub.voxel_size_z,5.5) + self.assertEqual(sub.center_x,0) + self.assertEqual(sub.center_y,0) + self.assertEqual(sub.center_z,-0.5*sub.voxel_size_z) + + def test_get_slice_vertical_centre(self): + non_default_dimension_labels = [ImageDimension["HORIZONTAL_X"], ImageDimension["CHANNEL"], ImageDimension["HORIZONTAL_Y"], + ImageDimension["VERTICAL"]] + self.ig.set_labels(non_default_dimension_labels) + sub = self.ig.get_slice(vertical = 'centre') + self.assertTrue( sub.shape == (2,5,3)) + self.assertEqual(sub.center_x,0) + self.assertEqual(sub.center_y,0) + self.assertEqual(sub.center_z,0) + + + def test_get_slice_horizontal_x(self): + non_default_dimension_labels = [ImageDimension["HORIZONTAL_X"], ImageDimension["CHANNEL"], ImageDimension["HORIZONTAL_Y"], + ImageDimension["VERTICAL"]] + self.ig.set_labels(non_default_dimension_labels) + sub = self.ig.get_slice(horizontal_x = 1) + self.assertTrue(sub.shape == (5,3,4)) + self.assertEqual(sub.center_x,0.5) + self.assertEqual(sub.center_y,0) + self.assertEqual(sub.center_z,0) + + + def test_get_slice_channel(self): + non_default_dimension_labels = [ImageDimension["HORIZONTAL_X"], ImageDimension["CHANNEL"], ImageDimension["HORIZONTAL_Y"], + ImageDimension["VERTICAL"]] + self.ig.set_labels(non_default_dimension_labels) + sub = self.ig.get_slice(channel = 1) + self.assertTrue(sub.shape == (2,3,4)) + self.assertTrue(sub.channels == 1) + + def test_get_slice_channel_centre(self): + non_default_dimension_labels = [ImageDimension["HORIZONTAL_X"], ImageDimension["CHANNEL"], ImageDimension["HORIZONTAL_Y"], + ImageDimension["VERTICAL"]] + self.ig.set_labels(non_default_dimension_labels) + sub = self.ig.get_slice(channel = 'centre') + self.assertTrue(sub.shape == (2,3,4)) + self.assertTrue(sub.channels == 1) + + def test_get_slice_horizontal_y(self): + non_default_dimension_labels = [ImageDimension["HORIZONTAL_X"], ImageDimension["HORIZONTAL_Y"]] + self.ig.set_labels(non_default_dimension_labels) + sub = self.ig.get_slice(horizontal_y = 0) + self.assertTrue( sub.shape == (2,)) + self.assertEqual(sub.center_x,0) + self.assertEqual(sub.center_y,-1) + self.assertEqual(sub.center_z,0) + + def test_get_slice_horizontal_x_and_horizontal_y(self): + sub = self.ig.get_slice(horizontal_x=0,horizontal_y=0) + self.assertTrue( sub.shape == (5,4)) + + # test get_centre_slice --------------------------------------------------------------------------------------------------------------- + + def test_get_centre_slice(self): + sub = self.ig.get_centre_slice() + sub2 = self.ig.get_slice(vertical='centre') + self.assertTrue( sub == sub2) + + # ------------------------------------------------------------------------------------------------------------------------------- + + def test_shape_change_with_new_labels(self): + new_dimension_labels = [ImageDimension["HORIZONTAL_Y"], ImageDimension["CHANNEL"], ImageDimension["VERTICAL"], ImageDimension["HORIZONTAL_X"]] + self.ig.set_labels(new_dimension_labels) + self.assertTrue( self.ig.shape == (3,5,4,2)) \ No newline at end of file diff --git a/Wrappers/Python/test/test_reconstructors.py b/Wrappers/Python/test/test_reconstructors.py index 572e4309c..9f6d5bedf 100644 --- a/Wrappers/Python/test/test_reconstructors.py +++ b/Wrappers/Python/test/test_reconstructors.py @@ -80,7 +80,6 @@ def setUp(self): .set_angles(angles)\ .set_panel((det_pix_x,det_pix_y), (pix_size,pix_size))\ .set_labels(['angle','vertical','horizontal']) - self.ig3D = self.ag3D.get_ImageGeometry() self.ad3D = self.ag3D.allocate('random', seed=3) self.ig3D = self.ag3D.get_ImageGeometry() @@ -164,7 +163,6 @@ def setUp(self): .set_angles(angles)\ .set_panel((det_pix_x,det_pix_y), (pix_size,pix_size))\ .set_labels(['angle','vertical','horizontal']) - self.ig3D = self.ag3D.get_ImageGeometry() self.ad3D = self.ag3D.allocate('random', seed=4) self.ig3D = self.ag3D.get_ImageGeometry() @@ -351,7 +349,6 @@ def setUp(self): .set_angles(angles)\ .set_panel((det_pix_x,det_pix_y), (pix_size,pix_size))\ .set_labels(['angle','vertical','horizontal']) - self.ig3D = self.ag3D.get_ImageGeometry() self.ad3D = self.ag3D.allocate('random', seed=5) self.ig3D = self.ag3D.get_ImageGeometry() @@ -454,7 +451,6 @@ def setUp(self): .set_angles(angles)\ .set_panel((det_pix_x,det_pix_y), (pix_size,pix_size))\ .set_labels(['angle','vertical','horizontal']) - self.ig3D = self.ag3D.get_ImageGeometry() self.ad3D = self.ag3D.allocate('random', seed=3) self.ig3D = self.ag3D.get_ImageGeometry() diff --git a/Wrappers/Python/test/test_subset.py b/Wrappers/Python/test/test_subset.py index 638e88231..bfd73a200 100644 --- a/Wrappers/Python/test/test_subset.py +++ b/Wrappers/Python/test/test_subset.py @@ -161,7 +161,7 @@ def test_DataContainer(self): self.assertEqual(data_new.shape,(5,)) numpy.testing.assert_array_equal(data_new.array, arr[1,1,3,:]) - def test_ImageData(self): + def test_ImageData(self): ig = ImageGeometry(voxel_num_x=5, voxel_num_y=4, voxel_num_z=3, channels=2, dimension_labels=['channel','vertical','horizontal_y','horizontal_x']) data = ig.allocate(None) data_new = data.get_slice(vertical=1) @@ -299,6 +299,12 @@ def test_ImageDataSubset5a(self): sub = data.get_slice(horizontal_y = 1) self.assertTrue( sub.shape == (2,)) + def test_ImageDataCentreSubset(self): + data = self.ig.allocate() + sub = data.get_centre_slice() + sub2 = data.get_slice(vertical='centre') + self.assertTrue( sub == sub2) + def test_ImageDataSubset1b(self): non_default_dimension_labels = [ImageDimension["HORIZONTAL_X"], ImageDimension["CHANNEL"], ImageDimension["HORIZONTAL_Y"], ImageDimension["VERTICAL"]] @@ -313,7 +319,6 @@ def test_ImageDataSubset1c(self): sub = data.get_slice(channel=0,horizontal_x=0,horizontal_y=0) self.assertTrue( sub.shape == (4,)) - def test_AcquisitionDataAllocate1a(self): data = self.ag.allocate() default_dimension_labels = [AcquisitionDimension["CHANNEL"] , @@ -429,3 +434,9 @@ def test_AcquisitionDataSubset1l(self): sub = data.get_slice(projection = 0) sub = sub.get_slice(channel = 0) self.assertTrue(sub.shape == (2, 20)) + + def test_AcquisitionDataCentreSubset(self): + data = self.ag.allocate() + sub = data.get_centre_slice() + sub2 = data.get_slice(vertical='centre') + self.assertTrue( sub == sub2) diff --git a/Wrappers/Python/test/utils_projectors.py b/Wrappers/Python/test/utils_projectors.py index 9d767f51e..91a5c1393 100644 --- a/Wrappers/Python/test/utils_projectors.py +++ b/Wrappers/Python/test/utils_projectors.py @@ -656,10 +656,17 @@ def test_FBP_roi(self): np.testing.assert_allclose(reco.as_array(), self.gold_roi, atol=self.tolerance_fbp_roi) if AcquisitionType.DIM3 & self.ag.dimension: + # Manually generated single slice Image Geometry FBP = self.FBP(self.ig_single_slice, self.ag, **self.FBP_args) reco = FBP(self.acq_data) np.testing.assert_allclose(reco.as_array(), self.gold_roi_single_slice, atol=self.tolerance_fbp_roi) + # CIL generated single slice Image Geometry + central_slice_ig = self.ig.get_slice(vertical='centre') + FBP = self.FBP(central_slice_ig, self.ag, **self.FBP_args) + reco = FBP(self.acq_data) + np.testing.assert_allclose(reco.as_array(), self.img_data.get_slice(vertical='centre').as_array(),atol=self.tolerance_fbp) + def test_input_arguments(self): # default image_geometry, named parameter acquisition_geometry