probeOptimalColors: compute optimal colors by ray tracing

Description

Consider a colorSpec object x with type equal to 'responsivity.material'. The set of all possible material reflectance functions (or transmittance functions) is convex, closed, and bounded (in any reasonable function space), and this implies that the set of all possible output responses from x is also convex, closed, and bounded. The latter set is called the object-color solid or Rösch Farbkörper for x. A color on the boundary of the object-color solid is called an optimal color. The special points W (the response to the perfect reflecting diffuser) and 0 are on the boundary of this set. The interior of the line segment of neutrals joining 0 to W is in the interior of the object-color solid. It is natural to parameterize this segment from 0 to 1 (from 0 to W).

A ray \(r\) that is based at a point on the interior of the neutral line segment must intersect the boundary of the object-color solid in a unique optimal color. The purpose of the function probeOptimalColors() is to compute that intersection point.

Currently the function only works if the number of spectra in x is 3 (e.g. RGB or XYZ).

Before colorSpec v 0.8-1 this function used a 2D root-finding method that could only find optimal colors whose spectra contain 0, 1, or 2 transitions. But starting with v0.8-1, we have switched to zonohedral representation of the object-color solid, which makes it possible to discover more than 2 transitions. The inspiration for this change is the article by Centore. To inspect these computed spectra, the argument spectral must be set to TRUE.

Usage

# S3 method for colorSpec
probeOptimalColors( x, gray, direction, aux=FALSE, spectral=FALSE, tol=1.e-6 )

Value

If argument spectral=FALSE, probeOptimalColors() returns a data.frame with a row for each traced ray. There are length(gray) * nrow(direction) rays. The columns in the output are:

gray: the graylevel defining the \(basepoint\) of the ray. \(basepoint = gray*W\)
direction: the \(direction\) of the ray
s: computed scalar so that \(basepoint + s*direction\) is optimal
optimal: the optimal color on the boundary; \(optimal = basepoint + s*direction\)
lambda: lambda.1 and lambda.2 at the 2 transitions, in nm. lambda.1 < lambda.2 => bandpass, and lambda.1 > lambda.2 => bandstop. It will happen that the optimal spectrum has more than 2 transitions; in this case both lambdas are set to NA.
dol: delta and omega - the Logvinenko parameters \((\delta,\omega)\) for optimal colors, plus lambda (\(\lambda\)) in nm. \(\omega\) is the reparameterization of \(\lambda\) ; see Logvinenko. If there are more than 2 transistions, these are set to NA.

If aux is TRUE, these auxiliary columns related to performance and diagnostics are added:

timetrace: time to trace the ray, in seconds
parallelograms: # of parallelograms in the (possibly compound) face. 1 means just a single parallelogram.
tested: # of parallelograms actually tested for ray intersection. This only has meaning for compound faces.
alpha: the 2 coordinates of the intersection point inside the parallelogram

If argument spectral=TRUE, probeOptimalColors() returns a colorSpec object with quantity 'reflectance'. This object contains the optimal spectra, and can be used to inspect the spectra with more than 2 transitions, which will happen. The above-mentioned data.frame can then be obtained by applying extradata() to the returned object.

If an individual ray could not be traced (which should be rare), the row contains NA in appropriate columns.

In case of global error, the function returns NULL.

Arguments

x: a colorSpec object with type equal to 'responsivity.material' and 3 spectra
gray: vector of numbers in the open interval (0,1) that define neutral grays on the line segment from black to white; this neutral gray point is the basepoint of a probe ray
direction: a numeric Nx3 matrix with directions of the probe rays in the rows, or a numeric vector that can be converted to such a matrix, by row.
aux: a logical that specifies whether to return extra performance and diagnostic data; see Details
spectral: if TRUE, the function returns a colorSpec object with the optimal spectra, see Value.
tol: error tolerance for the intersection of probe and object-color boundary

WARNING

The preprocessing calculation of the zonohedron dominates the total time. And this time goes up rapidly with the number of wavelengths. We recommend using a wavelength step of 5nm, as in the Examples. For best results, batch a lot of rays into a single function call and then process the output.
Moreover, the preprocessing time is dominated by the partitioning of the compound faces into parallelograms. This is made worse by spectral responses with little overlap, as in scanner.ACES. In these cases, try a larger step size, and then reduce. Optimizing these compound faces is a possible topic for the future.

Details

Each gray level and each direction defines a ray. So the total number of rays traced is length(gray) * nrow(direction). The 3 responsivities are regarded not as continuous functions, but as step functions. This implies that the color solid is a zonohedron. In the preprocessing phase the zonohedral representation is calculated. The faces of the zonohedron are either parallelograms, or compound faces that can be partitioned into parallelograms. The centers of all these parallelograms are computed, along with their normals and plane constants.
This representation of the color solid is very strict regarding the 2-transition assumption. During use, one can count on there being some spectra with more than two transitions. Forcing the best 2-transition spectrum is a possible topic for the future.

References

Centore, Paul. A zonohedral approach to optimal colours. Color Research & Application. Vol. 38. No. 2. pp. 110-119. April 2013.

Logvinenko, A. D. An object-color space. Journal of Vision. 9(11):5, 1-23, (2009).
https://jov.arvojournals.org/article.aspx?articleid=2203976. doi:10.1167/9.11.5.

Schrödinger, E. (1920). Theorie der Pigmente von grösster Leuchtkraft. Annalen der Physik. 62, 603-622.

West, G. and M. H. Brill. Conditions under which Schrödinger object colors are optimal. Journal of the Optical Society of America. 73. pp. 1223-1225. 1983.

Examples

Run this code

wave    = seq(400,700,by=5)
D50.eye = product( D50.5nm, 'material', xyz1931.1nm, wavelength=wave )
probeOptimalColors( D50.eye, c(0.2,0.5,0.9), c(1,2,1, -1,-2,-1) )

##    gray direction.1 direction.2 direction.3         s  optimal.1  optimal.2
##  1  0.2           1           2           1 32.306207  52.533143  85.612065
##  2  0.2          -1          -2          -1  8.608798  11.618138   3.782055
##  3  0.5           1           2           1 20.993144  71.560483  94.485416
##  4  0.5          -1          -2          -1 20.993144  29.574196  10.512842
##  5  0.9           1           2           1  4.333700  95.354911 103.165832
##  6  0.9          -1          -2          -1 35.621938  55.399273  23.254556

##     optimal.3 lambda.1 lambda.2    dol.delta    dol.omega   dol.lambda
##  1  49.616046 451.8013 598.9589   0.63409966   0.48287469 536.97618091
##  2   8.701041 636.3031 429.4659   0.08458527   0.99624955 674.30015903
##  3  64.267740 441.9105 615.0822   0.78101041   0.49048222 538.73234859
##  4  22.281453 615.0822 441.9105   0.21898959   0.99048222 662.20606601
##  5  82.227974 422.9191 648.7404   0.95800430   0.49825407 540.49590064
##  6  42.272337 593.2415 455.2425   0.42035428   0.97962398 650.57382749


# create a 0-1 spectrum with 2 transitions
rectspec = rectangularMaterial( lambda=c(579.8697,613.7544), alpha=1, wave=wave )

# compute the corresponding color XYZ
XYZ = product( rectspec, D50.eye )
XYZ
##                             X        Y          Z
##  BP_[579.87,613.754] 33.42026 21.96895 0.02979764

# trace a ray from middle gray through XYZ
white.XYZ   = product( neutralMaterial(1,wave=wave), D50.eye )
direction   = XYZ - white.XYZ/2

res = probeOptimalColors( D50.eye, 0.5, direction, aux=FALSE )
res$s         
##  1.00004   the ray has gone past the original color to the boundary

res$optimal
##              X        Y          Z
##  [1,] 33.41958 21.96774 0.02808178

res$lambda    
##  NA NA     because there are more than 2 transitions in the true optimal

# since s=1.00004 > 1,
# XYZ is actually in the interior of the color solid, and not on the boundary.
# The boundary is a little-bit further along the ray,
# and the corresponding spectrum has more than 2 transitions.

Run the code above in your browser using DataLab