From 95fdb61488202942bd8496c34c2615643ccc84fd Mon Sep 17 00:00:00 2001 From: Paweł Redman Date: Tue, 9 Jun 2020 12:13:44 +0200 Subject: 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. --- jakob_hanika.py | 50 ++++++++++++++++++++++++++++++-------------------- test_colorchecker.py | 38 ++++++++++++++++++++++++-------------- 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) -- cgit