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)
)

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])

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))

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')
  )

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)

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)

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)

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

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 Additional Examples

2.1 Tinkering with plotColorMixture

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

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

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

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

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

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

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

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

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

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

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')

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

Top 3 mixture candidates

plotColorMixture(c('5B 5/10', '5Y 8/8'), label.cex = 0.65, showMixedSpec = TRUE, mixingMethod = 'reference', n = 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')

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, length.out = 6)
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))

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))

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))

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.46 NA dE00 exact
5YR 4/3 0.01 0.46 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'))

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')

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')

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')

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)

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'))

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)

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)

knitr::kable(t(quantile(d.adj)), digits = 3)
0% 25% 50% 75% 100%
0.001 0.013 0.026 0.063 0.253
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)

##
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())

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.2.