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.

  • 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
  • 1
    @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. – Niki Estner Aug 28 '13 at 16:57
up vote 13 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]]]
  • Nice use of Clip, with padding a default value could easily be achieved. – ssch Aug 29 '13 at 0:05

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]

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

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

img = Import[""]


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];
                                          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"]

jellyfish after remapping

Your Answer


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.