Rolling ball background estimation using a binary 3D volume

rollingBallBackgroundViaBinary3D.m

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