Take the 2-minute tour ×

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:

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.7s 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.7s for one color channel, so it's even slower for 3 channels.

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

share|improve this question
Very good question! Ages time ago I ended up generating scripts for mathmap with Mathematica to do similar things. Quite horribly, too... –  Yves Klett Aug 28 '13 at 12:37
Doesn't ImageTransformation do the same thing? "ImageTransformation[image,function] gives an image in which each pixel at position {x,y} corresponds to the position function[{x,y}] in image." –  Vitaliy Kaurov Aug 28 '13 at 16:40
@VitaliyKaurov: Basically, the function I'm looking for is to ImageTransformation what ListPlot is to Plot: Does the same thing, but with an array of data instead of a function to generate the data. –  nikie Aug 28 '13 at 16:57
add comment

2 Answers 2

up vote 6 down vote accepted

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]]]
share|improve this answer
Nice use of Clip, with padding a default value could easily be achieved. –  ssch Aug 29 '13 at 0:05
add comment

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}},
     nx, ny, nc,
     indx, indy
    {nx, ny, nc} = Dimensions[data];
    newdata = Table[0., {nx}, {ny}, {nc}];
     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}];

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

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

share|improve this answer
add comment

Your Answer


