# Using ListInterpolation to get a vector-valued function

I wanted to construct a ColorFunction from an image, so I tried

im = Image@
ReliefPlot[
Table[Sin[Abs@x - Abs@y], {x, -π, π, π/100}, {y, -π, π, π/100}],

cf = ListInterpolation[ImageData[im], {{0, 1}, {0, 1}}]


Of course, this doesn't work, because ImageData is an n * m * 3 list. However, I don't want to interpolate the third dimension, the resulting interpolating function should return a vector. I expected to find a relevant example in the Possible issues section, but alas, there is none.

The following approaches work:

cf = ListInterpolation[
Apply[RGBColor, ImageData@im, {2}], {{0, 1}, {0, 1}},
InterpolationOrder -> 0]


But then there's really no interpolation (interpolation order zero). A better approach:

cf = With[{f =
ListInterpolation[ImageData@im, {{0, 1}, {0, 1}, {1, 3}},
InterpolationOrder -> {3, 3, 1}]},
(RGBColor @@ {f[##, 1], f[##, 2], f[##, 3]}) &]


But I'm not sure if that's the best way.

So the general question is "how can one do a ListInterpolation of a list of vectors?" and a follow-up to it is "how can I make my approach cleaner and minimize any unnecessary overhead?"

For the general case (i.e. not building a ColorFunction, but simply interpolating an array of vectors) I prefer to not have to know anything about the structure of my list beforehand, except that it is a non-ragged n+1 dimensional array and I want a function of n arguments that returns a vector.

To make it absolutely clear, a vector-valued InterpolationFunction is quite possible, and indeed, the following is probably the ideal result for me:

With[{dim = (Most@Dimensions@ImageData@im - 1)},
Interpolation[
Flatten[MapIndexed[{(#2 - 1)/dim, #1} &, ImageData@im, {2}],
Length[dim] - 1]]]


What I don't like here, is that I'm forced to restructure the whole list and get the correct dimensions and all that, when ListInterpolation can do almost everything for me (except for vector values).

• Can you note in your question that Interpolation does support interpolating arbitrary dimensional vectors? Just to get this out of the way. – Szabolcs Feb 10 '16 at 11:46
• You could use Interpolation instead of ListInterpolation, but settling for three ListInterpolations might be easier actually. Personally I would go for three ListInterpolations. – Szabolcs Feb 10 '16 at 11:48
• @AlexeyPopkov what I meant is a 2d array of colors, between which I want to interpolate: colorarr = {{Red, Green, Blue}, {Yellow, Cyan, Magenta},{Gray, Black,White}};f = Function[{x, y}, Blend[Blend[#, x] & /@ colorarr, y]];RegionPlot[True, {x, 0, 1}, {y, 0, 1}, ColorFunction -> f] – LLlAMnYP Feb 22 '16 at 10:23
• @AlexeyPopkov data = colorarr = Map[RGBColor, ImageData[ExampleData[{"TestImage", "Lena"}]], {2}]; then define f as nested Blends as above. The interpolation is forced to re-evaluate at every call to f and is extremely slow. Interpolate the entire table of channel values first, as I show in my answer, and the RegionPlot is rendered almost instantly. – LLlAMnYP Feb 22 '16 at 10:34
• @AlexeyPopkov The problem is slightly different. Blend re-evaluates, because the Function has attributes HoldAll. The problem is, that Blend does not support a symboli second argument, that would allow to code that nested construct in a faster manner. I think, it would not be too hard to roll a myBlend that does support symbolic arguments. In fact, most of the work is already done in this QA. – LLlAMnYP Feb 22 '16 at 10:56

I came up with a pretty clean way of doing this. Interpolation and ListInterpolation support interpolation of symbolic values. I use this to hide the fact that the data is an n+1 dimensional array.

Module[{a},
Evaluate[Array[a, Most@Dimensions@data]] = data;
Block[{a}, ListInterpolation[Array[a, Most@Dimensions@data]]]]


Explanation:

• To avoid unnecessary garbage I localize a with Module.

• Then I place my data into the DownValues of a.

• At this point ListInterpolation[Array[a, Most@Dimensions@data]] would see an n+1-dimensional list with scalar output. However...

• ... if I Block a, then this is simply an n-dimensional list with symbolic output.

• Inside the block an InterpolationFunction is created, which is a normal Mathematica expression, with a[i,j] explicitly in the arguments.

• As soon as the Block is exited, the a[i,j] are replaced with the vector values and that rolls fine with InterpolationFunction.

EDIT

The code above works if the data range is not specified. Today I noticed, that if the data range is specified, the InterpolatingFunction returned by ListInterpolation inside the Block contains a[1., 1.]... (that is, real indices) instead of a[1, 1]... (integer indices).

This undesirable behavior is in itself worth another question.

Over there AlexeyPopkov suggests a clean workaround, by giving a the attribute NHoldAll:

Module[{a},
Evaluate[Array[a, Most@Dimensions@data]] = data;
a[idx__Real] := a @@ (Floor /@ {idx});
Block[{a}, SetAttributes[a, NHoldAll];
ListInterpolation[Array[a, Most@Dimensions@data], datarange]]]


EDIT 06.05.16

JasonB observes quite bad performance for this approach, and I'm, so far, not sure of the exact cause. But a workaround was quick enough to find: simply construct the appropriate arguments for InterpolatingFunction directly.

quickIntFunc[data_] :=
InterpolatingFunction[
{{1, 1}, Most@Dimensions@data} // Transpose,
{5, 3, 0, Most@Dimensions@data, {4, 4}, 0, 0, 0, 0, Automatic, {}, {}, False},
Range /@ (Most@Dimensions@data),
Map[List] /@ data,
{Automatic, Automatic}]

• For some reason this method is very slow sometimes, slower than building a list of interpolation functions. I tried to use this to simplify a code where the array dimensions are the rows and columns of an image and the vector is the RGB colors. Here is the example code, it takes a long time to build up the interpolating function using your method. – Jason B. May 6 '16 at 9:59
• @JasonB yes, interpolation of symbolic stuff takes a while, I'm seeing that. Subsequent replacement of a[x,y] with actual rgb values also takes a long time. And the interpolation of a list of {{{x1,y1},{r1,g1,b1}}...} is also atrociously long. I'll see if things could be sped up, as I don't see, why a vector-valued function absolutely must be 20 times slower to construct. – LLlAMnYP May 6 '16 at 10:27
• @JasonB see latest edit, this should suit your needs and performs 10x faster on my machine, (though a bit slower, than constructing 3 one-dimensional-output interpolations). – LLlAMnYP May 6 '16 at 10:55
• I will check it out, although I think I found a pretty decent solution to the problem, going to fill in the details on a new answer – Jason B. May 6 '16 at 10:56

This will build an interpolating function off of a n-dimensional array, that, when given n-1 arguments, returns a list whose length is the last value of Dimensions[array].

arrayInterpolation[array_] :=
Through@*(ListInterpolation[#] & /@
Transpose[array, RotateLeft@Range@(Depth@array - 1)]);


We can test it out,

testarray = RandomReal[1, {4, 4, 4, 4, 7}];
testfunc = arrayInterpolation[testarray];

testarray[[1, 2, 1, 2]]
testfunc[1, 2, 1, 2]
(* {0.862862, 0.566754, 0.311709, 0.734644, 0.982298, \
0.516959, 0.200609} *)
(* {0.862862, 0.566754, 0.311709, 0.734644, 0.982298, \
0.516959, 0.200609} *)


So it reproduces the data exactly, and it also interpolates:

ListPlot[{testarray[[1, 2, 1, 2]], testarray[[1, 2, 1, 3]],
testfunc[1, 2, 1, 2.5]},
PlotLegends -> {"vector 1", "vector 2", "interpolation"}]


• If someone knows a better way to make a single function that returns a list of values, when you have a list of functions, like a reverse of Map that'd be great. (Function[func, func[##]] /@ listoffunctions &) just feels hacky. – Jason B. Feb 10 '16 at 12:50
• I think, that's Through – LLlAMnYP Feb 10 '16 at 12:52
• Through@*cf1 instead of what you had in the last line of the module works just fine. Also fixed a typo - you had testarray in the module but localized array. – LLlAMnYP Feb 10 '16 at 12:55
• Thanks, also thanks for catching the typo – Jason B. Feb 10 '16 at 12:57
• no problem. Personally I prefer to use the composed form Through@*cf1, rather than nesting Function[Through[cf1[##]]], because Function has attributes HoldAll and every call to it will introduce extra evaluation, but I think, in this case there's no difference. +1 in any case. – LLlAMnYP Feb 10 '16 at 13:01