EDIT: I've found a bug in the code; Now the result looks a lot better, too.
This isn't perfect, but it's a start.
The first step is to find the "chain elements" (what are those anyway? I'm guessing cells?)
The chain elements have a distinctive scale, so I can filter them out easily using a median filter:
Subtracting the median filtered image from the original removes the background:
img = Import["https://i.stack.imgur.com/gGTro.jpg"];
(* I apply the median filter to a downsampled image - the result is almost the
same, but it's much much faster this way *)
diff = ImageAdjust[
ImageDifference[
ImageResize[MedianFilter[ImageResize[img, 250], 10],
ImageDimensions[img][[1]]], img]]
I don't want to work with RGB images, so I simply take the mean of the three channels. (ColorConvert[..., "Grayscale"] doesn't make much sense here, since these aren't "natural" colors, where the green channel contains most of the brightness).
meanImg = Total[ImageData[diff], {-1}]/3;
To find the centers of the chain links, I use a laplacian of gaussian filter, with a filter size that's approximately as large as the chain links themselves.
σ = 25;
log = Image[LaplacianGaussianFilter[meanImg, σ]*-σ^2]
I'm looking for local maxima in this image that are brighter than some threshold:
maxima = ImageMultiply[MaxDetect[log], Binarize[log, 0.5]];
pts = ComponentMeasurements[maxima, "Centroid"][[All, 2]];
This finds the maxima, but the list contains some "duplicates", where the same object contains multiple maxima in the LoG image. I can remove those relatively easily by using Mathematica's MeanShift function to "cluster" points that are close together, then removing duplicate points:
pts = Union[MeanShift[pts, 25, MaxIterations -> 10]];
In addition to the chain element centers, I'd like to have an estimate of the direction at each point. I can get a good estimate using calculus (I calculate the angle of the major eigenvector of the Hessian matrix at each point):
Clear[gaussianDerivative];
Do[gaussianDerivative[i, 2 - i] =
GaussianFilter[meanImg, 2 σ, {i, 2 - i}], {i, 0, 2}];
angles = Image[
ArcTan @@
Eigenvectors[{{m[0, 2], m[1, 1]}, {m[1, 1], m[2, 0]}}][[1]] /.
m -> gaussianDerivative];
getDirectionAt =
Function[pt,
With[{α = π/2 - ImageValue[angles, pt]}, {Cos[α],
Sin[α]}]];
Show[img,
Graphics[{Red, Thick,
Line[{# - 15 getDirectionAt[#], # + 15 getDirectionAt[#]}] & /@
pts}]]
Now the idea is to :
- pick a point at random
- pick it's nearest neighbor
- estimate where the next point in the chain would be - if there's a point close to that location, add it to the list
- rinse and repeat until no more points are found
This is the code:
nf = Nearest[pts];
continueChain = Function[chain,
Module[{expectedNextPoint, nearestNextPoint, direction},
(
direction = chain[[-1]] - chain[[-2]];
expectedNextPoint = chain[[-1]] + direction;
nearestNextPoint = nf[expectedNextPoint][[1]];
If[
Norm[nearestNextPoint - expectedNextPoint] <
Norm[direction]*0.75 &&
Abs[
getDirectionAt[chain[[-1]]].getDirectionAt[
nearestNextPoint]] > 0.5 &&
Norm[direction] < 100
, Append[chain, nearestNextPoint], chain]
)]];
findRandomChain[] := Module[{chain},
(
chain = nf[RandomChoice[pts], 2];
chain = FixedPoint[continueChain, chain];
chain = FixedPoint[continueChain, Reverse[chain]]
)]
This function finds one random chain. If I call it often enough, I'll (almost always) get all of the potential chains:
allChains =
SortBy[Union[Table[findRandomChain[], {1000}]], -Length[#] &];
The only problem is that some of those chains "overlap", i.e. some points are members of more than one chain. I think a good solution to this problem would be to apply a mixture of gaussians/expectation maximization algorithm to decide which point belongs to which chain. But it's getting late here, so I'll just go with a simple greedy algorithm:
- Start with the longest chain
- remove the points from the chain that are already used in some other chain
- If the chain still contains at least two or more points, keep it
- Remove all points in this chain from the set of "unused" points
- repeat for all chains
This is it:
removeChainPoints = Function[{unusedPoints, chain},
With[{modifiedChain = Select[chain, MemberQ[unusedPoints, #] &]},
(
If[Length[modifiedChain] >= 2, Sow[modifiedChain]];
Complement[unusedPoints, chain]
)]];
noOverlap = Reap[Fold[removeChainPoints, pts, allChains]][[2, 1]];
Now, noOverlap contains a set of non-overlapping chains:
totalChainLength = Total[Norm /@ Differences[#]] &;
colorFn = ColorData[3];
Show[
ColorConvert[img, "Grayscale"],
Graphics[
MapThread[
Module[{direction, normal},
(
direction = #1[[1]] - #1[[-1]];
direction = direction/Norm[direction];
normal = {-direction[[2]], direction[[1]]};
{
colorFn[#2],
Thick,
Line[#1],
Dotted,
Line[{#1[[1]], #1[[1]] + 50*normal, #1[[-1]] +
50*normal, #1[[-1]]}],
Text[totalChainLength[#1], Mean[#1] + 50*normal]
}
)] &,
{noOverlap, Range[Length[noOverlap]]}]]]
I didn't spend much time optimizing the thresholds in continueChain, so you could try to improve those. But you could probably come up with a smarter algorithm than simply adding one nearest point after the other. There's probably a well-known graph algorithm that will find globally optimal solutions in no time. (FindCurvePath almost does what I want, but I couldn't find a way to tweak the algorithm, so it groups far too many points into one curve.)
This is the result for the other image:
