Skip to main content

graphs and networks - Getting cell edges from a picture


I would like to write a program that, given a picture of an epithelium (2D cell array), for example



enter image description here


automatically detects the cell edges and returns the corresponding lattice. Illustratively,


enter image description here


Naturally, if such program is based on contrast and colour detection, the original picture might need to be edited so that cell membrane contrasts enough with cell interior. Furthermore, unlike the above sketch, I would need such polygons to be convex (maybe too tricky?).


Now, I know this might be a lot to ask, so as a first step I would like to know whether there are already inbuilt functions or packages that might help doing this type of image processing (maybe some neural network implementation?). Once edges are obtained as Line-type objects, for example, building the graph or mesh from them shouldn't be that hard.


Just as a reference, I would be interested in building something along the lines of Packing Analyzer.


Edit: Following Lukas Lang's answer below, it seems the presented code doesn't recognise images with "more evident" edges, like the image


enter image description here


or even a similar picture to the first


enter image description here



Image sources: 1, 2 and 3. Might have to do with the way the image is processed via preprocImg or the mergedCells function. Any ideas?



Answer



Here is an approach based on WatershedComponents and MorphologicalGraph. Some of the steps feel a bit over-complicated, so feel free to point out any improvements.


The end result is a Graph expression describing the cell walls:


enter image description here


Here is the code with some intermediate results:


Get the original image:


img = Import["https://i.stack.imgur.com/elbTN.png"]

enter image description here



Do some blurring & sharpening, followed by an extraction of the red color channel. The goal of this step is to get an image with the cell walls as visible as possible.


preprocImg = First@ColorSeparate@Sharpen[#, 5] &@Blur[img, 3]

enter image description here


The next step is the call to WatershedComponents. Unfortunately, I didn't manage to preprocess the image enough to get perfect results, so we have to postprocess them instead.


wsComponents = 
WatershedComponents[preprocImg, Method -> {"MinimumSaliency", .65}];
wsComponents // Colorize

enter image description here



As can be seen, some of the cells are split into multiple pieces. The idea of the next step is to exploit the fact that the cells are all convex. First, we compute the convex hulls of the individual components:


cellMeshes = Map[
ConvexHullMesh@*
Map[{#2, -#} & // Apply](*
convert from image coordinates to plot coordinates *)
]@
Values@GroupBy[First -> Last]@(* group positions by component *)
Catenate@
MapIndexed[List,
wsComponents, {2}](* add position to component indices *);

Show@cellMeshes

enter image description here


We can now merge those that overlap by some amount (I compare to the "reduced area", in analogy to the reduced mass from physics):


mergedCells =
Graph[(* create graph where overlapping cells are connected *)
cellMeshes,
If[(* check whether overlap is big enough *)
Area@RegionIntersection@##*(1/Area@# + 1/Area@#2) > 0.35,
UndirectedEdge@##,

Nothing
] & @@@ Subsets[cellMeshes, {2}](* look at all cell pairs *)
] // Map[RegionUnion]@*
ConnectedComponents(* merge overlapping cells *);
Show@mergedCells

enter image description here


Now we are almost done - we convert the result back into an image, so that we can finally use MorphologicalGraph. For this, we apply some styling to the regions and rasterize:


procImg = Region[(* apply cell styling *)
#,

BaseStyle -> {EdgeForm@{White, Thick}, FaceForm@Black}
] & /@ mergedCells //
Show[#, PlotRangePadding -> 0, ImageMargins -> 0] & //(*
remove image border *)
Rasterize[#, ImageSize -> ImageDimensions@img] & //
Binarize //
ImagePad[ImageCrop@#, BorderDimensions@#] &(* make border black *)

enter image description here


Now we are at the finish line - a call to MorphologicalGraph and some nice presentation is all that's needed now:



MorphologicalGraph[
#,
EdgeStyle -> Directive[Thick, Red],
VertexStyle -> Blue,
VertexSize -> 2,
Prolog -> Inset[img, {0, 0}, {0, 0}, ImageDimensions@img]
] &@procImg

enter image description here


Notes



The key difficulty with this approach is to get preprocImg to be sufficiently "nice" for WatershedComponents to work. For the three images in the question, the following three approaches seem to work:


img = Import["https://i.stack.imgur.com/elbTN.png"]
preprocImg = First@ColorSeparate@Sharpen[#, 5] &@Blur[img, 3]
wsComponents = WatershedComponents[preprocImg, Method -> {"MinimumSaliency", 0.65}];
Row@{img, preprocImg, wsComponents // Colorize}

enter image description here


img = Import["https://i.stack.imgur.com/5RPz5.png"]
preprocImg = ColorNegate@First@ColorSeparate@Sharpen[#, 5] &@Blur[img, 3]
wsComponents = WatershedComponents[preprocImg, Method -> {"MinimumSaliency", 0.65}];

Row@{img, preprocImg, wsComponents // Colorize}

enter image description here


img = Import["https://i.stack.imgur.com/dgz9H.jpg"]
preprocImg =
ColorNegate[20 (#2 - #)*#3] & @@ ColorSeparate@Sharpen[#, 3] &@
Blur[img, 10]
wsComponents = WatershedComponents[preprocImg, Method -> {"MinimumSaliency", 0.45}];

enter image description here




  • As can be seen, each image requires a different approach - unfortunately I couldn't get it to work with a single one yet

  • In the end preprocImg needs to be bright between the cells and dark inside the cells. For the first and second image, this is pretty straightforward using the brightness of the image. (Note that the image needs to be inverted in the second case) For the third image, I had to do some math on the color channels to get a meaningful result.

  • The blur radius is increased in the third case to smooth out the bright and dark areas.

  • The "MinimumSaliency" parameter of WatershedComponents can be used to control the number of cell "candidates" in wscomponents - the best value will depend on the contrast of preprocimg among other things.

  • The components in wscomponents need to resolve the individual cells - in the remaining steps, components are only merged, never split. Too many components on the other hand make the post-processing slow and unreliable (since the overlap criterion doesn't work anymore)


Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...