Sorry for the late reply to your question. When I went to examine the code, I realized that I couldn’t easily tell what it was doing, either. I eventually unraveled it far enough to figure out that it’s doing two things, filling holes and updating for differences between the preliminary and edited primary objects.
The filling holes code (Holey* images and such) fills any background region entirely enclosed by the same secondary object. Holes that tough two secondary objects or the edge are not filled.
The preliminary/edited code is there because when we identify primary objects (often nuclei), they are often filtered out by size or if they touch the border of the image, and removes from the measurements downstream. However, when identifying secondary objects, we still want to use the possibly-removed nuclei for identifying boundaries. So the other code here is responsible for finding out which primary objects have been removed due to size and border filtering, and removing secondary objects that correspond.
I have rewritten this section of code, and taken the hole filling code out into a separate subfunction. The new code looks like:
%%% STEP 3: We used the PrelimPrimaryLabelMatrixImage as the
%%% source for primary objects, but that label-matrix is built
%%% before small/large objects and objects touching the
%%% boundary are removed. We need to filter the label matrix
%%% from propagate to make the labels match, and remove any
%%% secondary objects that correspnd to size- or
%%% boundary-filtered primaries.
%%% Map preliminary labels to edited labels based on maximum
%%% overlap from prelim to edited. We can probably assume
%%% that no borders are adjusted during editing (i.e., any
%%% changes from Prelim to Edited only involves removing
%%% entire objects), but this is safer.
%%% (add one so that zeros are remapped correctly.)
PrelimToEditedHist = sparse(EditedPrimaryLabelMatrixImage(:) + 1, PrelimPrimaryLabelMatrixImage(:) + 1, 1);
[ignore, PrelimToEditedRemap] = sort(PrelimToEditedHist, 1);
PrelimToEditedRemap = PrelimToEditedRemap(end, :) - 1;
%%% make sure zeros map to zeros (note the off-by-one for the
%%% index because Matlab doesn't do 0-indexing).
PrelimToEditedRemap(1) = 0;
EditedLabelMatrixImage = PrelimToEditedRemap(PropagatedImage + 1);
%%% STEP 4:
%%% Fill holes (any contiguous, all-0 regions that are
%%% surrounded by a single value).
FinalLabelMatrixImage = CPfill_holes(EditedLabelMatrixImage);
Hopefully this is more clear.
Here is the CPfill_holes code:
function FilledLabelMatrix = CPfill_holes(LabelMatrix)
%%% function FilledLabelMatrix = CPfill_holes(LabelMatrix)
%%% Fill holes (0-values) in a label matrix, as defined by a
%%% 4-connected region. Holes are filled only if surrounded by a
%%% single value.
MaximaImage = ordfilt2(LabelMatrix, 5, [0 1 0; 1 1 1 ; 0 1 0]);
ZerosMaxima = MaximaImage .* (LabelMatrix == 0);
%%% Likewise, for every zero pixel, find its smallest adjacent
%%% label. A little trickier, since we need to ignore zeros.
%%% replace 0s with a unique label
UniqueLabel = max(LabelMatrix(:)) + 1;
NoZeros = LabelMatrix;
NoZeros(NoZeros == 0) = UniqueLabel;
MinimaImage = ordfilt2(NoZeros, 1, [0 1 0; 1 1 1 ; 0 1 0]);
%%% replace the unique label with zeros
MinimaImage(MinimaImage == UniqueLabel) = 0;
ZerosMinima = MinimaImage .* (LabelMatrix == 0);
%%% Find the zero regions, removing any that touch the border.
ZeroRegions = CPclearborder(bwlabel(LabelMatrix == 0, 4));
%%% Boundaries of the zero regions are those with a nonzero MaximaImage
ZeroBoundaries = ((ZerosMaxima ~= 0) & (ZeroRegions ~= 0));
%%% Now, build a map from zero region labels to object labels, based
%%% on ZerosMaxima
ZeroLocations = find(ZeroBoundaries);
%%% Check for no holes and return
if length(ZeroLocations) == 0,
FilledLabelMatrix = LabelMatrix;
LocationsZerosAndMaxima = sparse(ZeroLocations, ZeroRegions(ZeroLocations), ZerosMaxima(ZeroLocations));
LZMaxSorted = sort(LocationsZerosAndMaxima);
ZeroRemapperMax = LZMaxSorted(end, :);
%%% Now the same for ZerosMinima, except to find the minimum, we have to reverse the order of labels.
Reverser = [0 max(ZerosMinima(:)):-1:1];
LocationsZerosAndMinima = sparse(ZeroLocations, ZeroRegions(ZeroLocations), Reverser(ZerosMinima(ZeroLocations) + 1));
LZMinSorted = sort(LocationsZerosAndMinima);
ZeroRemapperMin = Reverser(LZMinSorted(end, :) + 1);
%%% Create the zero remapper
ZeroRemapper = ZeroRemapperMax;
%%% Anywhere that disagrees, set the remapper to 0
ZeroRemapper(ZeroRemapperMax ~= ZeroRemapperMin) = 0;
%%% Pad for zeros in the ZeroRegions image (i.e., nonzero regions in
%%% the LabelMatrix)
ZeroRemapper = [0 ZeroRemapper];
%%% Finally, remap the ZeroRegions, and combine them with the
%%% LabelMatrix to fill the holes.
ZerosReplaced = ZeroRemapper(ZeroRegions + 1);
FilledLabelMatrix = full(max(LabelMatrix, ZerosReplaced));