# Image segmentation by pixel value

I'm looking for a way to segment a grayscale image (or a 2D array) by its pixel values. We look at two neighboring pixels, if their values are close within a certain threshold, then we say they are connected. How to find all the connected segments in the image?

Here are some three examples,

example 1:

threshold = 1;
data = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};


give

{{1, 1, 1}, {1, 1, 1}, {1, 1, 1}}


This is because 1 is connected to 2, 2 is connected to 3, etc. And all the pixels are connected. So there is only 1 component in the image.

example 2:

threshold = 0;
data = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};


give

{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}


This is because all pixels are isolated and thus we have 9 components.

example 3:

threshold = 2;
data = {{5, 0, 2, 3, 10}, {9, 10, 10, 2, 2}, {10, 10, 3, 3, 3}, {8, 3, 1, 7,
5}, {4, 5, 0, 5, 8}};


give

{{1, 2, 2, 2, 3}, {4, 4, 4, 2, 2}, {4, 4, 2, 2, 2}, {4, 2, 2, 2,
2}, {2, 2, 2, 2, 2}}


This question is a spectral version of this general question: Is there a way to supply a custom function to specify the connectivity in a function like MorphologicalComponents? For example, something like:

MorphologicalComponents[data, Abs[#1 - #2] <= threshold &]

• You could use a pair of difference filters (x and y), threshold and add the result to get the boundaries between the components, invert and use MorphologicalComponents. But you'll need some post-processing to fill the borders between the components Jul 8, 2017 at 6:56

Update: If you want 4-neighborhood, you can use MorphologicalComponents to do most of the work, which is fast and easy to implement (that was my original attempt, see below). But I don't think this can be adapted for 8-neighborhood. For 8-neighborhood, I would implement the standard 2-pass connected component labeling algorithm (this might be what MorphologicalComponents does internally). The idea is:

• in the first pass, scan pixels from top to bottom
• if there's a connected neighbor pixel left above, assign that label to the current pixel
• if there are more than two different labels, store the label pair for later
• if there's no connected neighbor, assign a new label

So, the pipeline looks something like this: So the first pass assigns the labels 4 and 5 to the pixels at the right, and stores the information that 4 and 2 should be in the same component, and that 5 and 1 should be in the same component. The second pass then only requires a table lookup for each pixel.

The first pass does most of the work, and returns a "preliminary" label matrix, and a set of label pairs that should go into the same component.

Sadly, there's no efficient compilable set-data structure, so I'm using a global Association variable (connectedIndices) which means a few MainEvaluate calls in the compiled function:

isConnected = With[{t = threshold}, Compile[{x, y}, Abs[x - y] <= t]];

firstPass = With[{connectedFn = isConnected},
Compile[{{pixels, _Real, 2}},
Module[{x, y, componentIndex, neighborOffsets, offset, index,
newIndex, componentCount = 0, w, h, parent, child, relabel},
{h, w} = Dimensions[pixels];
componentIndex = ConstantArray[0, {h, w}];
neighborOffsets = {{-1, -1}, {-1, 0}, {-1, 1}, {0, -1}};
(* scan every pixel *)
connectedIndices = <||>;
Do[
(
index = 0;
Do[
(* find connected neighbors above and to the left *)
If[y + offset[] >= 1 && x + offset[] >= 1 &&
y + offset[] <= h && x + offset[] <= w &&
connectedFn[pixels[[y + offset[], x + offset[]]],
pixels[[y, x]]],
newIndex = componentIndex[[y + offset[], x + offset[]]];
If[index != 0 && index != newIndex,
(* more than one label found in neighborhood:
newIndex and index really are the same component -
save information for second pass *)
AssociateTo[connectedIndices, MinMax[{index, newIndex}] -> 1];
];
index = newIndex;
], {offset, neighborOffsets}];
If[index == 0,
index = ++componentCount];
componentIndex[[y, x]] = index;
),
{y, 1, Length[pixels]},
{x, 1, Length[pixels[]]}
];
Return[componentIndex]
], {{relabelConnectedComponents[__], _Integer, 1}},
CompilationOptions -> {"InlineCompiledFunctions" -> True},
CompilationTarget -> "C"]]


The second pass is straightforward: We can use graph functions to get a relabeling-lookup table from the connected label pairs:

Clear[relabelConnectedComponents]
relabelConnectedComponents[connectedIndices_, componentCount_] :=
Module[{comp, relabel, count},
comp = SortBy[
ConnectedComponents[
Graph[Range[componentCount],
UndirectedEdge @@@ connectedIndices]], Min];
relabel = ConstantArray[0, componentCount];
Do[
Do[relabel[[oldIndex]] = index, {oldIndex, comp[[index]]}],
{index, Length[comp]}];
relabel]


(note that this graph only contains one node per component not per pixel, so this is much much faster than using ConnectedComponents on the pixels directly)

And apply that lookup table:

Clear[secondPass]
secondPass[componentIndex_] := Module[{relabel},
(* second pass - relabel connected components *)
If[Length[Keys@connectedIndices] == 0, componentIndex,
relabel =
relabelConnectedComponents[Keys[connectedIndices],
Max[componentIndex]];
relabel[[#]] & /@ componentIndex
]];


Usage: secondPass@firstPass[data]

This seems to be about 5-6 times slower than the built-in MorphologicalComponents, which means it's about 50% slower than the method below, which works on upsample data. But it should work for any neighborhood.

(Original answer, using MorphologicalComponents)

First, two small utility functions: upsample duplicates every value in a list, showGrid displays a matrix of numbers:

upsample = Riffle[#, #] &;

showGrid[d_, opt___] :=
Grid[d /. n_?NumberQ :> Item[n, Background -> ColorData[n]], opt,
ItemSize -> {1.5, 2},
Dividers -> {{{{True, False}}, -1 -> True}, {{{True, False}}, -1 ->
True}}]


You could upsample the data (duplicate every row and column):

dataUpsampled = upsample /@ upsample[data];
showGrid[dataUpsampled] Then you take the differences between rows, columns and across corners, and compare the differences to the threshold:

padAndThreshold =
Dimensions[dataUpsampled]] &;
dataUpsampled[[2 ;;, 2 ;;]] - dataUpsampled[[;; -2, ;; -2]]];
dataUpsampled[[2 ;;, ;; -2]] - dataUpsampled[[;; -2, 2 ;;]]];


and combine them:

boundaries = UnitStep[-(dx + dy + dxy + dyx)];
showGrid[boundaries] now you can use MorphologicalComponents on this array:

showGrid[MorphologicalComponents[boundaries]] To get the original sized result, simply remove the even rows and columns:

comp = MorphologicalComponents[boundaries][[;; ;; 2, ;; ;; 2]];
showGrid[comp, Dividers -> All] • I always don't konw why the most right-bottom number is 2..
– yode
Jul 8, 2017 at 8:14
• @yode: Because he's using 8-neighborhood and you're using 4-neighborhood. Jul 8, 2017 at 8:16
• Thanks,changed.
– yode
Jul 8, 2017 at 8:30
• Is it possible to achieve the same without duplicating the image size (with Image as input)? Jul 9, 2017 at 12:44
• @AlexeyPopkov: Good question. All answers so far multiply the data (think of the original pixels as nodes and the three "extra" pixels as nodes connecting them). You could write your own flood-fill or component labeling algorithm of course, but I'm not sure if that's going to be faster than using an inbuilt function on larger data. Jul 9, 2017 at 13:01

threshold = 2;
coor = Catenate[Array[{##} &, Dimensions[data]]];
Rotate[g = EdgeDelete[NearestNeighborGraph[coor,DistanceFunction ->
ChessboardDistance], _?(Abs[Subtract @@ (# /. rule)] > threshold &)], -Pi/2] MatrixForm[SparseArray[Catenate[Tuples /@ MapIndexed[Rule, ConnectedComponents[g]]]]]


$\left( \begin{array}{ccccc} 4 & 1 & 1 & 1 & 3 \\ 2 & 2 & 2 & 1 & 1 \\ 2 & 2 & 1 & 1 & 1 \\ 2 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ \end{array} \right)$

• Thanks for the great answer! But it seems that your solution scales poorly with the size of data. For instance, data = RandomInteger[20, {100, 100}]; takes about 60 seconds for your solution, while it takes about 4 seconds for my dfs solution, and takes about 0.01 second for @nikie's solution. Jul 10, 2017 at 6:41

Here is a solution based on yode's code, but it uses depth-first search which runs in linear time with the size of data.

N0 = First@Dimensions[data];
coor = Catenate[
CoordinateBoundingBoxArray[{{1, 1}, Dimensions[data]}]];
g = NearestNeighborGraph[coor, DistanceFunction -> ChessboardDistance];
vls = VertexList[g];
els = {};

isConnected = (If[
Abs[data[[Sequence @@ #[]]] - data[[Sequence @@ #[]]]] <=
threshold, AppendTo[els, UndirectedEdge[#[], #[]]]];
Pause[0.1]; &);

Dynamic[connect =
Graph[vls, els, VertexLabels -> "Name",
VertexCoordinates -> (coor /. {i_, j_} -> {j, N0 - i})]]

DepthFirstScan[g, {"FrontierEdge" -> isConnected,
"BackEdge" -> isConnected, "ForwardEdge" -> isConnected,
"CrossEdge" -> isConnected}]


This shows the action of the process • Doesn't ConnectedComponents have O(n) running time, too? (MorphologicalComponents probably too) If I had to guess, I'd say Mathematica's ConnectedComponents is probably implemented using depth-first search, too. Jul 8, 2017 at 20:47
• @nikie You are right. I just saw that yode edited his answer. He was using RelationGraph before, which runs at O(n^2). Thanks for your answer by the way. Jul 8, 2017 at 20:55