# Surface Closing in Breast Segmentation from MRI Images
**URL:** https://discourse.itk.org/t/surface-closing-in-breast-segmentation-from-mri-images/6403
**Category:** Algorithms
**Created:** [January 21, 2024, 7:28am UTC](https://discourse.itk.org/t/surface-closing-in-breast-segmentation-from-mri-images/6403 "2024-01-21T07:28:02Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@Miri\_Trope](https://discourse.itk.org/u/Miri_Trope)
#### Post date: [January 21, 2024, 7:28am UTC](https://discourse.itk.org/t/surface-closing-in-breast-segmentation-from-mri-images/6403/1 "2024-01-21T07:28:02Z")
Hello everyone, I’m currently facing an issue in obtaining a closed surface for the segmentation output of the attached MRI breast image. The algorithm’s objective is to segment the breast (body). As you can see, the surface lacks smoothness due to the points not being connected. I attempted using a smoothing algorithm with high sigmas, but there was minimal improvement. Here’s the code responsible for segmenting the breast body.
```auto
import SimpleITK as sitk
# Load the MRI image
image_path = '16_breast_image.nrrd'
image = sitk.ReadImage(image_path)
# Create an empty segmentation image with the same size and metadata
segmentation_image = sitk.Image(image.GetSize(), sitk.sitkUInt8)
segmentation_image.CopyInformation(image)
# Define a threshold to distinguish between object and background
threshold = 40 # You can adjust this threshold based on your image characteristics
# Iterate through the image to find the first non-zero voxel
for k in range(segmentation_image.GetDepth()):
for j in range(segmentation_image.GetHeight()):
for i in range(segmentation_image.GetWidth()):
intensity = image.GetPixel(i, j, k)
# Check if the voxel corresponds to an object (not background)
if intensity > threshold:
# Set the corresponding voxel in the segmentation image to 255 (white)
segmentation_image.SetPixel(i, j, k, 1)
break
# Iterate through the image to find the first non-zero voxel in the reverse direction
for k in range(segmentation_image.GetDepth() - 1, -1, -1):
for j in range(segmentation_image.GetHeight() - 1, -1, -1):
for i in range(segmentation_image.GetWidth() - 1, -1, -1):
intensity = image.GetPixel(i, j, k)
# Check if the voxel corresponds to an object (not background)
if intensity > threshold:
# Set the corresponding voxel in the reverse segmentation image to 255 (white)
segmentation_image.SetPixel(i, j, k, 1)
break
```
Any comments you may have would be highly appreciated.
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [January 21, 2024, 4:29pm UTC](https://discourse.itk.org/t/surface-closing-in-breast-segmentation-from-mri-images/6403/2 "2024-01-21T16:29:54Z")
Hello @Miri_Trope,
Looking at your code, you are performing simple thresholding. In Python you should pretty much avoid iterating whenever possible. In this case, it’s easy to threshold the whole image, then smooth the result using the [AntiAliasBinaryImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1AntiAliasBinaryImageFilter.html#details) followed by extraction of the zero level set using the [ZeroCrossingImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1ZeroCrossingImageFilter.html).
You may want to consider using an automatically detected threshold, if possible (e.g. the [OtsuThresholdImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1OtsuThresholdImageFilter.html), other methods are available, see [this notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/90a70c1840474599c6362dabd594bbe2bc744913/Python/35_Segmentation_Shape_Analysis.ipynb#L107), search for `threshold_filters`). See code below:
```auto
import SimpleITK as sitk
image_path = '16_breast_image.nrrd'
image = sitk.ReadImage(image_path)
# segment using manually set threshold
threshold = 40
segmentation = image > threshold
#segment using Otsu threshold, assumes bi-modal distribution,
#foreground-background content
#foreground_value=1
#background_value=0
#segmentation = sitk.OtsuThreshold(image, background_value, foreground_value)
# smooth the image and then extract the zero level set
level_set_image = sitk.AntiAliasBinary(segmentation, maximumRMSError = 0.07, numberOfIterations = 1000)
smoothed_segmentation = sitk.ZeroCrossing(level_set_image)
```
---
### Author:  [@Miri\_Trope](https://discourse.itk.org/u/Miri_Trope)
#### Post date: [January 21, 2024, 7:04pm UTC](https://discourse.itk.org/t/surface-closing-in-breast-segmentation-from-mri-images/6403/3 "2024-01-21T19:04:40Z")
> [@zivy](#):
>
> ```auto
> segmentation = image > threshold
>
> #segment using Otsu threshold, assumes bi-modal distribution,
> #foreground-background content
> #foreground_value=1
> #background_value=0
> #segmentation = sitk.OtsuThreshold(image, background_value, foreground_value)
>
> # smooth the image and then extract the zero level set
> level_set_image = sitk.AntiAliasBinary(segmentation, maximumRMSError = 0.07, numberOfIterations = 1000)
> smoothed_segmentation = sitk.ZeroCrossing(level_set_image)
>
> ```
Dear Ziv,
I appreciate your feedback. The code I’ve written incorporates an algorithm that identifies and marks the initial pixel encountered within the image using raster order. This approach ensures the segmentation of the outer body itself. Please note that this method differs from a conventional threshold filter.
In any case, I’ve included the result image generated by the algorithm. Despite attempting smoothing and closing operations, I have been unsuccessful in achieving a fully closed surface due to numerous unconnected points that are missing.
Thank you,
Miri
[breast\_body\_segmentation\_output.nii.gz](https://discourse.itk.org/uploads/short-url/vRdDWgbw5KVLfT7mkI7jk86Y1Mc.gz) (155.5 KB)
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [January 21, 2024, 11:50pm UTC](https://discourse.itk.org/t/surface-closing-in-breast-segmentation-from-mri-images/6403/4 "2024-01-21T23:50:42Z")
Hello @Miri_Trope,
What you want to do is not trivial (i.e. not a single filter). Possibly try Slicer’s [surface wrap solidify](https://github.com/sebastianandress/Slicer-SurfaceWrapSolidify/tree/master) functionality. Try it via Slicer’s GUI and if this solves your issue, look at the code to get the details.