Refactoring; use primed c's everywhere

3 files changed, 48 insertions, 35 deletions

diff --git a/jakob_hanika.py b/jakob_hanika.py
index f881afd..f9629d4 100644
--- a/jakob_hanika.py
+++ b/jakob_hanika.py
@@ -6,28 +6,40 @@ from colour.difference import *
 
 # The same wavelength grid is used throughout
 wvl = np.arange(360, 830, 1)
-
-# Map wavelenghts from 360--830 nm to 0-1
-def remap(wvl):
-	return (wvl - 360) / (830 - 360)
+wvlp = (wvl - 360) / (830 - 360)
 
 # This is the model of spectral reflectivity described in the article.
-def model(wvl, cc):
-	x = cc[0] * wvl ** 2 + cc[1] * wvl + cc[2]
-	return 1 / 2 + x / (2 * np.sqrt(1 + x ** 2))
+def model(wvlp, ccp):
+	yy = ccp[0] * wvlp ** 2 + ccp[1] * wvlp + ccp[2]
+	return 1 / 2 + yy / (2 * np.sqrt(1 + yy ** 2))
+
+# Create a SpectralDistribution using given coefficients
+def model_sd(ccp):
+	# 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")
 
 # 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, ill_sd, ill_xy):
+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).
-	ev = model((wvl - 360) / (830 - 360), cc)
-	sd = SpectralDistribution(ev, wvl)
+	sd = model_sd(ccp)
 	Lab = XYZ_to_Lab(sd_to_XYZ(sd, illuminant=ill_sd) / 100, ill_xy)
 	return delta_E_CIE1976(target, Lab)
 
 
 
+# Map primed (dimensionless) coefficients to normal
+# FIXME: is this even useful for anything?
+def cc_from_primed(ccp):
+	return np.array([
+		ccp[0] / 220900,
+		ccp[1] / 470 - 36/11045 * ccp[0],
+		ccp[2] - 36/47 * ccp[1] + 1296/2209 * ccp[0]
+	])
+
+
 # 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).
@@ -35,39 +47,37 @@ def cb_basinhopping(x, f, accept):
 	return f < 0.1
 
 # Finds the parameters for Jakob and Hanika's model
-def jakob_hanika(target, ill_sd, ill_xy, x0=(0, 0, 0)):
+def jakob_hanika(target, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True):
 	# 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
 	opt = optimize.minimize(
-		error_function, x0, (target, ill_sd, ill_xy),
-		method="L-BFGS-B", options={"disp": True, "ftol": 1e-5}
+		error_function, ccp0, (target, ill_sd, ill_xy),
+		bounds=[[-300, 300]] * 3,
+		options={"disp": verbose}
 	)
-	print(opt)
+	if verbose:
+		print(opt)
 
 	error = error_function(opt.x, target, ill_sd, ill_xy)
-	print("Delta E is %g" % error)
+	if verbose:
+		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")
+		if verbose:
+			print("Error too large, trying global optimization")
 		opt = optimize.basinhopping(
 			lambda cc: error_function(cc, target, ill_sd, ill_xy),
-			x0, disp=True, callback=cb_basinhopping
+			x0, disp=verbose, callback=cb_basinhopping
 		)
-		print(opt)
+		if verbose:
+			print(opt)
 		error = error_function(opt.x, target, ill_sd, ill_xy)
-		print("Global delta E is %g" % 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")
+		if verbose:
+			print("Global delta E is %g" % error)
 
-	return cc, matched_sd, error
+	return opt.x, error
diff --git a/test_colorchecker.py b/test_colorchecker.py
index c934518..0368058 100644
--- a/test_colorchecker.py
+++ b/test_colorchecker.py
@@ -2,7 +2,7 @@ import numpy as np
 from colour import *
 from colour.plotting import *
 from matplotlib import pyplot as plt
-from jakob_hanika import jakob_hanika
+from jakob_hanika import jakob_hanika, model_sd
 
 
 
@@ -13,7 +13,7 @@ ill_sd = SpectralDistribution(ILLUMINANTS_SDS[illuminant])
 
 
 # Makes a comparison plot with SDs and swatches
-def plot_comparison(target, matched_sd, label, error, ill_sd):
+def plot_comparison(target, matched_sd, label, error, ill_sd, show=True):
 	if type(target) is SpectralDistribution:
 		target_XYZ = sd_to_XYZ(target, illuminant=ill_sd) / 100
 	else:
@@ -33,7 +33,8 @@ def plot_comparison(target, matched_sd, label, error, ill_sd):
       
 	axes = plt.subplot(2, 1, 2)
 	plt.title("ΔE = %g" % error)
-	plot_multi_colour_swatches([target_swatch, matched_swatch], axes=axes)
+	plot_multi_colour_swatches([target_swatch, matched_swatch],
+	                           standalone=show, axes=axes)
 
 if __name__ == "__main__":
 	# This demo goes through SDs in a color checker
@@ -42,6 +43,7 @@ if __name__ == "__main__":
 		target = XYZ_to_Lab(XYZ, ill_xy)
 
 		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)
+		ccp, error = jakob_hanika(target, ill_sd, ill_xy)
+		matched_sd = model_sd(ccp)
 
 		plot_comparison(sd, matched_sd, name, error, ill_sd)
diff --git a/test_example.py b/test_example.py
index 412cece..f497bbf 100644
--- a/test_example.py
+++ b/test_example.py
@@ -3,7 +3,7 @@ from colour import *
 from colour.plotting import *
 from colour.difference import *
 from matplotlib import pyplot as plt
-from jakob_hanika import jakob_hanika, model, remap, wvl
+from jakob_hanika import jakob_hanika, model_sd
 from test_colorchecker import plot_comparison
 
 
@@ -23,12 +23,13 @@ if __name__ == "__main__":
 	target = XYZ_to_Lab(XYZ, ill_xy)
 
 	print("Target: X=%g, Y=%g, Z=%g, L=%g, a=%g, b=%g" % (*XYZ, *target))
-	cc, found_sd, error = jakob_hanika(target, ill_sd, ill_xy)
+	ccp, error = jakob_hanika(target, ill_sd, ill_xy)
 
-	reference_sd = SpectralDistribution(model(remap(wvl), cc_ref), wvl)
+	reference_sd = model_sd(cc_ref)
 	reference_XYZ = sd_to_XYZ(reference_sd, illuminant=ill_sd)
 	reference_Lab = XYZ_to_Lab(reference_XYZ, ill_xy)
 
+	found_sd = model_sd(ccp)
 	found_XYZ = sd_to_XYZ(found_sd, illuminant=ill_sd)
 	found_Lab = XYZ_to_Lab(found_XYZ, ill_xy)