Rewrite the error function w/ analytic derivatives

This greatly speeds up the code and improves convergence. Something's slightly off with the code, though, and the output colors are slightly skewed.

3 files changed, 127 insertions, 15 deletions

diff --git a/gsoc_common.py b/gsoc_common.py
index b028e45..5d4aea0 100644
--- a/gsoc_common.py
+++ b/gsoc_common.py
@@ -9,8 +9,8 @@ D65_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"]["D65"]
 D65 = SpectralDistribution(ILLUMINANT_SDS["D65"])
 
 # The same wavelength grid is used throughout
-wvl = np.arange(360, 830, 5)
-wvlp = (wvl - 360) / (830 - 360)
+wvl = SpectralShape(360, 830, 5)
+wvlp = (wvl.range() - 360) / (830 - 360)
 
 # This is the model of spectral reflectivity described in the article.
 def model(wvlp, ccp):
@@ -21,27 +21,87 @@ def model(wvlp, 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
-	grid = wvlp if primed else wvl
-	return SpectralDistribution(model(grid, ccp), wvl, name="Model")
+	grid = wvlp if primed else wvl.range()
+	return SpectralDistribution(model(grid, ccp), wvl.range(), name="Model")
 
 # The goal is to minimize the color difference between a given distrbution
 # and the one computed from the model above.
-def error_function(ccp, 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).
-	sd = model_sd(ccp)
-	Lab = XYZ_to_Lab(sd_to_XYZ(sd, illuminant=ill_sd) / 100, ill_xy)
-	return delta_E_CIE1976(target, Lab)
+# This function also calculates the first derivatives with respect to c's.
+def error_function(c, target, wvl, cmfs, illuminant, illuminant_XYZ):
+	U = c[0] * wvl**2 + c[1] * wvl + c[2]
+	t1 = np.sqrt(1 + U**2)
+	R = 1 / 2 + U / (2 * t1)
 
+	t2 = 1 / (2 * t1) - U**2 / (2 * t1**3)
+	dR_dc = [wvl**2 * t2, wvl * t2, t2]
 
+	E = illuminant.values * R / 100
+	dE_dc = illuminant.values * dR_dc / 100
+
+	XYZ = np.empty(3)
+	dXYZ_dc = np.empty((3, 3))
+
+	dlambda = cmfs.wavelengths[1] - cmfs.wavelengths[0]
+	for i in range(3):
+		XYZ[i] = np.dot(E, cmfs.values[:, i]) * dlambda
+		for j in range(3):
+			dXYZ_dc[i, j] = np.dot(dE_dc[j], cmfs.values[:, i]) * dlambda
+
+	# FIXME: this isn't the full CIE 1976 lightness function
+	f = (XYZ / illuminant_XYZ)**(1/3)
+
+	# FIXME: this can be vectorized
+	df_dc = np.empty((3, 3))
+	for i in range(3):
+		for j in range(3):
+			df_dc[i, j] = 1 / (3 * illuminant_XYZ[i]**(1/3)
+			                   * XYZ[i]**(2/3)) * dXYZ_dc[i, j]
+
+	Lab = np.array([
+		116 * f[1] - 16,
+		500 * (f[0] - f[1]),
+		200 * (f[1] - f[2])
+	])
+
+	dLab_dc = np.array([
+		116 * df_dc[1],
+		500 * (df_dc[0] - df_dc[1]),
+		200 * (df_dc[1] - df_dc[2])
+	])
+
+	error = np.sqrt(np.sum((Lab - target)**2))
+
+	derror_dc = np.zeros(3)
+	for i in range(3):
+		for j in range(3):
+			derror_dc[i] += dLab_dc[j, i] * (Lab[j] - target[j])
+		derror_dc[i] /= error
+
+	#print("c=[%.3g, %.3g, %.3g], XYZ=[%.3g, %.3g, %.3g], Lab=[%.3g, %.3g, %.3g], error=%g" % (*c, *XYZ, *Lab, error))
+
+	return error, derror_dc
 
 # Finds the parameters for Jakob and Hanika's model
 def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_hard=True):
+	ill_sd = ill_sd.align(wvl)
+	ill_XYZ = sd_to_XYZ(ill_sd) / 100
+
 	def do_optimize(XYZ, ccp0):
 		Lab = XYZ_to_Lab(XYZ, ill_xy)
+
+		def fun(*args):
+			error, derror_dc = error_function(*args)
+			return error
+
+		def jac(*args):
+			error, derror_dc = error_function(*args)
+			return derror_dc
+
+		cmfs = STANDARD_OBSERVER_CMFS["CIE 1931 2 Degree Standard Observer"].align(wvl)
+
 		opt = optimize.minimize(
-			error_function, ccp0, (Lab, ill_sd, ill_xy),
-			method="Nelder-Mead", options={"disp": verbose}
+			fun, ccp0, (Lab, wvlp, cmfs, ill_sd, ill_XYZ), jac=jac,
+			method="L-BFGS-B", options={"disp": verbose}
 		)
 		if verbose:
 			print(opt)
@@ -61,7 +121,7 @@ def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_h
 	good_ccp = (2.1276356, -1.07293026, -0.29583292) # FIXME: valid only for D65
 
 	divisions = 3
-	while divisions < 10:
+	while divisions < 30:
 		if verbose:
 			print("Trying with %d divisions" % divisions)
 
@@ -76,7 +136,7 @@ def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_h
 				print("Intermediate step %d/%d, XYZ=%s with ccp0=%s" %
 				      (i + 1, divisions, intermediate_XYZ, ccp0))
 			opt = do_optimize(intermediate_XYZ, ccp0)
-			if opt.fun > 0.1:
+			if opt.fun > 1e-3:
 				if verbose:
 					print("WARNING: intermediate optimization failed")
 				break
diff --git a/test_diff.py b/test_diff.py
new file mode 100644
index 0000000..b634ca6
--- /dev/null
+++ b/test_diff.py
@@ -0,0 +1,51 @@
+import numpy as np
+from scipy.optimize import minimize
+from colour import *
+from colour.colorimetry import STANDARD_OBSERVER_CMFS, ILLUMINANT_SDS
+from colour.models import eotf_inverse_sRGB, sRGB_to_XYZ
+from matplotlib import pyplot as plt
+from gsoc_common import plot_comparison, error_function, model_sd, D65_xy
+
+shape = SpectralShape(360, 830, 1)
+cmfs = STANDARD_OBSERVER_CMFS["CIE 1931 2 Degree Standard Observer"].align(shape)
+
+illuminant = SpectralDistribution(ILLUMINANT_SDS["D65"]).align(shape)
+illuminant_XYZ = sd_to_XYZ(illuminant) / 100
+wvl = np.linspace(0, 1, len(shape.range()))
+
+target = np.array([50, -20, 30]) # Some arbitrary Lab coordinates
+xs = np.linspace(-10, 10, 500)
+h = xs[1] - xs[0]
+
+# This test checks if derivatives are calculated correctly by comparing them
+# to finite differences.
+for c_index in range(3):
+	errors = np.empty(len(xs))
+	derrors = np.empty(len(xs))
+
+	for i, x in enumerate(xs):
+		c = np.array([1.0, 1, 1])
+		c[c_index] = x
+
+		error, derror_dc = error_function(
+			c, target, wvl, cmfs, illuminant, illuminant_XYZ
+		)
+
+		errors[i] = error
+		derrors[i] = derror_dc[c_index]
+
+
+	plt.subplot(2, 3, 1 + c_index)
+	plt.xlabel("c%d" % c_index)
+	plt.ylabel("ΔE")
+	plt.plot(xs, errors)
+
+	plt.subplot(2, 3, 4 + c_index)
+	plt.xlabel("c%d" % c_index)
+	plt.ylabel("dΔE/dc%d" % c_index)
+
+	plt.plot(xs, derrors, "k-")
+	plt.plot(xs[:-1] + h / 2, np.diff(errors) / h, "r:")
+
+
+plt.show()
diff --git a/test_example.py b/test_example.py
index 6be982d..80b07c7 100644
--- a/test_example.py
+++ b/test_example.py
@@ -13,7 +13,8 @@ cc_ref = np.array([ 18.70184886, -13.19804478, 2.12180137])
 
 
 if __name__ == "__main__":
-	XYZ = sRGB_to_XYZ(eotf_inverse_sRGB(RGB_ref), D65_xy)
+	#XYZ = sRGB_to_XYZ(eotf_inverse_sRGB(RGB_ref), D65_xy)
+	XYZ = np.array([1/3, 1/4, 1/3])
 	Lab = XYZ_to_Lab(XYZ, D65_xy)
 
 	print("Target: X=%g, Y=%g, Z=%g, L=%g, a=%g, b=%g" % (*XYZ, *Lab))