Image Geometry get_slice returns IG with offset centre

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:

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):
         '''

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):

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 = []

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)))

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))

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()

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)

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