Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Rolling ball background estimation using a binary 3D volume
function [bg,image,binary3D] = rollingBallBackgroundViaBinary3D(image,ballRadius,smooth)
%rollingBallBackgroundViaBinary3D Calculate rolling ball background
%
% INPUT
% image - 2D image of positive integers
% ballRadius - (optional) radius of rolling ball, default: 5
% smooth - (optional) radius for Gaussian smoothing
%
% OUTPUT
% bg - rolling ball estimated background
% image - smoothed and rounded image
% binary3D - binary 3D volume
%
% If no output is requested a montage is shown from left to right of
% 1. smoothed and rounded image
% 2. rolling ball background
% 3. rolling ball background subtracted from smoothed nad rounded image
%
% References
% 1. SR Sternberg. Biomedical Image Processing. IEEE Computer 1983.
% https://doi.ieeecomputersociety.org/10.1109/MC.1983.1654163
% 2. SR Sternberg. Grayscale Morphology. CVGIP 1986
% https://doi.org/10.1016/0734-189X(86)90004-6
% 3. R Adams. “Radial Decomposition of Discs and Spheres”.
% CVGIP: Graphical Models and Image Processing,
% vol. 55, no. 5, September 1993, pp. 325-332.
% https://doi.org/10.1006/cgip.1993.1024
%
% AUTHORS
% Mark Kittisopikul - Dec 2019
if nargin < 2
% Default ball radius
ballRadius = 5;
end
if nargin > 2 && smooth > 0
% Do Gaussian smoothing if specified
image = imgaussfilt(image,smooth);
end
% Image must consist of positive integers
image = round(image);
% clijx.maximumOfAllPixels
binary3DSize = [size(image) double(max(image(:)))];
% Get linear indices for the top of the stack from xyz coordinates
[xc,yc] = meshgrid(1:size(image,2),1:size(image,1));
ind = sub2ind(binary3DSize,yc,xc,image);
% Create a binary 3D volume where the intensity values are stacks of 1s
binary3D = false(binary3DSize);
binary3D(ind) = true;
% clijx.binaryAnd
binary3D = cumsum(binary3D,3,'reverse');
% The structure element is then just a sphere
sphereStructureElement = strel('sphere',ballRadius);
% Pad array in all dimensions using the ball radius with 1s
binary3D = padarray(binary3D,ones(1,3)*ballRadius,true);
% Do erosion and then dilation using the sphere
% clijx.erodeSphere
% clijx.dilateSphere
binary3D = imopen(binary3D,sphereStructureElement);
% Crop the padding off
firstElement = ballRadius+1;
% clijx.crop 3D???
binary3D = binary3D(firstElement:end-ballRadius, ...
firstElement:end-ballRadius, ...
firstElement:end-ballRadius);
% The result is the remaining number of binary voxels per xy position
% clijx.sumZProjection
bg = sum(binary3D,3);
% If no output is requested, show a montage of the image, bg, and tophat
if nargout == 0
figure;
imshow([image bg double(image)-bg],[]);
end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.