Ok, here's the remainder of the algorithm to get to nanometers.
I couldn't find an implementation of this ANYWHERE, so I winged it.
The goal is to get the wavelength of the dominant color in the spectral data. You've got the color now (x, y) and the white point (wX, wY) so you draw a line starting at the data color (x, y) and passing through the white point (wX, wY) and determine where it intersects the boundary of the so-called "spectral locus" or the space where the colors exist.
Turns out its actually easy to draw the boundary line of the spectral locus by just iterating over the wavelengths we care about and using the same color space functions as before, except using 1.0 as the transmission value. This time, rather than pooping out the exact coordinate of the color that was matched, it poops out the edge point of the CIE diagram along which the wavelengths are labelled in nanometers for reference (see CIE 1931 diagram for example).
So my strategy was to just intersect a line segment with the current and previous point along the curve as I iterated through the wavelengths. As soon as something intersected, I had my wavelength;
NOTE: precision could be improved greatly as before by interpolating between wavelengths based on the exact intersection point's distance from the start of the wavelength range (ie. 500 and 501). I also used a 2d line segment vs. line segment function that I found on the interwebs (and cited) to perform the test. I believe the accuracy once I fix a few things should be plenty for our application once it's properly dialed in.
Seems to work. Haven't verified results yet, but gave me a reasonably yellow color for Kaolik so I'm posting it as-is.
I'll post everything on a git hub soon so that hopefully others can improve it and I can con someone in color science to check my math.
Here's the code tack this onto the previous code I posted, where I got x2d and y2d and the white point:
double diffX = x2d - xWhitePt;
double diffY = y2d - yWhitePt;
double normX = diffX / sqrt(diffX * diffX + diffY * diffY);
double normY = diffY / sqrt(diffX * diffX + diffY * diffY);
// create line segment to collide against spectral locus below
double segmentLength = 2.0;
double segX = x2d + normX * segmentLength;
double segY = y2d + normY * segmentLength;
// Iterate over spectral locus and collide with line segment S
// as given by [(x2d, y2d), (segX, segY)]
Point2d segStart = { x2d, y2d };
Point2d segEnd = { segX, segY };
double lambdaX = 0.0;
double lambdaY = 0.0;
double lambdaZ = 0.0;
double lambdaN = 0.0;
double dominantWavelength = 0.0;
Point2d pointOnLocus = { 0.0, 0.0 };
Point2d lastPoint = { 0.0, 0.0 };
for (double wave = 360.0; wave <= 830.0; wave += 1.0)
{
double transmittance = 1.0;
double refIllum = refIlluminant(wave);
CMFData matching = cmf(wave);
lambdaX = matching.m_xBar * transmittance * refIllum;
lambdaY = matching.m_yBar * transmittance * refIllum;
lambdaZ = matching.m_zBar * transmittance * refIllum;
lambdaN = matching.m_yBar * refIllum;
double locusN = sumN;
// Tristimulus values
double locusX = (1.0 / locusN) * lambdaX;
double locusY = (1.0 / locusN) * lambdaY;
double locusZ = (1.0 / locusN) * lambdaZ;
double locus_x2d = locusX / (locusX + locusY + locusZ);
double locus_y2d = locusY / (locusX + locusY + locusZ);
double locus_z = locusZ / (locusX + locusY + locusZ);
Point2d point;
point.m_x = locus_x2d;
point.m_y = locus_y2d;
// Intersect [lastPoint, point] with [segStart, segEnd]
if (lineSegmentIntersection(
lastPoint.m_x, lastPoint.m_y,
point.m_x, point.m_x,
segStart.m_x, segStart.m_y,
segEnd.m_x, segEnd.m_y,
&pointOnLocus.m_x, &pointOnLocus.m_y))
{
dominantWavelength = wave;
}
lastPoint = point;
}
if (dominantWavelength > 0.0)
{
char s[1024];
sprintf_s(s, 1024, "DomWavelen: %.2lf nm, Locus Coord (%.4lf, %.4lf)\n", dominantWavelength, pointOnLocus.m_x, pointOnLocus.m_y);
OutputDebugStringA(s);
}
