Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

Any ideas? Thanks in advance.

share|improve this question
    
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 at 9:27
    
Related: How to change ColorFunction after plotting - also Jens, but a different person... –  Jens Apr 17 at 4:52

3 Answers 3

up vote 1 down vote accepted

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[
     ColorData["Physical"]~Join~ColorData["Gradients"], ColorScheme] &&
     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} ]
 ]

enter image description here

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]

enter image description here

share|improve this answer

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

DataPaclets`ColorDataDump`gradientSchemeMain

and "RedBlueTones" corresponds to the 14th position.

DataPaclets`ColorDataDump`gradientSchemeMain[[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[DataPaclets`ColorDataDump`gradientSchemeMain[[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[DataPaclets`ColorDataDump`gradientSchemeMain, grad][[1, 1]];
  knots = DataPaclets`ColorDataDump`gradientSchemeMain[[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 = DataPaclets`ColorDataDump`gradientSchemeMain[[14, 5]]},
  Labeled[
    ListLinePlot[
      Transpose[{Range[0, 1, .1], #}]& /@ Transpose[List @@@ knots],
      PlotStyle -> {Red, Blue, Green}
    ],
    "RedBlueTones"
  ]
]

enter image description here

share|improve this answer
2  
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 at 17:21
1  
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 at 17:54
    
Thank you for spotting the issue. One could adapt via Reduce I suppose. –  Chip Hurst Apr 16 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 at 12:36
    
Mind you, my solution also rounds it up to whatever resolution you pick. –  gpap Apr 17 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]]
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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