# Align two images with respect to sperm head

My biologist friend has to align sperm images according to their heads. He usually does that manually by rotating each image. The problem is that there are so many images and it takes a lot of time. So, he asked me if there is any way to automate this?

To keep it simple, I have used only two images.

I tried using ImageAlign but I don't see the first image in the result. Is there any other way I can align these images w.r.t their heads?

{im1, im2} = {Import["https://i.stack.imgur.com/ijKgl.png"],
Import["https://i.stack.imgur.com/bi0MB.png"]};
ImageAlign[im1, im2, TransformationClass -> "Similarity"]


Edit:

I have added some more images which can be found here.

• Good question. 1. Exactly do you mean by "align the images"? Is it just the sperm-tail orientation or you really need to have the heads overlap when stacking the images? 2. "To keep it simple, I have used only two images." Actually, it might be beneficial to have a large collection of images to test the procedures people would come up with. Such a collection can be made by making variations to the two images your provided. (I might post code for that...) – Anton Antonov Jun 28 '17 at 14:42
• @AntonAntonov Thank you. He actually wants to align the heads, so as to study the behavior of tail motion. The images I have uploaded are from video1 of elifesciences.org/articles/05161. I could not convert the video into images as MMA says unsupported type of video file. Because of which I cropped two images from the video and uploaded here. To convert the video into Images I tried this: Import["elife-05161-media1.mp4", {"AVI", "ImageList"}]. If I can export the video into images, I will definitely put it here. – Anjan Kumar Jun 28 '17 at 14:58
• @AntonAntonov I have added a batch of images. Please check the edit. – Anjan Kumar Jun 28 '17 at 15:28
• I have a somewhat promising result, but not a complete solution. I will post an extended comment later today. (If someone does not come with solution...) – Anton Antonov Jun 28 '17 at 21:40
• You need to explain what you mean by "align the heads". For example, draw a red line in a few of the images by hand, starting at the center of the "head", pointing in the direction you expect the output images to be oriented. – Niki Estner Jun 29 '17 at 9:46

My answer is going to find the sperm, extract the coordinates of the midline and fit a polynomial to it to get the orientation near the head. The answer got a bit longer than I expected when I started, but I think each of the individual steps is simple enough and you'll need to do something similar anyway if you want to analyze the sperm's movement.

Ok, first step, load all the images and pick one to try each of the steps:

imgs = Import[
"https://www.dropbox.com/s/1tjw0agq8ap3a8n/elife.tiff?dl=1"];
img = imgs[[1]];


To get a clean binarization of the sperm, I'll use the Hessian matrix, i.e. the 2nd order derivatives at each point:

Do[derivative[i, 2 - i] =
GaussianFilter[ImageData[img][[All, All, 1]], 10, {i, 2 - i}], {i,
0, 2}];


To be precise, I'm going to look at the eigenvalues of the Hessian. If you think of the brightness value as a 2d function, it will look "ridge-like" (or "canyon-like") where the sperm is, i.e. it will have one large (absolute) eigenvalue, the other eigenvalue will be close to 0. At the ends and near the head, both (absolute) eigenvalues will be large, and everywhere else both eigenvalues will be small.

So, let's get the Hessian eigenvalues for each pixel. These are the eigenvalues for a generic symmetric matrix:

es = Eigensystem[{{m[2, 0], m[1, 1]}, {m[1, 1], m[0, 2]}}]


To prevent negative numbers under the square root in case of numerical errors, I'll replace the Sqrt with a clipped square root:

es = es /. Sqrt[x_] :> Sqrt[Clip[x, {0, Infinity}]];


And apply the resulting expression to the image derivatives:

eVals = es[[1]] /. m -> derivative;
eVals = eVals/Max[Abs[eVals]];

Image /@ Rescale[eVals]


As expected, the first eigenvalue is large (negative) and the second eigenvalue is close to zero at the sperm's tail. So we can binarize using these two conditions:

bin = Image[
UnitStep[-0.1 - eVals[[1]]] * UnitStep[0.15 - Abs[eVals[[2]]]]]


To make the following steps simpler, I'll fill all the holes and only select the largest connected component (that should be the sperm):

binSimple = SelectComponents[FillingTransform[bin], "Area", 1, Greater]


Thinning will thin this binary mask to a 1-pixel wide line:

HighlightImage[img, Thinning[binSimple]]


And we can get the points along that line using PixelValuePositions:

whitePts = PixelValuePositions[Thinning[binSimple], 1];


Unfortunately, these points are ordered from top to bottom, and we would like them ordered from head to tail of the sperm.

In simple cases, if the graph is just a simple path, you can use FindCurvePath or FindShortestTour to order these points. But in practice, Thinning will often produce tree, where we want a path that ignores the branches. So we want the diameter of the neighborhood graph, i.e. the shortest path between the two points that are the farthest apart:

graph = NearestNeighborGraph[whitePts, {All, 1.5}];

findDiameter[(g_)?UndirectedGraphQ] :=
Module[{d = GraphDistanceMatrix[g], u, v, pos},
pos = First[Position[d, Max[d]]]; {u, v} = VertexList[g][[pos]];
FindShortestPath[g, u, v]]

pts = findDiameter[graph];


Now pts contains the single longest line along the sperm's tail. But we don't know if pts[1] is the head or the tail. I'll simply assume that the head is always brighter than the tail (it was in your images), and reverse the sequence so pts[1] is the head:

endBrightness =
ImageValue[GaussianFilter[img, 10], pts[[{1, -1}]]][[All, 1]];
If[endBrightness[[1]] < endBrightness[[2]], pts = Reverse[pts]];


Now we can take the points nearest to the head:

headPts = Take[pts, UpTo[100]];



and fit a polynomial to them: the polynomial will map arc lengths to complex values, so I can fit the x and y coordinates in a single call:

arcLength = Accumulate[Prepend[Norm /@ Differences[N[headPts]], 0]];

fit = Fit[Transpose[{Rescale[arcLength], N[headPts . {1, I}]}],
u^Range[0, 4], {u}]


Now fit is a polynomial approximation to the sperm's tail near the head. The head is at u==0, so we get the head's position:

headPos = ReIm[fit /. u -> 0]


and the direction at the head's position:

dir = Normalize[ReIm[D[fit, u] /. u -> 0]]

Show[img,
ParametricPlot[ReIm[fit], {u, 0, 1}, PlotStyle -> {Dashed, Red}],


The neat thing here is that we use 100 points to fit 5 coefficients, so a lot of noise is averaged out.

To align the images, I'll choose a target position for the head:

targetPos = ImageDimensions[img]*{0.5, 0.2};


and rotate/translate the sperm so the head is at the position, looking down:

transform =
RotationTransform[ArcTan @@ dir - 90*Degree, headPos] @*

Show[ImageTransformation[img, transform, PlotRange -> Full,
DataRange -> Full],
Arrow[{targetPos, targetPos + {0, 50}}], Dashed,
Line[InverseFunction[transform] /@ pts]}]]


Whew. That was a lot of code. Let's put it all in a single function:

es = Eigensystem[{{m[2, 0], m[1, 1]}, {m[1, 1], m[0, 2]}}];
es = es /. Sqrt[x_] :> Sqrt[Clip[x, {0, Infinity}]];

findDiameter[(g_)?UndirectedGraphQ] :=
Module[{d = GraphDistanceMatrix[g], u, v, pos},
pos = First[Position[d, Max[d]]]; {u, v} = VertexList[g][[pos]];
FindShortestPath[g, u, v]]
align[img_] := (
Do[derivative[i, 2 - i] =
GaussianFilter[ImageData[img][[All, All, 1]],
10, {i, 2 - i}], {i, 0, 2}];
eVals = es[[1]] /. m -> derivative;
eVals = eVals/Max[Abs[eVals]];
bin = Image[
UnitStep[-0.1 - eVals[[1]]]*UnitStep[0.15 - Abs[eVals[[2]]]]];
binSimple =
SelectComponents[FillingTransform[bin], "Area", 1, Greater];
whitePts = PixelValuePositions[Thinning[binSimple], 1];
graph = NearestNeighborGraph[whitePts, {All, 1.5}];
pts = findDiameter[graph];
endBrightness =
ImageValue[GaussianFilter[img, 10], pts[[{1, -1}]]][[All, 1]];
If[endBrightness[[1]] < endBrightness[[2]], pts = Reverse[pts]];
arcLength =
fit = Fit[Transpose[{Rescale[arcLength], N[headPts . {1, I}]}],
u^Range[0, 4], {u}];
headPos = ReIm[fit /. u -> 0];
dir = Normalize[ReIm[D[fit, u] /. u -> 0]];
targetPos = ImageDimensions[img]*{0.5, 0.2};
transform =
RotationTransform[ArcTan @@ dir - 90*Degree, headPos] @*
Show[ImageTransformation[img, transform, PlotRange -> Full,
DataRange -> Full],
Arrow[{targetPos, targetPos + {0, 50}}],
Dashed, Line[InverseFunction[transform] /@ pts]}]]
)


and apply it to all the other images:

frames = Monitor[Table[align[imgs[[n]]], {n, Length[imgs]}], n];

ListAnimate[frames]


(As always, I've played with each of the filter sizes, thresholds etc. to get good results for this set of images. You might have to adjust them for different images. The binarization is probably the most fragile link of the whole chain, so you might want to play with different ideas here, too.)

• This is great -- instructive and exhaustive. (+1 of course). As for "the answer got a bit longer than I expected when I started [...]" I got the sense of that too with my experiments and gave up... – Anton Antonov Jun 29 '17 at 19:06
• @nikie Thank you for your answer. Will need some time to go through it. – Anjan Kumar Jun 30 '17 at 0:08
• @nikie Thanks again for your wonderful answer. Can you explain why you did not want to use Binarize[](an in-built function) to get the binary image? Is there some inherent disadvantage to it? – Anjan Kumar Jun 30 '17 at 10:33
• @AnjanKumar: I have tried Binarize on the raw image (it's always good to try in-built simple functions first), but it was hard to get the head and the tail as two different objects. That's why I did the Hessian eigenvalue thing. You could of course use Binarize in place of UnitStep, on the eigenvalue array - that was just convenience, because the data was already in list format at that point. – Niki Estner Jun 30 '17 at 14:02

Not an answer but an extended comment...

This is what I tried/did:

1. I cropped one of the images around the head;

2. experimented with aligning the images using different methods, and chose one that gives good results, ("KAZE");

3. tried to compute the angles of rotation using matching of keypoints, but I think some outliers make the results not that good.

### Definition of AngleOfImageAlignment

Clear[AngleOfImageAlignment]
AngleOfImageAlignment[i1_?ImageQ, i2_?ImageQ] :=
Block[{matches, a1, a2},
matches = ImageCorrespondingPoints[i1, i2, Method -> "KAZE"];
a1 = ArcTan @@@ (# - 1/2*ImageDimensions@i1 & /@ matches[[1]]);
a2 = ArcTan @@@ (# - 1/2*ImageDimensions@i2 & /@ matches[[2]]);
RootMeanSquare[Mod[a1 - a2, 2 Pi]]
];

• Thank you for your effort. It looks like a nice result. Will the results improve if we crop the small head a little bit more? – Anjan Kumar Jun 29 '17 at 13:10
• @AnjanKumar Yeah, I think that is worth experimenting with. I only experimented with two head crops. – Anton Antonov Jun 29 '17 at 13:23
• Is AngleOfImageAlignment[] a built-in function? I don't have that function in my MMA. – Anjan Kumar Jun 29 '17 at 13:42
• @AnjanKumar You can find a similar function in the help browser, ImageCorrespondingPoints>>Applications. (I included it in the answer.) – Anton Antonov Jun 29 '17 at 14:36
• Thank you for adding the definition. – Anjan Kumar Jun 30 '17 at 4:08