Given a light responder (e.g. an eye or a camera), two light spectra that produce the same response from the responder are called metamers for that responder. Similarly, given a material responder (e.g. a scanner), two reflectance spectra that produce the same response from the responder are called metamers for that responder.
For a given responder and response,
there are typically infinitely many metamers.
The set of all of them is often called the metameric suite.
The goal of the function invert() is to calculate a "good" metamer
in the "suite".
Koenderink calls this topic inverse colorimetry.
In the case that the estimated spectrum is a reflectance spectrum,
the topic is often called reflectance estimation or reflectance recovery,
see Bianco.
The centroid method, which is the default and the featured method in this package, computes the centroid of the set of all the metamers (if any). The centroid is computed in an infinite-dimensional context and is expounded further in Davis.
The Hawkyard method, see Hawkyard and Bianco, has been around a long time. The centroid and Hawkyard methods have similarities, e.g. both are low-dimensional with the number of variables equal to the number of responses (usually 3). The Hawkyard method is very fast, but has a key problem, see below.
The Transformed Least Slope Squared (TLSS) method was developed
by Scott Burns, see References.
This is my name for it, not his.
What I call TLLS is actually is a combination of Burns' LHTSS and LLSS methods;
the one that invert() chooses depends on type(x), see below.
Both of these are high-dimensional,
with the number of variables equal to #(responses) + #(wavelengths).
The first argument to invert() is the responder x,
and the second is the matrix response of responses (e.g. XYZs or RGBs).
The goal is to return a "good" spectrum for each response so that:
product( invert(x,response), x ) \(~\cong~\) response |
The error is returned as column estim.precis, see below.
First consider the case where x has type type='responsivity.material'.
The goal is to compute a reflectance spectra.
All the methods will fail if the response is on the object-color boundary
(an optimal color) or outside the boundary.
They may also fail if the response is inside the object-color
solid (the Rösch Farbkörper) and very close to the boundary.
The centroid method solves a non-linear system that contains a
Langevin-function-based squashing function, see Davis for details.
When successful it always returns a feasible spectrum
with small estim.precis.
The Hawkyard method is linear and very fast,
but in raw form it may return a non-feasible reflectance spectrum.
In this case invert() simply clamps to the interval [0,1] and so
estim.precis can be large.
The TLSS method solves a non-linear system that contains
the squashing function \((\tanh(z) + 1)/2\), see Burns for details.
When successful it always returns a feasible spectrum
with small estim.precis.
Now consider the case where x has type='responsivity.light'.
The goal is to compute the spectrum of a light source.
All the methods will fail if chromaticity of the response is on the boundary
of the inverted-U (assuming x models the human eye) or outside the boundary.
They may also fail if the response is inside the inverted-U
and very close to the boundary.
The centroid method works on a relatively small range of chromaticities;
it will fail if the response is too far from the response to Illuminant E.
See Davis for the details.
When successful it always returns an everywhere positive spectrum
with small estim.precis.
This method has the desirable property that if the response is multiplied by
a positive number, the computed spectrum is multiplied by that same number.
The Hawkyard method does not work in this case.
The TLSS method solves a non-linear system that contains
the squashing function \(\exp(z)\), see Burns for the details.
When successful it always returns an everywhere positive spectrum
with small estim.precis.
This method succeeds on a larger set of chromaticities than the centroid method.
It also has the desirable scaling multiplication property mentioned above.
The centroid and Hawkyard methods have an equalization option,
which is controlled by the argument alpha
and is enabled by default, see below.
When enabled, if the response comes from a constant spectrum
(a perfectly neutral gray material, or a multiple of Illuminant E),
then the computed spectrum is that same constant spectrum (up to numerical precision).
I call this the neutral-exact property.
Equalization is a complicated mechanism, for details see Davis.
For the TLSS method, the neutral-exact property is intrinsic,
and alpha is ignored.
NOTE:
If the responder has only one output channel (e.g. a monochrome camera)
and equalization is enabled,
then all responses are inverted to a constant spectrum.
This may or may not be desirable.
# S3 method for colorSpec
invert( x, response, method='centroid', alpha=1 )If type(x)='responsivity.material' it returns a colorSpec object
with type = 'material' (quantity = 'reflectance').
If type(x)='responsivity.light' it returns a colorSpec object
with type = 'light'
(quantity='energy' or quantity='photons'
depending on quantity(x)).
In either case, the returned object has organization = 'df.row'
and the extradata is a data.frame with these columns:
the input matrix of desired responses
the difference between the desired response and actual response. It is the mean of the absolute value of the differences.
See rootSolve::multiroot()
the time to compute the spectrum, in msec.
When method='Hawkyard', all N spectra are computed at once,
so all N spectra are assigned the same mean time.
the number of iterations that were required to find the relevant root.
This is not present when method='Hawkyard'.
a logical indicating whether the reflectance was clamped to [0,1]. This is present only when method='Hawkyard'.
If a response could not be estimated,
the row contains NA in appropriate columns,
and a warning is issued.
In case of global error, the function returns NULL.
a colorSpec object with type(x) =
'responsivity.material' or 'responsivity.light' and M responsivities.
The wavelengths must be regular (equidistant).
a numeric NxM matrix, or a numeric vector that can be converted
to such matrix, by row. The N responses are contained in the rows.
The rownames(response) are copied to the output specnames.
either 'centroid' or 'Hawkyard' or 'TLSS'.
'Hawkyard' is only valid when
type(x) is 'responsivity.material'.
Matching is partial and case-insensitive.
a vector of M weighting coefficients,
or a single number that is replicated to length M.
When method='centroid', alpha is used for equalizing
the responsivities, which is recommended.
For alpha to be valid, the linear combination of the M responsitivies,
with coefficients alpha, must be positive.
To disable equalization (not recommended) and use the original responsivities,
set alpha=NULL.
Similarly, when method='Hawkyard', alpha is used for equalizing
the responsivities, which is also recommended.
When method='TLSS', alpha is ignored.
If type(x)='responsivity.light' the centroid method may fail
(not converge) if the response is too far from that of Illuminant E.
For method='centroid' the function calls the non-linear root-finder
rootSolve::multiroot(),
which is general purpose and "full Newton".
For method='Hawkyard' the function solves a linear system by
inverting a small matrix (#[responses] x #[responses]).
The spectra are then clamped to [0,1].
For method='TLSS' the function solves a constrained least-squares problem
using Lagrange multipliers.
A critical point is found using a "full Newton" iteration.
The original MATLAB code is posted at Burns,
and was ported from MATLAB to R with only trivial changes.
When computing a reflectance spectrum, the Hawkyard method is used for the
initial guess, after little extra clamping.
This improved guess cuts the number of iterations substantially,
and the extra computation time is negligible.
Davis, Glenn. A Centroid for Sections of a Cube in a Function Space, with Application to Colorimetry. https://arxiv.org/abs/1811.00990. [math.FA]. 2018.
Bianco, Simone. Reflectance spectra recovery from tristimulus values by adaptive estimation with metameric shape correction. vol. 27, no 8. Journal of the Optical Society of America A. pages 1868-1877. 2010 https://opg.optica.org/josaa/abstract.cfm?uri=josaa-27-8-1868.
Burns, Scott A. Generating Reflectance Curves from sRGB Triplets. http://scottburns.us/reflectance-curves-from-srgb/.
Hawkyard, C. J. Synthetic reflectance curves by additive mixing. Journal of the Society of Dyers and Colourists. vol. 109. no. 10. Blackwell Publishing Ltd. pp. 323-329. 1993.
Koenderink, J.J. Color for the Sciences. MIT Press. 2010.
type(),
quantity(),
organization(),
specnames(),
product(),
extradata(),
rootSolve::multiroot(),
vignette Estimating a Spectrum from its Response
wave = 400:700
E.eye = product( illuminantE(1,wave), "material", xyz1931.1nm, wavelength=wave )
path = system.file( 'extdata/targets/CC_Avg30_spectrum_CGATS.txt', package='colorSpec' )
MacbethCC = readSpectra( path, wavelength=wave )
XYZ = product( MacbethCC, E.eye, wavelength=wave )
est.eq = invert( E.eye, XYZ, method='centroid', alpha=1 )
extra = extradata(est.eq)
range(extra$estim.precis) # prints 0.000000e+00 3.191741e-08
Run the code above in your browser using DataLab