8
$\begingroup$

I have a bunch of small images that I need to register with larger images with sub-pixel resolution (ie. based on image interpolation )

This works.. seeking advice to improve performance.

Create a sample image which is offset by a fraction of a pixel from a larger image:

 i0 = ColorConvert[ExampleData[{"TestImage", "Lena"}], "GrayScale"];
 shift = {-0.2, 0.3};
 w = 100;
 h = 200;
 x0 = {200, 300}
 small = ImageTake[ImageTransformation[i0, (# + shift) &, DataRange -> Full], 
     x0[[1]] + {1, h}, x0[[2]] + {1, w}];

ImageAlign nicely (quickly) locates the sub image in the host to within a pixel:

 foundregion = Transpose[{First@#, Last@# } &@Position[  Map[Norm,
     ImageData[ImageAlign[i0, small]], {2}], 
           x_ /; x > 0, {2}]] + {{1, 0}, {0, -1}} ;
 nadd = 3;
 localmatch = ImageTake[i0, Sequence @@ ( foundregion + nadd {{-1, 1}, {-1, 1}})];

now we have two small images, the left is 3 pix larger to avoid edge effects when we do interpolation.

 GraphicsRow[ {localmatch , small }]

enter image description here

here is where the time issue is, brute force image difference..

 err[x_?NumericQ, y_?NumericQ] :=
 Flatten[(ImageData[
    ImageSubtract[ ImageTake[ ImageTransformation[localmatch,
          (# + {x, y}) & , DataRange -> Full] , nadd + {1, h}, 
           nadd + {1, w}] , small ]]), 1] // Norm;

  (r = NMinimize[ err[x, y], {{x, -2, 2}, {y, -2, 2}}] ) // Timing

{26.722971, {0., {x -> 0.800107, y -> 0.300113}}}

  { ({x, y} /. Last@r ) -  foundregion[[All, 1]] + x0 , shift  } 

verify we recover the input:

{{-0.199893, 0.300113}, {-0.2, 0.3}}

Any thoughts on speeding it up?

edit

After making @nikie's changes (dramitcally faster) and playing with it a bit I found the way I was extracting the region of the aligned image was not reliable. this does the job well:

 foundregion = #[[Ordering[#][[{1, -1}]]]] & /@ 
                  (Transpose@Position[  Map[Norm,
                    ImageData[ImageAlign[i0, small]], {2}], x_ /; x > 0, {2}]);

However this is now my bottleneck. Frustrating, ImageAlign itself is very fast, but there seems to be no simple way to retrieve the actual displacement it has computed.

$\endgroup$
6
  • 4
    $\begingroup$ The simplest way is to replace (# + {x, y}) & with TranslationTransform[{x, y}] - that makes it about 10x faster on my PC (I'm guessing ImageTransformation is optimized for TransformationFunctions). Replacing NMinimize with FindMinimum makes it 10x faster, again. You can probably make it faster if you carefully think about how ImageTransformation interpolates between pixels (there might even be a closed form solution?), but that's a lot more work for 0.2 seconds. $\endgroup$ – Niki Estner Oct 25 '14 at 9:21
  • $\begingroup$ @nikie Thanks! you should post that as an answer. BTW FindMinimum produces a slightly less accurate result for some reason, but for present purpose that's not an issue. $\endgroup$ – george2079 Oct 27 '14 at 14:56
  • 1
    $\begingroup$ FindMinimum is of course only looking for a local minimum, not global. That might be the reason. $\endgroup$ – dr.blochwave Oct 27 '14 at 18:36
  • 1
    $\begingroup$ To extract the transformation directly from ImageAlign, use the secret option ImageAlign[i0, small, "FindGeometricTransformOutputQ" -> True] $\endgroup$ – Simon Woods Oct 27 '14 at 21:48
  • $\begingroup$ What about splitting off the answer part as a self-answer? $\endgroup$ – Yves Klett Oct 28 '14 at 9:53
3
$\begingroup$

So this doesn't sit unanswered, here is revised code incorporating @nikies comment: I also fixed some bugs (unrelated to performance) in the original.

 i0 = ExampleData[{"TestImage", "Lena"}];
 shift =  {-10.1, 3.6}; (* {row,col}  *)
 w = 80;
 h = 100;
 x0 = {220, 230};(*{row,col}*)
 small = ImageTake[ImageTransformation[i0,
           TranslationTransform[{1, -1} Reverse@ shift],
           DataRange -> Full], x0[[1]] + {1, h}, x0[[2]] + {1, w}]
 subregister[i0_, small_] := Module[{
      wh = ImageDimensions[small],
      foundregion, err, nadd = 2, rng, result, localmatch },
   subregister::badalign = "invalid initial alignment result";
   subregister::fiterror = "fit error too large";
   Catch[
     foundregion = #[[Ordering[#][[{1, -1}]]]] & /@ 
       Transpose@Position[  Map[Norm,
           ImageData[ImageAlign[i0, small]], {2}], x_ /; x > 0, {2}];
     If[ ! Subtract @@@ foundregion == -Reverse@wh ,
        Message[subregister::badalign]; Throw[Null]];
     rng = nadd + {1, #} & /@ wh // Reverse;
     localmatch = 
       ImageTake[i0, Sequence @@ ( foundregion + nadd {{-1, 1}, {-1, 1}})];
    err[x_?NumericQ, y_?NumericQ] :=Norm@Flatten[
       ImageData@
           ImageSubtract[ ImageTake[ ImageTransformation[localmatch,
             TranslationTransform[{x, y}] , DataRange -> Full] , 
               Sequence @@ rng] , small ], 1];
    result = FindMinimum[ err[x, y], {{x, -2, 2}, {y, -2, 2}}] ;
    If[result[[1]] > .1 , Message[subregister::fiterror]; 
         Throw[{{0, 0}, result[[1]]}]];
    { ({-y, x } /. Last@result ) + foundregion[[All, 1]] - {1, 1} , First@result} ]]
    subregister[i0, small] // Timing
    x0 + shift (* expected result *)

{0.780005, {{209.899, 233.6}, 0.0146732}}

{209.9, 233.6}

Note this isn't as robust as i'd like because ImageAlign occasionally returns a very poor fit (v.9). Just looking at the docs ImageAlign looks to be significantly revised for v10..

V10 update: using FindGeometricTransformOutputQ per comments

 subregister[i0_, small_] := {
      ImageDimensions[i0][[2]] -
      ImageDimensions[small][[2]] - #[[2]], #[[1]]}      &@
         ImageAlign[i0, small, 
            "FindGeometricTransformOutputQ" -> True][[2]][{0, 0}] 
 subregister[i0, small] // AbsoluteTiming

{0.000566317, {209.893, 233.596}}

x0 + shift (expected result)

{209.9, 233.6}

Just a little faster. (note a slight loss of accuracy for this example though)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.