Is it possible to extract edges shared between two labeled objects?

For the past couple of weeks I have been trying to implement a cell-tracking and lineage mapping algorithm using which biological cells can be tracked in a time-lapse image. The code can be found at GitHub. It is not necessary to check the GitHub code for the question at hand.

One thing I wish to implement on top of tracking cells is to be able to label edges shared between two cells (or simply track individual edges as well).

As an example below: is it possible to get all the edges shared between two cells which I can label as {1,2} ... where 1 and 2 are indices of cells?

I know that ComponentMeasurements perhaps cannot be used to get edges of the shared cells. A close but wonderful approach is mentioned by @Alexey Popkov for a seemingly unrelated answer:

But I do not know if or how can I use centroids of the cells and Alexey's code to label edges based on the two cells sharing it.

The matrix from which a Colorize image is generated is created as segmentImage[img] (a different image is posted below):

img =

Clear[segmentImage];
Module[{seg, areas, indexMaxarea, maxArea},
areas = ComponentMeasurements[seg, "Area"];
{indexMaxarea, maxArea} = First@MaximalBy[areas, Last] /. Rule -> List;
If[maxArea > threshCellsize,
ArrayComponents[seg, Length@areas, indexMaxarea -> 0], seg] ~ Dilation~1
];

• also i think this must be incorporated in ComponentMeasurements in the future. It is perhaps a fundamental property which cannot be ignored. – Ali Hashmi Jun 8 '17 at 22:39
• A possible scheme that i could think of today is to choose a random point on each of the edges and take the union of the labels found within a search radius from that random point. If one of the labels is zero then the edge belongs to the non-zero label. If the edge has two labels as neighbours then we get the edge labelled as i,j – Ali Hashmi Jun 9 '17 at 9:21
• How do you wish to label the edges? Via Tooltip (what requires conversion into a Graphics object) or in some other way? – Alexey Popkov Jun 9 '17 at 10:40
• Also I should point out that in the sense of ComponentMeasurements the term "edge" isn't well-defined. In this and previous question you implicitly use your own definition of "edge" between adjacent regions. But you still haven't defined what you actually mean by the term "edge", especially how "edge" is related to the regions themselves. This makes your question very ambiguous and unclear. – Alexey Popkov Jun 9 '17 at 11:42
• Please formulate in clear words answers at least to the following questions: 1) From where "edges" come? 2) How "edge" is related to the regions? 3) How wide the "edge" is expected to be? 4) What information is primary/initial for you: "edge" or regions? In other words, do you generate regions from "edges" or "edges" from regions? I understand that there may be tautology in these questions but it reflects ambiguity and confusion of the original question. – Alexey Popkov Jun 9 '17 at 11:47

Solution: labeling cells and edges between them

Here is an efficient and straightforward solution for this problem. Note that all the code is for Mathematica 8.

First, you don't need Dilation, for obtaining the segmentation it is sufficient to specify CornerNeighbors -> False (this also is much more efficient):

img = Import["http://i.stack.imgur.com/2a2j6.png"];

cellM = MorphologicalComponents[ColorNegate@img, CornerNeighbors -> False];


Second, you could use MorphologicalTransform as a workaround for the ImageFilter bug as I suggest in this answer (what also is much more efficient):

edges = MorphologicalTransform[img, If[#[[2, 2]] == 1 && Total[#, 2] == 3, 1, 0] &];


(also consider applying Thinning first as I recommend in the "UPDATE 2" section here).

The label matrix for the edges:

edgesM = MorphologicalComponents[edges];


Note that cellM and edgesM has common labels what isn't appropriate. In order to make all the labels unique, we can add Max[cellM] to every non-zero label in edgesM:

nCells = Max@cellM;
edgesM = Map[If[# != 0, # + nCells, #] &, edgesM, {2}];


I use Map instead of Replace here in order to avoid unpacking of the matrix edgesM. An alternative to the addition is to multiply edgesM by a sufficiently large number (which will guarantee obtaining a set of labels not intersecting with cellM). In order to keep the labels informative we could use as a multiplication factor a power of 10:

factor = 10^Ceiling[Log10[nCells]];
edgesM = edgesM*factor;


This method is also more than 10 times faster than the previous.

Now we simply add the two matrices for obtaining the complete matrix of labels:

fullM = cellM + edgesM;


Obtaining the labels of neighborhood components is straightforward:

neighbours = ComponentMeasurements[fullM, "ExteriorNeighbors"];


If you wish to ignore components that are connected to the border, simply add "BorderComponents" -> False:

neighbours = ComponentMeasurements[fullM, "ExteriorNeighbors", "BorderComponents" -> False];


List of neighborhood cells for every edge:

• for the addition method (see above):

Select[neighbours, First@# > nCells &]

• for the multiplication method:

Select[neighbours, First@# >= factor  &]


How to visualize the result

RGB colorspace

Here is a way to visualize the matrix fullM assigning darker colors to the edges and lighter to the cells. The implementation is based on the code of ImageColorOperationsDumphashcolor (the default coloring function of Colorize). We have 3 times more edges than cells, so I split the full intensity range (in this case 0 .. 255) into two unequal parts: channel values from 0 to 170 are for edges, and values from 171 to 255 are for cells. The following implementation doesn't unpack fullM and works both for addition and multiplication methods shown above, the second argument is threshold value which will differ for those methods (nCells and factor - 1 correspondingly):

hashColor =
Compile[{{i, _Integer}, {thr, _Integer}},
Which[i == 0, {0, 0, 0},
i <= thr, 171 + IntegerDigits[Hash[i], 84, 3],
True, IntegerDigits[Hash[i], 170, 3]]/255.,
RuntimeAttributes -> {Listable}, "RuntimeOptions" -> {"Speed"}]


Here is how it can be used (addition method):

Image[hashColor[fullM, nCells]]


Or equivalently (but slower):

Colorize[fullM, ColorFunctionScaling -> False,
ColorFunction -> (RGBColor[hashColor[#, nCells]] &)]


Another (and potentially better) way to implement the colorizing function may go through IntegerPartitions, RandomChoice and SeedRandom (for reproducibility).

LAB colorspace

And here is similar approach but using the "LAB" colorspace which provides explicit control over lightness. According to the Documentation page for LABColor,

RGBColor approximately corresponds to l between 0 and 1, a between -0.8 and 0.94, and b between -1.13 and 0.94.

So I rescale l, a and b accordingly and split the lightness values into two diapasons: from 0.8 to 1 it is for the cells, and from 0.1 to 0.7 it is for the edges. The gaps between 0 and 0.1 and between 0.7 and 0.8 aren't used in order to make colors more easily distinguishable:

hashLABColor = Inactivate@Compile[{{i, _Integer}, {thr, _Integer}},
Module[{L, a, b},
{L, a, b} = IntegerDigits[Hash[i], 1000, 3];
Which[i == 0, {0, 0, 0},
i <= thr, {rl1, ra, rb},
True, {rl2, ra, rb}]],
RuntimeAttributes -> {Listable}, RuntimeOptions -> {"Speed"}];
hashLABColor = Block[{L, a, b}, Activate@With[{
rl1 = Rescale[L, {0, 1000}, {.8, 1}],
rl2 = Rescale[L, {0, 1000}, {.1, .7}],
ra = Rescale[a, {0, 1000}, {-.8, .94}],
rb = Rescale[b, {0, 1000}, {-1.13, .94}]}, Evaluate@hashLABColor]];


It can be used similarly to the previous:

Image[hashLABColor[fullM, nCells], ColorSpace -> "LAB"]


The colors here seem lesser diverse than the ones obtained within RGB colorspace in the previous section. The reason is probably that the algorithm selects too many LAB colors which are out of RGB gamut and hence look very close or identical when displayed with a monitor working withing sRGB colorspace.

• i need the Dilation in my code but for another reason (its a long story and i will tend not to digress). I found your way to be definitely fast. So i am accepting and upvoting your answer. I am using version 11.1 but i adopted latter part of your code and it works like a charm !! – Ali Hashmi Jun 11 '17 at 19:04
• never knew ExteriorNeighbors can come handy – Ali Hashmi Jun 11 '17 at 19:07
• thanks ! i have given credits to you here: github.com/alihashmiii/Lineage-Mapper--a-cell-tracking-method/… – Ali Hashmi Jun 13 '17 at 9:54
• @AliHashmi Thanks, but when giving credits it is good idea to include a link to the corresponding post (no one knows what is "the faster approach proposed by Alexey Popkov"). It will also help possible future contributors to understand the code. – Alexey Popkov Jun 13 '17 at 10:01
• you are right. will add the link to the post – Ali Hashmi Jun 13 '17 at 11:03

Although slow than what i would have hoped for but it does the job

(* borrowed from Alexey Popkov's answer: creates edges *)
createEdge[img_] := Module[{imagetemp},
imagetemp = ImageFilter[If[#[[3, 3]] == 1 && Total[#[[2 ;; -2, 2 ;; -2]], 2] == 3, 1,
0] &, img, 2];
MorphologicalComponents[imagetemp]
]

(* nearest function to determine which cell label is closer to edge *)
nearestCellFunc[matdim_, labeledMat_] := Module[{table},
table = Tuples[{Range[First@matdim], Range[Last@matdim]}];
]

(* discard edges that are 3 pixels or less *)
edgeLabels[edges_, thresh_: 3] := Keys@DeleteCases[
ComponentMeasurements[edges, "Count"],
HoldPattern[_ -> x_ /; x <= thresh]];

(* to relate the edge label to the cell label(s) sharing it *)

createEdgeAssoc[labeledMat_, img_] :=
Reap[Module[{nf, edges, pos, elabels, positioninfo, positionordering,
nf = nearestCellFunc[Dimensions@labeledMat, labeledMat]; (* nearest function
to tell which cell label is closer *)
edges = createEdge[img]; (* create edges *)
elabels = edgeLabels[edges]; (* take edges larger than 3 pixels *)
(* masks for edges are extracted *)
Do[
positioninfo = Position[masks[[i]], 1]; (* take positions of non-zero pixels *)
positionordering = Last@FindShortestTour[positioninfo]; (*shortest
path to reorder white pixel positions*)
pos = positioninfo[[Round@(Length[positionordering]/2)]]; (* get middle
pixel on the edge *)
Sow[i -> DeleteCases[Union@nf[pos, {All, 3}], 0] /. {x_Integer} :> x],
{i, elabels}]] (* find which cell labels are around the middle pixel of the edge  *)
;][[2, 1]]


using img

createEdgeAssoc[segmentImage[img], img]

(* {1 -> 1, 2 -> 2, 3 -> 4, 4 -> 3, 5 -> 5, 6 -> 6, <<122>>, 135 -> {41, 45},
137 -> 43, 138 -> {44, 45}, 140 -> 42, 141 -> 45, 142 -> 44}*)


the key is the edge label and the value is the cell label(s) surrounding the edge.

• The slowness is probably due to FindShortestTour – Ali Hashmi Jun 10 '17 at 19:06