Distance map process is very memory expensive

Engineering
1

Hi,
I am using SimpleITK to compute distance map of a volume of 3Go on my machine of 16Go of RAM. But i am getting crash saying “Failed to allocate memory for image”. I tested the maurer and daniels functions. Is there something to do about this?
Thanks

2

Hello @Limz,

This is a “feature” not a bug. For a binary segmentation, the pixel type is 8-bit unsigned integer, for the distance map the pixel type is 32-bit float, so four times larger. Thus for a 3GB segmentation the corresponding distance map takes 12GB. Just these two images are already stressing your machine’s memory.

As images are only an approximation of a continuous world, you can possibly work with a coarser approximation.

Resample the segmentation to a coarser grid and work with that. Note that in ITK/SimpleITK these images will still occupy the same physical space, they just have fewer samples. How aggressively you resample depends on the accuracy required by the application. Below is code that illustrates the concepts so that you can play with it and decide if using a coarser approximation maintains the desired accuracy while addressing the memory constraints.

import SimpleITK as sitk

segmentation = sitk.Image([512, 512], sitk.sitkUInt8)
segmentation[80:120, 80:200] = 1
segmentation[360:440, 280:400] = 1

distance_map = sitk.SignedMaurerDistanceMap(
    segmentation, squaredDistance=False, useImageSpacing=True
)

# Code below downsamples the segmentation by a factor of 2 in all dimensions
# Reduces memory to 0.25 of original size in 2D and 0.125 in 3D
# Downsampling by different factors along each of the dimensions is possible.
# In microscopy it is not uncommon to only downsample in the x-y plane and keep
# the Z dimension as is (original Z spacing is very large compared to X and Y spacing):
# new_size = [segmentation.GetWidth()//2, segmentation.GetHeight()//2, segmentation.GetDepth()]


new_size = [sz // 2 for sz in segmentation.GetSize()]
new_spacing = [
    ((original_sz - 1) * original_spc) / (new_sz - 1)
    for original_sz, original_spc, new_sz in zip(
        segmentation.GetSize(), segmentation.GetSpacing(), new_size
    )
]
downsampled_segmentation = sitk.Resample(
    segmentation,
    new_size,
    sitk.Transform(),
    sitk.sitkNearestNeighbor,
    segmentation.GetOrigin(),
    new_spacing,
    segmentation.GetDirection(),
)

downsampled_distance_map = sitk.SignedMaurerDistanceMap(
    downsampled_segmentation, squaredDistance=False, useImageSpacing=True
)

# Difference between original and upsampled version of coarser distance map, using the
# default linear interpolator in the Resample function. Computationally more complex interpolators
# are available and may provide more accurate results (required accuracy is application dependent)
distance_map_diff = distance_map - sitk.Resample(downsampled_distance_map, segmentation)

sitk.Show(distance_map_diff)
3

Thanks @zivy , yes indeed , that’s due to the data type from the distance map. Okay then, I will try for the resampling.