I’m not sure if this will be useful, but I tested constructing reflectance curves from a set of precalculated curves. I saw that in Add spectral blending of mypaint (!1249) · Merge requests · Graphics / Krita · GitLab David Revoy mentioned it was slow. So, hopefully this will be useful in that regard (although that has been closed). The precalculated spectral curves were created using Ronald van Wijnen’s method. The method for constructing the curves is pretty simple. It is a linear interpolation between the two closest colors in the dictionary containing the precalculated curves. This result is then interpolated with a white reflectance curve. The curve is then scaled, which is the final approximate curve.
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from IPython.display import display
SIZE = 38
GAMMA = 2.4
EPSILON = 0.00000001
SPD_C = [0.96853629, 0.96855103, 0.96859338, 0.96877345, 0.96942204, 0.97143709, 0.97541862, 0.98074186, 0.98580992, 0.98971194, 0.99238027, 0.99409844, 0.995172, 0.99576545, 0.99593552, 0.99564041, 0.99464769, 0.99229579, 0.98638762, 0.96829712, 0.89228016, 0.53740239, 0.15360445, 0.05705719, 0.03126539, 0.02205445, 0.01802271, 0.0161346, 0.01520947, 0.01475977, 0.01454263, 0.01444459, 0.01439897, 0.0143762, 0.01436343, 0.01435687, 0.0143537, 0.01435408]
SPD_M = [0.51567122, 0.5401552, 0.62645502, 0.75595012, 0.92826996, 0.97223624, 0.98616174, 0.98955255, 0.98676237, 0.97312575, 0.91944277, 0.32564851, 0.13820628, 0.05015143, 0.02912336, 0.02421691, 0.02660696, 0.03407586, 0.04835936, 0.0001172, 0.00008554, 0.85267882, 0.93188793, 0.94810268, 0.94200977, 0.91478045, 0.87065445, 0.78827548, 0.65738359, 0.59909403, 0.56817268, 0.54031997, 0.52110241, 0.51041094, 0.50526577, 0.5025508, 0.50126452, 0.50083021]
SPD_Y = [0.02055257, 0.02059936, 0.02062723, 0.02073387, 0.02114202, 0.02233154, 0.02556857, 0.03330189, 0.05185294, 0.10087639, 0.24000413, 0.53589066, 0.79874659, 0.91186529, 0.95399623, 0.97137099, 0.97939505, 0.98345207, 0.98553736, 0.98648905, 0.98674535, 0.98657555, 0.98611877, 0.98559942, 0.98507063, 0.98460039, 0.98425301, 0.98403909, 0.98388535, 0.98376116, 0.98368246, 0.98365023, 0.98361309, 0.98357259, 0.98353856, 0.98351247, 0.98350101, 0.98350852]
SPD_R = [0.03147571, 0.03146636, 0.03140624, 0.03119611, 0.03053888, 0.02856855, 0.02459485, 0.0192952, 0.01423112, 0.01033111, 0.00765876, 0.00593693, 0.00485616, 0.00426186, 0.00409039, 0.00438375, 0.00537525, 0.00772962, 0.0136612, 0.03181352, 0.10791525, 0.46249516, 0.84604333, 0.94275572, 0.96860996, 0.97783966, 0.98187757, 0.98377315, 0.98470202, 0.98515481, 0.98537114, 0.98546685, 0.98550011, 0.98551031, 0.98550741, 0.98551323, 0.98551563, 0.98551547]
SPD_G = [0.49108579, 0.46944057, 0.4016578, 0.2449042, 0.0682688, 0.02732883, 0.013606, 0.01000187, 0.01284127, 0.02636635, 0.07058713, 0.70421692, 0.85473994, 0.95081565, 0.9717037, 0.97651888, 0.97429245, 0.97012917, 0.9425863, 0.99989207, 0.99989891, 0.13823139, 0.06968113, 0.05628787, 0.06111561, 0.08987709, 0.13656016, 0.22169624, 0.32176956, 0.36157329, 0.4836192, 0.46488579, 0.47440306, 0.4857699, 0.49267971, 0.49625685, 0.49807754, 0.49889859]
SPD_B = [0.97901834, 0.97901649, 0.97901118, 0.97892146, 0.97858555, 0.97743705, 0.97428075, 0.96663223, 0.94822893, 0.89937713, 0.76070164, 0.4642044, 0.20123039, 0.08808402, 0.04592894, 0.02860373, 0.02060067, 0.01656701, 0.01451549, 0.01357964, 0.01331243, 0.01347661, 0.01387181, 0.01435472, 0.01479836, 0.0151525, 0.01540513, 0.01557233, 0.0156571, 0.01571025, 0.01571916, 0.01572133, 0.01572502, 0.01571717, 0.01571905, 0.01571059, 0.01569728, 0.0157002]
CIE_CMF_X = [0.00006469, 0.00021941, 0.00112057, 0.00376661, 0.01188055, 0.02328644, 0.03455942, 0.03722379, 0.03241838, 0.02123321, 0.01049099, 0.00329584, 0.00050704, 0.00094867, 0.00627372, 0.01686462, 0.02868965, 0.04267481, 0.05625475, 0.0694704, 0.08305315, 0.0861261, 0.09046614, 0.08500387, 0.07090667, 0.05062889, 0.03547396, 0.02146821, 0.01251646, 0.00680458, 0.00346457, 0.00149761, 0.0007697, 0.00040737, 0.00016901, 0.00009522, 0.00004903, 0.00002]
CIE_CMF_Y = [0.00000184, 0.00000621, 0.00003101, 0.00010475, 0.00035364, 0.00095147, 0.00228226, 0.00420733, 0.0066888, 0.0098884, 0.01524945, 0.02141831, 0.03342293, 0.05131001, 0.07040208, 0.08783871, 0.09424905, 0.09795667, 0.09415219, 0.08678102, 0.07885653, 0.0635267, 0.05374142, 0.04264606, 0.03161735, 0.02088521, 0.01386011, 0.00810264, 0.0046301, 0.00249138, 0.0012593, 0.00054165, 0.00027795, 0.00014711, 0.00006103, 0.00003439, 0.00001771, 0.00000722]
CIE_CMF_Z = [0.00030502, 0.00103681, 0.00531314, 0.01795439, 0.05707758, 0.11365162, 0.17335873, 0.19620658, 0.18608237, 0.13995048, 0.08917453, 0.04789621, 0.02814563, 0.01613766, 0.0077591, 0.00429615, 0.00200551, 0.00086147, 0.00036904, 0.00019143, 0.00014956, 0.00009231, 0.00006813, 0.00002883, 0.00001577, 0.00000394, 0.00000158, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
XYZ_RGB = [[3.24306333, -1.53837619, -0.49893282], [-0.96896309, 1.87542451, 0.04154303], [0.05568392, -0.20417438, 1.05799454]]
def uncompand(x):
return x / 12.92 if x < 0.04045 else ((x + 0.055) / 1.055) ** GAMMA
def compand(x):
return x * 12.92 if x < 0.0031308 else 1.055 * x ** (1.0 / GAMMA) - 0.055
def srgb_to_linear(srgb):
r = uncompand(srgb[0] / 255)
g = uncompand(srgb[1] / 255)
b = uncompand(srgb[2] / 255)
return [r, g, b]
def linear_to_srgb(lrgb):
r = compand(lrgb[0])
g = compand(lrgb[1])
b = compand(lrgb[2])
return [round(clamp(r, 0, 1) * 255), round(clamp(g, 0, 1) * 255), round(clamp(b, 0, 1) * 255)]
def reflectance_to_xyz(R):
x = dotproduct(R, CIE_CMF_X)
y = dotproduct(R, CIE_CMF_Y)
z = dotproduct(R, CIE_CMF_Z)
return [x, y, z]
def xyz_to_srgb(xyz):
r = dotproduct(XYZ_RGB[0], xyz)
g = dotproduct(XYZ_RGB[1], xyz)
b = dotproduct(XYZ_RGB[2], xyz)
return linear_to_srgb([r, g, b])
def spectral_upsampling(lrgb):
w = min(min(lrgb[0], lrgb[1]), lrgb[2])
lrgb = [lrgb[0] - w, lrgb[1] - w, lrgb[2] - w]
c = min(lrgb[1], lrgb[2])
m = min(lrgb[0], lrgb[2])
y = min(lrgb[0], lrgb[1])
r = max(0, min(lrgb[0] - lrgb[2], lrgb[0] - lrgb[1]))
g = max(0, min(lrgb[1] - lrgb[2], lrgb[1] - lrgb[0]))
b = max(0, min(lrgb[2] - lrgb[1], lrgb[2] - lrgb[0]))
return [w, c, m, y, r, g, b]
def linear_to_reflectance(lrgb):
weights = spectral_upsampling(lrgb)
R = [0] * SIZE
for i in range(SIZE):
R[i] = (
max(EPSILON,
weights[0]
+ weights[1] * SPD_C[i]
+ weights[2] * SPD_M[i]
+ weights[3] * SPD_Y[i]
+ weights[4] * SPD_R[i]
+ weights[5] * SPD_G[i]
+ weights[6] * SPD_B[i]
)
)
return R
def dotproduct(a, b):
return sum(x * y for x, y in zip(a, b))
def clamp(value, min_value, max_value):
return min(max(value, min_value), max_value)
def hex_to_rgb_255(hex_color):
hex_color = hex_color.lstrip('#')
return tuple(int(hex_color[i:i+2], 16) for i in (0, 2, 4))
class ColorLUT:
def __init__(self):
self.anchors = [
{'name': 'Red', 'hex': '#ff0000', 'rgb': (255,0,0)},
{'name': 'Orange', 'hex': '#ff8000', 'rgb': (255,128,0)},
{'name': 'Yellow', 'hex': '#ffff00', 'rgb': (255,255,0)},
{'name': 'Chartreuse','hex': '#80ff00', 'rgb': (128,255,0)},
{'name': 'Green', 'hex': '#00ff00', 'rgb': (0,255,0)},
{'name': 'Spring', 'hex': '#00ff80', 'rgb': (0,255,128)},
{'name': 'Cyan', 'hex': '#00ffff', 'rgb': (0,255,255)},
{'name': 'Azure', 'hex': '#0080ff', 'rgb': (0,128,255)},
{'name': 'Blue', 'hex': '#0000ff', 'rgb': (0,0,255)},
{'name': 'Violet', 'hex': '#8000ff', 'rgb': (128,0,255)},
{'name': 'Magenta', 'hex': '#ff00ff', 'rgb': (255,0,255)},
{'name': 'Rose', 'hex': '#ff0080', 'rgb': (255,0,128)},
]
for anchor in self.anchors:
lrgb = srgb_to_linear(anchor['rgb'])
anchor['R'] = linear_to_reflectance(lrgb)
white_lrgb = srgb_to_linear((255, 255, 255))
self.white_anchor = {
'name': 'White',
'rgb': (255, 255, 255),
'R': linear_to_reflectance(white_lrgb)
}
def get_neighbors_scale_and_whiteness(self, target_hex):
t_rgb_255 = hex_to_rgb_255(target_hex)
t_rgb_norm = np.array([x/255.0 for x in t_rgb_255])
max_val = np.max(t_rgb_norm)
if max_val == 0:
return self.anchors[0], self.anchors[0], 0, 0.0, 0.0, self.white_anchor
shell_rgb = t_rgb_norm / max_val
whiteness = np.min(shell_rgb)
if whiteness >= 1.0:
chroma_rgb = shell_rgb
else:
chroma_rgb = (shell_rgb - whiteness) / (1 - whiteness)
distances = []
for i, anchor in enumerate(self.anchors):
a_rgb = np.array([x/255.0 for x in anchor['rgb']])
dist = np.linalg.norm(chroma_rgb - a_rgb)
distances.append((dist, i))
distances.sort(key=lambda x: x[0])
n1 = self.anchors[distances[0][1]]
n2 = self.anchors[distances[1][1]]
v1 = np.array([x/255.0 for x in n1['rgb']])
v2 = np.array([x/255.0 for x in n2['rgb']])
vec_line = v2 - v1
vec_point = chroma_rgb - v1
len_sq = np.dot(vec_line, vec_line)
t = 0 if len_sq == 0 else np.dot(vec_point, vec_line) / len_sq
t = np.clip(t, 0, 1)
lin_t_rgb = srgb_to_linear(t_rgb_255)
lin_max = max(lin_t_rgb)
scale = lin_max
return n1, n2, t, scale, whiteness, self.white_anchor
class LUTValidatorUI:
def __init__(self):
self.lut = ColorLUT()
self.color_picker = widgets.ColorPicker(description='Target', value='#d3d1bb')
self.dampener = widgets.FloatSlider(
value=2.0, min=1.0, max=5.0, step=0.1,
description='White Dampening'
)
self.output = widgets.Output()
self.color_picker.observe(self.update, names='value')
self.dampener.observe(self.update, names='value')
self.update(None)
def update(self, change):
target_hex = self.color_picker.value
damp_factor = self.dampener.value
t_rgb_255 = hex_to_rgb_255(target_hex)
t_lrgb = srgb_to_linear(t_rgb_255)
R_true = linear_to_reflectance(t_lrgb)
n1, n2, t, scale, whiteness, white_anchor = self.lut.get_neighbors_scale_and_whiteness(target_hex)
w_adjusted = whiteness ** damp_factor
R1 = n1['R']
R2 = n2['R']
R_w = white_anchor['R']
R_approx = []
for i in range(SIZE):
val_hue = (1 - t) * R1[i] + t * R2[i]
val_desat = (w_adjusted * R_w[i]) + ((1 - w_adjusted) * val_hue)
R_approx.append(val_desat * scale)
R_true_arr = np.array(R_true)
R_approx_arr = np.array(R_approx)
rmse = np.sqrt(np.mean((R_true_arr - R_approx_arr)**2))
xyz = reflectance_to_xyz(R_approx)
rgb_approx_list = xyz_to_srgb(xyz)
hex_approx = '#{:02x}{:02x}{:02x}'.format(*rgb_approx_list)
wl = np.linspace(380, 750, SIZE)
with self.output:
self.output.clear_output(wait=True)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
ax1.set_title(f"Target: {target_hex}\nLinear White: {whiteness:.2f} -> Adjusted: {w_adjusted:.2f}")
ax1.plot(wl, R_true_arr, color='black', linewidth=4, alpha=0.3, label="True")
ax1.plot(wl, R_approx_arr, color=target_hex, linestyle='--', linewidth=2, label="Approx")
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-0.05, 1.05)
ax2.axis('off')
ax2.set_title(f"Accuracy Check\nRMSE: {rmse:.5f}")
ax2.add_patch(plt.Rectangle((0, 0.5), 1, 0.5, color=target_hex))
ax2.add_patch(plt.Rectangle((0, 0), 1, 0.5, color=hex_approx))
ax2.text(0.5, 0.75, "Ground Truth", ha='center', va='center', color='white')
ax2.text(0.5, 0.25, "Adjusted LUT", ha='center', va='center', color='white')
ax2.set_xlim(0, 1); ax2.set_ylim(0, 1)
plt.tight_layout()
plt.show()
def display(self):
display(widgets.VBox([
self.color_picker,
self.dampener,
self.output
]))
app = LUTValidatorUI()
app.display()
Below is the optimized version of the code above. I used Gemini 3.0 to do it; therefore, I figured including the initial code would be best. The main difference I noticed between the methods is the interpolation value ‘t’ is different. The optimized version uses properties of the color wheel to find the closest reference colors quicker than comparing the color to all the other references to find the closest two.
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from IPython.display import display
from numba import njit, prange
# ==========================================
# 1. DATA & CONSTANTS
# ==========================================
SIZE = 38
GAMMA = 2.4
EPSILON = 1e-8
# Raw SPDs for Ground Truth Calculation
_SPD_C = [0.9685, 0.9686, 0.9686, 0.9688, 0.9694, 0.9714, 0.9754, 0.9807, 0.9858, 0.9897, 0.9924, 0.9941, 0.9952, 0.9958, 0.9959, 0.9956, 0.9946, 0.9923, 0.9864, 0.9683, 0.8923, 0.5374, 0.1536, 0.0571, 0.0313, 0.0221, 0.0180, 0.0161, 0.0152, 0.0148, 0.0145, 0.0144, 0.0144, 0.0144, 0.0144, 0.0144, 0.0144, 0.0144]
_SPD_M = [0.5157, 0.5402, 0.6265, 0.7560, 0.9283, 0.9722, 0.9862, 0.9896, 0.9868, 0.9731, 0.9194, 0.3256, 0.1382, 0.0502, 0.0291, 0.0242, 0.0266, 0.0341, 0.0484, 0.0001, 0.0001, 0.8527, 0.9319, 0.9481, 0.9420, 0.9148, 0.8707, 0.7883, 0.6574, 0.5991, 0.5682, 0.5403, 0.5211, 0.5104, 0.5053, 0.5026, 0.5013, 0.5008]
_SPD_Y = [0.0206, 0.0206, 0.0206, 0.0207, 0.0211, 0.0223, 0.0256, 0.0333, 0.0519, 0.1009, 0.2400, 0.5359, 0.7987, 0.9119, 0.9540, 0.9714, 0.9794, 0.9835, 0.9855, 0.9865, 0.9867, 0.9866, 0.9861, 0.9856, 0.9851, 0.9846, 0.9843, 0.9840, 0.9839, 0.9838, 0.9837, 0.9837, 0.9836, 0.9836, 0.9835, 0.9835, 0.9835, 0.9835]
_SPD_R = [0.0315, 0.0315, 0.0314, 0.0312, 0.0305, 0.0286, 0.0246, 0.0193, 0.0142, 0.0103, 0.0077, 0.0059, 0.0049, 0.0043, 0.0041, 0.0044, 0.0054, 0.0077, 0.0137, 0.0318, 0.1079, 0.4625, 0.8460, 0.9428, 0.9686, 0.9778, 0.9819, 0.9838, 0.9847, 0.9852, 0.9854, 0.9855, 0.9855, 0.9855, 0.9855, 0.9855, 0.9855, 0.9855]
_SPD_G = [0.4911, 0.4694, 0.4017, 0.2449, 0.0683, 0.0273, 0.0136, 0.0100, 0.0128, 0.0264, 0.0706, 0.7042, 0.8547, 0.9508, 0.9717, 0.9765, 0.9743, 0.9701, 0.9426, 0.9999, 0.9999, 0.1382, 0.0697, 0.0563, 0.0611, 0.0899, 0.1366, 0.2217, 0.3218, 0.3616, 0.4836, 0.4649, 0.4744, 0.4858, 0.4927, 0.4963, 0.4981, 0.4989]
_SPD_B = [0.9790, 0.9790, 0.9790, 0.9789, 0.9786, 0.9774, 0.9743, 0.9666, 0.9482, 0.8994, 0.7607, 0.4642, 0.2012, 0.0881, 0.0459, 0.0286, 0.0206, 0.0166, 0.0145, 0.0136, 0.0133, 0.0135, 0.0139, 0.0144, 0.0148, 0.0152, 0.0154, 0.0156, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157]
SPD_MATRIX = np.ascontiguousarray(np.array([_SPD_C, _SPD_M, _SPD_Y, _SPD_R, _SPD_G, _SPD_B], dtype=np.float64))
# CMF for Visualization conversion
CMF_MATRIX = np.ascontiguousarray(np.array([
[0.00006469, 0.00021941, 0.00112057, 0.00376661, 0.01188055, 0.02328644, 0.03455942, 0.03722379, 0.03241838, 0.02123321, 0.01049099, 0.00329584, 0.00050704, 0.00094867, 0.00627372, 0.01686462, 0.02868965, 0.04267481, 0.05625475, 0.0694704, 0.08305315, 0.0861261, 0.09046614, 0.08500387, 0.07090667, 0.05062889, 0.03547396, 0.02146821, 0.01251646, 0.00680458, 0.00346457, 0.00149761, 0.0007697, 0.00040737, 0.00016901, 0.00009522, 0.00004903, 0.00002],
[0.00000184, 0.00000621, 0.00003101, 0.00010475, 0.00035364, 0.00095147, 0.00228226, 0.00420733, 0.0066888, 0.0098884, 0.01524945, 0.02141831, 0.03342293, 0.05131001, 0.07040208, 0.08783871, 0.09424905, 0.09795667, 0.09415219, 0.08678102, 0.07885653, 0.0635267, 0.05374142, 0.04264606, 0.03161735, 0.02088521, 0.01386011, 0.00810264, 0.0046301, 0.00249138, 0.0012593, 0.00054165, 0.00027795, 0.00014711, 0.00006103, 0.00003439, 0.00001771, 0.00000722],
[0.00030502, 0.00103681, 0.00531314, 0.01795439, 0.05707758, 0.11365162, 0.17335873, 0.19620658, 0.18608237, 0.13995048, 0.08917453, 0.04789621, 0.02814563, 0.01613766, 0.0077591, 0.00429615, 0.00200551, 0.00086147, 0.00036904, 0.00019143, 0.00014956, 0.00009231, 0.00006813, 0.00002883, 0.00001577, 0.00000394, 0.00000158, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
], dtype=np.float64))
XYZ_RGB_MATRIX = np.ascontiguousarray(np.array([
[3.24306333, -1.53837619, -0.49893282],
[-0.96896309, 1.87542451, 0.04154303],
[0.05568392, -0.20417438, 1.05799454]
], dtype=np.float64))
# ==========================================
# 2. NUMBA JIT KERNELS (The Optimized Parts)
# ==========================================
@njit(fastmath=True)
def _srgb_to_linear_one(x):
if x < 0.04045:
return x / 12.92
else:
return ((x + 0.055) / 1.055) ** GAMMA
@njit(parallel=True, fastmath=True)
def get_spectra_batch_parallel(rgb_array, anchor_curves, white_curve, damp_factor):
"""
Highly Optimized Parallel Batch Processor
Input: rgb_array (N, 3) 0-255
Output: results (N, 38)
"""
n_pixels = rgb_array.shape[0]
results = np.empty((n_pixels, SIZE), dtype=np.float64)
for i in prange(n_pixels):
r = rgb_array[i, 0] / 255.0
g = rgb_array[i, 1] / 255.0
b = rgb_array[i, 2] / 255.0
max_val = max(r, max(g, b))
if max_val < 1e-9:
results[i, :] = 0.0
continue
lin_r = _srgb_to_linear_one(r)
lin_g = _srgb_to_linear_one(g)
lin_b = _srgb_to_linear_one(b)
scale = max(lin_r, max(lin_g, lin_b))
# Shell Projection
s_r = r / max_val
s_g = g / max_val
s_b = b / max_val
whiteness = min(s_r, min(s_g, s_b))
# Hue Calculation
shell_chroma = 1.0 - whiteness
hue = 0.0
if shell_chroma > 1e-9:
if s_r == 1.0:
hue = (s_g - s_b) / shell_chroma
elif s_g == 1.0:
hue = 2.0 + (s_b - s_r) / shell_chroma
else:
hue = 4.0 + (s_r - s_g) / shell_chroma
hue *= 60.0
if hue < 0: hue += 360.0
# O(1) Neighbor Indexing
segment = hue / 30.0
idx1 = int(segment) % 12
idx2 = (idx1 + 1) % 12
t = segment - int(segment)
w_adj = whiteness ** damp_factor
c1 = anchor_curves[idx1]
c2 = anchor_curves[idx2]
# Unroll for write
for k in range(SIZE):
hue_mix = (1.0 - t) * c1[k] + t * c2[k]
val = (w_adj * white_curve[k]) + ((1.0 - w_adj) * hue_mix)
results[i, k] = val * scale
return results
@njit(fastmath=True)
def linear_to_reflectance_jit(lrgb):
"""Ground Truth Calculator"""
w = np.min(lrgb)
rem = lrgb - w
c = min(rem[1], rem[2])
m = min(rem[0], rem[2])
y = min(rem[0], rem[1])
r = max(0.0, min(rem[0] - rem[2], rem[0] - rem[1]))
g = max(0.0, min(rem[1] - rem[2], rem[1] - rem[0]))
b = max(0.0, min(rem[2] - rem[1], rem[2] - rem[0]))
R = np.empty(SIZE, dtype=np.float64)
for i in range(SIZE):
sum_spd = (c * SPD_MATRIX[0, i] + m * SPD_MATRIX[1, i] + y * SPD_MATRIX[2, i] +
r * SPD_MATRIX[3, i] + g * SPD_MATRIX[4, i] + b * SPD_MATRIX[5, i])
R[i] = max(EPSILON, w + sum_spd)
return R
def hex_to_rgb_arr(hex_color):
h = hex_color.lstrip('#')
return np.array([int(h[0:2], 16), int(h[2:4], 16), int(h[4:6], 16)], dtype=np.float64)
# ==========================================
# 3. CLASS WRAPPER
# ==========================================
class ParallelLUT:
def __init__(self):
# 1. Prepare Anchors
anchors_raw = [(255,0,0), (255,128,0), (255,255,0), (128,255,0), (0,255,0), (0,255,128),
(0,255,255), (0,128,255), (0,0,255), (128,0,255), (255,0,255), (255,0,128)]
# Calculate curves for anchors
curves = []
for rgb in anchors_raw:
arr = np.array(rgb, dtype=np.float64) / 255.0
lrgb = np.array([_srgb_to_linear_one(x) for x in arr])
curves.append(linear_to_reflectance_jit(lrgb))
self.anchor_curves = np.array(curves, dtype=np.float64)
# Calculate curve for white
w_lrgb = np.array([_srgb_to_linear_one(1.0)]*3)
self.white_curve = linear_to_reflectance_jit(w_lrgb)
def get_curve_for_ui(self, hex_color, damp):
# UI passes one color, but the kernel expects a batch (N, 3)
# We wrap the single color in a 2D array: [[r, g, b]]
rgb_in = hex_to_rgb_arr(hex_color).reshape(1, 3)
# Call the Parallel Kernel
result_batch = get_spectra_batch_parallel(rgb_in, self.anchor_curves, self.white_curve, damp)
# Return the first (and only) row
return result_batch[0]
# ==========================================
# 4. VISUALIZATION UI
# ==========================================
class CheckOptimizedUI:
def __init__(self):
self.lut = ParallelLUT()
# Warmup the JIT kernel with dummy data so the UI doesn't lag on first click
dummy = np.array([[255, 255, 255]], dtype=np.float64)
print("Warming up JIT...")
_ = get_spectra_batch_parallel(dummy, self.lut.anchor_curves, self.lut.white_curve, 2.0)
print("Ready.")
self.color_picker = widgets.ColorPicker(description='Target', value='#d3d1bb')
self.dampener = widgets.FloatSlider(value=2.0, min=1.0, max=5.0, description='Dampening')
self.output = widgets.Output()
self.color_picker.observe(self.update, names='value')
self.dampener.observe(self.update, names='value')
self.update(None)
def update(self, change):
hex_c = self.color_picker.value
damp = self.dampener.value
# 1. Ground Truth
arr = hex_to_rgb_arr(hex_c) / 255.0
lrgb = np.array([_srgb_to_linear_one(x) for x in arr])
R_true = linear_to_reflectance_jit(lrgb)
# 2. Optimized LUT (The Parallel Kernel)
R_approx = self.lut.get_curve_for_ui(hex_c, damp)
# 3. RMSE
rmse = np.linalg.norm(R_true - R_approx) / np.sqrt(SIZE)
# 4. Plot
wl = np.linspace(380, 750, SIZE)
with self.output:
self.output.clear_output(wait=True)
plt.figure(figsize=(10, 5))
plt.title(f"Optimization Check: Is Parallel Batch Logic Correct?\nTarget: {hex_c} | RMSE: {rmse:.5f}")
plt.plot(wl, R_true, 'k-', linewidth=4, alpha=0.3, label="Ground Truth (Smits)")
plt.plot(wl, R_approx, '--', color=hex_c, linewidth=2, label="Parallel Batch Approx")
plt.legend()
plt.grid(True, alpha=0.3)
plt.ylim(-0.05, 1.05)
plt.show()
def display(self):
display(widgets.VBox([self.color_picker, self.dampener, self.output]))
app = CheckOptimizedUI()
app.display()
Here is a benchmark for going from sRGB → reflectance → mix with KM → xyz → sRGB for an 1920x1080 image.
Initializing Benchmark for 1920x1080 (2,073,600 pixels)…
Generating random image data…
Compiling Kernels…
— STARTING TIMING —
- Converted Image 1 to Reflectance: 0.6438s
- Converted Image 2 to Reflectance: 0.6040s
- Mixed (KM) and Rendered to RGB: 0.4204s
TOTAL PIPELINE TIME: 1.6683s
This is mixing all ~4 million colors.
Since this method only needs the two reference curves, whiteness and scale values, I figure each pixel on the canvas could store these values along with the RGBA values to further improve speed. As the reference curves are stored in a dictionary, you just have to store the keys which are strings. The whiteness and scale values are just floats.
All of the spectral mixing merge requests are closed, but hopefully this could still be useful.