# Implementing fast linear search for perfect vertical image alignment

After experimenting with image alignment aprroaches presented here and here and some other approaches involving ImageDistance and NMinimize, I realized that the most robust approach is the direct search, and it is probably can be implemented effficiently using Compile or FunctionCompile (but my experience with them is still very limited).

Suppose we have two overlapping screenshots with equal dimensions (and equal number of channels), but originally saved in a lossy format like JPG:

imgs = Import /@ {"https://i.imgur.com/sXrjTom.jpg", "https://i.imgur.com/gWvFAod.jpg"};
{i1, i2} = ImageTake[#, {230, 2210}] & /@ imgs


We need to find the correct vertical alignment for them. Unfortunately, ImageAlign is not only slow, but also gives unstable results, as I describe in this answer. So the best what we can do is to implement our own alignment algorithm.

The basic algorithm is fairly simple. Supposing that the images are successive, and the second image continuates the first:

1. Get ImageDimensions and convert our images into ImageData:

{width, height} = ImageDimensions[i1];
id1 = ImageData[i1, "Byte", Interleaving -> True];
id2 = ImageData[i2, "Byte", Interleaving -> True];

2. Define a distance finction depending on the alignment:

distance[align_] :=
Total@TakeLargest[Abs[Flatten[id1[[-align ;;]] - id2[[;; align]]]], 100]

3. Calculate the distances for all alignments starting from some minimal value:

distances = Monitor[Table[{align, distance[align]}, {align, 20, height}], align];

4. Take the alignment corresponding to the minimal distance as the correct alignment:

MinimalBy[distances, Last]

{{613, 2653}}

ListLinePlot[distances, PlotRange -> All]


Another obvious variation would be calculating the distances until we find the distance lesser than some predefined threshold:

thr = 5*width;
dist = Infinity;
align = 19;
While[dist > thr && align <= height - 1,
align += 1;
dist = distance[align]];
{align, dist}

{613, 2653}


What would be the most efficient way to implement both these approaches as compiled functions?

• I don't think that compiling would help much. I'd rather suggest to use ParallelTable. Another strategy might be to first blur i1 and i2, e.g. with GaussianFilter. The would smear out the peak in the distance plot such that it might become detectable by Table[{align, distance[align]}, {align, 20, height,step}] with larger stepsize step, e.g. step = 10. If you have found the best align for the blurred data, you can refine it by search between align-step and align+step. It is also possible to do more than one level of refinement... Commented Apr 16, 2022 at 13:42
• Another, similar strategy could be to Downsample or ArrayResample the images, derive a first guess for align from the downsamples, and then do a similar refinement. Commented Apr 16, 2022 at 13:44

Another, similar strategy could be to downsample the images by employing ArrayResample, derive a first guess for align from the downsamples, and then do a refinement. Somewhat like the following:

resample[a_, shrinkfactor_Integer] := Module[{m, n, l},
{m, n, l} = Dimensions[a];
ArrayResample[
a, {Quotient[m, shrinkfactor], Quotient[n, shrinkfactor], l},
"Bin"]
];

findAlignment[a_, b_] := MinimalBy[
ParallelTable[
{align,
Total[
TakeLargest[Flatten[Abs[a[[-align ;;]] - b[[;; align]]]],
100]]},
{align, 20, Min[Length[a], Length[b]]}
],
Last
][[1, 1]];

refineAlignment[a_, b_, guess_, shrinkfactor_] := MinimalBy[
ParallelTable[
{align,
Total[TakeLargest[Flatten[Abs[a[[-align ;;]] - b[[;; align]]]], 100]]},
{
align,
Max[1, shrinkfactor (guess - 1)],
Min[shrinkfactor (guess + 1), Min[Length[a], Length[b]]]
}
],
Last
][[1, 1]];

shrinkfactor = 4;

guess = findAlignment[resample[ImageData[i1], shrinkfactor],
resample[ImageData[i2], shrinkfactor]]; // AbsoluteTiming // First
align = refineAlignment[ImageData[i1], ImageData[i2], guess,
shrinkfactor] // AbsoluteTiming


0.423536

{0.391173, 613}

This is about 14 times as fast as calling findAlignment[ImageData[i1], ImageData[i2]]. Nonetheless, it still finds the same result.

• Thank you, works good at my image set. However, I think that it still can be fooled by parasite minimums which in some cases become almost equal to the true alignment. This is a consequence of the fact that with ArrayResample we loose the most remarkable feature of the algorithm which makes it robust: the true minimum is always very sharp and very low as compared to any parasite minimum. ArrayResample not only introduces additional parasite minimas, but also makes the true minimum compatible with the latter. Commented Apr 16, 2022 at 15:02
• Yes, of course. One can remedy this to some extent by taking not just the minimum alignment from the coarsened images, but, say, the best 10 alignments. Commented Apr 16, 2022 at 15:05
• Probably you meant Min[Length[a], Length[b]], not Max[Length[a], Length[b]]. Commented Apr 16, 2022 at 15:19
• Good point! Thanks! =) Commented Apr 16, 2022 at 16:51

Here is a way to speed up the search for the perfect alignment without loosing the most remarkable feature of the algorithm given in the question, which makes it robust: the true minimum is always very sharp and very low as compared to any parasite minimum.

1. It keeps arrays packed.
2. Even without ParallelTable, it uses all available processor cores for the computation.

We can safely downsample both images horizontally, but not vertically. And we also can safely limit the height of the window for computing the distance:

{width, height} = ImageDimensions[i1];
id1 = ImageData[i1, "Byte", Interleaving -> True];
id2 = ImageData[i2, "Byte", Interleaving -> True];

Clear[distance]
distance[id1_, id2_, align_, maxWindow_, vstep_] :=
Total@TakeLargest[
Abs[Flatten[
id1[[-align ;; Min[-1, maxWindow - align - 1], ;; ;; vstep]] -
id2[[;; Min[align, maxWindow], ;; ;; vstep]]]], 100]


Now

distances =
Monitor[Table[{align, distance[id1, id2, align, 200, 10]}, {align, 20, height}],
align]; // AbsoluteTiming // First
MinimalBy[distances, Last]

1.54873

{{613, 203}}

ListLinePlot[distances, PlotRange -> All]


Even without ParallelTable, the minimum is found in 1.5 seconds, and it is as sharp and low, as originally!

Update

It looks like that we get much more stable results with low windows heights when instead of Total[TakeLargest[dists, 100]] we use Max[Max[dists] - Min[dists], Total[TakeLargest[dists, 100]]/100.]:

Clear[distance]
distance[id1_, id2_, align_, maxWindow_, vstep_] :=
Module[{dists =
Abs[Flatten[
id1[[-align ;; Min[-1, maxWindow - align - 1], ;; ;; vstep]] -
id2[[;; Min[align, maxWindow], ;; ;; vstep]]]]}, {align,
Max[Max[dists] - Min[dists], Total[TakeLargest[dists, 100]]/100.]}]

distances =
Monitor[Table[distance[id1, id2, align, 200, 10], {align, 1, height}], align]; //
AbsoluteTiming // First
MinimalBy[distances, Last]

0.997083

{{613, 3}}

ListLinePlot[distances, PlotRange -> All]


Update 2

Significant speed gains along with significantly less memory consumption can be achieved by storing the intermediate image data as a NumericArray and reversing the data of the first image (to take lines from the beginning, not from the end):

height = ImageDimensions[i1][[2]];
vstep = 10;
id1 = ImageInternalImageData[i1, Interleaving -> True,
DataReversed -> True][[All, ;; ;; vstep]]
id2 = ImageInternalImageData[i2, Interleaving -> True][[All, ;; ;; vstep]]


Clear[distance]
distance[id1_, id2_, align_, maxWindow_] :=
Module[{dists =
Abs[Flatten[
Normal@id1[[align ;; Max[1, align - maxWindow + 1] ;; -1]] -
Normal@id2[[;; Min[align, maxWindow]]]]]},
{align, Max[Max[dists] - Min[dists], Total[TakeLargest[dists, 100]]/100.]}]

distances = Monitor[Table[distance[id1, id2, align, 200], {align, 1, height}], align]; //
AbsoluteTiming // First
MinimalBy[distances, Last]

0.74263

{{613, 3}}


A complete implementation can be found in this answer.

Here is a method with two caveats.

(1) It is not quite linear but rather it's softly linear (it has a logarithmic factor)..

(2) While the result will show that it works quite well, it was failing miserably when I tried it last night. I can think of a few reasons for that. In any case, I had to restart my machine today and now I get a viable result.

The rough idea is to apply a correlation to the images and see where there is a spike. We expect this to happen at a location where they overlap.

Here are some refinements (which might or might not be necessary).

(1) The image is white-dominated. This corresponds to a value of 1 in the image arrays for one of the three scales. So it can contribute unduly to the correlation values. To remove this possibility I reverse the images, that is, subtract their values from 1.

(2) I flatten the resulting arrays. I do this because I know what to make of the spike in the result. Very likely this can be done without flattening.

imgs = Import /@ {"https://i.imgur.com/sXrjTom.jpg",
"https://i.imgur.com/gWvFAod.jpg"};
{i1, i2} = ImageTake[#, {230, 2210}] & /@ imgs;
reverseImageArrays = Map[1 - ImageData[#] &, {i1, i2}];
flatarrays = Map[Flatten, reverseImageArrays];


I want to see where the second image shows a solid overlap with the first, so I use the second as first argument to ListCorrelate below. The third argument of 1 forces a cyclic use of the second argument.

lcorr = ListCorrelate[flatarrays[[2]], flatarrays[[1]], 1];


Find the location of the max value. For purposes of comparison, also find the min and the mean. The idea is that at overlaps we expect a spike and elsewhere we expect a lower fluctuation from multiplying uncorrelated values, many of which are zeros in one or both factors.

Mean[lcorr]
MinMax[lcorr]
maxpos = First[FirstPosition[lcorr, Max[lcorr]]]

(* Out[59]= 25498.3

Out[60]= {10112.9, 60888.6}

Out[61]= 4432321 *)


Now find the row where this maxes out. One thing is quite encouraging: this next computation actually gives an integer (as opposed to rational) value.

(maxpos - 1)/
Apply[Times, Rest[Dimensions[reverseImageArrays[[1]]]]]

(* Out[53]= 1368 *)


(Why do I subtract 1 in the numerator? It had to be an increment of either a -1, 0, or 1, and that's the value that gave the integer result. Which is to say, there was a bit of voodoo in that.)

Now we can check this result.

Image[1 - reverseImageArrays[[1]][[1368 ;;]]]


Image[1 -
reverseImageArrays[[2]][[1 ;;
Length[reverseImageArrays[[2]]] - 1368]]]


Visually this seems to be about as good as can be.

• Thank you very much, very interesting! However, for the first image the span should be 1368 + 1 ;; for the exact match. Commented Apr 17, 2022 at 18:13
• Further to your comment, the correct row is found by (maxpos - 1)/ Apply[Times, Rest[Dimensions[reverseImageArrays[[1]]]]]. I now see why this works out: just consider the case of two identical images, where we want the result to be 1. (You maybe knew this already, but I had missed it.) Commented Apr 17, 2022 at 21:38
• Thank you again for the explanation, now it is clear: (maxpos - 1)/ Apply[Times, Rest[Dimensions[reverseImageArrays[[1]]]]] gives the number of rows from the start to the max value. In the case of identical images it is zero, and the max value is at the first element (zero shift). This explains why we must subtract 1 from maxpos: maxpos = 1 corresponds to zero shift. Commented Apr 18, 2022 at 4:15
• I've found that using the generalized correlation syntax it is possible (at least partially) to implement the algorithm from my answer here using ListCorrelate: lst=Max/@Abs[ListCorrelate[ImageData[i2,"Byte"][[;;200,;;;;21]],ImageData[i1,"Byte"][[;;,;;;;21]],{1,-1},{},Subtract,List]];Position[lst,Min[lst]] (returns {{1369}}). But unfortunately it is about 10 times slower than the implementation given in my answer... Commented Apr 18, 2022 at 13:55
• Unfortunately, it unpacks even arrays of machine reals, just try: On["Packing"];ListCorrelate[ImageData[i2][[;;200]],ImageData[i1],{1,-1},{},Subtract,List];. Also look at processor and memory usage: only one core is used, and it takes huge amount of memory (4 GB), and finally crashes the kernel... Commented Apr 19, 2022 at 1:28