Upload 7 files

Model_Class.py

@@ -0,0 +1,109 @@
+import pytorch_lightning as pl
+import torch
+import torch.nn as nn
+import utils
+from torchvision.models import resnet50
+import torch
+from monai.transforms import (
+    Compose, Resize, ResizeWithPadOrCrop, 
+)
+from pytorch_grad_cam import GradCAM
+import matplotlib.colors as mcolors
+import matplotlib.pyplot as plt
+import numpy as np
+from PIL import Image
+from io import BytesIO
+
+class ResNet(pl.LightningModule):
+    def __init__(self):
+        super().__init__()
+        self.save_hyperparameters()
+
+
+        backbone = resnet50()
+        num_input_channel = 1
+        layer = backbone.conv1
+        new_layer = nn.Conv2d(
+            in_channels=num_input_channel,
+            out_channels=layer.out_channels,
+            kernel_size=layer.kernel_size,
+            stride=layer.stride,
+            padding=layer.padding,
+            bias=layer.bias,
+        )
+        new_layer.weight = nn.Parameter(layer.weight.sum(dim=1, keepdim=True))
+
+        backbone.conv1 = new_layer
+
+
+        backbone.fc = nn.Sequential(
+            nn.Linear(2048, 1024),
+            nn.ReLU(),
+            nn.BatchNorm1d(1024),
+            nn.Dropout(0),
+            nn.Linear(1024, 2),
+        )
+
+        self.model = backbone
+
+    def forward(self, x):
+        out = self.model(x)
+        return out
+
+
+val_transforms_416x628 = Compose(
+    [
+        utils.CustomCLAHE(),  
+        Resize(spatial_size=628, mode="bilinear", align_corners=True, size_mode="longest"),
+        ResizeWithPadOrCrop(spatial_size=(416, 628)),
+    ]
+)
+
+device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+checkpoint = torch.load("classification_model.ckpt", map_location=torch.device('cpu'))
+model = ResNet().to(device)
+model.load_state_dict(checkpoint["state_dict"])
+model.eval()
+
+
+def load_and_classify_image(image_path):
+    image = val_transforms_416x628(image_path)
+    image = image.unsqueeze(0).to(device)
+
+    with torch.no_grad():
+        prediction = model(image)
+        prediction = torch.nn.functional.softmax(prediction, dim=1).squeeze(0)
+        return prediction.to('cpu'), image.to('cpu')
+    
+
+def make_GradCAM(image):
+
+    model.eval()
+    target_layers = [model.model.layer4[-1]]
+
+    arr = image.numpy().squeeze()
+    cam = GradCAM(model=model, target_layers=target_layers)
+    targets = None
+    grayscale_cam = cam(
+        input_tensor=image,
+        targets=targets,
+        aug_smooth=False,
+        eigen_smooth=True,
+    )
+    grayscale_cam = grayscale_cam.squeeze()
+
+    jet = plt.colormaps.get_cmap("inferno")
+    newcolors = jet(np.linspace(0, 1, 256))
+    newcolors[0, :3] = 0
+    new_jet = mcolors.ListedColormap(newcolors)
+  
+    plt.figure(figsize=(10, 10))
+    plt.imshow(arr, cmap='gray')
+    plt.imshow(grayscale_cam, cmap=new_jet, alpha=0.5)
+    plt.axis('off')
+    buffer2 = BytesIO()
+    plt.savefig(buffer2, format='png', bbox_inches='tight', pad_inches=0)
+    buffer2.seek(0)
+    gradcam_image = np.array(Image.open(buffer2)).squeeze()
+    
+    return gradcam_image

Model_Seg.py

@@ -0,0 +1,99 @@
+import torch
+import numpy as np
+from monai.transforms import Compose, LoadImage, EnsureChannelFirst, Lambda, Resize, NormalizeIntensity, GaussianSmooth, ScaleIntensity, AsDiscrete, KeepLargestConnectedComponent, Invert, Rotate90, SaveImage, Transform
+from monai.inferers import SlidingWindowInferer
+from monai.networks.nets import UNet
+
+class RgbaToGrayscale(Transform):
+    def __call__(self, x):
+        # squeeze last dimension, to ensure C, H, W format
+        x = x.squeeze(-1)
+        # Ensure the tensor is 3D (channels, height, width)
+        if x.ndim != 3:
+            raise ValueError(f"Input tensor must be 3D. Shape: {x.shape}")
+        
+        # Check the number of channels
+        if x.shape[0] == 4:  # Assuming RGBA
+            rgb_weights = torch.tensor([0.2989, 0.5870, 0.1140], device=x.device)
+            # Apply weights to RGB channels, output should retain one channel dimension
+            grayscale = torch.einsum('cwh,c->wh', x[:3, :, :], rgb_weights).unsqueeze(0)
+        elif x.shape[0] == 3:  # Assuming RGB
+            rgb_weights = torch.tensor([0.2989, 0.5870, 0.1140], device=x.device)
+            grayscale = torch.einsum('cwh,c->wh', x, rgb_weights).unsqueeze(0)
+        elif x.shape[0] == 1:  # Already grayscale
+            grayscale = x
+        else:
+            raise ValueError(f"Unsupported channel number: {x.shape[0]}")
+        return grayscale
+
+    def inverse(self, x):
+        # Simply return the input as the output
+        return x
+    
+model = UNet(
+    spatial_dims=2,
+    in_channels=1,
+    out_channels=4,
+    channels=[64, 128, 256, 512],
+    strides=[2, 2, 2],
+    num_res_units=3
+)
+
+device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
+
+checkpoint_path = 'segmentation_model.pt'
+checkpoint = torch.load(checkpoint_path, map_location=device)  
+assert model.state_dict().keys() == checkpoint['network'].keys(), "Model and checkpoint keys do not match"
+
+model.load_state_dict(checkpoint['network'])
+model.eval()
+
+# Define transforms for preprocessing
+pre_transforms = Compose([
+    LoadImage(image_only=True),
+    EnsureChannelFirst(),
+    RgbaToGrayscale(),  # Convert RGBA to grayscale
+    Resize(spatial_size=(768, 768)), 
+    Lambda(func=lambda x: x.squeeze(-1)),  # Adjust if the input image has an extra unwanted dimension
+    NormalizeIntensity(),
+    GaussianSmooth(sigma=0.1),
+    ScaleIntensity(minv=-1, maxv=1)
+])
+
+
+
+# Define transforms for postprocessing
+post_transforms = Compose([
+    AsDiscrete(argmax=True, to_onehot=4),
+    KeepLargestConnectedComponent(),
+    AsDiscrete(argmax=True),
+    Invert(pre_transforms),
+    #SaveImage(output_dir='./', output_postfix='seg', output_ext='.nii', resample=False)
+])
+
+
+
+def load_and_segment_image(input_image_path):
+        
+    
+    image_tensor = pre_transforms(input_image_path)
+    image_tensor = image_tensor.unsqueeze(0).to(device)
+
+    # Inference using SlidingWindowInferer
+    inferer = SlidingWindowInferer(roi_size=(512, 512), sw_batch_size=16, overlap=0.75)
+    with torch.no_grad():
+        outputs = inferer(image_tensor, model)
+
+
+    outputs = outputs.squeeze(0)
+
+    processed_outputs = post_transforms(outputs)
+
+    # rotate 
+    rotate = Rotate90(spatial_axes=(0, 1), k=3)
+    processed_outputs = rotate(processed_outputs).to('cpu')  
+
+    output_array = processed_outputs.squeeze().detach().numpy().astype(np.uint8)
+
+
+    return output_array

app.py

@@ -0,0 +1,133 @@
+import gradio as gr
+import utils
+import Model_Class 
+import Model_Seg
+
+import SimpleITK as sitk
+import torch
+from numpy import uint8
+import spaces
+image_base64 = utils.image_to_base64("anatomy_aware_pipeline.png")
+article_html = f"<img src='data:image/png;base64,{image_base64}' alt='Anatomical pipeline illustration' style='width:100%;'>"
+
+description_markdown = """
+- This tool combines a U-Net Segmentation Model with a ResNet-50 for Classification.
+- **Usage:** Just drag a pelvic x-ray into the box and hit run.
+- **Process:** The input image will be segmented and cropped to the SIJ before classification.
+- **Please Note:** This tool is intended for research purposes only.
+- **Privacy:** This tool runs completely locally, ensuring data privacy.
+"""
+
+css = """
+
+h1 {
+    text-align: center;
+    display:block;
+}
+.markdown-block {
+    background-color: #0b0f1a; /* Light gray background */
+    color: black;             /* Black text */
+    padding: 10px;            /* Padding around the text */
+    border-radius: 5px;       /* Rounded corners */
+    box-shadow: 0 0 10px rgba(11,15,26,1);
+    display: inline-flex;       /* Use inline-flex to shrink to content size */
+    flex-direction: column;
+    justify-content: center;    /* Vertically center content */
+    align-items: center;        /* Horizontally center items within */
+    margin: auto;               /* Center the block */
+}
+
+.markdown-block ul, .markdown-block ol {
+    background-color: #1e2936;
+    border-radius: 5px;
+    padding: 10px;
+    box-shadow: 0 0 10px rgba(0,0,0,0.3);
+    padding-left: 20px;         /* Adjust padding for bullet alignment */
+    text-align: left;           /* Ensure text within list is left-aligned */
+    list-style-position: inside;/* Ensures bullets/numbers are inside the content flow */
+}
+
+footer {
+    display:none !important
+}
+"""
+
+@spaces.GPU
+def predict_image(input_image, input_file):
+
+    if input_image is not None:
+        image_path = input_image
+
+    elif input_file is not None:
+        image_path = input_file
+    
+    else:
+        return None , None , "Please input an image before pressing run" , None , None
+
+    image_mask = Model_Seg.load_and_segment_image(image_path)
+
+    overlay_image_np, original_image_np = utils.overlay_mask(image_path, image_mask)
+
+    image_mask_im = sitk.GetImageFromArray(image_mask[None, :, :].astype(uint8))
+    image_im = sitk.GetImageFromArray(original_image_np[None, :, :].astype(uint8))
+    cropped_boxed_im, _ = utils.mask_and_crop(image_im, image_mask_im)
+
+    cropped_boxed_array = sitk.GetArrayFromImage(cropped_boxed_im)
+    cropped_boxed_array_disp = cropped_boxed_array.squeeze()
+    cropped_boxed_tensor = torch.Tensor(cropped_boxed_array)
+    prediction, image_transformed = Model_Class.load_and_classify_image(cropped_boxed_tensor)
+
+
+    gradcam = Model_Class.make_GradCAM(image_transformed)
+    
+    nr_axSpA_prob = float(prediction[0].item())
+    r_axSpA_prob = float(prediction[1].item())
+
+    # Decision based on the threshold
+    considered = "be considered r-axSpA" if r_axSpA_prob > 0.59 else "not be considered r-axSpA"
+
+    explanation = f"According to the pre-determined cut-off threshold of 0.59, the image should  {considered}. This Tool is for research purposes only."
+
+    pred_dict = {"nr-axSpA": nr_axSpA_prob, "r-axSpA": r_axSpA_prob}
+
+    return  overlay_image_np, pred_dict, explanation, gradcam, cropped_boxed_array_disp
+
+
+
+
+with gr.Blocks(css=css, title="Anatomy Aware axSpA") as iface:
+
+    gr.Markdown("# Anatomy-Aware Image Classification for radiographic axSpA")
+    gr.Markdown(description_markdown, elem_classes="markdown-block")
+
+    with gr.Row():
+        with gr.Column():
+
+            with gr.Tab("PNG/JPG"):
+                input_image = gr.Image(type='filepath', label="Upload an X-ray Image")
+
+            with gr.Tab("NIfTI/DICOM"):
+                input_file = gr.File(type='filepath', label="Upload an X-ray Image")
+
+            with gr.Row():
+                submit_button = gr.Button("Run", variant="primary")
+                clear_button = gr.ClearButton()
+
+        with gr.Column():
+            overlay_image_np = gr.Image(label="Segmentation Mask")
+
+            pred_dict = gr.Label(label="Prediction")
+            explanation= gr.Textbox(label="Classification Decision")
+
+            with gr.Accordion("Additional Information", open=False):
+                gradcam = gr.Image(label="GradCAM")
+                cropped_boxed_array_disp = gr.Image(label="Bounding Box")
+    
+    submit_button.click(predict_image, inputs = [input_image, input_file], outputs=[overlay_image_np, pred_dict, explanation, gradcam, cropped_boxed_array_disp])
+    clear_button.add([input_image,overlay_image_np, pred_dict, explanation, gradcam, cropped_boxed_array_disp])
+    gr.HTML(article_html)
+
+
+if __name__ == "__main__":
+    iface.queue()
+    iface.launch(server_name='0.0.0.0', server_port=8080)

requirements_small.txt

@@ -0,0 +1,12 @@
+gradio==4.29.0
+spaces==0.28.3
+numpy==1.22.2
+torch==1.13.0
+torchvision==0.14.0
+scikit-image==0.19.3
+pytorch-lightning==1.8.6
+monai-weekly==1.2.dev2320
+simpleitk==2.2.1
+nibabel==5.2.1
+itk==5.3.0
+grad-cam==1.4.6

utils.py

@@ -0,0 +1,353 @@
+from monai.transforms import Transform
+import torch
+import skimage
+import torch
+import SimpleITK as sitk
+import numpy as np
+from PIL import Image
+from io import BytesIO
+import matplotlib.pyplot as plt
+import SimpleITK as sitk
+from matplotlib.colors import ListedColormap
+import base64
+import numpy as np
+from cv2 import dilate
+from scipy.ndimage import label
+
+def image_to_base64(image_path):
+    with open(image_path, "rb") as image_file:
+        return base64.b64encode(image_file.read()).decode('utf-8')
+
+class CustomCLAHE(Transform):
+    """Implements Contrast-Limited Adaptive Histogram Equalization (CLAHE) as a custom transform, as described by Qiu et al.
+
+    Attributes:
+        p1 (float): Weighting factor, determines degree of of contour enhacement. Default is 0.6.
+        p2 (None or int): Kernel size for adaptive histogram. Default is None.
+        p3 (float): Clip limit for histogram equalization. Default is 0.01.
+
+    """
+
+    def __init__(self, p1=0.6, p2=None, p3=0.01):
+        self.p1 = p1
+        self.p2 = p2
+        self.p3 = p3
+
+    def __call__(self, data):
+        """Apply the CLAHE algorithm to input data.
+
+        Args:
+            data (Union[dict, np.ndarray]): Input data. Could be a dictionary containing the image or the image array itself.
+
+        Returns:
+            torch.Tensor: Transformed data.
+        """
+        
+        if isinstance(data, dict):
+            im = data["image"]
+
+        else:
+            im = data
+        im = im.numpy()
+        
+
+        # remove the first dimension
+        im = im[0]
+        im = im[None, :, :]
+        #im = np.expand_dims(im, axis=0)
+        im = skimage.exposure.rescale_intensity(im, in_range="image", out_range=(0, 1))
+        im_noi = skimage.filters.median(im)
+        im_fil = im_noi - self.p1 * skimage.filters.gaussian(im_noi, sigma=1)
+        im_fil = skimage.exposure.rescale_intensity(im_fil, in_range="image", out_range=(0, 1))
+        im_ce = skimage.exposure.equalize_adapthist(im_fil, kernel_size=self.p2, clip_limit=self.p3)
+        if isinstance(data, dict):
+            data["image"] = torch.Tensor(im_ce)
+        else:
+            data = torch.Tensor(im_ce)
+            
+        return data
+
+
+
+def custom_colormap():
+
+    cdict = [(0, 0, 0, 0),    # Class 0 - fully transparent (background)
+             (0, 1, 0, 0.5),  # Class 1 - Green with 50% transparency 
+             (1, 0, 0, 0.5),  # Class 2 - Red with 50% transparency 
+             (1, 1, 0, 0.5)]  # Class 3 - Yellow with 50% transparency 
+    cmap = ListedColormap(cdict)
+    return cmap
+
+def read_image(image_path):
+    try:
+        original_image = Image.open(image_path).convert('L')
+        original_image_np = np.array(original_image)
+        return original_image_np.squeeze()
+
+    except Exception as e:
+        try :
+            original_image = sitk.ReadImage(image_path)
+            original_image_np = sitk.GetArrayFromImage(original_image)
+            return original_image_np.squeeze()
+        except Exception as e:
+            print("Failed Loading the Image: ", e)
+            return None
+
+def overlay_mask(image_path, image_mask):
+    original_image_np = read_image(image_path).squeeze().astype(np.uint8)
+
+    #adjust mask intensities for display
+    image_mask_disp = image_mask
+    plt.figure(figsize=(10, 10))
+    plt.imshow(original_image_np, cmap='gray')
+
+    plt.imshow(image_mask_disp, cmap=custom_colormap(), alpha=0.5)
+    plt.axis('off')
+
+    # Save the overlay to a buffer
+    buffer = BytesIO()
+    plt.savefig(buffer, format='png', bbox_inches='tight', pad_inches=0)
+    buffer.seek(0)
+    overlay_image_np = np.array(Image.open(buffer))
+    return overlay_image_np, original_image_np
+
+
+def bounding_box_mask(image, label):
+    """Generates a bounding box mask around a labeled region in an image
+
+    Args:
+        image (SimpleITK.Image): The input image.
+        label (SimpleITK.Image): The labeled image containing the region of interest.
+
+    Returns:
+        SimpleITK.Image: An image containing the with the bounding box mask applied with the
+        same spacing as the original image.
+
+    Note:
+        This function assumes that the input image and label are SimpleITK.Image objects.
+        The returned bounding box mask is a binary image where pixels inside the bounding box
+        are set to 1 and others are set to 0.
+    """
+    # get original spacing
+    original_spacing = image.GetSpacing()
+
+    # convert image and label to arrays
+    image_array = sitk.GetArrayFromImage(image)
+    image_array = np.squeeze(image_array)
+    label_array = sitk.GetArrayFromImage(label)
+    label_array = np.squeeze(label_array)
+
+    # determine corners of the bounding box
+    first_nonzero_row_index = np.nonzero(np.any(label_array != 0, axis=1))[0][0]
+    last_nonzero_row_index = np.max(np.nonzero(np.any(label_array != 0, axis=1)))
+    first_nonzero_column_index = np.nonzero(np.any(label_array != 0, axis=0))[0][0]
+    last_nonzero_column_index = np.max(np.nonzero(np.any(label_array != 0, axis=0)))
+
+    top_left_corner = (first_nonzero_row_index, first_nonzero_column_index)
+    bottom_right_corner = (last_nonzero_row_index, last_nonzero_column_index)
+
+    # define the bounding box as an array mask
+    bounding_box_array = label_array.copy()
+    bounding_box_array[
+        top_left_corner[0] : bottom_right_corner[0] + 1,
+        top_left_corner[1] : bottom_right_corner[1] + 1,
+    ] = 1
+    
+    # add channel dimension
+    bounding_box_array = bounding_box_array[None, ...].astype(np.uint8)
+
+    # get Image from Array Mask and apply original spacing
+    bounding_box_image = sitk.GetImageFromArray(bounding_box_array)
+    bounding_box_image.SetSpacing(original_spacing)
+    return bounding_box_image
+
+
+def threshold_based_crop(image):
+    """
+    Use Otsu's threshold estimator to separate background and foreground. In medical imaging the background is
+    usually air. Then crop the image using the foreground's axis aligned bounding box.
+    Args:
+        image (SimpleITK image): An image where the anatomy and background intensities form a
+                                 bi-modal distribution
+                                 (the assumption underlying Otsu's method.)
+    Return:
+        Cropped image based on foreground's axis aligned bounding box.
+    """
+
+    inside_value = 0
+    outside_value = 255
+    label_shape_filter = sitk.LabelShapeStatisticsImageFilter()
+    # uncomment for debugging
+    #sitk.WriteImage(image, "./image.png")
+    label_shape_filter.Execute(sitk.OtsuThreshold(image, inside_value, outside_value))
+    bounding_box = label_shape_filter.GetBoundingBox(outside_value)
+    return sitk.RegionOfInterest(
+        image,
+        bounding_box[int(len(bounding_box) / 2) :],
+        bounding_box[0 : int(len(bounding_box) / 2)],
+    )
+
+def creat_SIJ_mask(image, input_label):
+    """
+    Create a mask for the sacroiliac joints (SIJ) from pelvis and sascrum segmentation mask
+
+    Args:
+        image (SimpleITK.Image): x-ray image.
+        input_label (SimpleITK.Image): Segmentation mask containing labels for sacrum, left- and right pelvis
+
+    Returns:
+        SimpleITK.Image: Mask of the SIJ
+
+    """
+    
+    original_spacing = image.GetSpacing()
+    # uncomment for debugging
+    #sitk.WriteImage(input_label, "./input_label.png")
+    mask_array = sitk.GetArrayFromImage(input_label).squeeze()
+    
+    sacrum_value = 1  
+    left_pelvis_value = 2  
+    right_pelvis_value = 3  
+    background_value = 0
+
+    
+    sacrum_mask = (mask_array == sacrum_value)
+
+    first_nonzero_column_index = np.nonzero(np.any(sacrum_mask != 0, axis=0))[0][0]
+    last_nonzero_column_index = np.max(np.nonzero(np.any(sacrum_mask != 0, axis=0)))
+    box_width=last_nonzero_column_index-first_nonzero_column_index
+
+    dilation_extent = int(np.round(0.05 * box_width))
+
+    dilated_sacrum_mask = dilate_mask(sacrum_mask, dilation_extent)
+
+    intersection_left = (dilated_sacrum_mask & (mask_array == left_pelvis_value))
+    if np.all(intersection_left == 0):
+        print("Warning: No left intersection")
+        left_pelvis_mask = (mask_array == 2)
+        intersection_left = create_median_height_array(left_pelvis_mask)
+        
+    intersection_left = keep_largest_component(intersection_left)
+    
+    intersection_right = (dilated_sacrum_mask & (mask_array == right_pelvis_value))
+    if np.all(intersection_right == 0):
+        print("Warning: No right intersection")
+        right_pelvis_mask = (mask_array == 3)
+        intersection_right = create_median_height_array(right_pelvis_mask)
+    intersection_right = keep_largest_component(intersection_right)
+    
+    intersection_mask = intersection_left +intersection_right
+    intersection_mask = intersection_mask[None, ...]
+                                    
+    instersection_mask_im = sitk.GetImageFromArray(intersection_mask)
+    instersection_mask_im.SetSpacing(original_spacing)
+    return instersection_mask_im
+
+def dilate_mask(mask, extent):
+    """
+    Keeps only the largest connected component in a binary segmentation mask.
+
+    Args:
+        mask (numpy.ndarray): A numpy array representing the binary segmentation mask, 
+                              with 1s indicating the label and 0s indicating the background.
+
+    Returns:
+        numpy.ndarray: A modified version of the input mask, where only the largest 
+                       connected component is retained, and other components are set to 0.
+
+    """
+    mask_uint8 = mask.astype(np.uint8)
+
+    kernel = np.ones((2*extent+1, 2*extent+1), np.uint8)
+    dilated_mask = dilate(mask_uint8, kernel, iterations=1)
+    return dilated_mask
+
+def mask_and_crop(image, input_label):
+    """
+    Performs masking and cropping operations on an image and its label.
+
+    Args:
+        image (SimpleITK.Image): The image to be processed.
+        label (SimpleITK.Image): The corresponding label image.
+
+    Returns:
+        tuple: A tuple containing two SimpleITK.Image objects.
+            - cropped_boxed_image: The image after applying bounding box masking and cropping.
+            - mask: The binary mask corresponding to the label after cropping.
+
+    Note:
+        This function relies on other functions: bounding_box_mask() and threshold_based_crop().
+    """
+    input_label = creat_SIJ_mask(image,input_label)
+    box_mask = bounding_box_mask(image, input_label)
+    
+    boxed_image = sitk.Mask(image, box_mask, maskingValue=0, outsideValue=0)
+    masked_image = sitk.Mask(image, input_label, maskingValue=0, outsideValue=0)
+
+    cropped_boxed_image = threshold_based_crop(boxed_image)
+    cropped_masked_image = threshold_based_crop(masked_image)
+
+    mask = np.squeeze(sitk.GetArrayFromImage(cropped_masked_image))
+    mask = np.where(mask > 0, 1, 0)
+    mask = sitk.GetImageFromArray(mask[None, ...])
+    return cropped_boxed_image, mask
+
+def create_median_height_array(mask):
+    """
+    Creates an array based on the median height of non-zero elements in each column of the input mask.
+
+    Args:
+        mask (numpy.ndarray): A binary mask with 1s representing the label and 0s the background.
+
+    Returns:
+        numpy.ndarray: A new binary mask array with columns filled based on the median height,
+                       or None if the input mask has no non-zero columns.
+                       
+    Note: 
+        This function is only used when there is no intersection between pelvis and sacrum, and creates an alternative
+        SIJ mask, that serves as an approximate replacement.
+    """
+    rows, cols = mask.shape
+    column_details = []
+
+    for col in range(cols):
+        column_data = mask[:, col]
+        non_zero_indices = np.nonzero(column_data)[0]
+        if non_zero_indices.size > 0:
+            height = non_zero_indices[-1] - non_zero_indices[0] + 1
+            start_idx = non_zero_indices[0]
+            column_details.append((height, start_idx, col))
+            
+    if not column_details:
+        return None  
+    median_height = round(np.median([h[0] for h in column_details]))
+    median_cols = [(col, start_idx) for height, start_idx, col in column_details if height == median_height]
+    new_array = np.zeros_like(mask, dtype=int)
+    for col, start_idx in median_cols:
+        start_col = max(0, col - 5)
+        end_col = min(cols, col + 5)
+        new_array[start_idx:start_idx + median_height, start_col:end_col] = 1
+    return new_array
+
+def keep_largest_component(mask):
+    """
+    Identifies and retains the largest connected component in a binary segmentation mask.
+
+    Args:
+        mask (numpy.ndarray): A binary mask with 1s representing the label and 0s the background.
+
+    Returns:
+        numpy.ndarray: The modified mask with only the largest connected component.
+    """
+    # Label the connected components
+    labeled_array, num_features = label(mask)
+
+    # If no features are found, return the original mask
+    if num_features <= 1:
+        return mask
+
+    # Find the largest connected component
+    largest_component = np.argmax(np.bincount(labeled_array.flat)[1:]) + 1
+
+    # Generate the mask for the largest component
+    return (labeled_array == largest_component).astype(mask.dtype)