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:
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:
Key assumptions include:
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
= ?.
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)
)
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')
)
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.
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)
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)
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')
# 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))
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))
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'))
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'))
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)
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 | 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)
##
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.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.