How can I efficiently "remap" an image?

Question

Many image processing libraries like OpenCV, Intel Performance Primitives or Octave have a useful function called "remap", that takes an image, an array with X coordinates and an array with Y coordinates, and returns an image that "transforms" the image by that geometric mapping. Or, as the IPP documentation puts it: "Pixel remapping is performed using pxMap and pyMap buffers to look-up the coordinates of the source image pixel that is written to the target destination image pixel:

dst_pixel[i, j] = src_pixel[pxMap[i, j], pyMap[i, j]]

The closest thing in Mathematica that I'm aware of is ImageTransformation, but that takes a function, not an array. Using arrays can often save a lot of time (for example, you can apply the same mapping-array to multiple image, and arithmetic operations on arrays are very fast).

The best I've come up with so far is to combine ImageTransformation and ListInterpolation to convert the arrays to functions:

img = ExampleData[{"TestImage", "Lena"}];
mapX = Table[i + j, {i, 500}, {j, 500}];
mapY = Table[i - j, {i, 500}, {j, 500}];

(Imagine some time-consuming operation here. I know that passing i + j, i - j as a function to ImageTransformation is probably much faster in this simple case, but that's not the point.)

{xFn, yFn} = 
  ListInterpolation[#, {{0, 1}, {0, 1}}, InterpolationOrder -> 1] & /@ {mapX, mapY};

Timing[ImageTransformation[img, {xFn @@ #, yFn @@ #} &, {500, 500}, 
  PlotRange -> {{0, 1}, {0, 1}}, DataRange -> Full]]

This works, but is extremely slow (4.7 s on my PC).

The second idea I had was to use ListInterpolation on the image directly (in this case, only on the red channel):

redFn = ListInterpolation[ImageData[img][[All, All, 1]], 
   InterpolationOrder -> 1];    
Timing[Image[redFn[mapY, mapX]]]

This takes 1.7 s for one color channel, so it's even slower for 3 channels.

For comparison: The IPP's remap function usually takes a few milliseconds.

Answer 1

Perhaps something like this may do. Convert the mapX and mapY arrays to a 1D list of index values corresponding to position in the flattened image data. A 1D list of positions can be used very quickly with Part, and the 2D image reconstructed using Partition

imagemap[img_, mapX_, mapY_] := Module[{a, b, mx, my, id, pix},
  {a, b} = ImageDimensions[img];
  mx = Clip[mapX, {1, a}];
  my = Clip[mapY, {1, b}];
  id = Flatten[ImageData[img], 1];
  pix = Flatten[Transpose[(my - 1) a + mx]];
  Image@Partition[id[[pix]], Length[mx]]]

Answer 2

Straightforward method without interpolation, takes 0.1s but it feels like there is lots of room for improvement. I used the definition with dst_pixel so the orientation is different from result in question.

I don't know of a quick way to have Extract or Part return a default value for indices that are out of range so I ended up using Compile instead:

cImageRemap = 
  Compile[{{data, _Real, 3}, {pxMap, _Integer, 2}, {pyMap, _Integer, 2}},
   Block[{
     nx, ny, nc,
     newdata,
     indx, indy
     },
    {nx, ny, nc} = Dimensions[data];
    newdata = Table[0., {nx}, {ny}, {nc}];
    Do[
     indx = pxMap[[ix, iy]];
     indy = pyMap[[ix, iy]];
     If[1 <= indx <= nx && 1 <= indy <= ny,
      newdata[[ix, iy]] = data[[indx, indy]]
      ];
     , {ix, nx}, {iy, ny}];
    newdata
    ]];

ImageRemap[img_Image, pxMap_?(MatrixQ[#, IntegerQ] &), pyMap_?(MatrixQ[#, IntegerQ] &)] /; 
  Dimensions[pxMap] == Dimensions[pyMap] == ImageDimensions[img] :=
 Image@cImageRemap[ImageData@img, pxMap, pyMap]

img = ExampleData[{"TestImage", "Lena"}];
mapX = Table[i + j, {i, #1}, {j, #2}] & @@ ImageDimensions[img];
mapY = Table[i - j, {i, #1}, {j, #2}] & @@ ImageDimensions[img];
Timing[ImageRemap[img, mapX, mapY];]
(* 0.1s *)

Edit: With CompilationTarget->"C" it takes 0.03s

Answer 3

First, here's an image from the docs which we'll use for testing:

img = Import["https://i.stack.imgur.com/bzkJM.png"]

Going with the definition given in the IPP docs, here is a remapping method based on the use of ImageValue[]:

Options[ImageRemap] = {Padding -> 0, Resampling -> "Bilinear"}; 
ImageRemap[img_Image, xm_?MatrixQ, ym_?MatrixQ, opts : OptionsPattern[]] /; 
Reverse[ImageDimensions[img]] == Dimensions[xm] == Dimensions[ym] :=
        Module[{h, w}, {w, h} = ImageDimensions[img];
               Image[Partition[ImageValue[img,
                                          Flatten[Transpose[{ym, h - xm} - 1/2,
                                                            {3, 1, 2}], 1], 
                                          DataRange -> Full, opts,
                                          Options[ImageRemap]], w]]]

I don't know if there is a better way of doing this, but this implementation at least inherits the Padding and Resampling options of ImageValue[], which may be needed for some mapping applications.

Let's try it out:

{w, h} = ImageDimensions[img];
mx = Array[Plus, {h, w}]; my = Array[BitXor, {h, w}];
ImageRemap[img, mx, my, Padding -> GrayLevel[3/4], Resampling -> "Lanczos"]