Some important fixes and general cleanup

I missed the 'illuminant' argument in sd_to_XYZ and effectively treated the reflectivities as emission spectra. This was probably very wrong and is now fixed.

2 files changed, 54 insertions, 34 deletions

diff --git a/jakob_hanika.py b/jakob_hanika.py
index 3b35670..f881afd 100644
--- a/jakob_hanika.py
+++ b/jakob_hanika.py
@@ -2,25 +2,28 @@ import numpy as np, scipy.optimize as optimize
 from colour import *
 from colour.difference import *
 
+
+
 # The same wavelength grid is used throughout
-wvl = np.arange(360, 830, 5)
+wvl = np.arange(360, 830, 1)
+
+# Map wavelenghts from 360--830 nm to 0-1
+def remap(wvl):
+	return (wvl - 360) / (830 - 360)
 
 # This is the model of spectral reflectivity described in the article.
 def model(wvl, cc):
-	# One thing they didn't mention in the text is that their polynomial
-	# is nondimensionalized (mapped from 360--800 nm to 0--1).
-	def remap(x):
-		return (x - 360) / (830 - 360)
-
-	x = cc[0] * remap(wvl) ** 2 + cc[1] * remap(wvl) + cc[2]
+	x = cc[0] * wvl ** 2 + cc[1] * wvl + cc[2]
 	return 1 / 2 + x / (2 * np.sqrt(1 + x ** 2))
 
 # The goal is to minimize the color difference between a given distrbution
 # and the one computed from the model above.
-def error_function(cc, target, illuminant):
-	ev = model(wvl, cc)
+def error_function(cc, target, ill_sd, ill_xy):
+	# One thing they didn't mention in the text is that their polynomial
+	# is nondimensionalized (mapped from 360--800 nm to 0--1).
+	ev = model((wvl - 360) / (830 - 360), cc)
 	sd = SpectralDistribution(ev, wvl)
-	Lab = XYZ_to_Lab(sd_to_XYZ(sd), illuminant)
+	Lab = XYZ_to_Lab(sd_to_XYZ(sd, illuminant=ill_sd) / 100, ill_xy)
 	return delta_E_CIE1976(target, Lab)
 
 
@@ -31,20 +34,20 @@ def error_function(cc, target, illuminant):
 def cb_basinhopping(x, f, accept):
 	return f < 0.1
 
-# Finds coefficients for Jakob and Hanika's function
-def jakob_hanika(target, illuminant):
+# Finds the parameters for Jakob and Hanika's model
+def jakob_hanika(target, ill_sd, ill_xy, x0=(0, 0, 0)):
 	# First, a conventional solver is called. For 'yellow green' this
 	# actually fails: gets stuck at a local minimum that's far away
 	# from the global one.
 	# FIXME: better parameters, a better x0, a better method?
-	# FIXME: stop iterating as soon as delta E is negligible 
+	# FIXME: stop iterating as soon as delta E is negligible
 	opt = optimize.minimize(
-		error_function, (0, 0, 0), (target, illuminant),
+		error_function, x0, (target, ill_sd, ill_xy),
 		method="L-BFGS-B", options={"disp": True, "ftol": 1e-5}
 	)
 	print(opt)
 
-	error = error_function(opt.x, target, illuminant)
+	error = error_function(opt.x, target, ill_sd, ill_xy)
 	print("Delta E is %g" % error)
 
 	# Basin hopping is far more likely to find the actual minimum we're
@@ -52,12 +55,19 @@ def jakob_hanika(target, illuminant):
 	if error > 0.1:
 		print("Error too large, trying global optimization")
 		opt = optimize.basinhopping(
-			lambda cc: error_function(cc, target, illuminant),
-			(0, 0, 0), disp=True, callback=cb_basinhopping
+			lambda cc: error_function(cc, target, ill_sd, ill_xy),
+			x0, disp=True, callback=cb_basinhopping
 		)
 		print(opt)
-		error = error_function(opt.x, target)
+		error = error_function(opt.x, target, ill_sd, ill_xy)
 		print("Global delta E is %g" % error)
 
-	matched_sd = SpectralDistribution(model(wvl, opt.x), wvl, name="Model")
-	return opt.x, matched_sd, error
+	# Map back to dimensional coefficients
+	cc = np.array([
+		opt.x[0] / 220900,
+		opt.x[1] / 470 - 36/11045 * opt.x[0],
+		opt.x[2] - 36/47 * opt.x[1] + 1296/2209 * opt.x[0]
+	])
+	matched_sd = SpectralDistribution(model(wvl, cc), wvl, name="Model")
+
+	return cc, matched_sd, error
diff --git a/test_colorchecker.py b/test_colorchecker.py
index 907322b..c934518 100644
--- a/test_colorchecker.py
+++ b/test_colorchecker.py
@@ -4,34 +4,44 @@ from colour.plotting import *
 from matplotlib import pyplot as plt
 from jakob_hanika import jakob_hanika
 
-D65 = ILLUMINANTS['CIE 1931 2 Degree Standard Observer']['D65']
+
+
+illuminant = "D65"
+ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant]
+ill_sd = SpectralDistribution(ILLUMINANTS_SDS[illuminant])
+
+
 
 # Makes a comparison plot with SDs and swatches
-def plot_comparison(target_sd, matched_sd, label, error):
-	target_XYZ = sd_to_XYZ(target_sd)
-	target_RGB = np.clip(XYZ_to_sRGB(target_XYZ / 100), 0, 1)
+def plot_comparison(target, matched_sd, label, error, ill_sd):
+	if type(target) is SpectralDistribution:
+		target_XYZ = sd_to_XYZ(target, illuminant=ill_sd) / 100
+	else:
+		target_XYZ = target
+	target_RGB = np.clip(XYZ_to_sRGB(target_XYZ), 0, 1)
 	target_swatch = ColourSwatch(label, target_RGB)
-	matched_XYZ = sd_to_XYZ(matched_sd)
-	matched_RGB = np.clip(XYZ_to_sRGB(matched_XYZ / 100), 0, 1)
+	matched_XYZ = sd_to_XYZ(matched_sd, illuminant=ill_sd) / 100
+	matched_RGB = np.clip(XYZ_to_sRGB(matched_XYZ), 0, 1)
 	matched_swatch = ColourSwatch("Model", matched_RGB)
       
 	axes = plt.subplot(2, 1, 1)
 	plt.title(label)
-	plot_multi_sds([target_sd, matched_sd], axes=axes, standalone=False)
+	if type(target) is SpectralDistribution:
+		plot_multi_sds([target, matched_sd], axes=axes, standalone=False)
+	else:
+		plot_single_sd(matched_sd, axes=axes, standalone=False)
       
 	axes = plt.subplot(2, 1, 2)
 	plt.title("ΔE = %g" % error)
 	plot_multi_colour_swatches([target_swatch, matched_swatch], axes=axes)
 
-
-
 if __name__ == "__main__":
 	# This demo goes through SDs in a color checker
 	for name, sd in COLOURCHECKERS_SDS['ColorChecker N Ohta'].items():
-		XYZ = sd_to_XYZ(sd)
-		target = XYZ_to_Lab(XYZ, D65)
+		XYZ = sd_to_XYZ(sd, illuminant=ill_sd) / 100
+		target = XYZ_to_Lab(XYZ, ill_xy)
 
-		print("The target is '%s' with L=%g, a=%g, b=%g" % (name, *target))
-		_, matched_sd, error = jakob_hanika(target, D65)
+		print("Color checker: The target is '%s' with L=%g, a=%g, b=%g" % (name, *target))
+		_, matched_sd, error = jakob_hanika(target, ill_sd, ill_xy)
 
-		plot_comparison(sd, matched_sd, name, error)
+		plot_comparison(sd, matched_sd, name, error, ill_sd)