33
$\begingroup$

Since ColorData["VisibleSpectrum"] is wrong, I would like to have a more accurate function to use.
Can this information be extracted from Mathematica itself?

$\endgroup$
35
$\begingroup$

Notice:

Simon Woods did just this months ago for an answer I missed:


It seems that it can. By spelunking ChromaticityPlot I found:

Image`ColorOperationsDump`$wavelengths

Image`ColorOperationsDump`tris

These are a list of wavelengths and their corresponding XYZ color values used by this plot command:

ChromaticityPlot["sRGB", Appearance -> {"VisibleSpectrum", "Wavelengths" -> True}]

enter image description here

We can therefore use them to generate a new color function:

ChromaticityPlot;  (* pre-load internals *)

newVisibleSpectrum =
  With[
    {colors =
      {Image`ColorOperationsDump`$wavelengths,
       XYZColor @@@ Image`ColorOperationsDump`tris}\[Transpose]},
    Blend[colors, #] &
  ];

A comparison with the old function:

ArrayPlot[
 {Range[385, 745]},
 ImageSize -> 400,
 AspectRatio -> 0.1,
 ColorFunctionScaling -> False,
 ColorFunction -> #
] & /@
  {"VisibleSpectrum", newVisibleSpectrum} // Column

enter image description here

589nm is now the bright sodium yellow that it should be:

Graphics[{newVisibleSpectrum @ 589, Disk[]}]

enter image description here

If you wish to integrate this into ColorData see:


As requested by J.M. red-green-blue plots for each function:

old = List @@@ Array[ColorData["VisibleSpectrum"], 361, 385];

new = List @@@ ColorConvert[Array[newVisibleSpectrum, 361, 385], "RGB"];

ListLinePlot[Transpose @ #,
   PlotStyle -> {Red, Green, Blue},
   DataRange -> {385, 745}
] & /@ {old, new}

enter image description here

Clipping occurs during conversion to screen RGB; the newVisibleSpectrum function actually produces unclipped XYZColor data. For example projected over gray:

newVSgray = 
  With[{colors = 
     Array[{#, Blend[{newVisibleSpectrum@#, ColorConvert[GrayLevel[.66], "XYZ"]}, 
         0.715]} &, 361, 385]}, Blend[colors, #] &];

ListLinePlot[
 List @@@ ColorConvert[Array[newVSgray, 361, 385], "RGB"] // Transpose, 
 PlotStyle -> {Red, Green, Blue}, DataRange -> {385, 745}, ImageSize -> 400]

enter image description here

Which corresponds to the plot:

ArrayPlot[{Range[385, 745]}, ImageSize -> 400, AspectRatio -> 0.1, 
 ColorFunctionScaling -> False, ColorFunction -> newVSgray, 
 Background -> GrayLevel[0.567]]

enter image description here

cf. "VisibleSpectrum" similarly over gray blended in XYZColor and RGBColor respectively:

enter image description here
enter image description here

Note that neither rendering of this spectrum has the vibrancy of newVisibleSpectrum.

$\endgroup$
  • 2
    $\begingroup$ It seems the internal functions are only available after you call ChromaticityPlot once. Is that normal? $\endgroup$ – Rahul Feb 3 '15 at 13:40
  • $\begingroup$ @Rahul Yes, it's common with many internal functions and I should have thought to mention it. $\endgroup$ – Mr.Wizard Feb 3 '15 at 14:12
  • 1
    $\begingroup$ +1 Nice. Though spelunking ChromaticityPlot is so last year :-) $\endgroup$ – Simon Woods Feb 3 '15 at 14:27
  • $\begingroup$ @Simon Gah! How'd I miss that? $\endgroup$ – Mr.Wizard Feb 3 '15 at 14:37
  • 1
    $\begingroup$ That gradient is very bad, and that makes me dubious about the black body one as well. I'll try to see where/how they are implemented and fix them. I should also read here more often... (btw if the color matching functions are such coveted material I may put them in a decent context or push for them to be in System altogether) $\endgroup$ – Batracos Aug 2 '15 at 9:19
7
$\begingroup$

(with many thanks to halirutan and kirma for their kind assistance)

Here's a different take. In this article, piecewise Gaussian functions that approximate the CIE color matching functions are presented. For this answer, instead of just taking the coefficients from the paper directly, I used their proposed model in FindFit[], using a 1 nm tabulation of the 1931 CMFs from here as the data the proposed model is being fitted to. Here are the resulting functions:

SetAttributes[{cieX, cieY, cieZ}, Listable];
cieX[λ_] := {0.36263412, 1.05593554, -0.2122116}.Exp[
            -MapThread[(λ - #1) Piecewise[{{#2, λ < #1}, {#3, True}}] &,
                       {{440.76797, 599.44051, 493.87482},
                        {0.066588512, 0.027633445, 0.048888466},
                        {0.020005468, 0.032068622, 0.039095442}}]^2/2]

cieY[λ_] := {0.82087906, 0.28579841}.Exp[
            -MapThread[(λ - #1) Piecewise[{{#2, λ < #1}, {#3, True}}] &,
                       {{568.78966, 530.87379},
                        {0.021328878, 0.061294726},
                        {0.024693511, 0.032178136}}]^2/2]

cieZ[λ_] := {1.21649968, 0.68126275}.Exp[
            -MapThread[(λ - #1) Piecewise[{{#2, λ < #1}, {#3, True}}] &,
                       {{436.96247, 459.03433},
                        {0.084453906, 0.038519512},
                        {0.027788061, 0.072502761}}]^2/2]

In version 10, one can of course directly do Through[XYZColor[cieX, cieY, cieZ][λ]]; for older versions, you will have to do the conversion to sRGB yourself, like so:

(* gamma correction *)
sRGBGamma = Function[x, With[{z = Abs[x]},
                             Sign[x] Piecewise[{{12.92 z, z <= 0.0031308}},
                                               1.055 z^(1/2.4) - 0.055]],
                     Listable];

myVisibleSpectrum[λ_] :=
     RGBColor @@ Clip[sRGBGamma[{{3.2404542, -1.5371385, -0.49853141},
                                 {-0.96926603, 1.8760108, 0.041556017},
                                 {0.055643431, -0.20402591, 1.0572252}} .
                                Through[{cieX, cieY, cieZ}[λ]]], {0, 1}]

Here are the corresponding gradient and component plots for myVisibleSpectrum[]: myVisibleSpectrum demonstration

Here's the result of coloring a Disk[] with myVisibleSpectrum[589]: sodium yellow?

The approximation looks good, at least with my eyes.

This replicates the internal data (as used in newVisibleSpectrum) quite well:

cieXYZ = Through[XYZColor[cieX, cieY, cieZ][#]] &;

{d1, d2} = Transpose /@ Table[List @@ fn[λ], {fn, {cieXYZ, newVisibleSpectrum}},
                              {λ, 385, 745}];

Show[ListLinePlot[d1, PlotStyle -> Directive[Thickness[0.01], Hue[0.55, 0.5, 1]],
                  PlotTheme -> "Pastel", InterpolationOrder -> 3],
     ListLinePlot[d2, PlotStyle -> Directive[Black, Dashed]]]

comparison of tabulated CMF with analytic fit


Alternatives

I had already discussed Bruton's algorithm in my other answer. I'm not sure how to convert wavelengths to the colors used in Sam's answer, but here is how to see them in Mathematica:

DensityPlot[x, {x, 0, 1}, {y, 0, 1/8}, AspectRatio -> Automatic, 
            ColorFunction -> Function[t, RGBColor @@ Haversine[2 π t - {π, 5 π/3, π/3}]]]

Sam's color gradient

$\endgroup$
  • $\begingroup$ Thanks for answering. It seems you've replicated the internal data quite well. I'll assume your gedanken mma isn't up to direct .PNG output yet so I included a plot to show this. $\endgroup$ – Mr.Wizard Jun 27 '15 at 2:56
  • $\begingroup$ No, I'm using a borrowed computer for now. :) (though it only has version 8). Thanks for putting in the plot! $\endgroup$ – J. M. will be back soon Jun 27 '15 at 3:14
5
$\begingroup$

Here's my mathematical rainbow spectrum plot.

I'm not sure how accurate it is, and I did not attempt to calibrate it with measured data; but it's the smoothest, best looking spectrum I've seen, so I thought I would post it here.

rainbow spectrum

Large images: http://sam.nipl.net/pix/spectrum.png http://sam.nipl.net/pix/spectrum2.png

Reasoning:

  • All colours in the rainbow spectrum should be shown with equal intensity.
  • In RGB, yellow(1,1,0) has greater intensity than red(1,0,0) or green(0,1,0).
  • So, we need a less intense yellow for our rainbow. I used cosines.
  • The rainbow is an octave of light, from 400THz-800Thz, or 750nm-375nm.
  • We can see 12 rainbow colours, like the 12 musical notes.
  • The refractive index varies with the wavelength of light.
  • The rainbow is a uniform plot over wavelength (not freq. or "note").
  • Several "rainbow echoes" are seen, due to resonance with higher octaves.

Here is a picture of rainbow echoes, these are due I suppose to harmonic resonance:

rainbow with echoes

Here is the code, it's not in Mathematica (sorry!):

#!/usr/local/bin/cz --

use b

Main()
    space(1900,100)
    origin(w_2, h_2)
    rainbow_init()

    for(x, 0.0, w)
        rainbow_colour(x/w)
        line(x, 0, x, h)
    paint()

    if args && cstr_eq(arg[0], "-s"):
        plot_12_tones()

num rb_red_angle, rb_green_angle, rb_blue_angle
num rb_red_power, rb_green_power, rb_blue_power

rainbow_init()
    rb_red_angle = deg2rad(0)
    rb_green_angle = deg2rad(120)
    rb_blue_angle = deg2rad(-120)

    rb_red_power = 1
    rb_green_power = 1
    rb_blue_power = 1

colour rainbow_colour(num h)
    num y = pow(2, h)-1
    num a = 2*pi*y
    num r = rb_red_power * (cos((a-rb_red_angle))+1)/2
    num g = rb_green_power * (cos((a-rb_green_angle))+1)/2
    num b = rb_blue_power * (cos((a-rb_blue_angle))+1)/2
    return rgb(r, g, b)

plot_12_tones()
    black()
    for(i, 0, 12+1):
        num freq = pow(2, i/12.0)
        num wl = 2 - 1/freq * 2
        num x = wl*w
        x = iclamp(x, 0, w-1)
        line(x, 0, x, 2)
        line(x, h, x, h-3)
$\endgroup$
  • 2
    $\begingroup$ Sam, welcome to Mathematica Stack Exchange. I appreciate you posting this despite what I am about to say: if you do not have Mathematica an/or are unfamiliar with its language it will be difficult for you to meaningfully participate here. If you choose to participate you should focus on clear English descriptions of algorithms or approaches rather than code in another language. $\endgroup$ – Mr.Wizard Feb 5 '15 at 17:29
  • 2
    $\begingroup$ Subjectively the yellow in your photo looks more intense than the one in your smooth spectrum; I take it your intent is not to replicate the photo but instead to produce a pleasingly smooth gradient, correct? $\endgroup$ – Mr.Wizard Feb 5 '15 at 17:30
  • 1
    $\begingroup$ Mr.Wizard, yes that's correct. Perhaps I'll have a go at translating it into Mathematica... $\endgroup$ – Sam Watkins Feb 5 '15 at 17:43
  • 2
    $\begingroup$ The "echoes" are called supernumerary bows. en.wikipedia.org/wiki/Rainbow#Supernumerary_rainbow $\endgroup$ – evanb Jun 17 '15 at 13:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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