How to remove CT bed/shadows in a CT image with python?

Question

I am working with 3D CT images and trying to remove the lines from the bed.

A slice from the original Image:

Following is my code to generate the mask:

segmentation = morphology.dilation(image_norm, np.ones((1, 1, 1)))
labels, label_nb = ndimage.label(segmentation)    
label_count = np.bincount(labels.ravel().astype(int))
label_count[0] = 0

mask = labels == label_count.argmax()
 
mask = morphology.dilation(mask, np.ones((40, 40, 40)))
mask = ndimage.morphology.binary_fill_holes(mask) 
mask = morphology.dilation(mask, np.ones((1, 1, 1)))

This results in the following image:

As you can see, in the above image the CT scan as distorted as well.

If I change: mask = morphology.dilation(mask, np.ones((40, 40, 40))) to mask = morphology.dilation(mask, np.ones((100, 100, 100))), the resulting image is as follows:

How can I remove only the two lines under the image without changing the image area? Any help is appreciated.

Answer 1

You've probably found another solution by now. Regardless, I've seen similar CT processing questions on SO, and figured it would be helpful to demonstrate a Scikit-Image solution. Here's the end result.

Here's the code to produce the above images.

from skimage import io, filters, color, morphology
import matplotlib.pyplot as plt
import numpy as np 

image = color.rgba2rgb(
    io.imread("ctimage.png")[9:-23,32:-9]
)
gray = color.rgb2gray(image)
tgray = gray > filters.threshold_otsu(gray)
keep_mask = morphology.remove_small_objects(tgray,min_size=463)
keep_mask = morphology.remove_small_holes(keep_mask)
maskedimg = np.einsum('ijk,ij->ijk',image,keep_mask)

fig,axes = plt.subplots(ncols=3)
image_list = [image,keep_mask,maskedimg]
title_list = ["Original","Mask","Imgage w/mask"]
for i,ax in enumerate(axes):
    ax.imshow(image_list[i])
    ax.set_title(title_list[i])
    ax.axis("off")
fig.tight_layout()

Notes on code

image = color.rgba2rgb(
    io.imread("ctimage.png")[9:-23,32:-9]
)
gray = color.rgb2gray(image)

The image saved as RGBA when I loaded it from SO. It needs to be in grayscale for use in the threshold function. Your image might already by in grayscale. Also, the downloaded image showed axis markings. That's why I've trimmed the image.

maskedimg = np.einsum('ijk,ij->ijk',image,keep_mask)

I wanted to apply keep_mask to every channel of the RGB image. The mask is a 2D array, and the image is a 3D array. I referenced this previous question in order to apply the mask to the image.