Appendix A - Munsell Renotation Details

For history prior to 1943, please consult [12], [13], and [14].

The OSA Subcommitte on the Spacing of the Munsell Colors issued their Final Report in 1943 [11]. There were 2734 published coordinates (what might be called aim points), for Munsell Values from 1 to 9. The work was based on visual interpolation and smoothing on xy graphs of spectrophotometrically measured physical samples from the Munsell Book of Color (1929). There were 421 samples from [7] and 561 samples from [5]. There were a few obvious transcription errors in this data, which will be discussed below. After 1943, the Munsell Book of Color used these aim points, which are called the renotations.

In 1956 Judd and Wyszecki extended the renotation to Values 0.2/, 0.4/, 0.6/, and 0.8/, in [6]. The Chromas are /1, /2, /3, /4, /6, and /8. However, Chroma /1 was obtained by averaging /2 with white, and /3 was obtained by averaging /2 with /4. If Chromas /1 and /3 are dropped there are 355 very dark samples that remain. There were now 2734 + 355 = 3089 renotations, with even Chroma, in 1956.

In 1958, at the Forty-Third Annual Meeting of the Optical Society of America, a 15-minute talk [19] was given on extrapolation of these renotations for a computer program in development. Since the abstract of [19] is so short, we think it is appropriate to quote it in full:

TA15. Extension of the Munsell Renotation System. J . C. SCHLETER, DEANE B. JUDD, AND H. J. KEEGAN, National Bureau of Standards, Washington, D. C. - In order to convert tristimulus data to Munsell renotations by using the tri-dimensional, linear, interpolation program developed at the NBS for the high-speed, digital IBM 704 computer, it was necessary to extend the published data of the Munsell re-notation system^1,2 by extrapolation so that, for each point in Munsell color space, the enclosing six-sided volume would be defined by its eight corners. The approximately 2000 corners to be defined by CIE chromaticity coordinates, x, y, are at a given value level the intersections of the hue and chroma loci that cross the MacAdam limit³ at that value level and at the next lower level. Tentative definitions were found by graphical extrapolation on plots at each of the 14 value levels similar to those published. A second set of tentative definitions were found by plotting x and y against Munsell value for each of the 40 Munsell hues and extrapolating to value 10/. The final definitions were found by adjustment of the two tentative sets to produce smooth loci in both types of plots. Plots at selected value levels will be shown illustrating the extended data. (15 min.)
¹ Newhall, Nickerson, and Judd, J. Opt. Soc. Am. 33, 385 (1943).
² D. B. Judd and G. Wyszecki, J. Opt. Soc. Am. 46, 281 (1956).
³ D. MacAdam, J. Opt. Soc. Am. 25, 361 (1935).

The “2000 corners” is the number of points in the extension; we will call these the Schleter points.

The next talk [16] was a preliminary report on the above-mentioned software program. The authors note that 1) the extended grid from [19] “consists of approximately 5000 points”, 2) the program stores these points on magnetic tape, and 3) an xyY \(\to\) HVC conversion takes “approximately 20 seconds per sample”.

In 1960 there was an updated report [17] on the program. Rheinboldt and Menard report that the 1956 grid of 3089 points was extended by 1907 points to yield 4996 points. So the exact number of Schleter points is 1907. We also learn that the memory of the NBS IBM 704 “had been increased to 8000 words” for “the second version of the code”. This made it possible to store the entire grid in memory. They do not report the improved conversion time per sample.

Skip now to contemporary times. The Program of Color Science at the Rochester Institute Technology has posted 2 files: real.dat and all.dat, see [18]. The first file real.dat is described as “those colors listed the original 1943 renotation article” [11]. There are 2734 samples in the file, so this is an exact match ! Paul Centore found that before June 2010 there were 5 samples missing, but this has now been corrected (in 2018). Note that real.dat does not include the very dark colors, even though they are inside the MacAdam limits.

These are all the Munsell data, including the extrapolated colors. Note that extrapolated colors are in some cases unreal. That is, some lie outsize the Macadam [sic] limits.

This file should be used for those performing multidimensional interpolation to/from Munsell data. You will need the unreal colors in order to completely encompass the real colors, which is required to do the interpolation when near the Macadam limits.

There are 4995 samples in all.dat, which differs from the count in [17] by only 1. Perhaps the latter count includes Illuminant C ? These two counts are so close that we are confident that the extended grid in the NBS software from 1958 survived, and made its way to the Program of Color Science at RIT, and then into this package. The above-quoted abstract is likely all we will ever know about the methods used to extrapolate the values in all.dat. Note that all.dat does include the very dark colors.

Now return to the NBS computer program in [17]. For the forward conversion HVC \(\to\) xy they first determine the two adjacent planes of constant V that the given point is between. In both of these two HC planes they use bilinear interpolation to determine xy. They then linearly interpolate these 2 xy vectors. So a single conversion requires a table lookup of 8 points in general. With this method, the plotted ovoids of constant chroma are 40-sided polygons, and not the smooth ovoids in the published plots by Dorothy Nickerson. For an example, which looks more like a spiderweb, see the Loci of Contant Hue and Chroma vignette.

To get smooth ovoids, the default method in munsellinterpol uses bicubic interpolation in the HC plane. The forward conversion is then class \(C^1\) in H and C (except at C=0). To get \(C^1\) continuity in V as well, the default conversion uses cubic interpolation for 4 planes of constant value - 2 above the input value and 2 below. This requires a table lookup of 64 points in general, and it requires extending the renotations even more ! In the terminology of mathematical morphology the required binary image is the dilation of the 4995 points with a 3x3x3 structuring element. The grid in munsellinterpol currently has 7606 points.

The extrapolation was done, again from the grid of 4995 points, using the function gam() in the package mgcv which is preinstalled with R itself. It is a thin-plate spline model. A separate model was computed for each of the Munsell values: 0.2,0.4,0.6,0,8,1,…10. Unfortunately, for V=1.5 (on almost any V near 1.5) the forward conversion HVC \(\to\) xy was badly non-injective for small Chroma in the area of the hue GY. It is because the extrapolated xy values at V=0.8 seem to diverge a lot from V=1 and V=2 (they went too far into halfplane x+y>1). Recall that a conversion at V=1.5 requires lookup at 4 planes V=0.8,1,2, and 3. Also note that the weight at V=0.8 is negative. Cubic interpolation has overshoot and undershoot - and linear interpolation (where weights are non-negative) does not.

So for dark colors V<1 we decided to ignore the extrapolations in [19], and re-extrapolate, using a mgcv::gam() model, from all real points with V \(\le\) 1. This extrapolation was satisfactory.

Appendix B - Munsell Renotation Discrepancies

There have been at least 4 versions of the renotation tables. In chronological order they are:

• original tables in [11]
• tables in the “Color Bible” [28]
• tables in ASTM standard [3]
• digital file all.dat at [18] mentioned in the previous Appendix

I have noticed 4 discrepancies, and there might be more. They appear to be either simple transcription errors, or adjustments to make the data smoother. This package uses all.dat, except for 2.5PB 10/2 noted below.

For samples 7.5BG 4/16 and 5BG 1/4 it looks like a simple transcription error in [28], which was corrected in [3] and [18]. For 10PB 8/4 it looks like [3] made a smoothing adjustment, which was carried over to [18]. For 2.5R 9/2 it looks like [18] made a smoothing adjustment.

Finally, in the extrapolated data in [18], we have made a small smoothing adjustment at 2.5PB 10/2 for this package. Before the adjustment, the forward mapping HVC \(\to\) xyY was not injective near N 10/, which caused the root-finding iteration to cycle. This was discovered during testing phase; more ambitious “vetting” software, which specifically searches for non-injectivity, might uncover more of these.

Appendix C - Inversion Details

Given an xy target point, the method in [17] searches for the exact quadrilateral that contains xy. It then solves a pair of quadratic equations to determine H and C. For bicubic interpolation, a similar explicit solution is not really feasible. The method in [2] starts with an initial approximation HC\(_0\), and modifies H and C alternately to converge to an HC that maps to xy, within a certain tolerance.

We decided to use a general root-solver from the package rootSolve that uses Newton-Raphson iteration. In our case, this iteration really does double-duty: the first iteration or two locate the quadrilateral that contains xy, and the remaining iterations polish the root. It typically takes about 4 or 5 total iterations. The theory of Newton’s Method only covers functions of class \(C^2\) (see [24]) and our function is at best \(C^1\). However, it is \(C^\infty\) on the quadrilateral interiors, so if the root is in the interior of a quadrilateral, and the iteration gets close enought to the root, then it typically stays inside that quadrilateral and converges. Even when a root is on the boundary of quadrilateral, we have found that the basin of attraction is still large enough to find it. The usual trouble is when the intial estimate is too far from the root, and then the iterations wanders outside the lookup table itself. Cycling is another possible problem - the iteration could jump around from quadrilateral to quadrilateral. But we have only seen this happen once, when the function was non-injective near white (mentioned at the end of the previous section). The quadrilaterals become small triangular slivers there, which causes more trouble near all neutrals.

Instead of HC itself, we decided to use a rectangular system: \[A = C \cos \left( \frac{\pi}{50} H \right) ~~~~~~~~ B = C \sin \left( \frac{\pi}{50} H \right)\] with the inverse: \[H=\left( \frac{50}{\pi} \right) \text{atan2}(B,A) ~~~~~~~~ C = \text{sqrt}( A^2 + B^2 )\]

\(A\) and \(B\) are analogous to \(a\) and \(b\) in \(Lab\) space. The symbols \(a''\) and \(b''\) are used for the same coordinate system in [4] page 207. So now we are trying to invert \(AB \to HC \to xy\). Each iteration requires an extra conversion \(AB \to HC\), but this is very fast.

For the initial estimate, the algorithm in [2] recommends: \[A_0 = a/5.5 ~~~~~~~~~ B_0 = b/5.5\] We use a cubic polynomial in \(a\) and \(b\) for better accuracy. A separate polynomial is computed for each V-plane, using R’s built-in function for linear models with no intercept:

stats::lm( A ~ polym(a,b,degree=3,raw=TRUE) + 0, ... )

and similarly for \(B\). We precompute the model, but only save the coefficients. For V < 1, a cubic polynomial in \(a\) and \(b\) was still not good enough, and it was found that a cubic in \(\Delta x\) and \(\Delta y\) worked better; where \(\Delta x = x - x_C\) and \(\Delta y = y - y_C\), and \(x_C,y_C\) are the chromaticity of Illuminant C. Evaluation of the polynomial takes much longer than division by 5.5, but it is still negligible compared to the iterations to come.

The root-solver is called like this:

rootSolve::multiroot( forwardfun, AB0, rtol=1.e-8, atol=1.e-6, ctol=0, xy_target=xy_target, ...)

The function forwardfun() performs the equivalent of MunsellToxyY() but does not call it explicitly. The argument AB0 is the initial estimate. It usually takes 4 or 5 iterations to converge with the indicated tolerances, which are in the xy plane.

Appendix D - Glossary

luminance factor

The luminance factor of an object is the quotient of the luminance of the object (for a specific illuminant and in a specific viewing direction) divided by that of the perfect reflecting diffuser. In Munsell work this quotient is denoted by Y and is a percentage in the interval [0,100]. For Lambertian surfaces, Y does not depend on the viewing direction.

uniform lightness scale

A uniform lightness scale is a reparameterization of luminance factor, so that equal numerical steps in lightness produce equal perceptual steps. The Munsell Value is an example of such a scale; in early investigations it was found that the luminance factor of the background affects the perception. Another popular scale is CIE Lightness (1976). For more on these, see the Uniform Lightness Scales vignette.

Session Information

R version 4.1.3 (2022-03-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=C                           LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] flextable_0.6.10      spacesXYZ_1.2-1       spacesRGB_1.5-0       munsellinterpol_3.0-0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.7           knitr_1.37           xml2_1.3.3           magrittr_2.0.1      
 [5] uuid_1.0-3           R6_2.5.1             rlang_0.4.12         fastmap_1.1.0       
 [9] stringr_1.4.0        highr_0.9            tools_4.1.3          data.table_1.14.2   
[13] rootSolve_1.8.2.3    xfun_0.28            htmltools_0.5.2      systemfonts_1.0.3   
[17] yaml_2.2.1           digest_0.6.28        zip_2.2.0            officer_0.4.1       
[21] microbenchmark_1.4.9 base64enc_0.1-3      evaluate_0.14        rmarkdown_2.11      
[25] stringi_1.7.5        compiler_4.1.3       gdtools_0.2.3