Code covered by the BSD License  

Highlights from
Fast morphological reconstruction of large logical masks.

image thumbnail

Fast morphological reconstruction of large logical masks.

by

 

MATLAB's imreconstruct is slow for large 3D logical masks. bwreconstruct is a faster replacement.

bwreconstruct(marker, mask, conn)
function maskOut = bwreconstruct(marker, mask, conn)
%BWRECONSTRUCT Fast morphological reconstruction of large logical masks.
%   IM = BWRECONSTRUCT(MARKER,MASK) performs morphological reconstruction
%   of the logical MARKER under the logical MASK.  MARKER and MASK are two
%   two binary images with the same size, resulting in a binary image IM of
%   equivalent size.
%
%   IM = BWRECONSTRUCT(MARKER,MASKCC) performs similarly, but uses the
%   connected component structure MASKCC. MASKCC is the output from
%   BWCONNCOMP(MASK), and when it is already available, using MASKCC is
%   faster than providing a MASK binary image.
%
%   IM = IMRECONSTRUCT(MARKER,MASK,CONN) performs morphological
%   reconstruction with the specified connectivity.  CONN may have the
%   following scalar values:
%
%       4     two-dimensional four-connected neighborhood
%       8     two-dimensional eight-connected neighborhood
%       6     three-dimensional six-connected neighborhood
%       18    three-dimensional 18-connected neighborhood
%       26    three-dimensional 26-connected neighborhood
%
%   Connectivity may be defined in a more general way (see IMRECONSTRUCT
%   for details).
%
%   Class support
%   -------------
%   MARKER and MASK must be logical arrays with the same dimension. If
%   MASKCC is used instead of MASK, it must relate to this same dimension.
%
%   Performance Note
%   ----------------
%   This function performs faster than IMRECONSTRUCT for large binary 
%   images of 3 or higher dimensions. For 2D or small images (sizes below
%   50-by-50-by-50, for example), IMRECONSTRUCT is still the faster option.
%
%   Example
%   ---------
%       mask = false(100,100,100);
%       mask(:,[3 50 75],:) = true;
%       marker = false(size(mask));
%       marker(45:55,45:55,50:55) = true;
%       tic, im1 = imreconstruct(marker,mask); toc
%       tic, im2 = bwreconstruct(marker,mask); toc
%       isequal(im1, im2)
%       figure, patch(isosurface(im2,0.5),'FaceColor','g','EdgeColor','none'), hold on
%       patch(isosurface(marker,0.5),'FaceColor','b','EdgeColor','none')
%       patch(isosurface(mask & ~im2,0.5),'FaceColor','r','EdgeColor','none')
%       axis image, view(3), camlight
%
%   See also IMRECONSTRUCT, BWSELECT.

% Sven Holcombe, 2013-04-20.

% Detect input as either BW mask, or BWCONNCOMP structure.
if islogical(mask)
    if nargin==3
        connCell = {conn};
    else
        connCell = {};
    end
    maskCC = bwconncomp(mask,connCell{:});
elseif isstruct(mask) && all(isfield(mask, {'ImageSize','PixelIdxList'}))
    maskCC = mask;
else
    error('bwreconstruct:badMask','Mask must either be a logical matrix or the result of bwconncomp()')
end

% Find all pixel indices of the marker
markerInds = find(marker);

% Ask which connected components intersec the marker
% Equivalent (but slower) code:
% CCkeepInds = cellfun(@(x)any(ismember(x,markerInds)), maskCC.PixelIdxList);
pxNums = cell2mat(arrayfun(@(i)ones(length(maskCC.PixelIdxList{i}),1)*i, (1:length(maskCC.PixelIdxList))', 'Un',0));
CCkeepInds = unique(pxNums(ismember(cat(1,maskCC.PixelIdxList{:}),markerInds)));

% Build result from the connected components which intersected the marker
maskOut = false(maskCC.ImageSize);
maskOut(cat(1,maskCC.PixelIdxList{CCkeepInds})) = true;

Contact us