# Speeding up ImageCorrespondingPoints

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

• For a start, you could try to downsample the images just for the detection of corresponding points. Moreover, I saw that ImageCorrespondingPoints has also an option TransformationClass. Seemingly, you do something redudantly here... – Henrik Schumacher Nov 18 '17 at 16:22
• Thanks, I missed that with the transformation class. I'll test whether it could speed up the whole process. However, I think simple downsampling is not the way to go, because to fit the images perfectly to each other I need the key points at the final resolution. Of course, I could adaptively crop the image to the area around the coursely found keypoint and then fine-adjust it at a higher resolution. But (1) I'm not convinced the speed will improve significantly and (2) using downsampled images prevents Mma from finding fine structures in those cases where the image contains only those. – Quit007 Nov 18 '17 at 19:00
• ImageCorrespondingPoints by default seems to use ImageKeyPoints and "RANSAC" form of FindGeometricTransform. Both are apt to be slow. In any case, are there services or libraries out there that can do this faster on images of your size? If not, I'm not sure there's really any reason this procedure should be fast. I think Mathematica's doing as much as is possible at the C level anyway. – b3m2a1 Nov 19 '17 at 19:47

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]] :=
trafoClass = OptionValue[TransformationClass];
numOfFeatures = OptionValue[MaxFeatures];
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.. :)

• I'm curious if you've also tried this with specialized image-processing software (e.g., align_image_stack in Hugin, or autoalign/merge in Photoshop) and, if so, how the computation time and results compare with what you've obtained with MMA. – theorist Aug 4 at 5:37
• @theorist At that time, I tried to find some software which would allow to do this efficiently and precisely. I don't remember exactly which software packages I tested. But I found them to be less precise in comparison with Mma. Because, if done carefully, Mma will give you a whole lot more control than "standard" software. Maybe there is some specialized vid software out there, but I lost interest with the speed up here. (here is one example: youtu.be/5cSXkuymV3Y). I believe YouTube still offers deshaking of vids. Depending on your problem the result might be ok. – Quit007 Aug 5 at 22:42