
function [seg, Am, mask] = mcleanupregions(seg, seRadius)
option = 2;
    % 1) Ensure every segment is distinct 
    [seg, maxlabel] = makeregionsdistinct(seg); 
    
    % 2) Perform a morphological opening on each segment, subtract the opening
    % from the orignal segment to obtain regions to be reassigned to
    % neighbouring segments.
    if seRadius
        se = circularstruct(seRadius);   % Accurate and not noticeably slower
                                         % if radius is small
%       se = strel('disk', seRadius, 4);  % Use approximated disk for speed
        mask = zeros(size(seg));

        if option == 1        
            for l = 1:maxlabel
                b = seg == l;
                mask = mask | (b - imopen(b,se));
            end
            
        else   % Rather than perform a morphological opening on every
               % individual region in sequence the following finds separate
               % lists of unconnected regions and performs openings on these.
               % Typically an image can be covered with only 5 or 6 lists of
               % unconnected regions.  Seems to be about 2X speed of option
               % 1. (I was hoping for more...)
            list = finddisconnected(seg);
            
            for n = 1:length(list)
                b = zeros(size(seg));
                for m = 1:length(list{n})
                    b = b | seg == list{n}(m);
                end

                mask = mask | (b - imopen(b,se));
            end
        end
        
        % Compute distance map on inverse of mask
        [~, idx] = bwdist(~mask);
        
        % Assign a label to every pixel in the masked area using the label of
        % the closest pixel not in the mask as computed by bwdist
        seg(mask) = seg(idx(mask));
    end
    
    % 3) As some regions will have been relabled, possibly broken into several
    % parts, or absorbed into others and no longer exist we ensure all regions
    % are distinct again, and renumber the regions so that they sequentially
    % increase from 1.  We also need to reconstruct the adjacency matrix to
    % reflect the changed number of regions and their relabeling.

    seg = makeregionsdistinct(seg);
    [seg, minLabel, maxLabel] = renumberregions(seg);
    Am = regionadjacency(seg);    
    