IdentifySecondary Module - Doubt

cellprofiler

#1

Hi

I’m a new user to cell profiler and I’ve been trying to understand how the ‘propagation routine(IdentifySecondary)’ works and the implementation of propagate in Matlab.

The code is well commented - in fact the propagate section of IdentifySecondary.m is neatly divided into four steps. However, I am not able to understand Steps 3 and 4.

The documentation says:-

STEP 3: Remove objects that are not desired, edited objects. The
edited primary object image is used rather than the preliminary one, so
that objects whose nuclei are on the edge of the image and who are
larger or smaller than the specified size are discarded.

I do not understand how this is being done in the Matlab code. A sparse matrix is generated(why?). Some filtering is done and various new images are created.
For eg:- LargestLabelImage, HoleyPrelimLabelMatrixImage,
MaximaImage,MinimaImage,ZeroRegionImage etc.

Could anyone please explain what these images are and what is happening in Step 3 of propagation in IdentifySecondary.m(roughly lines 375-475)?

I’d be grateful for any pointers/information/links you can provide me, and I would be able to appreciate the module better.

Thanks for your time
Aravind


#2

Hi,

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;
    return;
end
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));

#3

Hi Thouis

Thanks for the reply. I am looking into the relevant files in the latest development version.

I tried my hand at implementing the IdentifySecPropagateSubFunction.cpp to make it work without MATLAB - using another library CImg which is publicly available. I just give the propagate function the 3 images(Original,Thresholded, and PrelimPrimarlyLabelMatrix Images) and the lambda factor as parameters and save the output file.

I just can’t get it to produce the same output given by the MATLAB-with Cell Profiler Developer version. Is there a difference between the cpp and the compiled mexmac versions provided?

My implementation gives an output which almost resembles the PrelimPrimaryLabelMatrixImage whereas the correct output looks a little bit like the ThresholdedImage. What could I be missing?

Also - it takes a lot longer for my C++ program to produce the result - Is matlab optimized in some way to execute C++ functions rapidly?

Thanks
Aravind


#4

I’m not sure what might be different. The .cpp code was used to generate the .mex file. If you’re getting different answers, you might go through step by step and compare the two algorithms. Also, the .cpp code doesn’t use anything special to run fast, it just gets two arrays from matlab and uses those directly.


#5

hi Thouis

thanks for the quick reply.

I havent been able to get CP(MATLAB licenses, version conflicts etc) do what I want - thats one of the reasons why I’m trying to implement things on my own.

Do you have any idea on what software I can use to track/debug c++ along with matlab other than Microsoft Visual Studio? any open source software?

Thanks
Aravind