1 Background

The Munsell system provides a quantitative (numeric) description of color that is approximately linear in terms of average human perception, and, conveniently isolates the independent qualities of color into hue, value, and chroma. Moving from one chip to the next represents a linear traversal through a continuous, 3D (hue, value, chroma) color space. The presentation (soil color book) and notation (e.g. 10YR 3/3) that we are familiar with is based on a discretization of that space into perceptually-compact units (the color chips) that offers sufficient contrast for field identification without an overwhelming number of possibilities. For example, our ability to perceive contrast in lightness (Munsell value) is roughly 2-4x greater than our ability to perceive contrast in chroma: hence the lack of chips representing /5, /7, /9, etc..

Given that the Munsell system is quantitative and roughly linear in terms of human color perception, it is entirely possible to:

  • treat value and chroma as numeric entities
  • interpolate between chips (hue, value, or chroma)
  • simulate mixtures of two or more colors in Munsell notation

I’d recommend SSSA special publication #31, any book on color science, or one of several excellent websites. Our definition of “soil color contrast” is based on the quantitative, linear nature of the Munsell system.

An accurate simulation of pigment mixtures (“subtractive” color mixtures) is incredibly complex due to factors that aren’t easily measured or controlled: pigment solubility, pigment particle size distribution, water content, substrate composition, and physical obstruction to name a few. That said, it is possible to simulate reasonable, subtractive color mixtures given a reference spectra library (350-800nm) and some assumptions about pigment qualities and lighting. For the purposes of estimating a mixture of soil colors (these are pigments after all) we can relax these assumptions and assume a standard light source. The only missing piece is the spectral library for all Munsell chips in our color books.

Thankfully, Scott Burns has outlined the entire process, and Paul Centore has provided a Munsell color chip reflectance spectra library. The estimation of a subtractive mixture of soil colors can proceed as follows:

  1. look-up the associated spectra for each color
  2. compute the weighted (area proportion) geometric mean of spectra
  3. either:
  • search for the closest matching spectra in the reference library, or
  • compute a new reluctance spectra by integration over standard observer and illuminant
  1. convert the final spectra (reference, or computed) to closest Munsell chip

Key assumptions include:

  • similar particle size distribution
  • similar mineralogy (i.e. pigmentation qualities)
  • similar water content
  • similar scattering coefficients

For the purposes of estimating (for example) a “mixed soil color within the top 18cm of soil” these assumptions are usually valid. Again, these are estimates that are ultimately “snapped” to the nearest chip and not do not need to approach the accuracy of paint-matching systems.

So what is a reasonable answer to the following soil color mixture question: 50% 5YR 2/2 + 50% 10YR 4/4 = ?.

1.1 Munsell Spectral Library

library(aqp)
library(lattice)
library(tactile)
library(sharpshootR)

# need this for mixingMethod = 'reference'
library(gower)

# need this for colorMixtureVenn()
library(venn)

# local copy of the Munsell chip spectral library
# c/o http://www.munsellcolourscienceforpainters.com/
# odd chroma spectra via interpolation
# see ?munsell.spectra for details
# try aqp:::.summarizeMunsellSpectraRanges()
data(munsell.spectra)

# all hues, limit to specific hue / chroma slice
x <- munsell.spectra[munsell.spectra$value  == 6 & munsell.spectra$chroma == 8, ]

# each Munsell chip has a 36-element spectra
# ranging from 380-730 nm
# table(x$munsell)

# spectra IDs
x$ID <- factor(x$munsell)
# create a color / chip
cols <- parseMunsell(as.character(levels(x$ID)))

# plot style
tps <- tactile.theme(superpose.line = list(col = cols, lwd = 2))

# final figure
xyplot(
  reflectance ~ wavelength, groups = ID, data = x, 
  par.settings = tps,
  main = 'Value 6 / Chroma 8',
  type = c('l', 'g'),
  ylab = 'Reflectance',
  xlab = 'Wavelength (nm)',
  scales = list(tick.number = 12),
  xlim = c(370, 740)
)

Value 6 / Chroma 8 Wavelength (nm) Reflectance 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 400 450 500 550 600 650 700

1.2 Implementation

Required packages, be sure to get the latest {aqp} and {sharpshootR} from GitHub.

library(aqp)
library(lattice)
library(latticeExtra)
library(tactile)
library(sharpshootR)
library(gower)
library(venn)
library(cluster)

Create a vector of colors to mix, and evaluate color contrast metrics.

colors <- c('5YR 2/2', '10YR 4/4')
colorContrast(m1 = colors[1], m2 = colors[2])
##        m1       m2 dH dV dC    dE00        cc
## 1 5YR 2/2 10YR 4/4  2  2  2 19.8106 Prominent

We can see that these colors are 2 hue, 2 value, and 2 chroma “chips” apart from each other. This evaluates to a prominent soil color contrast class and \(\Delta{E_{00}}\) of about 20.

A picture is far more intuitive.

par(bg = 'black', fg = 'white', mar = c(0, 0, 0, 0))
colorContrastPlot(m1 = colors[1], m2 = colors[2])

m1 m2 5YR 2/2 10YR 4/4 Prominent Δ E 00 19.8

Simulate a mixture via mixMunsell. See the online documentation for details on the method and relevant citations.

# just the closest match
mixed <- mixMunsell(colors, mixingMethod = 'reference')

# top 5 matches, decreasing order
mixed.5 <- mixMunsell(colors, mixingMethod = 'reference', n = 5)

Left to right: original colors-to-mix and final mixture candidate (10YR 3/3). This is the Munsell chip associated with the reference spectra that is closest to the weighted geometric mean of the source spectra.

par(bg = 'black', fg = 'white', mar = c(0, 0, 0, 0))
soilPalette(parseMunsell(c(colors, mixed$munsell)), lab = c(colors, mixed$munsell))

5YR 2/2 10YR 4/4 10YR 3/3

Evaluate the top 5 mixture candidates (ranked by spectral distance) using metrics of color contrast. Note that color contrast does not always correlate with spectral distance.

par(bg = 'black', fg = 'white')
colorContrastPlot(
  m1 = rep(mixed$munsell, times = nrow(mixed.5)),
  m2 = mixed.5$munsell,
  labels = c('closest match', 'top 5 matches')
  )

closest match top 5 matches 10YR 3/3 10YR 3/3 10YR 3/3 10YR 3/3 10YR 3/3 10YR 3/3 7.5YR 3/3 2.5YR 3/2 2.5Y 3/3 7.5YR 3/2 Faint Faint Prominent Faint Faint Δ E 00 0 Δ E 00 3.31 Δ E 00 8.29 Δ E 00 3.42 Δ E 00 4.52

1.3 Results

Demonstrate reference spectra for the source colors and final mixture. The dashed line is the weighted geometric mean spectra.

plotColorMixture(colors, mixingMethod = 'reference', showMixedSpec = TRUE)

Wavelength (nm) Reflectance 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 400 450 500 550 600 650 700 750 5YR 2/2 50% 10YR 4/4 50% 10YR 3/3 #1 color 1 color 2 mix #1

50% 5YR 2/2 + 50% 10YR 4/4 = 10YR 3/3.

What do you think, is the mixture plausible? Experiments with real soil material + Nix Pro sensor measurements, visible-spectrum spectrophotometer pending.

1.3.1 Mollic Epipedon Mixing Rule

Consider the physical mixture of soil material 0 to 18cm, with two horizons spanning 0-12cm, and 12-18cm. Assume that the soil material is of roughly the same soil texture and has similar reflective properties.

60% 10YR 2/2 + 40% 10Y 6/6 = 5Y 3/3

# horizon thickness (cm)
w <- c(12, 8)

# moist soil colors
m <- c('10YR 2/2', '10Y 6/6')

plotColorMixture(m, w, mixingMethod = 'exact', label.cex = 0.65, overlapFix = TRUE)

Wavelength (nm) Reflectance 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 400 450 500 550 600 650 700 750 10YR 2/2 60% 10Y 6/6 40% 5Y 3/3 #1 color 1 color 2 mix #1

1.3.2 Experimenting with Color Mixtures

The colorMixtureVenn function from the sharpshootR package is a handy way to view simulated mixtures of 2-7 colors.

colorMixtureVenn(colors)

5YR 2/2 10YR 4/4 7.5YR 3/3

chips <- c('10YR 8/1', '2.5YR 3/6', '10YR 2/2')
colorMixtureVenn(chips)

10YR 8/1 2.5YR 3/6 10YR 2/2 2.5YR 5/5 10YR 4/2 5YR 2.5/4 5YR 4/3

chips <- c('2.5YR 8/2', '2.5B 3/6', '2.5G 4/4', '2.5R 4/8', '2.5Y 6/10')
colorMixtureVenn(chips, ellipse = TRUE)

2.5YR 8/2 2.5B 3/6 2.5G 4/4 2.5R 4/8 2.5Y 6/10 10BG 4/1 7.5GY 5/3 2.5Y 5/2 10YR 5/4 2.5GY 4/2 5BG 5/4 7.5PB 4/1 7.5GY 5/3 2.5YR 5/2 10Y 6/4 5YR 6/6 10BG 3/2 10GY 4/4 5Y 4/2 10YR 4/4 2.5B 5/4 10GY 6/2 2.5R 6/6 2.5Y 7/7 5BG 3/11 5PB 3/2 7.5GY 4/4 2.5YR 4/2 10Y 5/5 5YR 5/7 10Y 4/2

2 Additional Examples

2.1 Tinkering with plotColorMixture

plotColorMixture(c('10YR 5/3', '10YR 3/2'))

Wavelength (nm) Reflectance 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 400 450 500 550 600 650 700 750 10YR 5/3 50% 10YR 3/2 50% 10YR 4/2 #1 color 1 color 2 mix #1

Compare “reference” vs. “exact” mixing methods.

plotColorMixture(c('10YR 6/2', '5YR 5/6'), mixingMethod = 'reference', showMixedSpec = TRUE)

Wavelength (nm) Reflectance 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 400 450 500 550 600 650 700 750 10YR 6/2 50% 5YR 5/6 50% 7.5R 5/5 #1 color 1 color 2 mix #1

plotColorMixture(c('10YR 6/2', '5YR 5/6'), mixingMethod = 'exact')

Wavelength (nm) Reflectance 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 400 450 500 550 600 650 700 750 10YR 6/2 50% 5YR 5/6 50% 5YR 6/4 #1 color 1 color 2 mix #1

# is this right? can the resulting spectra be "higher" than the source?
plotColorMixture(c('10YR 6/2', '5YR 5/6'), w = c(2,1))

Wavelength (nm) Reflectance 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 400 450 500 550 600 650 700 750 10YR 6/2 67% 5YR 5/6 33% 7.5YR 6/4 #1 color 1 color 2 mix #1

# label collision
plotColorMixture(c('10YR 5/3', '10YR 3/2', '5R 2/2'))

Wavelength (nm) Reflectance 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 400 450 500 550 600 650 700 750 10YR 5/3 33% 10YR 3/2 33% 5R 2/2 33% 5YR 3/2 #1 color 1 color 2 color 3 mix #1

plotColorMixture(c('10YR 5/3', '10YR 3/2', '5R 2/2'), label.cex = 0.65)

Wavelength (nm) Reflectance 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 400 450 500 550 600 650 700 750 10YR 5/3 33% 10YR 3/2 33% 5R 2/2 33% 5YR 3/2 #1 color 1 color 2 color 3 mix #1

plotColorMixture(c('10YR 4/6', '2.5Y 6/2', '5Y 2/2'), w = c(1, 1, 2))

Wavelength (nm) Reflectance 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 400 450 500 550 600 650 700 750 10YR 4/6 25% 2.5Y 6/2 25% 5Y 2/2 50% 2.5Y 3/3 #1 color 1 color 2 color 3 mix #1

plotColorMixture(c('5G 6/5', '5R 5/4'), w = c(1, 2))

Wavelength (nm) Reflectance 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 400 450 500 550 600 650 700 750 5G 6/5 33% 5R 5/4 67% 7.5YR 5/2 #1 color 1 color 2 mix #1

plotColorMixture(c('5G 6/5', '5R 5/4'), w = c(3, 1), label.cex = 0.65)

Wavelength (nm) Reflectance 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 400 450 500 550 600 650 700 750 5G 6/5 75% 5R 5/4 25% 2.5G 6/3 #1 color 1 color 2 mix #1

Check: “blue” + “yellow” = “green”, yes. Note weighted geometric mean spectra (dotted line): closest chip isn’t a perfect fit but close enough.

plotColorMixture(c('5B 5/10', '5Y 8/8'), label.cex = 0.65, showMixedSpec = TRUE, mixingMethod = 'reference')

Wavelength (nm) Reflectance 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 400 450 500 550 600 650 700 750 5B 5/10 50% 5Y 8/8 50% 7.5G 6/6 #1 color 1 color 2 mix #1

plotColorMixture(c('5B 5/10', '5Y 8/8'), label.cex = 0.65, mixingMethod = 'exact')

Wavelength (nm) Reflectance 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 400 450 500 550 600 650 700 750 5B 5/10 50% 5Y 8/8 50% 7.5G 6/5 #1 color 1 color 2 mix #1

Top 3 mixture candidates

plotColorMixture(c('5B 5/10', '5Y 8/8'), label.cex = 0.65, showMixedSpec = TRUE, mixingMethod = 'reference', n = 3)

Wavelength (nm) Reflectance 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 400 450 500 550 600 650 700 750 5B 5/10 50% 5Y 8/8 50% 7.5G 6/6 #1 7.5G 6/5 #2 7.5G 6/7 #3 color 1 color 2 mix #1 mix #2 mix #3

Interesting mixture of three colors, closest Munsell chip of mixture is not in spectral library.

plotColorMixture(c('5Y 8/10', '5B 4/10', '5R 4/8'), label.cex = 0.65, showMixedSpec = TRUE, mixingMethod = 'exact')

Wavelength (nm) Reflectance 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 400 450 500 550 600 650 700 750 5Y 8/10 33% 5B 4/10 33% 5R 4/8 33% 2.5GY 4/2 #1 color 1 color 2 color 3 mix #1

2.2 Simulation of continuous mixtures

# m: two colors
# s: sequence of weights for the first color
.continuousMixture <- function(m, s) {
  z <- lapply(s, function(i) {
  m <- mixMunsell(m, w = c(1-i, i), mixingMethod = 'exact')
  return(m$munsell)
})

z <- do.call('c', z)
return(z)
}

Blue + yellow.

s <- seq(0, 1, by = 0.25)
z <- .continuousMixture(c('5B 5/10', '5Y 8/8'), s)

par(bg = 'black', fg = 'white', mar = c(0,0,0,0))
soilPalette(parseMunsell(z), sprintf('%s (%s / %s))', z, 1 - s, s))

5B 5/10 (1 / 0)) 7.5BG 5/16 (0.75 / 0.25)) 7.5G 6/5 (0.5 / 0.5)) 5GY 7/5 (0.25 / 0.75)) 5Y 8/8 (0 / 1))

Near-white + dark brown.

s <- seq(0, 1, length.out = 6)
z <- .continuousMixture(c('10YR 9/1', '10YR 2/1'), s)

par(bg = 'black', fg = 'white', mar = c(0,0,0,0))
soilPalette(parseMunsell(z), sprintf('%s (%s / %s))', z, 1 - s, s))

7.5YR 9/1 (1 / 0)) 10YR 7/1 (0.8 / 0.2)) 10YR 5/1 (0.6 / 0.4)) 10YR 4/1 (0.4 / 0.6)) 10YR 3/1 (0.2 / 0.8)) 10YR 2/1 (0 / 1))

Cobalt blue + grey.

s <- seq(0, 1, length.out = 6)
z <- .continuousMixture(c('7.5B 6/10', '7.5B 8/2'), s)

par(bg = 'black', fg = 'white', mar = c(0,0,0,0))
soilPalette(parseMunsell(z), sprintf('%s (%s / %s))', z, 1 - s, s))

7.5B 6/10 (1 / 0)) 2.5B 6/15 (0.8 / 0.2)) 7.5B 7/8 (0.6 / 0.4)) 7.5B 7/6 (0.4 / 0.6)) 7.5B 8/4 (0.2 / 0.8)) 5B 8.5/2 (0 / 1))

2.3 Interpreting Spectral (Gower) Distances

Should the mixture of two 10YR colors always result in a 10YR mixture candidate? Depends on the mixing method

m <- c('10YR 8/6', '10YR 2/2')

knitr::kable(
  rbind(
    mixMunsell(m, mixingMethod = 'exact'),
    mixMunsell(m, mixingMethod = 'reference'),
    mixMunsell(m, mixingMethod = 'estimate')
  ),
  digits = 2
)
munsell distance scaledDistance distanceMetric mixingMethod
10YR 4/3 3.50 NA dE00 exact
7.5YR 4/3 0.01 0.28 Gower reference
10YR 5/4 1.21 NA dE00 estimate

The weighted geometric mean of the two spectra seems plausible:

plotColorMixture(c('10YR 8/6', '10YR 2/2'))

Wavelength (nm) Reflectance 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 400 450 500 550 600 650 700 750 10YR 8/6 50% 10YR 2/2 50% 10YR 4/3 #1 color 1 color 2 mix #1

Evaluate the top 20 mixture candidates in terms of spectral distance and \(\Delta{E_{00}}\) relative to the most likely match.

m <- c('10YR 8/6', '10YR 2/2')
x <- mixMunsell(m, n = 20, mixingMethod = 'reference')

z <- expand.grid(m1 = x$munsell[1], m2 = x$munsell)
d <- colorContrast(z[[1]], z[[2]])

n <- nrow(x)
s <- 1:n
cols <- parseMunsell(x$munsell)

Prepare some figures to assist with the exploration.

par(mar = c(4.5, 4.5, 2, 1), bg = 'black', fg = 'white')

plot(s, x$distance, type = 'n', ylab = 'Spectral Distance', xlab = 'Match Candidate', col.axis = 'white', col.lab = 'white', axes = FALSE)
grid()
lines(s, x$distance, type = 's')
points(s, x$distance, pch = 15, col = cols, cex= 4)
text(s, x$distance, gsub(' ', '\n', x$munsell), cex = 0.5, col = invertLabelColor(cols))
axis(side = 1, at = s, col.axis = 'white')
axis(side = 2, col.axis = 'white', las = 1, cex.axis = 0.75)
title(sprintf("Mixture: %s + %s", m[1], m[2]), col.main = 'white')

Match Candidate Spectral Distance 7.5YR 4/3 5YR 4/3 5YR 4/4 2.5YR 4/4 7.5YR 4/4 10R 4/3 2.5YR 4/3 2.5Y 4/2 7.5Y 4/2 10YR 4/3 10YR 4/4 5Y 4/2 2.5Y 4/3 7.5YR 4/2 10R 4/4 10YR 4/2 5YR 4/2 2.5YR 4/2 10YR 4/5 7.5R 4/3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 0.007 0.008 0.009 0.010 0.011 0.012 0.013 Mixture: 10YR 8/6 + 10YR 2/2

par(mar = c(4.5, 4.5, 2, 1), bg = 'black', fg = 'white')

plot(s, d$dE00, type = 'n', ylab = expression(~Delta*E['00']), xlab = 'Match Candidate', col.axis = 'white', col.lab = 'white', axes = FALSE)
grid()
lines(s, d$dE00, type = 's')
points(s, d$dE00, pch = 15, col = cols, cex= 4)
text(s, d$dE00, gsub(' ', '\n', x$munsell), cex = 0.5, col = invertLabelColor(cols))
axis(side = 1, at = s, col.axis = 'white')
axis(side = 2, at = 0:15, col.axis = 'white', las = 1, cex.axis = 0.75)
title(sprintf("Mixture: %s + %s", m[1], m[2]), col.main = 'white')

Match Candidate Δ E 00 7.5YR 4/3 5YR 4/3 5YR 4/4 2.5YR 4/4 7.5YR 4/4 10R 4/3 2.5YR 4/3 2.5Y 4/2 7.5Y 4/2 10YR 4/3 10YR 4/4 5Y 4/2 2.5Y 4/3 7.5YR 4/2 10R 4/4 10YR 4/2 5YR 4/2 2.5YR 4/2 10YR 4/5 7.5R 4/3 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 0 1 2 3 4 5 6 7 8 9 10 11 12 Mixture: 10YR 8/6 + 10YR 2/2

The relationship between spectral distance and \(\Delta{E_{00}}\) is complex.

par(mar = c(4.5, 4.5, 2, 1), bg = 'black', fg = 'white')

plot(x$distance, d$dE00, type = 'n', xlab = 'Spectral Distance', ylab = expression(~Delta*E['00']), col.axis = 'white', col.lab = 'white', axes = FALSE)
grid()
points(x$distance, d$dE00, pch = 15, col = cols, cex= 4.5)
text(x$distance, d$dE00, gsub(' ', '\n', x$munsell), cex = 0.66, col = invertLabelColor(cols))
axis(side = 1, col.axis = 'white')
axis(side = 2, at = 0:15, col.axis = 'white', las = 1, cex.axis = 0.75)
title(sprintf("Mixture: %s + %s", m[1], m[2]), col.main = 'white')

Spectral Distance Δ E 00 7.5YR 4/3 5YR 4/3 5YR 4/4 2.5YR 4/4 7.5YR 4/4 10R 4/3 2.5YR 4/3 2.5Y 4/2 7.5Y 4/2 10YR 4/3 10YR 4/4 5Y 4/2 2.5Y 4/3 7.5YR 4/2 10R 4/4 10YR 4/2 5YR 4/2 2.5YR 4/2 10YR 4/5 7.5R 4/3 0.007 0.008 0.009 0.010 0.011 0.012 0.013 0 1 2 3 4 5 6 7 8 9 10 11 12 Mixture: 10YR 8/6 + 10YR 2/2

Try again, this time using the entire spectral library.

m <- c('10YR 8/6', '10YR 2/2')

# return distances to entire spectral library
x <- mixMunsell(m, n = 2486, mixingMethod = 'reference')

# do not compare 1st mixture candidate to itself (dE00 = 0)
z <- expand.grid(m1 = x$munsell[1], m2 = x$munsell[-1])
d <- colorContrast(z[[1]], z[[2]])

# color vector, note that we are skipping the first color
cols <- parseMunsell(x$munsell[-1])

par(mar = c(4.5, 4.5, 2, 1), bg = 'black', fg = 'white')

# skip first mixture candidate
plot(x$distance[-1], d$dE00, xlab = 'Spectral Distance', ylab = expression(~Delta*E['00']), col.axis = 'white', col.lab = 'white', las = 1, type = 'n', ylim = c(0, 50))

grid()

# skip first mixture candidate
points(x$distance[-1], d$dE00, pch = 15, col = cols, cex = 0.5)
title(sprintf("Mixture Candidates: %s + %s", m[1], m[2]), col.main = 'white')

# add first mixture candidate
points(x$distance[1], d$dE00[1], bg = cols[1], pch = 22)

0.00 0.05 0.10 0.15 0.20 0 10 20 30 40 50 Spectral Distance Δ E 00 Mixture Candidates: 10YR 8/6 + 10YR 2/2

Reasonable thresholds on spectral distance between mixture and mixture candidates.

# assuming threshold of dE00
reasonable.matches <- d$m2[which(d$dE00 <= 5)]

# find these colors in the mixture candidates
idx <- which(x$munsell %in% reasonable.matches)

# combine mixture candidates + contrast metrics
z <- merge(x[idx, ], d, by.x = 'munsell', by.y = 'm2', sort = FALSE, incomparables = FALSE)

# display
colorContrastPlot(z$m1, z$munsell, labels = c('First Mixture Canidate', 'Subsequent Candidates'))

First Mixture Canidate Subsequent Candidates 7.5YR 4/3 7.5YR 4/3 7.5YR 4/3 7.5YR 4/3 7.5YR 4/3 7.5YR 4/3 5YR 4/3 5YR 4/4 7.5YR 4/4 10YR 4/3 7.5YR 4/2 5YR 4/2 Faint Faint Faint Faint Faint Faint Δ E 00 2.73 Δ E 00 3.95 Δ E 00 3.17 Δ E 00 3.48 Δ E 00 3.91 Δ E 00 4.55

2.3.1 Interpreting Spectral (Gower) Distance

Note that these distances are always within 0-1.

library(cluster)

# spectral library
data("munsell.spectra.wide")

# remove 'wavelength' column and transpose
z <- t(munsell.spectra.wide[, -1])

# all pair-wise distances, using Gower's distance metric
d <- daisy(z, metric = 'gower', stand = FALSE)

# quick plot
par(mfcol = c(1, 2), mar = c(4.5, 4, 3, 1))
hist(as.vector(d), main = 'All Spectra', axes = FALSE, xlab = 'Spectral Distance (Gower)', ylab = 'Frequency', breaks = 50)
box()
axis(1, at = seq(0, 1, by = 0.1), cex.axis = 0.75)
axis(2, las = 1, cex.axis = 0.75)

hist(as.vector(d)[as.vector(d) < 0.1], main = 'Spectra d < 0.1', axes = FALSE, xlab = 'Spectral Distance (Gower)', ylab = '', breaks = 50)
box()
axis(1, at = seq(0, 0.1, by = 0.01), cex.axis = 0.75)
axis(2, las = 1, cex.axis = 0.75)

All Spectra Spectral Distance (Gower) Frequency 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0 50000 100000 150000 200000 Spectra d < 0.1 Spectral Distance (Gower) 0.00 0.02 0.04 0.06 0.08 0.10 0 10000 20000 30000 40000

2.3.2 Spectral Differences: Adjacent Chips

Use spectra differences between adjacent chips to determine a reasonable threshold.

# Munsell LUT
data(munsell)

# start with a copy
x <- munsell
x$munsell <- sprintf('%s %s/%s', x$hue, x$value, x$chroma)

# only chips with available spectra
# 2442 chips
idx <- which(x$munsell %in% dimnames(z)[[1]])
x <- x[idx, ]

# create all pair-wise combinations
# ~ 2,980,461 combinations!
x.comb <- combn(x$munsell, 2)

# convert to 2-column data.frame
x.comb <- data.frame(t(x.comb), stringsAsFactors = FALSE)

head(x.comb)
##        X1        X2
## 1 10B 2/1   10B 2/2
## 2 10B 2/1   10B 2/3
## 3 10B 2/1   10B 2/4
## 4 10B 2/1   10B 2/5
## 5 10B 2/1   10B 2/6
## 6 10B 2/1 10B 2.5/1
## TODO: make this a new function
## TODO: the following takes a while

# borrowed from aqp::colorContrast
m1.pieces <- parseMunsell(x.comb[[1]], convertColors = FALSE)
m2.pieces <- parseMunsell(x.comb[[2]], convertColors = FALSE)

# convert to value and chroma to numeric
m1.pieces[[2]] <- as.numeric(m1.pieces[[2]])
m1.pieces[[3]] <- as.numeric(m1.pieces[[3]])
m2.pieces[[2]] <- as.numeric(m2.pieces[[2]])
m2.pieces[[3]] <- as.numeric(m2.pieces[[3]])

# difference in number of hue chips, clock-wise, as specified in:
# https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/ref/?cid=nrcs142p2_053569
dH <- abs(huePosition(m1.pieces[[1]]) - huePosition(m2.pieces[[1]]))
# difference in number of value chips
dV <- abs(m1.pieces[[2]] - m2.pieces[[2]])
# difference in number of chroma chips
dC <- abs(m1.pieces[[3]] - m2.pieces[[3]])

## ^^^^ need a new function

# encode combinations
x.comb$code <- sprintf("%s-%s-%s", dH, dV, dC)

# select adjacent chips (single hue, value, or chroma)

# first-pass, delta = 1
# ~ 17,060 combinations
idx <- which(dH %in% c(0,1) & dV %in% c(0,1) & dC %in% c(0,1))
x.comb <- x.comb[idx, ]

# second-pass, single chip changes only
# ~ 4,378 combinations
idx <- which(x.comb$code %in% c('1-0-0', '0-1-0', '0-0-1'))
x.comb <- x.comb[idx, ]
# expand full distance matrix to upper/diag/lower
m <- as.matrix(d)

# vector of pair-wise distances for adjacent chips only
d.adj <- sapply(1:nrow(x.comb), FUN = function(i) {
  # distance by pair
  m[x.comb$X1[i], x.comb$X2[i]]
})

par(mfcol = c(1, 2), mar = c(4.5, 4, 3, 1))
hist(as.vector(d), main = 'All Spectra', axes = FALSE, xlab = 'Spectral Distance (Gower)', ylab = 'Frequency', breaks = 50)
box()
axis(1, at = seq(0, 1, by = 0.1), cex.axis = 0.75)
axis(2, las = 1, cex.axis = 0.75)

hist(d.adj, main = 'Adjacent Munsell Chips', axes = FALSE, xlab = 'Spectral Distance (Gower)', ylab = 'Frequency', breaks = 50)
box()
axis(1, at = seq(0, 1, by = 0.1), cex.axis = 0.75)
axis(2, las = 1, cex.axis = 0.75)

All Spectra Spectral Distance (Gower) Frequency 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0 50000 100000 150000 200000 Adjacent Munsell Chips Spectral Distance (Gower) Frequency 0.0 0.1 0.2 0 500 1000 1500

knitr::kable(t(quantile(d.adj)), digits = 3)
0% 25% 50% 75% 100%
0 0.008 0.018 0.043 0.216
idx.hue <- which(x.comb$code == '1-0-0')
idx.value <- which(x.comb$code == '0-1-0')
idx.chroma <- which(x.comb$code == '0-0-1')

# evaluate by hue
d.adj.hue <- sapply(idx.hue, FUN = function(i) {
  m[x.comb$X1[i], x.comb$X2[i]]
})

# evaluate by value
d.adj.value <- sapply(idx.value, FUN = function(i) {
  m[x.comb$X1[i], x.comb$X2[i]]
})

# evaluate by chroma
d.adj.chroma <- sapply(idx.chroma, FUN = function(i) {
  m[x.comb$X1[i], x.comb$X2[i]]
})

par(mfcol = c(1, 3), mar = c(4.5, 4, 3, 1))

hist(d.adj.value, main = 'Adjacent Munsell Chips (Value)', axes = FALSE, xlab = 'Spectral Distance (Gower)', ylab = 'Frequency', breaks = 50, xlim = c(0, 0.3))
box()
axis(1, at = seq(0, 0.5, by = 0.05), cex.axis = 0.75)
axis(2, las = 1, cex.axis = 0.75)

hist(d.adj.hue, main = 'Adjacent Munsell Chips (Hue)', axes = FALSE, xlab = 'Spectral Distance (Gower)', ylab = 'Frequency', breaks = 50, xlim = c(0, 0.3))
box()
axis(1, at = seq(0, 0.5, by = 0.05), cex.axis = 0.75)
axis(2, las = 1, cex.axis = 0.75)

hist(d.adj.chroma, main = 'Adjacent Munsell Chips (Chroma)', axes = FALSE, xlab = 'Spectral Distance (Gower)', ylab = 'Frequency', breaks = 50, xlim = c(0, 0.3))
box()
axis(1, at = seq(0, 0.5, by = 0.05), cex.axis = 0.75)
axis(2, las = 1, cex.axis = 0.75)

Adjacent Munsell Chips (Value) Spectral Distance (Gower) Frequency 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0 50 100 150 Adjacent Munsell Chips (Hue) Spectral Distance (Gower) Frequency 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0 100 200 300 400 500 600 Adjacent Munsell Chips (Chroma) Spectral Distance (Gower) Frequency 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0 200 400 600 800

##
adj.spec.dist <- make.groups(
  all = data.frame(x = d.adj),
  hue = data.frame(x = d.adj.hue),
  value = data.frame(x = d.adj.value),
  chroma = data.frame(x = d.adj.chroma)
)

bwplot(which ~ x, data = adj.spec.dist, par.settings = tactile.theme())

x all hue value chroma 0.00 0.05 0.10 0.15 0.20

2.3.3 Comparison with dE000

x.comb$code <- factor(x.comb$code)
cc <- colorContrast(x.comb$X1, x.comb$X2)


xyplot(
  cc$dE00 ~ d.adj,
  groups = x.comb$code,
  main = 'Adjacent Munsell Chips',
  scales = list(
    x = list(log = 10)
  ),
  par.settings = tactile.theme(
    background = list(col = 'black'),
    axis.line = list(col = 'white'),
    axis.text = list(col = 'white'),
    add.text = list(col = 'white'),
    par.main.text = list(col = 'white'),
    par.xlab.text = list(col = 'white'),
    par.ylab.text = list(col = 'white'),
    superpose.symbol = list(pch = 15, cex = 1, alpha = 0.5)
  ),
  xscale.components = xscale.components.log10.3,
  xlab = 'Spectral Distance', ylab = expression(~Delta*E['00']),
  auto.key = list(columns = 3, lines = FALSE, points = TRUE, cex = 1),
  panel = function(...) {
    panel.grid(h = -1, v = -1, lty = 3)
    panel.xyplot(...)
  }
)


This document is based on aqp version 2.0.4.