Adapt Jakob and Hanika's interpolation code

6 files changed, 62 insertions, 2 deletions

diff --git a/data/aces2065_1.coeff b/data/aces2065_1.coeff
new file mode 100644
index 0000000..25614e7
--- /dev/null
+++ b/data/aces2065_1.coeff
diff --git a/data/prophotorgb.coeff b/data/prophotorgb.coeff
new file mode 100644
index 0000000..2c15813
--- /dev/null
+++ b/data/prophotorgb.coeff
diff --git a/data/rec2020.coeff b/data/rec2020.coeff
new file mode 100644
index 0000000..7ce1269
--- /dev/null
+++ b/data/rec2020.coeff
diff --git a/data/srgb.coeff b/data/srgb.coeff
new file mode 100644
index 0000000..a0b2ba2
--- /dev/null
+++ b/data/srgb.coeff
diff --git a/jakob_hanika.py b/jakob_hanika.py
index af55c94..22b87fc 100644
--- a/jakob_hanika.py
+++ b/jakob_hanika.py
@@ -15,10 +15,11 @@ def model(wvlp, ccp):
 	return 1 / 2 + yy / (2 * np.sqrt(1 + yy ** 2))
 
 # Create a SpectralDistribution using given coefficients
-def model_sd(ccp):
+def model_sd(ccp, primed=True):
 	# FIXME: don't hardcode the wavelength grid; there should be a way
 	#        of creating a SpectralDistribution from the function alone
-	return SpectralDistribution(model(wvlp, ccp), wvl, name="Model")
+	grid = wvlp if primed else wvl
+	return SpectralDistribution(model(grid, ccp), wvl, name="Model")
 
 # The goal is to minimize the color difference between a given distrbution
 # and the one computed from the model above.
diff --git a/test_coeff.py b/test_coeff.py
new file mode 100644
index 0000000..e53c6f1
--- /dev/null
+++ b/test_coeff.py
@@ -0,0 +1,59 @@
+import numpy as np, struct
+from scipy.interpolate import RegularGridInterpolator
+from colour import *
+from colour.plotting import *
+from colour.models import eotf_inverse_sRGB
+from jakob_hanika import jakob_hanika, model_sd
+from test_colorchecker import plot_comparison
+
+
+
+illuminant = "D65"
+ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant]
+ill_sd = SpectralDistribution(ILLUMINANTS_SDS[illuminant])
+
+
+
+class JakobHanikaInterpolator:
+	def __init__(self):
+		pass
+
+	def read_file(self, path):
+		with open(path, "rb") as fd:
+			if fd.read(4).decode("ISO-8859-1") != "SPEC":
+				raise ValueError("bad magic number")
+
+			self.res = struct.unpack("i", fd.read(4))[0]
+			self.scale = np.fromfile(fd, count=self.res, dtype=np.float32) # FIXME: default dtype
+			coeffs = np.fromfile(fd, count=3*self.res**3*3, dtype=np.float32) # FIXME: default dtype
+			coeffs = coeffs.reshape(3, self.res, self.res, self.res, 3)
+
+		t = np.linspace(0, 1, self.res)
+		axes = [self.scale, t, t]
+		self.cubes = [RegularGridInterpolator(axes, coeffs[j, :, :, :, :], bounds_error=False) for j in range(3)]
+
+	def __call__(self, RGB):
+		RGB = np.asarray(RGB, dtype=np.float) # FIXME: default dtype
+		vmax = RGB.max()
+		imax = RGB.argmax()
+		cube = self.cubes[imax]
+		chroma = RGB / (vmax + 1e-10)
+		chroma = np.array([
+			vmax,
+			chroma[..., (imax + 2) % 3],
+			chroma[..., (imax + 1) % 3]
+		])
+		return cube(chroma).squeeze()
+
+
+
+if __name__ == "__main__":
+	jh = JakobHanikaInterpolator()
+	jh.read_file("data/srgb.coeff")
+
+	linear_RGB = np.array([1, 0, 0]) # A random test color
+	RGB = eotf_inverse_sRGB(linear_RGB)
+	XYZ = sRGB_to_XYZ(RGB)
+	cc = jh(linear_RGB)
+
+	plot_comparison(XYZ, model_sd(cc, primed=False), "Model", 123, ill_sd)