Inverting a real-valued index grid

Question

OpenCV's remap() uses a real-valued index grid to sample a grid of values from an image using bilinear interpolation, and returns the grid of samples as a new image.

To be precise, let:

A = an image  X = a grid of real-valued X coords into the image.  Y = a grid of real-valued Y coords into the image. B = remap(A, X, Y) 

Then for all pixel coordinates i, j,

B[i, j] = A(X[i, j], Y[i, j])  

Where the round-braces notation A(x, y) denotes using bilinear interpolation to solve for the pixel value of image A using float-valued coords x and y.

My question is: given an index grid X, Y, how can I generate an "inverse grid" X^-1, Y^-1 such that:

X(X^-1[i, j], Y^-1[i, j]) = i Y(X^-1[i, j], Y^-1[i, j]) = j 

And

X^-1(X[i, j], Y[i, j]) = i Y^-1(X[i, j], Y[i, j]) = j 

For all integer pixel coordinates i, j?

FWIW, the image and index maps X and Y are the same shape. However, there is no a priori structure to the index maps X and Y. For example, they're not necessarily affine or rigid transforms. They may even be uninvertible, e.g. if X, Y maps multiple pixels in A to the same exact pixel coordinate in B. I'm looking for ideas for a method that will find a reasonable inverse map if one exists.

The solution need not be OpenCV-based, as I'm not using OpenCV, but another library that has a remap() implementation. While any suggestions are welcome, I'm particularly keen on something that's "mathematically correct", i.e. if my map M is perfectly invertible, the method should find the perfect inverse, within some small margin of machine precision.

Answer 1

Well I just had to solve this remap inversion problem myself and I'll outline my solution.

Given X, Y for the remap() function that does the following:

B[i, j] = A(X[i, j], Y[i, j])    

I computed Xinv, Yinv that can be used by the remap() function to invert the process:

A[x, y] = B(Xinv[x,y],Yinv[x,y]) 

First I build a KD-Tree for the 2D point set {(X[i,j],Y[i,j]} so I can efficiently find the N nearest neighbors to a given point (x,y). I use Euclidian distance for my distance metric. I found a great C++ header lib for KD-Trees on GitHub.

Then I loop thru all the (x,y) values in A's grid and find the N = 5 nearest neighbors {(X[i_k,j_k],Y[i_k,j_k]) | k = 0 .. N-1} in my point set.

• If distance d_k == 0 for some k then Xinv[x,y] = i_k and Yinv[x,y] = j_k, otherwise...

• Use Inverse Distance Weighting (IDW) to compute an interpolated value:

  • let weight w_k = 1 / pow(d_k, p) (I use p = 2)
  • Xinv[x,y] = (sum_k w_k * i_k)/(sum_k w_k)
  • Yinv[x,y] = (sum_k w_k * j_k)/(sum_k w_k)

Note that if B is a W x H image then X and Y are W x H arrays of floats. If A is a w x h image then Xinv and Yinv are w x h arrays for floats. It is important that you are consistent with image and map sizing.

Works like a charm! My first version I tried brute forcing the search and I never even waited for it to finish. I switched to a KD-Tree then I started to get reasonable run times. I f I ever get time I would like to add this to OpenCV.

The second image below is use remap() to remove the lens distortion from the first image. The third image is a result of inverting the process.