I found my error! Two errors, in fact. What I get now does indeed look like it should. Conversion from spectrum to color seems to work too, pretty much. The only thing I'm still not entirely sure about is how to properly normalize the color I get. If I simply divide by the number of samples taken I get some slightly colored grey rather than white. From what I gathered I could also normalize it such that the Y value ends up at 1 which, if I understood right, means the color I get will be at the maximum brightness of the color system (though I might still get some hue since it isn't "perfect white")
Perhaps this is usable as is to you either way. It's fairly straight forward code. If you need me to explain anything, just ask away. I'd love to see full-spectrum Koch gears.
from math import exp, floor, ceil
import random
import numpy as np
from matplotlib import pyplot as plt
import csv
# a bunch of utility functions to get to proper units
def to_eV(v):
return v * 4.135668
def to_nm(v):
return 299.8/v
def from_nm(v):
return 299.8/v
def from_K(T):
return T * 1.e-4
def from_THz(v):
return v * 1.e-3
C = 77.766 # prefactor such that peak is 1
K = 4.79924 # prefactor if units are in PHz and 10^4 K
# desired range for v is .4 to .789
# desired range for T is 7.98e-2 to 2.
# typical range would be .5 to .7, commonly .65
def Planck(v, T): # the distribution of spectral frequencies for a given temperature
vT = v / T
return C * vT * vT * vT / (exp(K * vT) - 1)
topPHz = from_nm(390)
botPHz = from_nm(830)
def MC(T): # generates samples from the distribution using Markov Chain Monte Carlo method
v0 = random.uniform(botPHz, topPHz)
I0 = Planck(v0, T)
while True:
v = random.uniform(botPHz, topPHz)
I = Planck(v, T)
a = I/I0
if a >= 1:
v0 = v
I0 = I
else:
p = random.random()
if p < a:
v0 = v
I0 = I
yield v0
# example histogram for some temperature
samples = 1000000
temp = 5777 # in Kelvin
data = [to_nm(v) for v, _ in zip(MC(from_K(temp)), range(samples))] # in nanometers, from 100k samples.
specxyz = []
with open('lin2012xyz10e_fine_7sf.csv', newline='') as inputfile:
for row in csv.reader(inputfile):
specxyz.append(row)
col = (0, 0, 0)
r, g, b = 0., 0., 0.
for d in data:
x = d*10 - 3900
i = floor(x)
y = x - i
vals = zip(specxyz[i], specxyz[i + 1])
ivals = []
for v in vals:
ivals.append(float(v[0])*(1-y) + float(v[1])*y)
r += ivals[1]
g += ivals[2]
b += ivals[3]
t = r + g + b
r /= g
b /= g
g = 1
def xyztorgb(x, y, z):
r = 3.2404542 * x - 1.5371385 * y - 0.4985314 * z
g = -0.969266 * x + 1.8760108 * y + 1.8760108 * z
b = 0.0556434 * x - 0.2040259 * y + 1.0572252 * z
return r, g, b
print('XYZ')
print(r)
print(g)
print(b)
print()
r, g, b = xyztorgb(r, g, b)
r = int(min(1, max(0, r)) * 255)
g = int(min(1, max(0, g)) * 255)
b = int(min(1, max(0, b)) * 255)
print('RGB')
print(r)
print(g)
print(b)
print()
binbot = ceil(min(data))
bintop = floor(max(data))
bincount = int(bintop-binbot)
bins = np.linspace(binbot,
bintop,
bincount)
plt.xlim([min(data), max(data)])
plt.hist(data, bins=bins, alpha=0.5)
plt.title('spectral distribution of ' + repr(temp) + 'K light source')
plt.xlabel('nm')
plt.ylabel('count')
plt.show()
It will output a plot of the sampled spectrum's histogram as well as the corresponding XYZ and RGB values (using a XYZ-to-RGB matrix I found somewhere)
The file "lin2012xyz10e_fine_7sf.csv" is just a csv I got from here:
http://cvrl.ioo.ucl.ac.uk/In particular, and this might be a slight error, not quite sure, I used
"10-deg XYZ CMFs transformed from the CIE (2006) 2-deg LMS cone fundamentals"
in
"New CIE XYZ (from LMS) functions (proposed)"
I picked (and in this code it's hardcoded) the 0.1nm option as a csv file.
What might be problematic about that is the "proposed": It's a new standard that's apparently not yet in common use. Presumably it's more accurate than the old version, however it might actually need a slightly different matrix to convert from XYZ to sRGB. Not sure.
And I chose the variant where I normalize my XYZ such that Y = 1 which might not actually be what you'd want to do.
My output for a 5777K lightsource ended up being
XYZ
0.9739410242962122
1
0.9990353465096308
RGB
255
255
231
Which agrees really well with an actual visible portion of a black body spectrum of this temperature.
EDIT:
So the way this works would be, you'd use spectral data to calculate the color of your gear (use lots of samples. The more samples, the closer to actual white your light source will be). You don't actually need to store the whole spectrum: Just convert each final ray as it hits the camera into XYZ colors using the lookup table from the page I linked above. In the XYZ color space, add everything up. If I understood right, it's conveniently linear, so adding up in spectral space and then converting to XYZ is the same as converting to XYZ and then adding up in that space.
Then, once your render is finished, first find your maximum Y value and normalize the entire image by that. That will mean that no point is actually brighter than this. If you are concerned with saturation-clipping, it might help to actually make the brightest point darker, so normalize it such that the maximum Y is only like .8 or something. Then, convert it all over to your desired color space, for instance, for screens, the sRGB color space. That then is the final image.
Oh and the way I set it up in the python script above, I actually interpolate wavelengths linearly to get an approximate color for wavelengths that aren't exactly .1nm. Presumably, at that resolution, you could also just round them. Would probably work just as well.