srgb_to_reflectance.py

# Generating Reflectance Curves from sRGB Triplets
# http://scottburns.us/reflectance-curves-from-srgb-10/

# Original Matlab code & article:
# http://scottburns.us/wp-content/uploads/2019/05/LHTSS-text-file.txt
# http://scottburns.us/reflectance-curves-from-srgb-10/

# Note: This file was transpiled from Matlab by GPT-4 with
# some additional human help by @mattdesl to fix bugs. There
# may still be crucial errors—notably, it cannot solve due to "Singular matrix" error.

# This file is licensed the same as the Matlab file:
# Licensed under a Creative Commons Attribution-ShareAlike 4.0 International
# License (http://creativecommons.org/licenses/by-sa/4.0/).

import numpy as np
import sys
np.set_printoptions(threshold=sys.maxsize)
from scipy.sparse import diags

def LHTSS(T, sRGB):
  """
  This is the Least Hyperbolic Tangent Slope Squared (LHTSS) algorithm for
  generating a "reasonable" reflectance curve from a given sRGB color triplet.
  The reflectance spans the wavelength range 380-730 nm in 10 nm increments.

  :param T: 3x36 matrix converting reflectance to D65-weighted linear rgb
  :param sRGB: 3 element vector of target D65 referenced sRGB values in 0-255 range
  :return rho: 36x1 vector of reconstructed reflectance values, all (0->1)
  """

  rho = np.zeros(36)

  # handle special case of (0, 0, 0)
  if np.all(sRGB == 0):
    rho = 0.0001 * np.ones(36)
    return rho

  # handle special case of (255, 255, 255)
  if np.all(sRGB == 255):
    rho = np.ones(36)
    return rho
  
  # 36x36 difference matrix for Jacobian
  D = diags([-2, 4, -2], [-1, 0, 1], shape=(36, 36)).toarray()
  D[0, 0] = 2
  D[35, 35] = 2

  # compute target linear rgb values
  sRGB = sRGB / 255  # convert to 0-1
  rgb = np.zeros(3)
  # remove gamma correction to get linear rgb
  for i in range(3):
    if sRGB[i] < 0.04045:
      rgb[i] = sRGB[i] / 12.92
    else:
      rgb[i] = ((sRGB[i] + 0.055) / 1.055) ** 2.4

  # initialize
  z = np.zeros(36)  # starting point all rho=1/2
  lambda_ = np.zeros(3)  # starting Lagrange mult
  maxit = 100  # max number of iterations
  ftol = 1.0e-8  # function solution tolerance
  count = 0  # iteration counter

  # with np.errstate(over='ignore'):
  # Newton's method iteration
  while count <= maxit:
    d0 = (np.tanh(z) + 1) / 2
    d1 = np.diag((1 - np.tanh(z)**2) / 2)
    d2 = np.diag(-np.tanh(z) * (1 - np.tanh(z)**2))

    F = np.concatenate([D @ z + d1 @ T.T @ lambda_, T @ d0 - rgb])  # 39x1 F vector

    # Calculate the four blocks of the J matrix
    top_left_block = D + np.diag(np.diag(d2) * (T.T @ lambda_))
    top_right_block = d1 @ T.T
    bottom_left_block = T @ d1
    bottom_right_block = np.zeros((3, 3))

    # Combine the blocks to create the 39x39 J matrix
    J = np.block([[top_left_block, top_right_block], [bottom_left_block, bottom_right_block]])
    
    delta = np.linalg.solve(J, -F)  # solve Newton system of equations J*delta = -F
    z += delta[:36]  # update z
    lambda_ += delta[36:]  # update lambda
    if np.all(np.abs(F) < ftol):
      # solution found
      rho = (np.tanh(z) + 1) / 2
      return rho
    count += 1
    print(count)
  print(f'No solution found in {maxit} iterations.')
  return None

T = np.array([
  [
    0.0000646936115727633,  0.000219415369171578,
      0.00112060228414359,   0.00376670730427686,
        0.0118808497572766,    0.0232870228938867,
        0.0345602796797156,    0.0372247180152918,
        0.0324191842208867,    0.0212337349018611,
        0.0104912522835777,   0.00329591973705558,
      0.000507047802540891,  0.000948697853868474,
      0.00627387448845597,    0.0168650445840847,
        0.0286903641895679,    0.0426758762490725,
        0.0562561504260008,    0.0694721289967602,
        0.0830552220141023,    0.0861282432155783,
        0.0904683927868683,    0.0850059839999687,
        0.0709084366392777,    0.0506301536932269,
        0.0354748461653679,    0.0214687454102844,
        0.0125167687669176,   0.00680475126078526,
      0.00346465215790157,   0.00149764708248624,
      0.000769719667700118,  0.000407378212832335,
      0.000169014616182123, 0.0000952268887534793
  ],
  [
    0.00000184433541764457, 0.0000062054782702308,
      0.0000310103776744139,  0.000104750996050908,
      0.000353649345357243,  0.000951495123526191,
        0.00228232006613489,   0.00420743392201395,
        0.00668896510747318,   0.00988864251316196,
          0.015249831581587,    0.0214188448516808,
        0.0334237633103485,    0.0513112925264347,
        0.0704038388936896,    0.0878408968669549,
        0.0942514030194481,    0.0979591120948518,
          0.094154532672617,    0.0867831869897857,
          0.078858499565938,    0.0635282861874625,
        0.0537427564004085,    0.0426471274206905,
        0.0316181374233466,    0.0208857265390802,
        0.0138604556350511,   0.00810284218307029,
        0.00463021767605804,     0.002491442109212,
        0.00125933475912608,  0.000541660024106255,
      0.000277959820700288,  0.000147111734433903,
      0.0000610342686915558, 0.0000343881801451621
  ],
  [
      0.000305024750978023,    0.00103683251144092,
        0.00531326877604233,     0.0179548401495523,
          0.057079004340659,       0.11365445199637,
          0.173363047597462,      0.196211466514214,
          0.186087009289904,      0.139953964010199,
        0.0891767523322851,     0.0478974052884572,
        0.0281463269981882,     0.0161380645679562,
        0.00775929533717298,    0.00429625546625385,
        0.00200555920471153,   0.000861492584272158,
      0.000369047917008248,   0.000191433500712763,
      0.000149559313956664,  0.0000923132295986905,
      0.0000681366166724671,  0.0000288270841412222,
      0.0000157675750930075, 0.00000394070233244055,
    0.00000158405207257727,                      0,
                          0,                      0,
                          0,                      0,
                          0,                      0,
                          0,                      0
  ]
])
sRGB = np.array([250, 92, 8])
print(LHTSS(T, sRGB))