Implementing fast linear search for perfect vertical image alignment

Question

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?

Answer 1

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.

Answer 2

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.

Additional advantages of this implementation:

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 = Image`InternalImageData[i1, Interleaving -> True, 
   DataReversed -> True][[All, ;; ;; vstep]]
id2 = Image`InternalImageData[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.

Answer 3

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.