I'm currently fighting with the spectrum->RGB color conversions, my algorithm seems to make an error somewhere, e.g. I get values >1 for some spectral responses.
Now there is a way to calculate an rgb value for normal incidence from index of refraction:
F(0°) = $(\frac{n-1}{n+1})^2$
Where n is the index of refraction (Real-Time Rendering 3rd, p.234 f.) And if $n$ is constant, you can use it for all RGB channels.
However, if I calculate the spectral value with my spectrum->RGB conversion and just use a constant $n$, I get different values for R, G and B.
Specifically, $n = 1.52$ would lead to $RGB_{single} = (0.04257999496;0.04257999496;0.04257999496)$ whereas my program calculates $RGB_{spectrum} = 0.054223;0.042692;0.041095$ (evaluated at $390nm - 780nm$ with $\Delta\lambda = 10nm$).
So my question now is: would a constant index of refraction really lead to a grayscale (i.e. all channels the same) value? Because the conversion functions depend on the CIE functions, taking the integral of those and clearly the integrals are different for $x, y, z$.
Also: since I seem to have a problem with my program in the first place: is there a step by step tutorial for creating the conversion function somewhere?
Edit 01.09.2017 10:25 After being asked to share my code, here goes (in a somewhat top-down approach): First I read the CIE functions from file and store them in a std::vector> with wavelength as key and function value as value. the vector just encapsulates x, y and z independently.
std::vector<std::map<float, float>> cieFunctions = readCieFunctions(ciePath);
Now cieFunctions.at(0) is the x function and so forth. I read in the spectral values basically the same way, into a map with key being the wavelength and value being the reflectivity. Here I take care to match the units, wavelength as $nm$ and reflectivity in values $\in [0, 1]$ (rather than $\in [0, 100]$%).
Having reached this, I call my calculateRgbValue(std::map<float, float>spectralValue) function, where I use the matrix for a D65 white point monitor. I calculate the coordinate, divide it by the integral of Y over the visible wavelength spectrum (constant taken from Physically Based Rendering book) and use that in a multiplcation with said matrix to get the RGB value:
QString SpectralToRGBConverter::calculateRgbValue(std::map<float, float> spectralValues)
{
QString result;
Vec3 rgb;
Mat3 conversionMatrix(
3.240479, -1.537150, -0.498535,
-0.969256, 1.875992, 0.041556,
0.055648, -0.204043, 1.057311);
float X = calcCieCoordinate(spectralValues, cieFunctions.at(0));
float Y = calcCieCoordinate(spectralValues, cieFunctions.at(1));
float Z = calcCieCoordinate(spectralValues, cieFunctions.at(2));
float yIntegral = 106.856895; //from Phyiscally Based Rendering 3rd, p. 325
Vec3 XYZ(X, Y, Z);
XYZ /= yIntegral;
rgb = conversionMatrix*XYZ;
result = (std::to_string(rgb.x) + ";" + std::to_string(rgb.y) + ";" + std::to_string(rgb.z)).c_str();
return result;
}
The idea of calculating the cie coordinate is best described with the mathematical formula: $coord = \sum_{i=0}^{n-1}x_i c_i \Delta\lambda$ where $n$ is the number of samples in the spectrum, $x_i$ is the $i-th$ cie function value for the wavelength of the $i-th$ sample and $c_i$ is the $i-th$ function value of the spectrum and $\Delta\lambda$ is the stride between two samples $\Delta\lambda = \lambda_{i+1} - \lambda_i$
float SpectralToRGBConverter::calcCieCoordinate(std::map<float, float> spectralValues, std::map<float, float> function)
{
//http://www.brucelindbloom.com/index.html?Eqn_Spect_to_XYZ.html
//Coordinate = Int_lambda function(lambda) * spectralValue(lambda) dlambda
//for discrete values:
//Coordinate = Sum_i function(i) * spectalValue(i) * Delta lambda
//where Delta lambda is the spacing between two measurements
float deltaLambda;// in nm
auto first = spectralValues.begin();
auto second = std::next(first, 1);
if (ui.radioButtonMicroMeter->isChecked())
{
deltaLambda = (second->first - first->first) * 1000;
}
else
{
deltaLambda = second->first - first->first;
}
float coordinate = 0;
for (auto pairWavelengthReflectance : spectralValues)
{
coordinate += getFunctionValue(function, pairWavelengthReflectance.first) * pairWavelengthReflectance.second * deltaLambda;
}
return coordinate;
}
getFunctionValue is basically a nearest neighbor sampling:
float SpectralToRGBConverter::getFunctionValue(std::map<float, float> function, float x)
{
float conversionFactorWaveLength;
if (ui.radioButtonMicroMeter->isChecked())
{
//if x is not in nm like function.first
conversionFactorWaveLength = 1000;
}
else
{
conversionFactorWaveLength = 1;
}
float conversionFactorReflectivityScale;
if (ui.radioButtonPercentage->isChecked())
{
//if the reflectivity is given in percentages [0, 100] rather than proportions [0, 1]
conversionFactorReflectivityScale = 0.01;
}
else
{
conversionFactorReflectivityScale = 1;
}
float xLookUp = x * conversionFactorWaveLength;
float previous = 0;
for (auto pairXY : function)
{ //maps are sorted
float current = pairXY.first;
if (current == xLookUp)
{
return function[current] * conversionFactorReflectivityScale;
}
else
if (current > xLookUp)
{
//get the closest value to x
float min;
if (xLookUp - previous < current - xLookUp)
{
min = previous;
}
else
{
min = current;
}
return function[min] * conversionFactorReflectivityScale;
}
previous = current;
}
return -1;
}
