diff options
Diffstat (limited to 'jakob_hanika.py')
| -rw-r--r-- | jakob_hanika.py | 50 |
1 files changed, 30 insertions, 20 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 |