Recover data points from RGBColor of known ColorScheme

In general I want to recover the z-values from an image of a DensityPlot. A simple example might just be

Flatten[Table[{x, y, x Sin[x y]}, {x, 0, 2 \[Pi], \[Pi]/50},
{y, 0, 2 \[Pi], \[Pi]/50}], 1]
ListDensityPlot[%, PlotRange -> All, ColorFunction ->"RedBlueTones",
Frame -> False,ImageSize -> 800, PlotRangePadding -> None]
Export["img.jpg", %]


I can get all the color data just by

pixelData = ImageData[Import["img.jpg"], "Real32"]


but when I want to get a z-value according to the RGBColor of any individual pixel within "pixelData" I'm lost. I guess it's rather straight forward using something like the inverse of

ColorData["RedBlueTones",z]


I also don't care about any range specifications since I can rescale to the desired range provided with the bar legend of the image of interest.

-
I think the the Export[] function makes this irreversible. However, if you can figure the inverse function of your colorfunction, you could extract the z values from a new densityplot of the function z[x_,y_]=InverseColorFunction[imagedata[[x,y]]] –  Coolwater Apr 16 '14 at 9:27
Related: How to change ColorFunction after plotting - also Jens, but a different person... –  Jens Apr 17 '14 at 4:52

With the caveat that this has an answer up to scaling, one can write a function that creates a NearestFunction object that converts from RGBColor triplets to z-values at the unit interval:

ClearAll[colorSchemeToUnitZ];
colorSchemeToUnitZ[ColorScheme_String,
resolution_: 100] /; (MemberQ[
resolution < 150) :=
colorSchemeToUnitZ[ColorScheme] = Module[
{rules},
rules = Table[(ColorData[ColorScheme][i] /. RGBColor -> List) -> {i}, {i,0, 1, 1/N@resolution}];
Nearest[rules]
]


Essentially, every color scheme defines a 1D path in a 3D color space i.e.

With[{
pts1 = Table[(ColorData["RedBlueTones"][i] /. RGBColor -> List), {i,
0, 1, .01}],
pts2 = Table[(ColorData["GreenBrownTerrain"][i] /.
RGBColor -> List), {i, 0, 1, .01}]},
ListPointPlot3D[{pts1, pts2},
ColorFunction -> (RGBColor@## &),
AxesLabel -> {"R", "G", "B"},
PlotRange -> {0, 1} ]
]


and the NearestFunction maps to a parameter along the path every RGB value of the color scheme. It is exactly what you wan: an inverse for the color function corresponding to the color scheme:

invData =
First@Flatten[
Map[colorSchemeToUnitZ["RedBlueTones"],
pixelData, {2}], {{3}, {2}}];


(note, I've only used pixelData here). Now

ArrayPlot[Transpose@invData]


-

After tracing through ColorFunction, I found the gradients live in the list

DataPacletsColorDataDumpgradientSchemeMain


and "RedBlueTones" corresponds to the 14th position.

DataPacletsColorDataDumpgradientSchemeMain[[14, 5]]
(* {RGBColor[0.450385,0.157961,0.217975],RGBColor[0.599449,0.262748,0.294618],
RGBColor[0.721701,0.434448,0.400225],RGBColor[0.813151,0.617722,0.507726],
RGBColor[0.865768,0.767491,0.623596],RGBColor[0.857126,0.848339,0.734867],
RGBColor[0.771923,0.848195,0.811697],RGBColor[0.61971,0.781131,0.831119],
RGBColor[0.433786,0.670834,0.793785],RGBColor[0.256859,0.523007,0.711644],
RGBColor[0.139681,0.311666,0.550652]} *)


After some more investigating, it appears ColorFunction just interpolates linearly between these:

ColorData["RedBlueTones"][.02]
(* RGBColor[0.4801978,0.17891839999999998,0.2333036] *)

Blend[DataPacletsColorDataDumpgradientSchemeMain[[14, 5, 1;;2]], .2]
(* RGBColor[0.4801978,0.17891839999999998,0.2333036] *)


So we can use this to write a function that inverts a gradient ColorFunction:

InverseColorData[grad_String, color_RGBColor] := Module[{pos, knots, step, interpdata, interp, sol, x, y, z},
pos = Position[DataPacletsColorDataDumpgradientSchemeMain, grad][[1, 1]];
knots = DataPacletsColorDataDumpgradientSchemeMain[[pos, 5]];
step = 1/(Length[knots]-1);

interpdata = Transpose[{Range[0., 1, step], #}]& /@ Transpose[List @@@ knots];
interp = MapThread[LineBetween, {interpdata, {x, y, z}}];

sol = MapThread[#3 /. Solve[#1 == #2 && 0 <= #3 <= 1, #3]&, {interp, List @@ color, {x, y, z}}];

First[Intersection[##, SameTest -> (Chop[#1 - #2] == 0&)]& @@ sol]
]

LineBetween[lis_, x_] := Module[{par = Partition[lis, 2, 1], x1, y1, x2, y2},
Piecewise[
(
{{x1, y1}, {x2, y2}} = #;
{(y2-y1)/(x2-x1)*(x - x1) + y1, x1 <= x <= x2}
)& /@ par
]
]


and now we test:

With[{rand = RandomReal[{0, 1}, 10]},
Transpose[{
rand,
InverseColorData["RedBlueTones", ColorData["RedBlueTones"][#]]& /@ rand
}]
]
(* {{0.244535,0.244535}, {0.470216,0.470216}, {0.65504,0.65504}, {0.729637,0.729637},
{0.205673,0.205673}, {0.664396,0.664396}, {0.363139,0.363139}, {0.671589,0.671589},
{0.329273,0.329273}, {0.119555,0.119555}} *)


Bonus:

Since ColorFunction linearly interpolates, we can easily visualize the channels of "RedBlueTones":

With[{knots = DataPacletsColorDataDumpgradientSchemeMain[[14, 5]]},
Labeled[
ListLinePlot[
Transpose[{Range[0, 1, .1], #}]& /@ Transpose[List @@@ knots],
PlotStyle -> {Red, Blue, Green}
],
"RedBlueTones"
]
]


-
Nice +1. Just to make some extra work it doesn't work quite right with some other color schemes eg. ColorData["AvocadoColors"] (random choice.). I'm guessing it needs just a little tweak to make more robust. –  george2079 Apr 16 '14 at 17:21
The issue is that color scheme has a region of a constant channel (Red==0 for z<1/4), so Solve fails. –  george2079 Apr 16 '14 at 17:54
Thank you for spotting the issue. One could adapt via Reduce I suppose. –  Chip Hurst Apr 16 '14 at 19:12
Works like a charm (at least for this color scheme) as long as the RGB channels are given in high precission. For (re)imported images, however, the RGB tripel is rounded at some point and exact solutions cannot be found. In this case, the solution proposed by gpap produces results more reliably. –  Jens Apr 17 '14 at 12:36
Mind you, my solution also rounds it up to whatever resolution you pick. –  gpap Apr 17 '14 at 14:02

You could do this:

pixelData = ImageData[Import["img.jpg"], "Real32"];
fit[u_, v_, w_] = Normal[LinearModelFit[Table[Append[List @@ Blend["RedBlueTones",
#] &[t], t], {t, 0, 1, 1/1000}], Flatten[Table[u^i v^j w^k, {i, 0, 2}, {j, 0, 2},
{k, 0, 2}]], {u, v, w}]];
Plot[fit @@ Function[Blend["RedBlueTones", #]][t] - t, {t, 0, 1}](*Error of fit*)
getZ[x_, y_] := fit @@ pixelData[[x, y]]

-