how to speedup subpixel registration

Question

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

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.

Answer 1

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)