1 files changed, 97 insertions, 0 deletions
diff --git a/jakob_hanika.py b/jakob_hanika.py
new file mode 100644
index 0000000..78803a1
--- /dev/null
+++ b/jakob_hanika.py
@@ -0,0 +1,97 @@
+import numpy as np, scipy.optimize as optimize
+from colour import *
+from colour.difference import *
+from colour.plotting import *
+from matplotlib import pyplot as plt
+
+
+
+# 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)
+	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_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)
+      
+	axes = plt.subplot(2, 1, 2)
+	plt.title("ΔE = %g" % error)
+	plot_multi_colour_swatches([target_swatch, matched_swatch], axes=axes)
+
+
+
+# The same illuminant is used throughout
+il = ILLUMINANTS['CIE 1931 2 Degree Standard Observer']['D65']
+
+# The same wavelength grid is used throughout
+wvl = np.arange(360, 830, 10)
+
+
+
+# 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]
+	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):
+	ev = model(wvl, cc)
+	sd = SpectralDistribution(ev, wvl)
+	Lab = XYZ_to_Lab(sd_to_XYZ(sd), il)
+	return delta_E_CIE1976(target, Lab)
+
+
+
+# This callback to scipy.optimize.basinhopping tells the solver to stop once
+# the error is small enough. The threshold was chosen arbitrarily, as a small
+# fraction of the JND (about 2.3 with this metric).
+def cb_basinhopping(x, f, accept):
+	return f < 0.1
+
+# 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, il)
+
+	print("The target is '%s' with L=%g, a=%g, b=%g" % (name, *target))
+
+	# 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 (instead of
+	#        going).
+	opt = optimize.minimize(
+		error_function, (0, 0, 0), target, method="L-BFGS-B",
+		options={"disp": True, "ftol": 1e-5}
+	)
+	print(opt)
+
+	error = error_function(opt.x, target)
+	print("Delta E is %g" % error)
+
+	# Basin hopping is far more likely to find the actual minimum we're
+	# looking for, but it's extremely slow in comparison.
+	if error > 0.1:
+		print("Error too large, trying global optimization")
+		opt = optimize.basinhopping(
+			lambda cc: error_function(cc, target),
+			(0, 0, 0), disp=True, callback=cb_basinhopping
+		)
+		print(opt)
+		error = error_function(opt.x, target)
+		print("Global delta E is %g" % error)
+
+	matched_sd = SpectralDistribution(model(wvl, opt.x), wvl, name="Model")
+	plot_comparison(sd, matched_sd, name, error)