diff --git a/src/readii/analyze/plot_correlation.py b/src/readii/analyze/plot_correlation.py
index 9cd800f..6c170d4 100644
--- a/src/readii/analyze/plot_correlation.py
+++ b/src/readii/analyze/plot_correlation.py
@@ -418,8 +418,8 @@ def plotCrossCorrHeatmap(correlation_matrix:pd.DataFrame,
     cross_corr_heatmap = plotCorrelationHeatmap(cross_corr, 
                                                 diagonal=False, 
                                                 cmap=cmap, 
-                                                xlabel=vertical_feature_name, 
-                                                ylabel=horizontal_feature_name,
+                                                xlabel=horizontal_feature_name, 
+                                                ylabel=vertical_feature_name,
                                                 title=f"{correlation_method.capitalize()} Cross Correlations", 
                                                 subtitle=f"{vertical_feature_name} vs {horizontal_feature_name}")
     
diff --git a/src/readii/process/images/__init__.py b/src/readii/process/images/__init__.py
new file mode 100644
index 0000000..2d55d48
--- /dev/null
+++ b/src/readii/process/images/__init__.py
@@ -0,0 +1,25 @@
+"""Module for processing and manipulating images."""
+
+from .crop import (
+    apply_bounding_box_limits,
+    check_bounding_box_single_dimension,
+    crop_image_to_mask,
+    crop_to_bounding_box,
+    crop_to_centroid,
+    crop_to_maxdim_cube,
+    find_bounding_box,
+    find_centroid,
+    validate_new_dimensions,
+)
+
+__all__ = [
+    "crop_image_to_mask",
+    "find_bounding_box",
+    "find_centroid",
+    "validate_new_dimensions",
+    "apply_bounding_box_limits",
+    "check_bounding_box_single_dimension",
+    "crop_to_maxdim_cube",
+    "crop_to_bounding_box",
+    "crop_to_centroid",
+]
\ No newline at end of file
diff --git a/src/readii/process/images/crop.py b/src/readii/process/images/crop.py
new file mode 100644
index 0000000..4469c45
--- /dev/null
+++ b/src/readii/process/images/crop.py
@@ -0,0 +1,433 @@
+
+from typing import Literal, Optional
+
+import numpy as np
+import SimpleITK as sitk
+from imgtools.ops.functional import resize
+from radiomics import imageoperations
+
+from readii.image_processing import getROIVoxelLabel
+from readii.utils import logger
+
+from readii.process.images.utils.bounding_box import Centroid, Coordinate, Size3D
+
+def validate_new_dimensions(image:sitk.Image,
+                            new_dimensions:tuple | int
+                            ) -> None:
+    """Validate that the input new dimensions are valid for the image.
+
+    Parameters
+    ----------
+    image : sitk.Image
+        Image to validate the new dimensions for.
+    new_dimensions : tuple or int
+        Tuple of values representing the new dimensions to validate, or a single integer representing the number of dimensions.
+
+    Raises
+    ------
+    ValueError
+        If the new dimensions are not valid for the image.
+    """
+    # Check that the number of dimensions in the new dimensions matches the number of dimensions in the image
+    if isinstance(new_dimensions, tuple):
+        if len(new_dimensions) != image.GetDimension():
+            msg = f"Number of dimensions in new_dimensions ({len(new_dimensions)}) does not match the number of dimensions in the image ({image.GetDimension()})."
+            logger.exception(msg)
+            raise ValueError(msg)
+
+    elif isinstance(new_dimensions, int):
+        if new_dimensions != image.GetDimension():
+            msg = f"Number of dimensions in new_dimensions ({new_dimensions}) does not match the number of dimensions in the image ({image.GetDimension()})."
+            logger.exception(msg)
+            raise ValueError(msg)
+    
+    else:
+        msg = "New dimensions must be a tuple of integers or a single integer."
+        logger.exception(msg)
+        raise ValueError(msg)
+
+
+
+def find_bounding_box(mask:sitk.Image,
+                      min_dim_size:int = 4
+                      ) -> tuple[Coordinate, Coordinate, Size3D]:
+    """Find the bounding box of a region of interest (ROI) in a given binary mask image.
+    
+    Parameters
+    ----------
+    mask : sitk.Image
+        Mask image to find the bounding box within.
+    min_dim_size : int, optional
+        Minimum size of the bounding box along each dimension. The default is 4.
+    
+    Returns
+    -------
+    bounding_box : tuple[Coordinate, Coordinate, Size3D]
+        Tuple containing the minimum and maximum coordinates for a bounding box, along with the size of the bounding box
+    """
+    # Convert the mask to a uint8 image
+    mask_uint = sitk.Cast(mask, sitk.sitkUInt8)
+    stats = sitk.LabelShapeStatisticsImageFilter()
+    stats.Execute(mask_uint)
+    # Get the bounding box starting/minimum coordinates and size
+    coord_x, coord_y, coord_z, size_x, size_y, size_z = stats.GetBoundingBox(1)
+    
+    # Ensure minimum size of 4 pixels along each dimension (requirement of cropping method)
+    size_x = max(size_x, min_dim_size)
+    size_y = max(size_y, min_dim_size)
+    size_z = max(size_z, min_dim_size)
+    
+    # Create an object to store the bounding box size in same manner as coordinates
+    bbox_size = Size3D(size_x, size_y, size_z)
+
+    # Create an object to store the minimum bounding box coordinate
+    bbox_min_coord = Coordinate(coord_x, coord_y, coord_z)
+    
+    # Calculate the maximum coordinate of the bounding box by adding the size to the starting coordinate
+    bbox_max_coord = bbox_min_coord + bbox_size
+
+    return bbox_min_coord, bbox_max_coord, bbox_size
+
+
+def check_bounding_box_single_dimension(bb_min_val:int, 
+                                        bb_max_val:int,
+                                        expected_dim:int, 
+                                        img_dim:int
+                                        ) -> tuple[int,int]:
+    """Check if minimum and maximum values for a single bounding box dimension fall within the same dimension in the image the bounding box was made for.
+    
+    Parameters
+    ----------
+    bb_min_val : int
+        Minimum value for the bounding box dimension.
+    bb_max_val : int
+        Maximum value for the bounding box dimension.
+    expected_dim : int
+        Expected dimension of the bounding box.
+    img_dim : int
+        Dimension of the image the bounding box was made for.
+    
+    Returns
+    -------
+    bb_min_val : int
+        Updated minimum value for the bounding box dimension.
+    bb_max_val : int
+        Updated maximum value for the bounding box dimension.
+
+    Examples
+    --------
+    >>> check_bounding_box_single_dimension(0, 10, 20, 30)
+    (0, 10)
+    >>> check_bounding_box_single_dimension(30, 40, 10, 30)
+    (20, 30)
+    >>> check_bounding_box_single_dimension(bb_x_min, bb_x_max, expected_dim_x, img_dim_x)
+    """
+    # Check if the minimum bounding box value is outside the image
+    if bb_min_val < 0:
+        # Set the minimum value to 0 (edge of image) and the max value to the minimum of the expected dimension or edge of image
+        bb_min_val, bb_max_val = 0, min(expected_dim, img_dim)
+    
+    # Check if the maximum bounding box value is outside the image
+    if bb_max_val > img_dim:
+        # Set the minimum value to the maximum of the image dimension or edge of image and the max value to the edge of image
+        bb_min_val, bb_max_val = max(0, img_dim - expected_dim), img_dim
+
+    return bb_min_val, bb_max_val
+
+
+
+def apply_bounding_box_limits(image:sitk.Image,
+                              bounding_box:tuple[int,int,int,int,int,int],
+                              expected_dimensions:tuple[int,int,int]
+                              ) -> tuple:
+    """Check that bounding box coordinates are within the image dimensions. If not, move bounding box to the edge of the image and expand to expected dimension.
+    
+    Parameters
+    ----------
+    image : sitk.Image
+        Image to check the bounding box coordinates against.
+    bounding_box : tuple[int,int,int,int,int,int]
+        Bounding box to check the coordinates of. Order is [min_x, max_x, min_y, max_y, min_z, max_z].
+    expected_dimensions : tuple[int,int,int]
+        Expected dimensions of the bounding box. Used if the bounding box needs to be shifted to the edge of the image.
+    
+    Returns
+    -------
+    min_x, min_y, min_z, max_x, max_y, max_z : tuple[int,int,int,int,int,int]
+        Updated bounding box coordinates. Order is [min_x, max_x, min_y, max_y, min_z, max_z].
+    """
+    # Get the size of the image to use to determine if crop dimensions are larger than the image
+    img_x, img_y, img_z = image.GetSize()
+
+    # Extract the bounding box coordinates
+    min_x, max_x, min_y, max_y, min_z, max_z = bounding_box
+
+    # Check each bounding box dimensions coordinates and move to image edge if not within image
+    min_x, max_x = check_bounding_box_single_dimension(min_x, max_x, expected_dimensions[0], img_x)
+    min_y, max_y = check_bounding_box_single_dimension(min_y, max_y, expected_dimensions[1], img_y)
+    min_z, max_z = check_bounding_box_single_dimension(min_z, max_z, expected_dimensions[2], img_z)
+
+    return min_x, max_x, min_y, max_y, min_z, max_z
+
+
+
+def find_centroid(mask:sitk.Image) -> np.ndarray:
+    """Find the centroid of a region of interest (ROI) in a given binary mask image.
+
+    Parameters
+    ----------
+    mask : sitk.Image
+        Mask image to find the centroid within.
+
+    Returns
+    -------
+    centroid : np.ndarray
+        Numpy array containing the coordinates of the ROI centroid.
+    """
+    # Convert the mask to a uint8 image
+    mask_uint = sitk.Cast(mask, sitk.sitkUInt8)
+    stats = sitk.LabelShapeStatisticsImageFilter()
+    stats.Execute(mask_uint)
+    # Get the centroid coordinates as a physical point in the mask
+    centroid_coords = stats.GetCentroid(1)
+    # Convert the physical point to an index in the mask array
+    centroid_x, centroid_y, centroid_z = mask.TransformPhysicalPointToIndex(centroid_coords)
+    # Convert to a Centroid object
+    centroid = Centroid(centroid_x, centroid_y, centroid_z)
+    return centroid
+
+
+
+def crop_to_centroid(image:sitk.Image,
+                     centroid:tuple,
+                     crop_dimensions:tuple,
+                     ) -> sitk.Image:
+    """Crop an image centered on the centroid with specified crop dimension. No resizing/resampling is performed.
+
+    Parameters
+    ----------
+    image : sitk.Image
+        Image to crop.
+    centroid : tuple
+        Tuple of integers representing the centroid of the image to crop. Must have the same number of dimensions as the image.
+    crop_dimensions : tuple
+        Tuple of integers representing the dimensions to crop the image to. Must have the same number of dimensions as the image.
+
+    Returns
+    -------
+    cropped_image : sitk.Image
+        Cropped image.
+    """
+    # Check that the number of dimensions in the crop dimensions matches the number of dimensions in the image
+    validate_new_dimensions(image, crop_dimensions)
+
+    # Check that the centroid dimensions match the image dimensions
+    validate_new_dimensions(image, centroid)
+
+    min_x = int(centroid[0] - crop_dimensions[0] // 2)
+    max_x = int(centroid[0] + crop_dimensions[0] // 2)
+    min_y = int(centroid[1] - crop_dimensions[1] // 2)
+    max_y = int(centroid[1] + crop_dimensions[1] // 2)
+    min_z = int(centroid[2] - crop_dimensions[2] // 2)
+    max_z = int(centroid[2] + crop_dimensions[2] // 2)
+
+    # Test if bounding box coordinates are within the image, move to image edge if not
+    min_x, max_x, min_y, max_y, min_z, max_z = apply_bounding_box_limits(image, 
+                                                                         bounding_box = (min_x, max_x, min_y, max_y, min_z, max_z),
+                                                                         expected_dimensions = crop_dimensions)
+
+    return image[min_x:max_x, min_y:max_y, min_z:max_z]
+
+
+
+def crop_to_bounding_box(image:sitk.Image,
+                         bounding_box:tuple[int,int,int,int,int,int],
+                         resize_dimensions:tuple[int,int,int]
+                         ) -> sitk.Image:
+    """Crop an image to a given bounding box and resize to a specified crop dimensions.
+
+    Parameters
+    ----------
+    image : sitk.Image
+        Image to crop.
+    bounding_box : tuple[int,int,int,int,int,int]
+        Bounding box to crop the image to. Order is [min_x, max_x, min_y, max_y, min_z, max_z].
+    resize_dimensions : tuple[int,int,int]
+        Dimensions to resize the image to.
+    
+    Returns
+    -------
+    cropped_image : sitk.Image
+        Cropped image.
+    """
+    # Check that the number of dimensions in the crop dimensions matches the number of dimensions in the image
+    validate_new_dimensions(image, resize_dimensions)
+
+    # Check that the number of bounding box dimensions match the image dimensions
+    validate_new_dimensions(image, int(len(bounding_box)/2))
+
+    # Current bounding box dimensions
+    current_image_dimensions = (max_x - min_x, max_y - min_y, max_z - min_z)
+    # bounding_box[1] - bounding_box[0], bounding_box[]
+
+    # Test if bounding box coordinates are within the image, move to image edge if not
+    min_x, max_x, min_y, max_y, min_z, max_z = apply_bounding_box_limits(image, 
+                                                                         bounding_box, 
+                                                                         expected_dimensions=current_image_dimensions)
+
+    # Crop image to the bounding box
+    img_crop = image[min_x:max_x, min_y:max_y, min_z:max_z]
+    # Resample the image to the new dimensions and spacing
+    img_crop = resize(img_crop, size = resize_dimensions)
+    return img_crop
+
+
+
+def crop_to_maxdim_cube(image:sitk.Image,
+                        bounding_box:tuple[int,int,int,int,int,int],
+                        resize_dimensions:tuple[int,int,int]
+                        ) -> sitk.Image:
+    """
+    Crop given image to a cube based on the max dim from a bounding box and resize to specified input size.
+
+    Parameters
+    ----------
+    image : sitk.Image
+        Image to crop.
+    bounding_box : tuple[int,int,int,int,int,int]
+        Bounding box to find maximum dimension from. The order is (min_x, min_y, min_z, max_x, max_y, max_z).
+    resize_dimensions : tuple[int,int,int]
+        Crop dimensions to resize the image to.
+
+    Returns
+    -------
+        sitk.Image: The cropped and resized image.
+    """
+    # Check that the number of dimensions in the crop dimensions matches the number of dimensions in the image
+    validate_new_dimensions(image, resize_dimensions)
+
+    # Check that the number of bounding box dimensions match the image dimensions
+    validate_new_dimensions(image, len(bounding_box)//2)
+
+    # Extract out the bounding box coordinates
+    min_x, max_x, min_y, max_y, min_z, max_z = bounding_box
+
+    # Get maximum dimension of bounding box
+    max_dim = max(max_x - min_x, max_y - min_y, max_z - min_z)
+    mean_x = int((max_x + min_x) // 2)
+    mean_y = int((max_y + min_y) // 2)
+    mean_z = int((max_z + min_z) // 2)
+
+    # define new bounding boxes based on the maximum dimension of ROI bounding box
+    min_x = int(mean_x - max_dim // 2)
+    max_x = int(mean_x + max_dim // 2)
+    min_y = int(mean_y - max_dim // 2)
+    max_y = int(mean_y + max_dim // 2)
+    min_z = int(mean_z - max_dim // 2)
+    max_z = int(mean_z + max_dim // 2)
+
+    # Test if bounding box coordinates are within the image, move to image edge if not
+    min_x, max_x, min_y, max_y, min_z, max_z = apply_bounding_box_limits(image, 
+                                                                         bounding_box = (min_x, max_x, min_y, max_y, min_z, max_z),
+                                                                         expected_dimensions = (max_dim, max_dim, max_dim))
+    # Crop image to the cube bounding box
+    img_crop = image[min_x:max_x, min_y:max_y, min_z:max_z]
+    # Resample the image to the new dimensions and spacing
+    img_crop = resize(img_crop, size=resize_dimensions)
+    return img_crop
+
+
+
+def crop_with_pyradiomics(image:sitk.Image,  
+                          mask:sitk.Image,  
+                          mask_label: int | None = None  
+                          ) -> tuple[sitk.Image, sitk.Image]:  
+    """Crop an image to a bounding box around a region of interest in the mask using PyRadiomics functions.  
+
+    Parameters
+    ----------
+    image : sitk.Image
+        Image to crop.
+    mask : sitk.Image
+        Mask to crop the image to.
+    mask_label : int, optional
+        Label of the region of interest to crop to in the mask. If not provided, will use the label of the first non-zero voxel in the mask.
+    
+    Returns
+    -------
+    image_crop : sitk.Image
+        Cropped image.
+    mask_crop : sitk.Image
+        Cropped mask.
+    """
+    # Get the label of the region of interest in the mask if not provided
+    if mask_label is None:
+        mask_label = getROIVoxelLabel(mask)
+    
+    # Check that CT and segmentation correspond, segmentationLabel is present, and dimensions match
+    bounding_box, corrected_mask = imageoperations.checkMask(image, mask, label=mask_label)
+
+    # Update the mask if correction was generated by checkMask
+    if corrected_mask:
+        mask = corrected_mask
+    
+    # Crop the image and mask to the bounding box
+    image_crop, mask_crop = imageoperations.cropToTumorMask(image, mask, bounding_box)
+
+    return image_crop, mask_crop
+
+
+
+def crop_image_to_mask(image:sitk.Image,
+                       mask:sitk.Image,
+                       crop_method:Literal["bounding_box", "centroid", "cube", "pyradiomics"],
+                       resize_dimensions:Optional[tuple[int,int,int]]
+                       ) -> tuple[sitk.Image, sitk.Image]:
+    """Crop an image and mask to an ROI in the mask and resize to a specified crop dimensions.
+
+    Parameters
+    ----------
+    image : sitk.Image
+        Image to crop.
+    mask : sitk.Image
+        Mask to crop the image to. Will also be cropped.
+    crop_method : str, optional
+        Method to use to crop the image to the mask. Must be one of "bounding_box", "centroid", "cube, or "pyradiomics".
+    resize_dimensions : tuple[int,int,int]
+        Dimensions to resize the image to.
+    
+    Returns
+    -------
+    cropped_image : sitk.Image
+        Cropped image.
+    cropped_mask : sitk.Image
+        Cropped mask.
+    """
+    if resize_dimensions is None and crop_method != "pyradiomics":  
+        msg = f"resize_dimensions is required for crop_method '{crop_method}'."  
+        raise ValueError(msg)
+
+    match crop_method:            
+        case "bounding_box":
+            bbox_coords = find_bounding_box(mask)
+            cropped_image = crop_to_bounding_box(image, bbox_coords, resize_dimensions)
+            cropped_mask = crop_to_bounding_box(mask, bbox_coords, resize_dimensions)
+        
+        case "centroid":
+            centroid = find_centroid(mask)
+            cropped_image = crop_to_centroid(image, centroid, resize_dimensions)
+            cropped_mask = crop_to_centroid(mask, centroid, resize_dimensions)
+        
+        case "cube":
+            bbox_coords = find_bounding_box(mask)
+            cropped_image = crop_to_maxdim_cube(image, bbox_coords, resize_dimensions)
+            cropped_mask = crop_to_maxdim_cube(mask, bbox_coords, resize_dimensions)
+        
+        case "pyradiomics":
+            cropped_image, cropped_mask = crop_with_pyradiomics(image, mask)
+        
+        case _:
+            msg = f"Invalid crop method: {crop_method}. Must be one of 'bounding_box', 'centroid', 'cube', or 'pyradiomics'."
+            raise ValueError(msg)
+    
+    return cropped_image, cropped_mask
\ No newline at end of file
diff --git a/src/readii/process/images/utils/bounding_box.py b/src/readii/process/images/utils/bounding_box.py
new file mode 100644
index 0000000..5e0aceb
--- /dev/null
+++ b/src/readii/process/images/utils/bounding_box.py
@@ -0,0 +1,63 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+
+@dataclass
+class Point3D:
+    """Represent a point in 3D space."""
+    x: int
+    y: int
+    z: int
+
+    @property
+    def as_tuple(self) -> tuple[int, int, int]:
+        """Return the point as a tuple."""
+        return self.x, self.y, self.z
+
+    def __add__(self, other: Point3D) -> Point3D:
+        """Add two points."""
+        return Point3D(x=self.x + other.x, y=self.y + other.y, z=self.z + other.z)
+
+    def __sub__(self, other: Point3D) -> Point3D:
+        """Subtract two points."""
+        return Point3D(x=self.x - other.x, y=self.y - other.y, z=self.z - other.z)
+
+
+@dataclass
+class Size3D(Point3D):
+    """Represent the size of a 3D object using its width, height, and depth."""
+
+    pass
+
+
+@dataclass
+class Coordinate(Point3D):
+    """Represent a coordinate in 3D space."""
+    
+    def __add__(self, other: Size3D) -> Coordinate:
+        """Add a size to a coordinate to get a second coordinate."""
+        return Coordinate(x=self.x + other.x, y=self.y + other.y, z=self.z + other.z)
+    
+    def __sub__(self, other: Size3D) -> Coordinate:
+        """Subtract a size from a coordinate to get a second coordinate."""
+        return Coordinate(x=self.x - other.x, y=self.y - other.y, z=self.z - other.z)
+    
+
+
+@dataclass
+class Centroid(Coordinate):
+    """Represent the centroid of a region in 3D space.
+
+    A centroid is simply a coordinate in 3D space that represents
+    the center of mass of a region in an image. It is represented
+    by its x, y, and z coordinates.
+
+    Attributes
+    ----------
+    x : int
+    y : int
+    z : int
+    """
+
+    pass
\ No newline at end of file
diff --git a/tests/process/images/test_crop.py b/tests/process/images/test_crop.py
new file mode 100644
index 0000000..101e081
--- /dev/null
+++ b/tests/process/images/test_crop.py
@@ -0,0 +1,177 @@
+import pytest
+import SimpleITK as sitk
+
+from readii.image_processing import loadDicomSITK, loadSegmentation
+from readii.process.images.crop import (
+    apply_bounding_box_limits,
+    check_bounding_box_single_dimension,
+    crop_image_to_mask,
+    crop_to_bounding_box,
+    crop_to_centroid,
+    crop_to_maxdim_cube,
+    find_bounding_box,
+    find_centroid,
+    validate_new_dimensions,
+)
+
+
+@pytest.fixture
+def nsclcCT():
+    return "tests/NSCLC_Radiogenomics/R01-001/09-06-1990-NA-CT_CHEST_ABD_PELVIS_WITH_CON-98785/3.000000-THORAX_1.0_B45f-95741"
+
+
+@pytest.fixture
+def nsclcSEG():
+    return "tests/NSCLC_Radiogenomics/R01-001/09-06-1990-NA-CT_CHEST_ABD_PELVIS_WITH_CON-98785/1000.000000-3D_Slicer_segmentation_result-67652/1-1.dcm"
+
+
+@pytest.fixture
+def lung4D_ct_path():
+    return "tests/4D-Lung/113_HM10395/11-26-1999-NA-p4-13296/1.000000-P4P113S303I10349 Gated 40.0B-29543"
+
+
+@pytest.fixture
+def lung4D_rt_path():
+    return "tests/4D-Lung/113_HM10395/11-26-1999-NA-p4-13296/1.000000-P4P113S303I10349 Gated 40.0B-47.35/1-1.dcm"
+
+
+@pytest.fixture
+def lung4D_image(lung4D_ct_path):
+    return loadDicomSITK(lung4D_ct_path)
+
+
+@pytest.fixture
+def lung4D_mask(lung4D_ct_path, lung4D_rt_path):
+    segDictionary = loadSegmentation(
+        lung4D_rt_path,
+        modality="RTSTRUCT",
+        baseImageDirPath=lung4D_ct_path,
+        roiNames="Tumor_c.*",
+    )
+    return segDictionary["Tumor_c40"]
+
+
+@pytest.mark.parametrize(
+    "crop_method, expected_size",
+    [
+        ("bounding_box", (50, 50, 50)),
+        ("centroid", (50, 50, 50)),
+        ("cube", (50, 50, 50)),
+        # ("pyradiomics", (22, 28, 14)),
+    ],
+)
+def test_crop_image_to_mask_methods(
+    lung4D_image,
+    lung4D_mask,
+    crop_method,
+    expected_size,
+    resize_dimensions=(50, 50, 50),
+):
+    """Test cropping image to mask with different methods"""
+    cropped_image, cropped_mask = crop_image_to_mask(
+        lung4D_image,
+        lung4D_mask,
+        crop_method,
+        resize_dimensions,
+    )
+    assert (
+        cropped_image.GetSize() == expected_size
+    ), f"Cropped image size is incorrect, expected {expected_size}, got {cropped_image.GetSize()}"
+    assert (
+        cropped_mask.GetSize() == expected_size
+    ), f"Cropped mask size is incorrect, expected {expected_size}, got {cropped_mask.GetSize()}"
+
+
+####################################################################################################
+# Parameterized tests for find_centroid and find_bounding_box
+
+
+@pytest.fixture
+def image_and_mask_with_roi(request, label_value: int = 1):
+    """
+    Fixture to create a 3D image and mask with a specified region of interest (ROI).
+
+    This fixture is used indirectly in parameterized tests via the `indirect` keyword
+    in `pytest.mark.parametrize`. The `request.param` provides the ROI coordinates, which
+    are used to define the region of interest in the mask. The ROI is specified as a
+    6-tuple (x_min, x_max, y_min, y_max, z_min, z_max).
+
+    """
+    # Create a 3D image
+    image = sitk.Image(100, 100, 100, sitk.sitkInt16)
+    mask = sitk.Image(100, 100, 100, sitk.sitkUInt8)
+
+    # Unpack ROI from the request parameter and apply it to the mask
+    roi = request.param
+    mask[roi[0] : roi[1], roi[2] : roi[3], roi[4] : roi[5]] = label_value
+
+    return image, mask
+
+
+def test_find_bounding_box_and_centroid_bad_label():
+    mask = sitk.Image(100, 100, 100, sitk.sitkUInt8)
+    mask[10:20, 10:20, 10:20] = 2
+
+    with pytest.raises(RuntimeError):
+        find_bounding_box(mask)
+
+    with pytest.raises(RuntimeError):
+        find_centroid(mask)
+
+
+# Test cases for find_bounding_box
+
+
+@pytest.mark.parametrize(
+    # First parameter passed indirectly via the fixture
+    # Second parameter is the expected bounding box
+    "image_and_mask_with_roi, expected_bbox",
+    [
+        #
+        # x_min, x_max, y_min, y_max, z_min, z_max
+        #
+        # Simple case: Perfect bounding box
+        ((85, 95, 70, 90, 60, 90), (85, 95, 70, 90, 60, 90)),
+        # Complex case: Non-standard bounding box dimensions
+        ((32, 68, 53, 77, 10, 48), (32, 68, 53, 77, 10, 48)),
+        # Single-plane ROI: ROI in only one slice
+        # since min_dim_size is 4, the ROI is expanded to 4 in if a side is too small
+        # x_max is expanded to 24
+        ((20, 21, 30, 60, 40, 80), (20, 24, 30, 60, 40, 80)),
+        # Minimum size ROI
+        # x_max is expanded to 49, y_max is expanded to 14, z_max is expanded to 9
+        ((45, 46, 10, 12, 5, 6), (45, 49, 10, 14, 5, 9)),
+    ],
+    indirect=["image_and_mask_with_roi"],  # Use the fixture indirectly
+)
+def test_find_bounding_box(image_and_mask_with_roi, expected_bbox):
+    _, mask = image_and_mask_with_roi
+    bounding_box = find_bounding_box(mask, min_dim_size=4)
+    assert (
+        bounding_box == expected_bbox
+    ), f"Bounding box is incorrect, expected {expected_bbox}, got {bounding_box}"
+
+
+# Test cases for find_centroid
+
+
+@pytest.mark.parametrize(
+    "image_and_mask_with_roi, expected_centroid",
+    [
+        # Simple case: Perfect centroid
+        ((85, 95, 70, 90, 60, 90), (90, 80, 75)),
+        # Complex case: Non-standard dimensions
+        ((32, 68, 53, 77, 10, 48), (50, 65, 29)),
+        # Single-plane ROI
+        ((20, 21, 30, 60, 40, 80), (20, 45, 60)),
+        # Minimum size ROI
+        ((45, 46, 10, 12, 5, 6), (45, 11, 5)),
+    ],
+    indirect=["image_and_mask_with_roi"],  # Use the fixture indirectly
+)
+def test_find_centroid(image_and_mask_with_roi, expected_centroid):
+    _, mask = image_and_mask_with_roi
+    centroid = find_centroid(mask)
+    assert (
+        centroid == expected_centroid
+    ), f"Centroid is incorrect, expected {expected_centroid}, got {centroid}"