Speeding up ImageCorrespondingPoints

Question

I have the following problem: I am taking free-handedly photographs from the same scene during e.g. hiking, because sometimes I don't want to lug around a tripod. The goal is e.g. to create later on with Mathematica an image with an artificial long term exposure (see here) or artificial "quasi"-aperture (see here). For that I create with Mma a binary Mask for the rigid regions BGMask and use this mask to find with ImageCorrespondingPoints the matching points between the images. (There is also ImageAlign but unfortunately it doesn't support a Mask, and I need that mask to get keypoints from rigid structures in the photographed scene.) After that I calculate the Transformation necessary to get the one Image transformed from one to the other by FindGeometricTransform. In this step I specify manually the expected transformation class, e.g. in case if the foreground or only the background of the scene is static. After that I export the image. The relevant code I'm using is given here:

RMA = ImageCorrespondingPoints[Image1, Image2, Masking -> BGMask, MaxFeatures -> 15];
tr = FindGeometricTransform[RMA[[1]], RMA[[2]], TransformationClass -> "Perspective"];
EntwackeltesBild = ImagePerspectiveTransformation[Image2, tr[[2]], DataRange -> Full, 
                     PlotRange -> Full];
Export[Name <> " Entwackelt.jpg", Show[EntwackeltesBild, 
                          ImageSize -> ImageDimensions[EntwackeltesBild][[1]]], 
    "CompressionLevel" -> 0];

While this method works well enough as you can see in the above image links, it takes unfortunately 1h to deshake an image, due to their size which is 21 MPixel. Given that sometimes I need to deshake up to 100 images this results in ~4 day long computations, for just a single photographic scene.

Thus my question, whether there is an efficient way to speed up the calculations besides using parallel-computing capabilities of Mma? Compiling isn't obviously intended for image processing functions (see SE: Mma compilable functions) and ImageCorrespondingPoints doesn't accept any Method-Options...

But of course there are 3 key components ImageCorrespondingPoints, FindGeometricTransform and ImagePerspectiveTransformation that might be slow, but somehow, as I am using only up to 15 key points to calculate the necessary transformation, I don't believe that FindGeometricTransform is the problem.

Edit 1: I will add later the corresponding time durations of each operation, but currently I am waiting for the images to get calculated...

Edit 2: I just saw on Mma-help that ImageCorrespondingPoints can take in the newest version also Options for Methods such as "AKAZE" "BRISK" or "ORB" and so on and obviously there is a lot of papers investigated the performance and stability of these methods (e.g. here (IEEE) here (PDF) or here (PDF)). Unfortunately, on my machine I got only Mathematica 10.2 so I can't test those features. I don't know where to look up which kind of algorithm Mma is using in my version.

Edit 3: Just to show what kind of images I am talking about (downsampled): and

Answer 1

Finally, I did some runtime testing with the above images in terms of size and using different transformation classes without using any masking. For the size-considerations I simply resized the whole image. This is the resulting graph:

I did the calculation with default setting up to the maximum image size of 21 MP, as you can see the other TranformatioClass-Settings converge more or less to this, so choosing a specific setting might make sense for later calculations, but it doesn't speed up the process. The green points is simple fit which fits the large image sizes relatively well with 10^-11.58 (numOfPxls)^2.1 for my system. I suppose Mma compares each feature from one image within the whole extents of the other, e.g. feature from the left lowest corner is compared in the image also if it matches the right above corner, in order to be thorough. But this "thorougness" is not really necessary in my case.

Following from this result it is clear what I'll need to do next: partition the image and get the image for each sub-image. (e.g. if I partition the image of 21 MP into 4 subimages with 5.05 MP then the computation will take for each one ~ 340 s, and in total 1360 s which is 4.5x faster than for the whole image). The speed-up will be even more effective when the image is partitioned into more parts.

Of course there are some limitations:

• the camera shake needs to be moderate, because otherwise scene features will get into different subimages and won't be recognized this way
• the more subdivisions are used for the main image, the less features are found due to the borders where the features don't exist in the corresponding subimages
• furthermore, there is a kind of a saturation where the computation time will not follow the above equation but stay more or less constant (see graph-> image sizes <5x10^5)
• last but not least when using a mask the mask will need to be also partitioned into the subimages.

Edit: new code and examples

Now to the code. First, we need a function that understands how far are the partitioned images from the base point, i.e. the origin {0,0}. The only thing which this function needs is the array which comes directly out of ImagePartition[]:

calcRelShiftCoord[partionedImg_List] := 
 Module[{partImgSizes, xVals, yVals},
  partImgSizes = Map[ImageDimensions, partionedImg, {2}];
  xVals = Accumulate[
Transpose[Prepend[Drop[partImgSizes[[-1]], -1], {0, 0}]][[1]]];
  yVals = Accumulate[
    Transpose[Prepend[Drop[partImgSizes[[;; , 1]], -1], {0, 0}]][[2]]];
  Table[{xVals[[j]], yVals[[i]]}, {i, Dimensions[partImgSizes][[2]], 
    1, -1}, {j, Dimensions[partImgSizes][[1]]}]
  ]

The main function takes both the large images, the number of subdivisions per image edge (e.g. if the number is 4 then the image would be subdivided into 16 images), and some options ImageCorrespondingPoints[] takes such as TransformationClass, Masking and MaxFeatures. I know there are other options but it is easy to include them if you need. The main computation is done in parallel to really benefit from Mma's capabilities to speed things up:

doDistributedParallelImageCorrespondingPoints[img1_, img2_, numOfLinearSubDivisions_, opts : OptionsPattern[ImageCorrespondingPoints]] := 
 Module[{maskingElement, trafoClass, numOfFeatures, img1Array, img2Array, maskArray, res, relShiftVectors, shiftedRelCoords, realCorrPnts},
  maskingElement = If[ImageQ[OptionValue[Masking]], OptionValue[Masking], Binarize[img1, 0]];
  trafoClass = OptionValue[TransformationClass];
  numOfFeatures = OptionValue[MaxFeatures];
  {img1Array, img2Array, maskArray} = ImagePartition[#, {Scaled[1/numOfLinearSubDivisions]}] & /@ {img1, img2, maskingElement};
  res = ParallelTable[ImageCorrespondingPoints[img1Array[[i, j]], img2Array[[i, j]], Masking -> maskArray[[i, j]], MaxFeatures -> numOfFeatures, 
          TransformationClass -> trafoClass], 
        {j, numOfLinearSubDivisions}, {i, numOfLinearSubDivisions}];
 relShiftVectors = calcRelShiftCoord[img1Array];
 shiftedRelCoords = Select[ParallelTable[
    If[Length[res[[i, j]][[1]]] != 0,
        Map[(relShiftVectors[[i, j]] + #) &, res[[i, j]], {2}], 
              Nothing], 
     {j, numOfLinearSubDivisions}, {i, numOfLinearSubDivisions}], Length[#] > 0 &];
  realCorrPnts = {Apply[Join[##] &, #[[1]] & /@ Flatten[shiftedRelCoords, 1]], 
                  Apply[Join[##] &, #[[2]] & /@ Flatten[shiftedRelCoords, 1]]};
 realCorrPnts
]

To see if it works all properly here are the images from the middle of the large image as well as the whole images:

and the overall image

calculated with the code

allPnts = doDistributedParallelImageCorrespondingPoints[bilder[[1]], bilder[[2]], 4, MaxFeatures -> 5];

The only thing which I am not totally satisfied with is that I can't calculate the keypoint strength, so that I keep only the strongest 15 points or so. Of course one could do an analysis with incremental decrease of a known key point strength and classify the found points by this, but as I wanted to speed up the code this doesn't make really sense. On the other side I can imagine some situations where a known keypoint strength helps to limit the errors and time for further calculations..

P.S.: The speed up is very efficient. With the mask the old method took on average 3h 11min to correct one image. Now it takes roughly 4.5 min.. :)