Major cleanup

- The interpolator code is no longer here (but in the Colour fork)
- Refactoring made the code shorter and cleaner
- All tests work properly again

8 files changed, 86 insertions, 153 deletions

diff --git a/jakob_hanika.py b/gsoc_common.py
index 22b87fc..b028e45 100644
--- a/jakob_hanika.py
+++ b/gsoc_common.py
@@ -2,8 +2,11 @@ import numpy as np, scipy.optimize as optimize
 from colour import *
 from colour.difference import delta_E_CIE1976
 from colour.colorimetry import *
+from colour.plotting import *
+from matplotlib import pyplot as plt
 
-
+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)
@@ -91,3 +94,27 @@ def jakob_hanika(target_XYZ, ill_sd, ill_xy, ccp0=(0, 0, 0), verbose=True, try_h
 
 	raise Exception("optimization failed for target_XYZ=%s, ccp0=%s" \
 	                % (target_XYZ, ccp0))
+
+# Makes a comparison plot with SDs and swatches
+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:
+		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, 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)
+	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],
+	                           standalone=show, axes=axes)
diff --git a/test_bad_convergence.py b/test_bad_convergence.py
index dadb21e..a47af5a 100644
--- a/test_bad_convergence.py
+++ b/test_bad_convergence.py
@@ -1,15 +1,6 @@
 import numpy as np
 from colour import *
-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(ILLUMINANT_SDS[illuminant])
-
-
+from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison
 
 # A set of inputs that make the Nelder-Mead solver work really hard.
 # This might not be interesting in future versions of jakob_hanika.py, once
@@ -19,6 +10,5 @@ if __name__ == "__main__":
 	XYZ = np.array([0.00263241, 0.00174905, 0.00092602])
 	ccp0 = np.array([5.11263599, -2.65344447, -2.37301856])
 
-	ccp, error = jakob_hanika(XYZ, ill_sd, ill_xy, ccp0)
-	if error > 0.1:
-		plot_comparison(XYZ, model_sd(ccp), "Target", error, ill_sd)
+	ccp, error = jakob_hanika(XYZ, D65, D65_xy, ccp0)
+	plot_comparison(XYZ, model_sd(ccp), "Target", error, D65)
diff --git a/test_coeff.py b/test_coeff.py
deleted file mode 100644
index c2a8970..0000000
--- a/test_coeff.py
+++ /dev/null
@@ -1,68 +0,0 @@
-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 colour.difference import delta_E_CIE1976
-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(ILLUMINANT_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 = [[0, 1, 2], self.scale, t, t]
-		self.cubes = RegularGridInterpolator(axes, coeffs[:, :, :, :, :], bounds_error=False)
-
-	def __call__(self, RGB):
-		RGB = np.asarray(RGB, dtype=np.float) # FIXME: default dtype
-		vmax = np.max(RGB, axis=-1)
-		imax = np.argmax(RGB, axis=-1)
-		chroma = RGB / (np.expand_dims(vmax, -1) + 1e-10) # Avoid division by zero
-		vmax = np.max(RGB, axis=-1)
-		v2 = np.take_along_axis(chroma, np.expand_dims((imax + 2) % 3, axis=-1), axis=-1).squeeze(axis=-1)
-		v3 = np.take_along_axis(chroma, np.expand_dims((imax + 1) % 3, axis=-1), axis=-1).squeeze(axis=-1)
-		coords = np.stack([imax, vmax, v2, v3], axis=-1)
-		return self.cubes(coords).squeeze()
-
-
-
-if __name__ == "__main__":
-	jh = JakobHanikaInterpolator()
-	jh.read_file("data/srgb.coeff")
-
-	RGBs = np.random.random((7, 6, 5, 4, 3))
-	ccs = jh(RGBs)
-
-	RGB = eotf_inverse_sRGB(RGBs[0, 0, 0, 0, :])
-	XYZ = sRGB_to_XYZ(RGB)
-	Lab = XYZ_to_Lab(XYZ, ill_xy)
-	cc = ccs[0, 0, 0, 0, :]
-
-	matched_sd = model_sd(cc, primed=False)
-	matched_XYZ = sd_to_XYZ(matched_sd, illuminant=ill_sd)
-	matched_Lab = XYZ_to_Lab(matched_XYZ, ill_xy)
-	error = delta_E_CIE1976(Lab, matched_Lab)
-
-	print(matched_sd)
-
-	plot_comparison(XYZ, matched_sd, "Model", error, ill_sd, ill_xy)
diff --git a/test_colorchecker.py b/test_colorchecker.py
index b0c9688..fc54d84 100644
--- a/test_colorchecker.py
+++ b/test_colorchecker.py
@@ -1,48 +1,15 @@
 import numpy as np
 from colour import *
-from colour.plotting import *
-from matplotlib import pyplot as plt
-from jakob_hanika import jakob_hanika, model_sd
-
-
-
-illuminant = "D65"
-ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant]
-ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant])
-
-
-
-# Makes a comparison plot with SDs and swatches
-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:
-		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, 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)
-	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],
-	                           standalone=show, axes=axes)
+from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison
 
+# This demo goes through SDs in a color checker
 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, illuminant=ill_sd) / 100
+	for name, sd in COLOURCHECKER_SDS['ColorChecker N Ohta'].items():
+		XYZ = sd_to_XYZ(sd, illuminant=D65) / 100
 
-		print("Color checker: The target is '%s' with X=%g, Y=%g, Z=%g" % (name, *XYZ))
-		ccp, error = jakob_hanika(XYZ, ill_sd, ill_xy)
+		print("Color checker: The target is '%s' with X=%g, Y=%g, Z=%g"
+		      % (name, *XYZ))
+		ccp, error = jakob_hanika(XYZ, D65, D65_xy)
 		matched_sd = model_sd(ccp)
 
-		plot_comparison(sd, matched_sd, name, error, ill_sd)
+		plot_comparison(sd, matched_sd, name, error, D65)
diff --git a/test_coverage.py b/test_coverage.py
index e73e015..90e3c79 100644
--- a/test_coverage.py
+++ b/test_coverage.py
@@ -1,14 +1,10 @@
 import numpy as np, os, multiprocessing
 from colour import *
 from colour.models import eotf_inverse_sRGB
-from jakob_hanika import jakob_hanika, model_sd
+from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison
 
 
 
-illuminant = "D65"
-ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant]
-ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant]) # Note: changed in dev (to ILLUMINANT_SD)
-
 # Resolution of the discretization cubes
 # Both should be 64, but this code is a bit too slow at the moment
 CHROMA_STEPS = 8
@@ -19,10 +15,10 @@ LIGHTNESS_STEPS = 8
 # Solve for a specific RGB color
 def optimize_RGB(linear_RGB, ccp0):
 	RGB = eotf_inverse_sRGB(linear_RGB)
-	target = sRGB_to_XYZ(RGB, ill_xy)
-	ccp, error = jakob_hanika(target, ill_sd, ill_xy, ccp0, verbose=False)
+	target = sRGB_to_XYZ(RGB, D65_xy)
+	ccp, error = jakob_hanika(target, D65, D65_xy, ccp0, verbose=False)
 	if error > 0.1:
-		print("WARNING: convergence failed with L=%s, starting from %s" \
+		print("WARNING: convergence failed with L=%s, starting from %s"
 		      % (linear_RGB, ccp0))
 
 	return ccp, error
@@ -90,7 +86,8 @@ if __name__ == "__main__":
 				dirname = "out/R%.3fG%.3fB%.3f" % (*RGB,)
 				total += 1
 
-				if os.path.exists(dirname) and len(os.listdir(dirname)) == LIGHTNESS_STEPS:
+				if os.path.exists(dirname) and \
+				   len(os.listdir(dirname)) == LIGHTNESS_STEPS:
 					done += 1
 					continue
 
diff --git a/test_example.py b/test_example.py
index f2e5345..6be982d 100644
--- a/test_example.py
+++ b/test_example.py
@@ -1,15 +1,11 @@
 import numpy as np
 from colour import *
 from colour.difference import delta_E_CIE1976
-from jakob_hanika import jakob_hanika, model_sd
-from test_colorchecker import plot_comparison
+from colour.models import eotf_inverse_sRGB
+from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison
 
 
 
-illuminant = "D65"
-ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant]
-ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant])
-
 # These numbers are taken from Jakob and Hanika's Jupyter notebook.
 RGB_ref = np.array([0.79264853,  0.4, 0.63703843]) # *linear* sRGB
 cc_ref = np.array([ 18.70184886, -13.19804478, 2.12180137])
@@ -17,21 +13,21 @@ cc_ref = np.array([ 18.70184886, -13.19804478, 2.12180137])
 
 
 if __name__ == "__main__":
-	XYZ = sRGB_to_XYZ(1.055 * RGB_ref ** (1/2.4) - 0.055, ill_xy)
-	Lab = XYZ_to_Lab(XYZ, ill_xy)
+	XYZ = sRGB_to_XYZ(eotf_inverse_sRGB(RGB_ref), D65_xy)
+	Lab = XYZ_to_Lab(XYZ, D65_xy)
 
 	print("Target: X=%g, Y=%g, Z=%g, L=%g, a=%g, b=%g" % (*XYZ, *Lab))
-	ccp, error = jakob_hanika(XYZ, ill_sd, ill_xy)
+	ccp, error = jakob_hanika(XYZ, D65, D65_xy)
 
 	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)
+	reference_XYZ = sd_to_XYZ(reference_sd, illuminant=D65)
+	reference_Lab = XYZ_to_Lab(reference_XYZ, D65_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)
+	found_XYZ = sd_to_XYZ(found_sd, illuminant=D65)
+	found_Lab = XYZ_to_Lab(found_XYZ, D65_xy)
 
 	print("Our results differ from the reference by ΔE = %g" \
 	      % delta_E_CIE1976(reference_Lab, found_Lab))
 
-	plot_comparison(XYZ, found_sd, "Reference", error, ill_sd)
+	plot_comparison(XYZ, found_sd, "Reference", error, D65)
diff --git a/test_integration.py b/test_integration.py
index 35c19ed..e6aa819 100644
--- a/test_integration.py
+++ b/test_integration.py
@@ -1,19 +1,16 @@
 import numpy as np
 from colour import *
+from colour.models import eotf_inverse_sRGB
 from colour.recovery import Jakob2019Interpolator
 from colour.plotting import plot_multi_sds
+from gsoc_common import D65_xy
 
 
 
-illuminant = "D65"
-ill_xy = ILLUMINANTS["CIE 1931 2 Degree Standard Observer"][illuminant]
-ill_sd = SpectralDistribution(ILLUMINANT_SDS[illuminant])
-
-
-
+# This is a basic test of the feature/jakob2019 branch
 if __name__ == "__main__":
 	RGB = np.array([0.79264853,  0.4, 0.63703843])
-	XYZ = sRGB_to_XYZ(1.055 * RGB ** (1/2.4) - 0.055, ill_xy)
+	XYZ = sRGB_to_XYZ(eotf_inverse_sRGB(RGB), D65_xy)
 	sd = XYZ_to_sd(XYZ, method="Jakob 2019", verbose=True)
 
 	interp = Jakob2019Interpolator()
diff --git a/test_interpolator.py b/test_interpolator.py
new file mode 100644
index 0000000..21a50c5
--- /dev/null
+++ b/test_interpolator.py
@@ -0,0 +1,27 @@
+import numpy as np, struct
+from colour import *
+from colour.difference import delta_E_CIE1976
+from colour.models import eotf_inverse_sRGB
+from colour.recovery import Jakob2019Interpolator
+from gsoc_common import D65, D65_xy, jakob_hanika, model_sd, plot_comparison
+
+# This script tests if the interpolator correctly handles multi-dimensional
+# inputs.
+if __name__ == "__main__":
+	interp = Jakob2019Interpolator()
+	interp.read_file("data/srgb.coeff")
+
+	RGBs = np.random.random((7, 6, 5, 4, 3))
+	ccs = interp.coeffs(RGBs)
+
+	RGB = eotf_inverse_sRGB(RGBs[0, 0, 0, 0, :])
+	XYZ = sRGB_to_XYZ(RGB)
+	Lab = XYZ_to_Lab(XYZ, D65_xy)
+	cc = ccs[0, 0, 0, 0, :]
+
+	matched_sd = model_sd(cc, primed=False)
+	matched_XYZ = sd_to_XYZ(matched_sd, illuminant=D65)
+	matched_Lab = XYZ_to_Lab(matched_XYZ, D65_xy)
+	error = delta_E_CIE1976(Lab, matched_Lab)
+
+	plot_comparison(XYZ, matched_sd, "Model", error, D65)