From 33584af20c1b52b78e4cd92abfbd4173435be423 Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Fri, 10 Oct 2025 15:11:06 +0000
Subject: [PATCH 01/11] First draft shrink volume code

---
 .../Python/cil/utilities/shrink_volume.py     | 210 ++++++++++++++++++
 1 file changed, 210 insertions(+)
 create mode 100644 Wrappers/Python/cil/utilities/shrink_volume.py

diff --git a/Wrappers/Python/cil/utilities/shrink_volume.py b/Wrappers/Python/cil/utilities/shrink_volume.py
new file mode 100644
index 0000000000..d3c0074fbf
--- /dev/null
+++ b/Wrappers/Python/cil/utilities/shrink_volume.py
@@ -0,0 +1,210 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.ndimage import label
+from skimage.filters import threshold_otsu
+
+from cil.processors import Binner, Slicer
+from cil.plugins.astra.processors import FBP
+
+import logging
+log = logging.getLogger(__name__)
+class VolumeShrinker(object):
+    """
+    Shrinks the reconstruction volume based on a supplied volume size or 
+    automatic detection of the region of interest using Otsu thresholding and 
+    connected components.   
+    """
+
+    def run(self, data, auto=True, threshold='Otsu', buffer=None, manual_limits=None):
+        """
+        Parameters
+        ----------
+        auto : bool, optional
+            If True, automatically detect and crop the reconstruction volume
+            If False, use manual_limits
+
+        threshold: string or float, optional
+            If automatically detecting the limits, specify the intensity threshold
+            between sample and background. By default use an Otsu filter. 
+
+        buffer: float, optional
+            Add a buffer around the automatically detected limits, expressed as 
+            a percentage of the axis size.
+
+        manual_limits : dict, optional
+            The limits {'axis_name1':(min, max), 'axis_name2':(min, max)}
+            The `key` being the axis name to apply the processor to, 
+            the `value` holding a tuple containing the min and max limits
+            or None, to specify no limit
+            Manual limits over-ride automatically detected limits
+        """
+
+        binning = min(int(np.ceil(data.geometry.config.panel.num_pixels[0] / 128)),16)
+        angle_binning = np.ceil(data.get_dimension_size('angle')/(data.get_dimension_size('horizontal')*(np.pi/2)))
+        roi = {
+                'horizontal': (None, None, binning),
+                'vertical': (None, None, binning),
+                'angle' : (None, None, angle_binning)
+            }
+        data_binned = Binner(roi)(data)
+
+        ag = data_binned.geometry
+        ig = ag.get_ImageGeometry()
+
+        fbp = FBP(ig, ag)
+        recon = fbp(data_binned)
+        recon.apply_circular_mask(0.9)
+
+        if auto:
+            bounds = self.reduce_reconstruction_volume(recon, binning, threshold, buffer)
+        else:
+            bounds = {}
+            for dim in recon.dimension_labels:
+                bounds[dim] = (0, recon.get_dimension_size(dim)*binning)
+
+        if manual_limits is not None:    
+            for dim, v in manual_limits.items():
+                if dim in recon.dimension_labels:
+                    if v is None:
+                        v = (0, recon.get_dimension_size(dim)*binning)
+                    elif v[0] is None:
+                        v[0] = 0
+                    elif v[1] is None:
+                        v[1] = recon.get_dimension_size(dim)*binning
+                    bounds[dim] = v
+                else:
+                    raise ValueError("dimension {} not recognised, must be one of {}".format(dim, recon.dimension_labels))
+
+        self.plot_with_bounds(recon, bounds, binning)
+
+        return self.update_ig(data.geometry.get_ImageGeometry(), bounds)
+
+    def update_ig(self, ig_unbinned, bounds):
+        ig = Slicer(roi={'horizontal_x':(bounds['horizontal_x'][0], bounds['horizontal_x'][1],1),
+                  'horizontal_y':(bounds['horizontal_y'][0], bounds['horizontal_y'][1], 1),
+                  'vertical':(bounds['vertical'][0], bounds['vertical'][1], 1)})(ig_unbinned)
+        return ig
+
+    def plot_with_bounds(self, recon, bounds, binning):
+        fig, axs = plt.subplots(nrows=1, ncols=recon.ndim, figsize=(14, 6))
+
+        dims = recon.dimension_labels
+        for i, dim in enumerate(dims):
+            ax = axs[i]
+
+            other_dims = [d for d in dims if d != dim]
+            y_dim, x_dim = other_dims
+            x_size = recon.get_dimension_size(x_dim)*binning
+            y_size = recon.get_dimension_size(y_dim)*binning
+
+            ax.imshow(recon.max(axis=dim).array, origin='lower', cmap='gray',
+                    extent=[0, x_size, 0, y_size])
+            
+            x_min, x_max = bounds[x_dim]
+            y_min, y_max = bounds[y_dim]
+
+            ax.plot([x_min, x_max], [y_min, y_min], '--r')
+            ax.plot([x_min, x_max], [y_max, y_max], '--r')
+            ax.plot([x_min, x_min], [y_min, y_max], '--r')
+            ax.plot([x_max, x_max], [y_min, y_max], '--r')
+
+            ax.set_xlabel(x_dim)
+            ax.set_ylabel(y_dim)
+            ax.set_title(f"Maximum values in direction: {dim}")
+
+    def reduce_reconstruction_volume(self, recon, binning, threshold, buffer):
+            
+        dims = recon.dimension_labels
+        all_bounds = {dim: [] for dim in dims}
+
+        for dim in dims:
+            arr = recon.max(axis=dim).array
+            mask, large_components_mask = self.otsu_large_components(arr, threshold)
+
+            x_indices = np.where(np.any(large_components_mask, axis=0))[0]
+            y_indices = np.where(np.any(large_components_mask, axis=1))[0]
+            x_min, x_max = x_indices[0], x_indices[-1]
+            y_min, y_max = y_indices[0], y_indices[-1]
+
+            axis = recon.get_dimension_axis(dim)
+            other_axes = [j for j in range(recon.ndim) if j != axis]
+
+            if buffer is not None:
+                y_full = recon.get_dimension_size(dims[other_axes[0]])
+                y_min_buffer = np.max([0, (y_min-y_full//buffer)])
+                y_max_buffer = np.min([y_full, y_max+(y_full//buffer)])
+
+                x_full = recon.get_dimension_size(dims[other_axes[1]])
+                x_min_buffer = np.max([0, (x_min-x_full//buffer)])
+                x_max_buffer = np.min([x_full, x_max+(x_full//buffer)])
+
+                all_bounds[dims[other_axes[0]]].append((y_min_buffer, y_max_buffer))
+                all_bounds[dims[other_axes[1]]].append((x_min_buffer, x_max_buffer))
+            else:
+                all_bounds[dims[other_axes[0]]].append((y_min, y_max))
+                all_bounds[dims[other_axes[1]]].append((x_min, x_max))
+
+        bounds = {}
+        for dim in dims:
+
+            mins = [b[0] for b in all_bounds[dim]]
+            maxs = [b[1] for b in all_bounds[dim]]
+            dim_min = np.min(mins)*binning
+            dim_max = np.max(maxs)*binning
+            
+            bounds[dim] = (dim_min, dim_max)
+
+            if log.isEnabledFor(logging.DEBUG):
+                print(f"{dim}: {bounds[dim][0]} to {bounds[dim][1]}")
+            
+        return bounds
+        
+    def otsu_large_components(self, arr, threshold):
+        
+
+        if isinstance(threshold, (int, float)):
+            thresh = threshold
+        elif isinstance(threshold, str) and threshold.lower() == 'otsu':
+            thresh = threshold_otsu(arr[arr > 0])
+        else:
+            raise ValueError(f"Threshold {threshold} not recognised, must be a number or 'Otsu'")
+        mask = arr > thresh
+
+
+        labeled_mask, num_features = label(mask)
+        component_sizes = np.bincount(labeled_mask.ravel())
+        min_size = 10
+
+        large_labels = np.where(component_sizes > min_size)[0]
+        large_labels = large_labels[large_labels != 0]  
+        large_components_mask = np.isin(labeled_mask, large_labels)
+
+        if log.isEnabledFor(logging.DEBUG):
+            fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(8, 2.5))
+
+            axes[0].imshow(arr, cmap=plt.cm.gray)
+            axes[0].set_title('Original')
+
+            axes[1].hist(arr.ravel(), bins=100)
+            axes[1].set_title('Histogram')
+            axes[1].axvline(thresh, color='r')
+
+            axes[2].imshow(mask, cmap=plt.cm.gray, extent=[axes[0].get_xlim()[0], axes[0].get_xlim()[1], axes[0].get_ylim()[0], axes[0].get_ylim()[1]])
+            axes[2].set_title('Thresholded')
+
+            axes[3].imshow(large_components_mask, cmap=plt.cm.gray, extent=[axes[0].get_xlim()[0], axes[0].get_xlim()[1], axes[0].get_ylim()[0], axes[0].get_ylim()[1]])
+            axes[3].set_title('Large components')
+
+            x_indices = np.where(np.any(large_components_mask, axis=0))[0]
+            y_indices = np.where(np.any(large_components_mask, axis=1))[0]
+            x_min, x_max = x_indices[0], x_indices[-1]
+            y_min, y_max = y_indices[0], y_indices[-1]
+
+            axes[3].plot([x_min, x_max], [y_min, y_min], '--r')
+            axes[3].plot([x_min, x_max], [y_max, y_max], '--r')
+            axes[3].plot([x_min, x_min], [y_min, y_max], '--r')
+            axes[3].plot([x_max, x_max], [y_min, y_max], '--r')
+
+            plt.tight_layout()
+
+        return mask, large_components_mask
\ No newline at end of file

From 10f42b7d4e92415643a9d547597785de8e9f14ce Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Mon, 27 Oct 2025 14:03:21 +0000
Subject: [PATCH 02/11] Make mask optional

---
 Wrappers/Python/cil/utilities/shrink_volume.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/Wrappers/Python/cil/utilities/shrink_volume.py b/Wrappers/Python/cil/utilities/shrink_volume.py
index d3c0074fbf..49e920317c 100644
--- a/Wrappers/Python/cil/utilities/shrink_volume.py
+++ b/Wrappers/Python/cil/utilities/shrink_volume.py
@@ -15,7 +15,7 @@ class VolumeShrinker(object):
     connected components.   
     """
 
-    def run(self, data, auto=True, threshold='Otsu', buffer=None, manual_limits=None):
+    def run(self, data, auto=True, threshold='Otsu', buffer=None, mask_radius=None, manual_limits=None):
         """
         Parameters
         ----------
@@ -31,6 +31,10 @@ def run(self, data, auto=True, threshold='Otsu', buffer=None, manual_limits=None
             Add a buffer around the automatically detected limits, expressed as 
             a percentage of the axis size.
 
+        mask_radius: float, optional
+            Radius of circular mask to apply on the reconstructed volume, before
+            automatically cropping the recontruction volume. Default is None.
+
         manual_limits : dict, optional
             The limits {'axis_name1':(min, max), 'axis_name2':(min, max)}
             The `key` being the axis name to apply the processor to, 
@@ -53,8 +57,9 @@ def run(self, data, auto=True, threshold='Otsu', buffer=None, manual_limits=None
 
         fbp = FBP(ig, ag)
         recon = fbp(data_binned)
-        recon.apply_circular_mask(0.9)
-
+        if mask_radius is not None:
+            recon.apply_circular_mask(mask_radius)
+            
         if auto:
             bounds = self.reduce_reconstruction_volume(recon, binning, threshold, buffer)
         else:

From 6dcec8b6ea7a6033f972232568044405fffebcf2 Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Tue, 9 Dec 2025 09:54:39 +0000
Subject: [PATCH 03/11] Update logic if auto=False

---
 .../Python/cil/utilities/shrink_volume.py     | 72 ++++++++++++-------
 1 file changed, 46 insertions(+), 26 deletions(-)

diff --git a/Wrappers/Python/cil/utilities/shrink_volume.py b/Wrappers/Python/cil/utilities/shrink_volume.py
index 49e920317c..f1b22ca03a 100644
--- a/Wrappers/Python/cil/utilities/shrink_volume.py
+++ b/Wrappers/Python/cil/utilities/shrink_volume.py
@@ -1,7 +1,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.ndimage import label
-from skimage.filters import threshold_otsu
+from skimage.filters import threshold_otsu, threshold_multiotsu
 
 from cil.processors import Binner, Slicer
 from cil.plugins.astra.processors import FBP
@@ -67,27 +67,32 @@ def run(self, data, auto=True, threshold='Otsu', buffer=None, mask_radius=None,
             for dim in recon.dimension_labels:
                 bounds[dim] = (0, recon.get_dimension_size(dim)*binning)
 
-        if manual_limits is not None:    
-            for dim, v in manual_limits.items():
-                if dim in recon.dimension_labels:
-                    if v is None:
-                        v = (0, recon.get_dimension_size(dim)*binning)
-                    elif v[0] is None:
-                        v[0] = 0
-                    elif v[1] is None:
-                        v[1] = recon.get_dimension_size(dim)*binning
-                    bounds[dim] = v
-                else:
-                    raise ValueError("dimension {} not recognised, must be one of {}".format(dim, recon.dimension_labels))
+            if manual_limits is not None:    
+                for dim, v in manual_limits.items():
+                    if dim in recon.dimension_labels:
+                        if v is None:
+                            v = (0, recon.get_dimension_size(dim)*binning)
+                        elif v[0] is None:
+                            v[0] = 0
+                        elif v[1] is None:
+                            v[1] = recon.get_dimension_size(dim)*binning
+                        bounds[dim] = v
+                    else:
+                        raise ValueError("dimension {} not recognised, must be one of {}".format(dim, recon.dimension_labels))
+            else:
+                bounds = None
 
         self.plot_with_bounds(recon, bounds, binning)
-
+        
         return self.update_ig(data.geometry.get_ImageGeometry(), bounds)
 
     def update_ig(self, ig_unbinned, bounds):
-        ig = Slicer(roi={'horizontal_x':(bounds['horizontal_x'][0], bounds['horizontal_x'][1],1),
-                  'horizontal_y':(bounds['horizontal_y'][0], bounds['horizontal_y'][1], 1),
-                  'vertical':(bounds['vertical'][0], bounds['vertical'][1], 1)})(ig_unbinned)
+        if bounds is None:
+            ig = ig_unbinned
+        else:
+            ig = Slicer(roi={'horizontal_x':(bounds['horizontal_x'][0], bounds['horizontal_x'][1],1),
+                    'horizontal_y':(bounds['horizontal_y'][0], bounds['horizontal_y'][1], 1),
+                    'vertical':(bounds['vertical'][0], bounds['vertical'][1], 1)})(ig_unbinned)
         return ig
 
     def plot_with_bounds(self, recon, bounds, binning):
@@ -105,13 +110,14 @@ def plot_with_bounds(self, recon, bounds, binning):
             ax.imshow(recon.max(axis=dim).array, origin='lower', cmap='gray',
                     extent=[0, x_size, 0, y_size])
             
-            x_min, x_max = bounds[x_dim]
-            y_min, y_max = bounds[y_dim]
+            if bounds is not None:
+                x_min, x_max = bounds[x_dim]
+                y_min, y_max = bounds[y_dim]
 
-            ax.plot([x_min, x_max], [y_min, y_min], '--r')
-            ax.plot([x_min, x_max], [y_max, y_max], '--r')
-            ax.plot([x_min, x_min], [y_min, y_max], '--r')
-            ax.plot([x_max, x_max], [y_min, y_max], '--r')
+                ax.plot([x_min, x_max], [y_min, y_min], '--r')
+                ax.plot([x_min, x_max], [y_max, y_max], '--r')
+                ax.plot([x_min, x_min], [y_min, y_max], '--r')
+                ax.plot([x_max, x_max], [y_min, y_max], '--r')
 
             ax.set_xlabel(x_dim)
             ax.set_ylabel(y_dim)
@@ -126,8 +132,8 @@ def reduce_reconstruction_volume(self, recon, binning, threshold, buffer):
             arr = recon.max(axis=dim).array
             mask, large_components_mask = self.otsu_large_components(arr, threshold)
 
-            x_indices = np.where(np.any(large_components_mask, axis=0))[0]
-            y_indices = np.where(np.any(large_components_mask, axis=1))[0]
+            x_indices = np.where(np.any(mask, axis=0))[0]
+            y_indices = np.where(np.any(mask, axis=1))[0]
             x_min, x_max = x_indices[0], x_indices[-1]
             y_min, y_max = y_indices[0], y_indices[-1]
 
@@ -170,7 +176,11 @@ def otsu_large_components(self, arr, threshold):
         if isinstance(threshold, (int, float)):
             thresh = threshold
         elif isinstance(threshold, str) and threshold.lower() == 'otsu':
-            thresh = threshold_otsu(arr[arr > 0])
+            # thresh = threshold_otsu(arr[arr > 0])
+            classes = 4
+            n_bins = 256
+            thresh = threshold_multiotsu(arr[arr>0], classes=classes, nbins=n_bins)
+            thresh = thresh[0]
         else:
             raise ValueError(f"Threshold {threshold} not recognised, must be a number or 'Otsu'")
         mask = arr > thresh
@@ -197,6 +207,16 @@ def otsu_large_components(self, arr, threshold):
             axes[2].imshow(mask, cmap=plt.cm.gray, extent=[axes[0].get_xlim()[0], axes[0].get_xlim()[1], axes[0].get_ylim()[0], axes[0].get_ylim()[1]])
             axes[2].set_title('Thresholded')
 
+            x_indices = np.where(np.any(mask, axis=0))[0]
+            y_indices = np.where(np.any(mask, axis=1))[0]
+            x_min, x_max = x_indices[0], x_indices[-1]
+            y_min, y_max = y_indices[0], y_indices[-1]
+
+            axes[2].plot([x_min, x_max], [y_min, y_min], '--r')
+            axes[2].plot([x_min, x_max], [y_max, y_max], '--r')
+            axes[2].plot([x_min, x_min], [y_min, y_max], '--r')
+            axes[2].plot([x_max, x_max], [y_min, y_max], '--r')
+
             axes[3].imshow(large_components_mask, cmap=plt.cm.gray, extent=[axes[0].get_xlim()[0], axes[0].get_xlim()[1], axes[0].get_ylim()[0], axes[0].get_ylim()[1]])
             axes[3].set_title('Large components')
 

From 69822bc76a2d833028f312381d68e4525d8f9e1d Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Tue, 10 Feb 2026 13:43:04 +0000
Subject: [PATCH 04/11] Rename functions

---
 Wrappers/Python/cil/utilities/shrink_volume.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/Wrappers/Python/cil/utilities/shrink_volume.py b/Wrappers/Python/cil/utilities/shrink_volume.py
index f1b22ca03a..e9fb1b4abb 100644
--- a/Wrappers/Python/cil/utilities/shrink_volume.py
+++ b/Wrappers/Python/cil/utilities/shrink_volume.py
@@ -33,7 +33,7 @@ def run(self, data, auto=True, threshold='Otsu', buffer=None, mask_radius=None,
 
         mask_radius: float, optional
             Radius of circular mask to apply on the reconstructed volume, before
-            automatically cropping the recontruction volume. Default is None.
+            automatically cropping the reconstruction volume. Default is None.
 
         manual_limits : dict, optional
             The limits {'axis_name1':(min, max), 'axis_name2':(min, max)}
@@ -130,7 +130,8 @@ def reduce_reconstruction_volume(self, recon, binning, threshold, buffer):
 
         for dim in dims:
             arr = recon.max(axis=dim).array
-            mask, large_components_mask = self.otsu_large_components(arr, threshold)
+            
+            mask, large_components_mask = self.threshold_large_components(arr, threshold)
 
             x_indices = np.where(np.any(mask, axis=0))[0]
             y_indices = np.where(np.any(mask, axis=1))[0]
@@ -170,7 +171,7 @@ def reduce_reconstruction_volume(self, recon, binning, threshold, buffer):
             
         return bounds
         
-    def otsu_large_components(self, arr, threshold):
+    def threshold_large_components(self, arr, threshold):
         
 
         if isinstance(threshold, (int, float)):
@@ -183,8 +184,8 @@ def otsu_large_components(self, arr, threshold):
             thresh = thresh[0]
         else:
             raise ValueError(f"Threshold {threshold} not recognised, must be a number or 'Otsu'")
-        mask = arr > thresh
 
+        mask = arr > thresh
 
         labeled_mask, num_features = label(mask)
         component_sizes = np.bincount(labeled_mask.ravel())

From e15f3db3dc17f8335a285b4c6549d7855cf320c0 Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Fri, 20 Feb 2026 17:34:37 +0000
Subject: [PATCH 05/11] Add tests

---
 .../Python/cil/utilities/shrink_volume.py     | 289 +++++++++++-------
 Wrappers/Python/test/test_shrink_volume.py    | 126 ++++++++
 2 files changed, 298 insertions(+), 117 deletions(-)
 create mode 100644 Wrappers/Python/test/test_shrink_volume.py

diff --git a/Wrappers/Python/cil/utilities/shrink_volume.py b/Wrappers/Python/cil/utilities/shrink_volume.py
index e9fb1b4abb..96f46aed16 100644
--- a/Wrappers/Python/cil/utilities/shrink_volume.py
+++ b/Wrappers/Python/cil/utilities/shrink_volume.py
@@ -1,90 +1,133 @@
+#  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
 import matplotlib.pyplot as plt
 from scipy.ndimage import label
-from skimage.filters import threshold_otsu, threshold_multiotsu
+from skimage.filters import threshold_multiotsu
+import importlib
+from cil.framework.labels import AcquisitionType
+
 
 from cil.processors import Binner, Slicer
-from cil.plugins.astra.processors import FBP
+
 
 import logging
 log = logging.getLogger(__name__)
 class VolumeShrinker(object):
     """
-    Shrinks the reconstruction volume based on a supplied volume size or 
-    automatic detection of the region of interest using Otsu thresholding and 
-    connected components.   
+    Shrinks the reconstruction volume from a dataset based on supplied volume 
+    limits, thresholding or automatic detection of the region of interest using 
+    Otsu thresholding.
+    
+    Parameters
+    ----------
+    data: AcquisitionData
+        The dataset to create a reduced reconstruction volume from.
+
+    recon_backend : {'tigre', 'astra'}
+        The backend to use for the reconstruction
+
     """
 
-    def run(self, data, auto=True, threshold='Otsu', buffer=None, mask_radius=None, manual_limits=None):
+    _supported_backends = ['astra', 'tigre']
+
+    def __init__(self, data, recon_backend='tigre'):
+         
+        self.data = data
+        self.recon_backend = recon_backend
+        self.recon_method = self._configure_recon(recon_backend)
+
+    def run(self, limits=None, preview=True, method='manual', **kwargs):
         """
         Parameters
         ----------
-        auto : bool, optional
-            If True, automatically detect and crop the reconstruction volume
-            If False, use manual_limits
-
-        threshold: string or float, optional
-            If automatically detecting the limits, specify the intensity threshold
-            between sample and background. By default use an Otsu filter. 
+        limits : dict, optional
+            ImageGeometry limits {'axis_name1':(min, max), 'axis_name2':(min, max)}
+            The `key` being the ImageGeometry axis name and `value` a tuple containing 
+            the min and max limits.
+            Default None, uses the full extent of the data
+
+        preview: bool, optional
+            If True, plots the maximum values from a binned reconstruction in each 
+            direction to check the ImageGeometry contains the full volume. 
+            Default is False.
+        
+        method : string, optional
+            If 'manual', use manual_limits
+            If 'threshold' crop the reconstruction volume based on a threshold
+            value between sample and background. Pass threshold as a kwarg.
+            If 'otsu', automatically detect and crop the reconstruction volume.
+            Default is 'manual'
+
+        **kwargs:
+        threshold: float, optional
+            If using 'threshold' method, specify the intensity threshold
+            between sample and background. Default is None.
 
         buffer: float, optional
-            Add a buffer around the automatically detected limits, expressed as 
-            a percentage of the axis size.
+            Buffer in pixels around the automatically detected limits. 
+            Default is None, no buffer added.
 
         mask_radius: float, optional
             Radius of circular mask to apply on the reconstructed volume, before
             automatically cropping the reconstruction volume. Default is None.
 
-        manual_limits : dict, optional
-            The limits {'axis_name1':(min, max), 'axis_name2':(min, max)}
-            The `key` being the axis name to apply the processor to, 
-            the `value` holding a tuple containing the min and max limits
-            or None, to specify no limit
-            Manual limits over-ride automatically detected limits
-        """
-
-        binning = min(int(np.ceil(data.geometry.config.panel.num_pixels[0] / 128)),16)
-        angle_binning = np.ceil(data.get_dimension_size('angle')/(data.get_dimension_size('horizontal')*(np.pi/2)))
-        roi = {
-                'horizontal': (None, None, binning),
-                'vertical': (None, None, binning),
-                'angle' : (None, None, angle_binning)
-            }
-        data_binned = Binner(roi)(data)
+        otsu_classes: int, optional
+            Number of material classes to use when automatically detecting the 
+            reconstruction volume using Otsu thresholding method. Default is 2.
 
-        ag = data_binned.geometry
-        ig = ag.get_ImageGeometry()
+        min_component_size: int, optional
+            Minimum size in pixels of connected components to keep when automatically 
+            cropping the reconstruction volume. Default is None.
+        """
 
-        fbp = FBP(ig, ag)
-        recon = fbp(data_binned)
-        if mask_radius is not None:
-            recon.apply_circular_mask(mask_radius)
-            
-        if auto:
-            bounds = self.reduce_reconstruction_volume(recon, binning, threshold, buffer)
-        else:
+        ig = self.data.geometry.get_ImageGeometry()
+        if method.lower() == 'manual':
             bounds = {}
-            for dim in recon.dimension_labels:
-                bounds[dim] = (0, recon.get_dimension_size(dim)*binning)
-
-            if manual_limits is not None:    
-                for dim, v in manual_limits.items():
-                    if dim in recon.dimension_labels:
+            for dim in ig.dimension_labels:
+                bounds[dim] = (0, ig.shape[ig.dimension_labels.index(dim)])
+            if limits is not None:
+                for dim, v in limits.items():
+                    if dim in ig.dimension_labels:
                         if v is None:
-                            v = (0, recon.get_dimension_size(dim)*binning)
+                            v = (0, ig.shape[ig.dimension_labels.index(dim)])
                         elif v[0] is None:
                             v[0] = 0
                         elif v[1] is None:
-                            v[1] = recon.get_dimension_size(dim)*binning
+                            v[1] = ig.shape[ig.dimension_labels.index(dim)]
                         bounds[dim] = v
                     else:
-                        raise ValueError("dimension {} not recognised, must be one of {}".format(dim, recon.dimension_labels))
-            else:
-                bounds = None
+                        raise ValueError("dimension {} not recognised, must be one of {}".format(dim, ig.dimension_labels))
+
+        elif method.lower() in ['threshold', 'otsu']:
+            mask_radius = kwargs.pop('mask_radius', None)
+            recon, binning = self.get_recon(mask_radius=mask_radius)
+            bounds = self.reduce_reconstruction_volume(recon, binning, method, kwargs)
+
+        if preview:
+            if method.lower() == 'manual':
+                mask_radius = kwargs.pop('mask_radius', None)
+                recon, binning = self.get_recon(mask_radius=mask_radius)
 
-        self.plot_with_bounds(recon, bounds, binning)
+            self.plot_with_bounds(recon, bounds, binning)
         
-        return self.update_ig(data.geometry.get_ImageGeometry(), bounds)
+        return self.update_ig(self.data.geometry.get_ImageGeometry(), bounds)
 
     def update_ig(self, ig_unbinned, bounds):
         if bounds is None:
@@ -94,6 +137,27 @@ def update_ig(self, ig_unbinned, bounds):
                     'horizontal_y':(bounds['horizontal_y'][0], bounds['horizontal_y'][1], 1),
                     'vertical':(bounds['vertical'][0], bounds['vertical'][1], 1)})(ig_unbinned)
         return ig
+    
+    def get_recon(self, mask_radius):
+        binning = min(int(np.ceil(self.data.geometry.config.panel.num_pixels[0] / 128)),16)
+        angle_binning = np.ceil(self.data.get_dimension_size('angle')/(self.data.get_dimension_size('horizontal')*(np.pi/2)))
+        roi = {
+                'horizontal': (None, None, binning),
+                'vertical': (None, None, binning),
+                'angle' : (None, None, angle_binning)
+            }
+        data_binned = Binner(roi)(self.data)
+
+        ag = data_binned.geometry
+        ig = ag.get_ImageGeometry()
+
+        fbp = self.recon_method(ig, ag)
+        recon = fbp(data_binned)
+        
+        if mask_radius is not None:
+            recon.apply_circular_mask(mask_radius)
+
+        return recon, binning
 
     def plot_with_bounds(self, recon, bounds, binning):
         fig, axs = plt.subplots(nrows=1, ncols=recon.ndim, figsize=(14, 6))
@@ -123,32 +187,65 @@ def plot_with_bounds(self, recon, bounds, binning):
             ax.set_ylabel(y_dim)
             ax.set_title(f"Maximum values in direction: {dim}")
 
-    def reduce_reconstruction_volume(self, recon, binning, threshold, buffer):
-            
+    def reduce_reconstruction_volume(self, recon, binning, method, kwargs):
+        threshold = kwargs.pop('threshold', None)
+        buffer = kwargs.pop('buffer', None)
+        min_component_size = kwargs.pop('min_component_size', None)
+        otsu_classes = kwargs.pop('otsu_classes', 2)   
+
         dims = recon.dimension_labels
         all_bounds = {dim: [] for dim in dims}
 
         for dim in dims:
             arr = recon.max(axis=dim).array
             
-            mask, large_components_mask = self.threshold_large_components(arr, threshold)
+            if method.lower() == 'threshold':
+                mask = arr > threshold
+
+            elif method.lower() == 'otsu':
+                n_bins = 256
+                threshold = threshold_multiotsu(arr[arr>0], classes=otsu_classes, nbins=n_bins)
+                threshold = threshold[0]
+                mask = arr > threshold
+
+            if min_component_size is not None:
+                mask = self.threshold_large_components(mask, min_component_size)
 
             x_indices = np.where(np.any(mask, axis=0))[0]
             y_indices = np.where(np.any(mask, axis=1))[0]
             x_min, x_max = x_indices[0], x_indices[-1]
             y_min, y_max = y_indices[0], y_indices[-1]
 
+            if log.isEnabledFor(logging.DEBUG):
+                fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(8, 2.5))
+
+                axes[0].imshow(arr, cmap=plt.cm.gray)
+                axes[0].set_title('Original')
+
+                axes[1].hist(arr.ravel(), bins=100)
+                axes[1].set_title('Histogram')
+                axes[1].axvline(threshold, color='r')
+
+                axes[2].imshow(mask, cmap=plt.cm.gray, extent=[axes[0].get_xlim()[0], axes[0].get_xlim()[1], axes[0].get_ylim()[0], axes[0].get_ylim()[1]])
+                axes[2].set_title('Thresholded')
+
+                axes[2].plot([x_min, x_max], [y_min, y_min], '--r')
+                axes[2].plot([x_min, x_max], [y_max, y_max], '--r')
+                axes[2].plot([x_min, x_min], [y_min, y_max], '--r')
+                axes[2].plot([x_max, x_max], [y_min, y_max], '--r')
+                plt.tight_layout()
+
             axis = recon.get_dimension_axis(dim)
             other_axes = [j for j in range(recon.ndim) if j != axis]
 
             if buffer is not None:
                 y_full = recon.get_dimension_size(dims[other_axes[0]])
-                y_min_buffer = np.max([0, (y_min-y_full//buffer)])
-                y_max_buffer = np.min([y_full, y_max+(y_full//buffer)])
+                y_min_buffer = np.max([0, y_min-(buffer//binning)])
+                y_max_buffer = np.min([y_full, y_max+(buffer//binning)])
 
                 x_full = recon.get_dimension_size(dims[other_axes[1]])
-                x_min_buffer = np.max([0, (x_min-x_full//buffer)])
-                x_max_buffer = np.min([x_full, x_max+(x_full//buffer)])
+                x_min_buffer = np.max([0, x_min-(buffer//binning)])
+                x_max_buffer = np.min([x_full, x_max+(buffer//binning)])
 
                 all_bounds[dims[other_axes[0]]].append((y_min_buffer, y_max_buffer))
                 all_bounds[dims[other_axes[1]]].append((x_min_buffer, x_max_buffer))
@@ -171,66 +268,24 @@ def reduce_reconstruction_volume(self, recon, binning, threshold, buffer):
             
         return bounds
         
-    def threshold_large_components(self, arr, threshold):
-        
-
-        if isinstance(threshold, (int, float)):
-            thresh = threshold
-        elif isinstance(threshold, str) and threshold.lower() == 'otsu':
-            # thresh = threshold_otsu(arr[arr > 0])
-            classes = 4
-            n_bins = 256
-            thresh = threshold_multiotsu(arr[arr>0], classes=classes, nbins=n_bins)
-            thresh = thresh[0]
-        else:
-            raise ValueError(f"Threshold {threshold} not recognised, must be a number or 'Otsu'")
-
-        mask = arr > thresh
-
-        labeled_mask, num_features = label(mask)
+    def threshold_large_components(self, mask, min_component_size):
+        labeled_mask, _ = label(mask)
         component_sizes = np.bincount(labeled_mask.ravel())
-        min_size = 10
+        min_component_size = 10
 
-        large_labels = np.where(component_sizes > min_size)[0]
+        large_labels = np.where(component_sizes > min_component_size)[0]
         large_labels = large_labels[large_labels != 0]  
         large_components_mask = np.isin(labeled_mask, large_labels)
 
-        if log.isEnabledFor(logging.DEBUG):
-            fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(8, 2.5))
-
-            axes[0].imshow(arr, cmap=plt.cm.gray)
-            axes[0].set_title('Original')
-
-            axes[1].hist(arr.ravel(), bins=100)
-            axes[1].set_title('Histogram')
-            axes[1].axvline(thresh, color='r')
-
-            axes[2].imshow(mask, cmap=plt.cm.gray, extent=[axes[0].get_xlim()[0], axes[0].get_xlim()[1], axes[0].get_ylim()[0], axes[0].get_ylim()[1]])
-            axes[2].set_title('Thresholded')
-
-            x_indices = np.where(np.any(mask, axis=0))[0]
-            y_indices = np.where(np.any(mask, axis=1))[0]
-            x_min, x_max = x_indices[0], x_indices[-1]
-            y_min, y_max = y_indices[0], y_indices[-1]
-
-            axes[2].plot([x_min, x_max], [y_min, y_min], '--r')
-            axes[2].plot([x_min, x_max], [y_max, y_max], '--r')
-            axes[2].plot([x_min, x_min], [y_min, y_max], '--r')
-            axes[2].plot([x_max, x_max], [y_min, y_max], '--r')
-
-            axes[3].imshow(large_components_mask, cmap=plt.cm.gray, extent=[axes[0].get_xlim()[0], axes[0].get_xlim()[1], axes[0].get_ylim()[0], axes[0].get_ylim()[1]])
-            axes[3].set_title('Large components')
-
-            x_indices = np.where(np.any(large_components_mask, axis=0))[0]
-            y_indices = np.where(np.any(large_components_mask, axis=1))[0]
-            x_min, x_max = x_indices[0], x_indices[-1]
-            y_min, y_max = y_indices[0], y_indices[-1]
-
-            axes[3].plot([x_min, x_max], [y_min, y_min], '--r')
-            axes[3].plot([x_min, x_max], [y_max, y_max], '--r')
-            axes[3].plot([x_min, x_min], [y_min, y_max], '--r')
-            axes[3].plot([x_max, x_max], [y_min, y_max], '--r')
+        return large_components_mask
+    
+    def _configure_recon(self, backend='tigre'):
+        """
+        Configures the recon for the right engine. Checks the geometry type and data order.
+        """
+        if backend not in self._supported_backends:
+            raise ValueError("Backend unsupported. Supported backends: {}".format(self._supported_backends))
 
-            plt.tight_layout()
+        module = importlib.import_module(f'cil.plugins.{backend}')
 
-        return mask, large_components_mask
\ No newline at end of file
+        return module.FBP
diff --git a/Wrappers/Python/test/test_shrink_volume.py b/Wrappers/Python/test/test_shrink_volume.py
new file mode 100644
index 0000000000..b2518b2221
--- /dev/null
+++ b/Wrappers/Python/test/test_shrink_volume.py
@@ -0,0 +1,126 @@
+#  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 unittest
+from unittest.mock import patch
+import numpy as np
+
+from cil.framework import ImageGeometry, ImageData
+from cil.framework.labels import AcquisitionType
+
+from cil.utilities.shrink_volume import VolumeShrinker
+from cil.utilities import dataexample
+
+from utils import has_astra, has_tigre, has_nvidia, initialise_tests
+
+initialise_tests()
+
+
+class TestVolumeShrinker(unittest.TestCase):
+
+    def setUp(self):
+        self.data_cone = dataexample.SIMULATED_CONE_BEAM_DATA.get()
+        self.data_parallel = dataexample.SYNCHROTRON_PARALLEL_BEAM_DATA.get()
+        arr = np.zeros((10, 10, 10), dtype=np.float32)
+        arr[2:8, 2:8, 2:8] = 0.2
+        arr[3:7, 3:7, 3:7] = 1
+        
+        self.test_recon = ImageData(arr, geometry=ImageGeometry(10,10,10))
+        
+
+    @unittest.skipUnless(has_tigre and has_nvidia, "TIGRE GPU not installed")
+    def test_init_tigre(self):
+        for data in [self.data_cone, self.data_parallel]:
+            data.reorder('tigre')
+            vs = VolumeShrinker(data, recon_backend='tigre')
+            self.assertEqual(vs.recon_backend, 'tigre')
+            self.assertEqual(vs.data, data)
+
+            data.reorder('astra')
+            with self.assertRaises(ValueError):
+                vs = VolumeShrinker(data, recon_backend='tigre')
+                vs.run()
+
+    @unittest.skipUnless(has_astra and has_nvidia, "Astra GPU not installed")
+    def test_init_astra(self):
+        for data in [self.data_cone, self.data_parallel]:
+            data.reorder('astra')
+            vs = VolumeShrinker(data, recon_backend='astra')
+            self.assertEqual(vs.recon_backend, 'astra')
+            self.assertEqual(vs.data, data)
+
+            data.reorder('tigre')
+            with self.assertRaises(ValueError):
+                vs = VolumeShrinker(data, recon_backend='astra')
+                vs.run()
+
+    @unittest.skipUnless(has_astra and has_nvidia, "Astra GPU not installed")
+    def test_run_manual_astra(self):
+
+        limits = {
+                'horizontal_x': (10, 50),
+                'horizontal_y': (20, 90),
+                'vertical': (2, 25)
+            }
+        
+        for data in [self.data_cone, self.data_parallel]:
+            data.reorder('astra')
+            vs = VolumeShrinker(data, recon_backend='astra')
+            new_ig = vs.run(limits=limits, preview=False, method='manual')
+
+            # expected sizes are stop - start
+            for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
+                self.assertEqual(new_ig.shape[new_ig.dimension_labels.index(dim)], limits[dim][1] - limits[dim][0])
+
+    @unittest.skipUnless(has_tigre and has_nvidia, "TIGRE GPU not installed")
+    def test_run_manual_astra(self):
+
+        limits = {
+                'horizontal_x': (10, 50),
+                'horizontal_y': (20, 90),
+                'vertical': (2, 25)
+            }
+        
+        for data in [self.data_cone, self.data_parallel]:
+            data.reorder('tigre')
+            vs = VolumeShrinker(data, recon_backend='tigre')
+            new_ig = vs.run(limits=limits, preview=False, method='manual')
+
+            # expected sizes are stop - start
+            for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
+                self.assertEqual(new_ig.shape[new_ig.dimension_labels.index(dim)], limits[dim][1] - limits[dim][0])
+
+    unittest.skipUnless(has_astra and has_nvidia, "Astra GPU not installed")
+    def test_reduce_reconstruction_volume_astra(self):
+        vs = VolumeShrinker(self.data_cone, recon_backend='astra')
+        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', kwargs={'threshold':0.5})
+        for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
+            self.assertEqual(bounds[dim], (3,6))
+
+        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', kwargs={'threshold':0.1})
+        for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
+            self.assertEqual(bounds[dim], (2,7))
+            
+        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='otsu', kwargs={})
+        for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
+            self.assertEqual(bounds[dim], (3,6))
+
+        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='otsu', kwargs={'buffer':1})
+        for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
+            self.assertEqual(bounds[dim], (2,7))
+

From e154b8b610487646eb702df589c0be26f4d45d05 Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Tue, 24 Feb 2026 10:57:56 +0000
Subject: [PATCH 06/11] Add colorbar in preview plots

---
 Wrappers/Python/cil/utilities/shrink_volume.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/Wrappers/Python/cil/utilities/shrink_volume.py b/Wrappers/Python/cil/utilities/shrink_volume.py
index 96f46aed16..9d79ae21e8 100644
--- a/Wrappers/Python/cil/utilities/shrink_volume.py
+++ b/Wrappers/Python/cil/utilities/shrink_volume.py
@@ -18,6 +18,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
 from scipy.ndimage import label
 from skimage.filters import threshold_multiotsu
 import importlib
@@ -171,9 +172,12 @@ def plot_with_bounds(self, recon, bounds, binning):
             x_size = recon.get_dimension_size(x_dim)*binning
             y_size = recon.get_dimension_size(y_dim)*binning
 
-            ax.imshow(recon.max(axis=dim).array, origin='lower', cmap='gray',
+            im = ax.imshow(recon.max(axis=dim).array, origin='lower', cmap='gray',
                     extent=[0, x_size, 0, y_size])
-            
+            divider = make_axes_locatable(ax)
+            cax = divider.append_axes('right', size='5%', pad=0.05)
+            fig.colorbar(im, cax=cax, orientation='vertical')
+
             if bounds is not None:
                 x_min, x_max = bounds[x_dim]
                 y_min, y_max = bounds[y_dim]
@@ -186,6 +190,7 @@ def plot_with_bounds(self, recon, bounds, binning):
             ax.set_xlabel(x_dim)
             ax.set_ylabel(y_dim)
             ax.set_title(f"Maximum values in direction: {dim}")
+        plt.tight_layout()
 
     def reduce_reconstruction_volume(self, recon, binning, method, kwargs):
         threshold = kwargs.pop('threshold', None)

From edc9b4445fb3d4e0debfff178b3142f3825f5cb0 Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Wed, 25 Feb 2026 14:32:19 +0000
Subject: [PATCH 07/11] Update CHANGELOG

---
 CHANGELOG.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 63e12e6484..1562852eda 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,8 @@
 * XX.X.X
   - New features:
     - LSQR algorithm added to the CIL algorithm class (#1975)
+    - Add `VolumeShrinker` tool to reduce the size of the reconstruction volume from an
+    `AcquisitionData` (#2221)
   - 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)

From 83002f668f53df354a46c8e4c3a909eab85c6a12 Mon Sep 17 00:00:00 2001
From: Hannah Robarts <77114597+hrobarts@users.noreply.github.com>
Date: Fri, 27 Feb 2026 08:35:45 +0000
Subject: [PATCH 08/11] Apply suggestions from code review

Co-authored-by: Laura Murgatroyd <60604372+lauramurgatroyd@users.noreply.github.com>
Signed-off-by: Hannah Robarts <77114597+hrobarts@users.noreply.github.com>
---
 Wrappers/Python/cil/utilities/shrink_volume.py | 2 +-
 Wrappers/Python/test/test_shrink_volume.py     | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/Wrappers/Python/cil/utilities/shrink_volume.py b/Wrappers/Python/cil/utilities/shrink_volume.py
index 9d79ae21e8..ee94aecf96 100644
--- a/Wrappers/Python/cil/utilities/shrink_volume.py
+++ b/Wrappers/Python/cil/utilities/shrink_volume.py
@@ -70,7 +70,7 @@ def run(self, limits=None, preview=True, method='manual', **kwargs):
             Default is False.
         
         method : string, optional
-            If 'manual', use manual_limits
+            If 'manual', use `limits`
             If 'threshold' crop the reconstruction volume based on a threshold
             value between sample and background. Pass threshold as a kwarg.
             If 'otsu', automatically detect and crop the reconstruction volume.
diff --git a/Wrappers/Python/test/test_shrink_volume.py b/Wrappers/Python/test/test_shrink_volume.py
index b2518b2221..5b549be49f 100644
--- a/Wrappers/Python/test/test_shrink_volume.py
+++ b/Wrappers/Python/test/test_shrink_volume.py
@@ -88,7 +88,7 @@ def test_run_manual_astra(self):
                 self.assertEqual(new_ig.shape[new_ig.dimension_labels.index(dim)], limits[dim][1] - limits[dim][0])
 
     @unittest.skipUnless(has_tigre and has_nvidia, "TIGRE GPU not installed")
-    def test_run_manual_astra(self):
+    def test_run_manual_tigre(self):
 
         limits = {
                 'horizontal_x': (10, 50),
@@ -108,7 +108,7 @@ def test_run_manual_astra(self):
     unittest.skipUnless(has_astra and has_nvidia, "Astra GPU not installed")
     def test_reduce_reconstruction_volume_astra(self):
         vs = VolumeShrinker(self.data_cone, recon_backend='astra')
-        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', kwargs={'threshold':0.5})
+        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', threshold=0.5)
         for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
             self.assertEqual(bounds[dim], (3,6))
 

From 6c322947cedb839497f05fe6dad3cd5fb1ca7742 Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Fri, 27 Feb 2026 12:12:49 +0000
Subject: [PATCH 09/11] Updates from review

---
 .../Python/cil/utilities/shrink_volume.py     | 122 ++++++++++++++----
 Wrappers/Python/test/test_shrink_volume.py    | 108 ++++++++++++++--
 2 files changed, 192 insertions(+), 38 deletions(-)

diff --git a/Wrappers/Python/cil/utilities/shrink_volume.py b/Wrappers/Python/cil/utilities/shrink_volume.py
index ee94aecf96..018f55e4ff 100644
--- a/Wrappers/Python/cil/utilities/shrink_volume.py
+++ b/Wrappers/Python/cil/utilities/shrink_volume.py
@@ -79,7 +79,7 @@ def run(self, limits=None, preview=True, method='manual', **kwargs):
         **kwargs:
         threshold: float, optional
             If using 'threshold' method, specify the intensity threshold
-            between sample and background. Default is None.
+            between sample and background in the reconstruction. Default is None.
 
         buffer: float, optional
             Buffer in pixels around the automatically detected limits. 
@@ -87,15 +87,46 @@ def run(self, limits=None, preview=True, method='manual', **kwargs):
 
         mask_radius: float, optional
             Radius of circular mask to apply on the reconstructed volume, before
-            automatically cropping the reconstruction volume. Default is None.
+            automatically cropping the reconstruction volume or displaying with 
+            preview. Default is None.
 
         otsu_classes: int, optional
             Number of material classes to use when automatically detecting the 
             reconstruction volume using Otsu thresholding method. Default is 2.
 
         min_component_size: int, optional
-            Minimum size in pixels of connected components to keep when automatically 
+            Minimum size in pixels of connected components to consider when automatically 
             cropping the reconstruction volume. Default is None.
+
+        Returns
+        -------
+        ig_reduced: ImageGeometry
+            The reduced size ImageGeometry
+
+        Example
+        -------
+        >>> vs = VolumeShrinker(data, recon_backend='astra')
+        >>> ig_reduced = vs.run(method='manual', limits={'horizontal_x':(10, 150)})
+        
+        Example
+        -------
+        >>> cil_log_level = logging.getLogger()
+        >>> cil_log_level.setLevel(logging.DEBUG)
+        >>> vs = VolumeShrinker(data, recon_backend='astra')
+        >>> ig_reduced = vs.run(method='threshold' threshold=0.9)
+
+        Example
+        -------
+        >>> cil_log_level = logging.getLogger()
+        >>> cil_log_level.setLevel(logging.DEBUG)
+        >>> vs = VolumeShrinker(data, recon_backend='astra')
+        >>> ig_reduced = vs.run(method='otsu', otsu_classes=3)
+
+        Note
+        ----
+        For the `otsu` and `threshold` methods, setting logging to 'debug' shows 
+        a plot of the threshold on a histogram of the data and the threshold mask 
+        on the reconstruction.
         """
 
         ig = self.data.geometry.get_ImageGeometry()
@@ -108,29 +139,42 @@ def run(self, limits=None, preview=True, method='manual', **kwargs):
                     if dim in ig.dimension_labels:
                         if v is None:
                             v = (0, ig.shape[ig.dimension_labels.index(dim)])
-                        elif v[0] is None:
-                            v[0] = 0
-                        elif v[1] is None:
-                            v[1] = ig.shape[ig.dimension_labels.index(dim)]
+                        else:
+                            if v[0] is None:
+                                v = list(v)
+                                v[0] = 0
+                                v = tuple(v)
+                            if v[1] is None:
+                                v = list(v)
+                                v[1] = ig.shape[ig.dimension_labels.index(dim)]
+                                v = tuple(v)
                         bounds[dim] = v
                     else:
                         raise ValueError("dimension {} not recognised, must be one of {}".format(dim, ig.dimension_labels))
 
         elif method.lower() in ['threshold', 'otsu']:
             mask_radius = kwargs.pop('mask_radius', None)
-            recon, binning = self.get_recon(mask_radius=mask_radius)
-            bounds = self.reduce_reconstruction_volume(recon, binning, method, kwargs)
+            recon, binning = self._get_recon(mask_radius=mask_radius)
+            threshold = kwargs.pop('threshold', None)
+            buffer = kwargs.pop('buffer', None)
+            min_component_size = kwargs.pop('min_component_size', None)
+            otsu_classes = kwargs.pop('otsu_classes', 2) 
+            bounds = self._reduce_reconstruction_volume(recon, binning, method, 
+                                                        threshold, buffer, min_component_size, otsu_classes)
 
         if preview:
             if method.lower() == 'manual':
                 mask_radius = kwargs.pop('mask_radius', None)
-                recon, binning = self.get_recon(mask_radius=mask_radius)
+                recon, binning = self._get_recon(mask_radius=mask_radius)
 
-            self.plot_with_bounds(recon, bounds, binning)
+            self._plot_with_bounds(recon, bounds, binning)
         
         return self.update_ig(self.data.geometry.get_ImageGeometry(), bounds)
 
     def update_ig(self, ig_unbinned, bounds):
+        """
+        Return a new unbinned ImageGeometry with the bounds applied
+        """
         if bounds is None:
             ig = ig_unbinned
         else:
@@ -139,7 +183,11 @@ def update_ig(self, ig_unbinned, bounds):
                     'vertical':(bounds['vertical'][0], bounds['vertical'][1], 1)})(ig_unbinned)
         return ig
     
-    def get_recon(self, mask_radius):
+    def _get_recon(self, mask_radius=None):
+        """
+        Gets a binned reconstruction from the dataset with an optional mask
+        Also returns the binning which has been used
+        """
         binning = min(int(np.ceil(self.data.geometry.config.panel.num_pixels[0] / 128)),16)
         angle_binning = np.ceil(self.data.get_dimension_size('angle')/(self.data.get_dimension_size('horizontal')*(np.pi/2)))
         roi = {
@@ -160,7 +208,12 @@ def get_recon(self, mask_radius):
 
         return recon, binning
 
-    def plot_with_bounds(self, recon, bounds, binning):
+    def _plot_with_bounds(self, recon, bounds, binning):
+        """
+        Plots the bounds on the maximum value in the reconstructed dataset along 
+        each direction.
+        """
+
         fig, axs = plt.subplots(nrows=1, ncols=recon.ndim, figsize=(14, 6))
 
         dims = recon.dimension_labels
@@ -192,11 +245,16 @@ def plot_with_bounds(self, recon, bounds, binning):
             ax.set_title(f"Maximum values in direction: {dim}")
         plt.tight_layout()
 
-    def reduce_reconstruction_volume(self, recon, binning, method, kwargs):
-        threshold = kwargs.pop('threshold', None)
-        buffer = kwargs.pop('buffer', None)
-        min_component_size = kwargs.pop('min_component_size', None)
-        otsu_classes = kwargs.pop('otsu_classes', 2)   
+    def _reduce_reconstruction_volume(self, recon, binning, method, 
+                                      threshold=None, buffer=None, min_component_size=None, otsu_classes=2):
+        """
+        Automatically finds the boundaries of the sample in a reconstructed
+        volume based on a threshold between sample and background. 
+        If method=`threshold`, the threshold must be provided.
+        If method=`otsu`, the threshold is calculated from an Otsu filter. The
+        number of `otsu_classes` can be passed. 
+        """
+          
 
         dims = recon.dimension_labels
         all_bounds = {dim: [] for dim in dims}
@@ -205,7 +263,12 @@ def reduce_reconstruction_volume(self, recon, binning, method, kwargs):
             arr = recon.max(axis=dim).array
             
             if method.lower() == 'threshold':
-                mask = arr > threshold
+                if threshold is not None:
+                    if threshold >= arr.max():
+                        raise ValueError("Threshold is greater than maximum value in dimension: no limits can be found. Try specifying a lower threshold.")
+                    mask = arr > threshold
+                else:
+                    raise ValueError("You must supply a threshold argument if method='threshold'")
 
             elif method.lower() == 'otsu':
                 n_bins = 256
@@ -214,7 +277,7 @@ def reduce_reconstruction_volume(self, recon, binning, method, kwargs):
                 mask = arr > threshold
 
             if min_component_size is not None:
-                mask = self.threshold_large_components(mask, min_component_size)
+                mask = self._threshold_large_components(mask, min_component_size)
 
             x_indices = np.where(np.any(mask, axis=0))[0]
             y_indices = np.where(np.any(mask, axis=1))[0]
@@ -225,19 +288,22 @@ def reduce_reconstruction_volume(self, recon, binning, method, kwargs):
                 fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(8, 2.5))
 
                 axes[0].imshow(arr, cmap=plt.cm.gray)
-                axes[0].set_title('Original')
+                axes[0].set_title('Maximum values')
 
                 axes[1].hist(arr.ravel(), bins=100)
-                axes[1].set_title('Histogram')
+                axes[1].set_title('Histogram of maximum values')
                 axes[1].axvline(threshold, color='r')
 
                 axes[2].imshow(mask, cmap=plt.cm.gray, extent=[axes[0].get_xlim()[0], axes[0].get_xlim()[1], axes[0].get_ylim()[0], axes[0].get_ylim()[1]])
-                axes[2].set_title('Thresholded')
+                axes[2].set_title('Calculated threshold mask')
 
                 axes[2].plot([x_min, x_max], [y_min, y_min], '--r')
                 axes[2].plot([x_min, x_max], [y_max, y_max], '--r')
                 axes[2].plot([x_min, x_min], [y_min, y_max], '--r')
                 axes[2].plot([x_max, x_max], [y_min, y_max], '--r')
+
+                plt.suptitle(dim)
+                
                 plt.tight_layout()
 
             axis = recon.get_dimension_axis(dim)
@@ -273,7 +339,13 @@ def reduce_reconstruction_volume(self, recon, binning, method, kwargs):
             
         return bounds
         
-    def threshold_large_components(self, mask, min_component_size):
+    def _threshold_large_components(self, mask, min_component_size):
+        """
+        Modify a threshold mask to ignore clusters of pixel values above the threshold 
+        if the cluster is below a specified `min_component_size` in pixels. This 
+        can be useful for noisy datasets where a few noisy pixels above the threshold 
+        may appear in the background and erroneously increase the bounds.
+        """
         labeled_mask, _ = label(mask)
         component_sizes = np.bincount(labeled_mask.ravel())
         min_component_size = 10
@@ -286,7 +358,7 @@ def threshold_large_components(self, mask, min_component_size):
     
     def _configure_recon(self, backend='tigre'):
         """
-        Configures the recon for the right engine. Checks the geometry type and data order.
+        Configures the recon for the right engine.
         """
         if backend not in self._supported_backends:
             raise ValueError("Backend unsupported. Supported backends: {}".format(self._supported_backends))
diff --git a/Wrappers/Python/test/test_shrink_volume.py b/Wrappers/Python/test/test_shrink_volume.py
index 5b549be49f..bee80ae35f 100644
--- a/Wrappers/Python/test/test_shrink_volume.py
+++ b/Wrappers/Python/test/test_shrink_volume.py
@@ -39,8 +39,12 @@ def setUp(self):
         arr = np.zeros((10, 10, 10), dtype=np.float32)
         arr[2:8, 2:8, 2:8] = 0.2
         arr[3:7, 3:7, 3:7] = 1
-        
         self.test_recon = ImageData(arr, geometry=ImageGeometry(10,10,10))
+
+        arr = np.zeros((10, 10, 10), dtype=np.float32)
+        arr[2:8, 2:8, 3:7] = 0.2
+        arr[3:7, 3:7, 4:6] = 1
+        self.test_recon_asymmetrical = ImageData(arr, geometry=ImageGeometry(10,10,10))
         
 
     @unittest.skipUnless(has_tigre and has_nvidia, "TIGRE GPU not installed")
@@ -80,12 +84,38 @@ def test_run_manual_astra(self):
         
         for data in [self.data_cone, self.data_parallel]:
             data.reorder('astra')
+            old_ig = data.geometry.get_ImageGeometry()
             vs = VolumeShrinker(data, recon_backend='astra')
             new_ig = vs.run(limits=limits, preview=False, method='manual')
-
-            # expected sizes are stop - start
+            # get the voxel size and centers that correspond to each dimension
+            voxel_sizes = {'horizontal_x': old_ig.voxel_size_x,'horizontal_y': old_ig.voxel_size_y,'vertical': old_ig.voxel_size_z}
+            centers = {'horizontal_x': new_ig.center_x, 'horizontal_y': new_ig.center_y, 'vertical': new_ig.center_z}
+            
             for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
-                self.assertEqual(new_ig.shape[new_ig.dimension_labels.index(dim)], limits[dim][1] - limits[dim][0])
+                # expected sizes are stop - start
+                new_size = limits[dim][1] - limits[dim][0]
+                ind = new_ig.dimension_labels.index(dim)
+                self.assertEqual(new_ig.shape[ind], new_size)
+
+                old_size = old_ig.shape[ind]
+
+                # expected center is (new center-old center)*voxel size
+                center = ((limits[dim][0]+(new_size/2))-(old_size/2))*voxel_sizes[dim]
+                self.assertEqual(center, centers[dim])
+
+
+        # test some non-sensical limits
+        limits = {
+                'horizontal_x': (10, 5),
+                'horizontal_y': (20, 9),
+                'vertical': (12, 2)
+            }
+
+        for data in [self.data_cone, self.data_parallel]:
+            data.reorder('astra')
+            vs = VolumeShrinker(data, recon_backend='astra')
+            with self.assertRaises(ValueError):
+                new_ig = vs.run(limits=limits, preview=False, method='manual')
 
     @unittest.skipUnless(has_tigre and has_nvidia, "TIGRE GPU not installed")
     def test_run_manual_tigre(self):
@@ -98,29 +128,81 @@ def test_run_manual_tigre(self):
         
         for data in [self.data_cone, self.data_parallel]:
             data.reorder('tigre')
+            old_ig = data.geometry.get_ImageGeometry()
             vs = VolumeShrinker(data, recon_backend='tigre')
             new_ig = vs.run(limits=limits, preview=False, method='manual')
 
-            # expected sizes are stop - start
+            # get the voxel size and centers that correspond to each dimension
+            voxel_sizes = {'horizontal_x': old_ig.voxel_size_x,'horizontal_y': old_ig.voxel_size_y,'vertical': old_ig.voxel_size_z}
+            centers = {'horizontal_x': new_ig.center_x, 'horizontal_y': new_ig.center_y, 'vertical': new_ig.center_z}
+            
             for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
-                self.assertEqual(new_ig.shape[new_ig.dimension_labels.index(dim)], limits[dim][1] - limits[dim][0])
+                # expected sizes are stop - start
+                new_size = limits[dim][1] - limits[dim][0]
+                ind = new_ig.dimension_labels.index(dim)
+                self.assertEqual(new_ig.shape[ind], new_size)
+
+                old_size = old_ig.shape[ind]
+
+                # expected center is (new center-old center)*voxel size
+                center = ((limits[dim][0]+(new_size/2))-(old_size/2))*voxel_sizes[dim]
+                self.assertEqual(center, centers[dim])
+
+        # test some non-sensical limits
+        limits = {
+                'horizontal_x': (10, 5),
+                'horizontal_y': (20, 9),
+                'vertical': (12, 2)
+            }
+        
+        for data in [self.data_cone, self.data_parallel]:
+            data.reorder('tigre')
+            vs = VolumeShrinker(data, recon_backend='tigre')
+            with self.assertRaises(ValueError):
+                new_ig = vs.run(limits=limits, preview=False, method='manual')
 
-    unittest.skipUnless(has_astra and has_nvidia, "Astra GPU not installed")
-    def test_reduce_reconstruction_volume_astra(self):
+    def test_reduce_reconstruction_volume(self):
         vs = VolumeShrinker(self.data_cone, recon_backend='astra')
-        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', threshold=0.5)
+
+        # test error with threshold higher than max value in volume
+        with self.assertRaises(ValueError):
+            bounds = vs._reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', threshold=500)
+
+        # test expected boundaries are found
+        bounds = vs._reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', threshold=0.5)
         for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
             self.assertEqual(bounds[dim], (3,6))
 
-        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', kwargs={'threshold':0.1})
+        bounds = vs._reduce_reconstruction_volume(self.test_recon, binning=1, method='threshold', threshold=0.1)
         for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
             self.assertEqual(bounds[dim], (2,7))
-            
-        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='otsu', kwargs={})
+
+        bounds = vs._reduce_reconstruction_volume(self.test_recon, binning=1, method='otsu')
         for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
             self.assertEqual(bounds[dim], (3,6))
 
-        bounds = vs.reduce_reconstruction_volume(self.test_recon, binning=1, method='otsu', kwargs={'buffer':1})
+        bounds = vs._reduce_reconstruction_volume(self.test_recon, binning=1, method='otsu', buffer=1)
         for dim in ['horizontal_x', 'horizontal_y', 'vertical']:
             self.assertEqual(bounds[dim], (2,7))
 
+        # test the asymmetrical volume
+        bounds = vs._reduce_reconstruction_volume(self.test_recon_asymmetrical, binning=1, method='threshold', threshold=0.5)
+        self.assertEqual(bounds['horizontal_x'], (4,5))
+        self.assertEqual(bounds['horizontal_y'], (3,6))
+        self.assertEqual(bounds['vertical'], (3,6))
+
+        bounds = vs._reduce_reconstruction_volume(self.test_recon_asymmetrical, binning=1, method='threshold', threshold=0.1)
+        self.assertEqual(bounds['horizontal_x'], (3,6))
+        self.assertEqual(bounds['horizontal_y'], (2,7))
+        self.assertEqual(bounds['vertical'], (2,7))
+
+        bounds = vs._reduce_reconstruction_volume(self.test_recon_asymmetrical, binning=1, method='otsu')
+        self.assertEqual(bounds['horizontal_x'], (4,5))
+        self.assertEqual(bounds['horizontal_y'], (3,6))
+        self.assertEqual(bounds['vertical'], (3,6))
+
+        bounds = vs._reduce_reconstruction_volume(self.test_recon_asymmetrical, binning=1, method='otsu', buffer=1)
+        self.assertEqual(bounds['horizontal_x'], (3,6))
+        self.assertEqual(bounds['horizontal_y'], (2,7))
+        self.assertEqual(bounds['vertical'], (2,7))
+

From ea98fcb5c2309ef70306289df971e617dbd28866 Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Fri, 27 Feb 2026 12:18:53 +0000
Subject: [PATCH 10/11] Add volumeshrinker to rst file

---
 docs/source/utilities.rst | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/docs/source/utilities.rst b/docs/source/utilities.rst
index 4e58d3dbdc..a9926e552f 100644
--- a/docs/source/utilities.rst
+++ b/docs/source/utilities.rst
@@ -136,4 +136,15 @@ link_islicer - link islicer objects by index
    :inherited-members:
 
 
+Shrink volume
+=============
+VolumeShrinker - Display 2D slices
+--------------------------
+
+.. autoclass:: cil.utilities.shink_volume.VolumeShrinker
+   :members:
+   :inherited-members:
+
+
+
 :ref:`Return Home <mastertoc>`

From 83df471f062598d1c62ead311cc86559590d18f0 Mon Sep 17 00:00:00 2001
From: hrobarts <hannah.robarts@stfc.ac.uk>
Date: Fri, 27 Feb 2026 13:38:19 +0000
Subject: [PATCH 11/11] Fix docs

---
 docs/source/utilities.rst | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/docs/source/utilities.rst b/docs/source/utilities.rst
index a9926e552f..79aa4105ef 100644
--- a/docs/source/utilities.rst
+++ b/docs/source/utilities.rst
@@ -138,13 +138,11 @@ link_islicer - link islicer objects by index
 
 Shrink volume
 =============
-VolumeShrinker - Display 2D slices
+VolumeShrinker - create a cropped reconstruction volume from image data
 --------------------------
 
-.. autoclass:: cil.utilities.shink_volume.VolumeShrinker
+.. autoclass:: cil.utilities.shrink_volume.VolumeShrinker
    :members:
    :inherited-members:
 
-
-
 :ref:`Return Home <mastertoc>`