Align two images with respect to sperm head

Question

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.

Answer 1

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]]; 

HighlightImage[img, {{Dashed, Line[pts]}, {Thick, Line[headPts]}}]

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}], 
 Graphics[{Green, Arrow[{headPos, headPos + 100*dir}]}]]

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] @* 
   TranslationTransform[headPos - targetPos]; 

Show[ImageTransformation[img, transform, PlotRange -> Full, 
  DataRange -> Full], 
   Graphics[{Red, Arrowheads[Medium], 
   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]]; 
  headPts = Take[pts, UpTo[100]]; 
    arcLength = 
   Accumulate[Prepend[Norm /@ Differences[N[headPts]], 0]]; 
  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] @* 
    TranslationTransform[headPos - targetPos]; 
    Show[ImageTransformation[img, transform, PlotRange -> Full, 
    DataRange -> Full], 
   Graphics[{Red, Arrowheads[Medium], 
     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.)

Answer 2

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]]
  ];