3D Morphological operation (imclose) taking too long time to run

I need to run imclose with a "sphere" shape or radius 50, on a 3D binary mask of size ~1K in each dimension. The function takes very long time to run. Is there a better way to do it?

% X = original image of size 600x1800x800

% Threshold image - global threshold
BW = imbinarize(X);

% Close mask with disk
radius = 50;
se = strel('sphere', radius);
BW = imclose(BW, se);

[Matt J]

The slow speed is certainly because of the large strel radius you're using, but recommendations for an alternative won't be possible withous seeing an image.

[Arash Mirhashemi]

Here is one slice of the 3D cube (CT scan of bone and cartilage), showing the binary mask overlayed on the original image. I want the holes to be filled, but ALSO, the dents on the borders to be filled. That is, if there is a half-hole on the border of the mask, I want it to be included as bone (e.g., some can be seen on bottom-left cut on the bone).

I am wondering if there is a cheat to run "imclose" faster. For example, as an approximate answer, to run it on a resized mask with down-sampled size, then resizing it back to the original size.

[Walter Roberson]

I wonder if you could take the fourier transform of the sphere and the fact that a number of morphological operations can be expressed as convolutions https://www.l3harrisgeospatial.com/docs/convolutionmorphologyfilters.html and that convolution can be expressed as operations on fourier transforms?

The disadvantage is that it might require calculations of bessel functions over space equivalent to the number of pixels.

Answer 1

[Sai Pavan]

Hello Arash,

I understand that you are looking for ways to reduce the amount of time taken to perform morphological operation 'imclose' on 3D binary masks.

Performing morphological operations like `imclose` on large 3D binary masks can indeed be computationally intensive, especially with a large structuring element like a sphere with a radius of 50. Please refer to a few strategies mentioned below to speed up the process, which may involve some trade-offs in terms of accuracy:

• Downsampling: As you suggested, one approach is to downsample the image, perform the morphological operation, and then upsample the result. This will reduce the computational load significantly but will result in a loss of precision.

% Downsample the image
scale_factor = 0.1;
BW_small = imresize3(BW, scale_factor, 'Method', 'nearest');
% Perform the closing operation on the smaller image
radius_small = round(radius * scale_factor);
se_small = strel('sphere', radius_small);
BW_small_closed = imclose(BW_small, se_small);
% Upsample the image back to the original size
BW_closed = imresize3(BW_small_closed, size(BW), 'Method', 'nearest');

Please refer to the below documentation to learn more about 'imresize3' function: https://www.mathworks.com/help/images/ref/imresize3.html

• Approximate Structuring Element: Instead of using a 'sphere', which is computationally expensive, you can approximate it with a 'cube' or 'ball' which may be faster. This will change the morphology of the result but can be much faster.

% Use an approximate structuring element
se_approx = strel('cube', 2 * radius + 1);
BW = imclose(BW, se_approx);

• Parallel Processing: If you have access to the Parallel Processing Toolbox, you can try to distribute the computation across multiple cores.

% Enable parallel processing
parpool; % or specify the number of workers if needed
% Perform the closing operation
BW = imclose(BW, se);

Please refer to the below documentation to learn more about 'parpool' function: https://www.mathworks.com/help/parallel-computing/parpool.html

• Block Processing: Divide the image into smaller blocks, perform `imclose` on each block, and then stitch the blocks back together. This is complex and may introduce artifacts at the block boundaries but can be effective for very large images.
• GPU Computing: If you have a compatible GPU and the Parallel Computing Toolbox, you can try to perform the operation on the GPU.

% Move the data to the GPU
BW_gpu = gpuArray(BW);
% Perform the closing operation on the GPU
BW_gpu = imclose(BW_gpu, se);
% Retrieve the result from the GPU
BW = gather(BW_gpu);

Please refer to the below documentation to learn more about:

1. 'gpuArray' function: https://www.mathworks.com/help/parallel-computing/gpuarray.html
2. 'gather' function: https://www.mathworks.com/help/parallel-computing/gpuarray.gather.html

Hope it helps!