Dental Shade Guides"
In munsellinterpol: Interpolate Munsell Renotation Data from Hue/Chroma to CIE/RGB

knitr::opts_chunk$set(echo = TRUE)
options( width=100 )

A Dental Shade Guide is a set of simulated teeth used to select prosthetic teeth by color. The simulated teeth are made of plastic or porcelain. Commercial shade guides have existed for almost a century. In 1933, Clark @CLARK1933 discussed the manufacture and use of porcelain shade guides based on the cylindrical color dimensions of hue, brilliance, and saturation, which correspond to Munsell's Hue, Value, and Chroma.

In @OBRIEN1989, the 24 teeth in a master Bioform shade guide were measured with a spectrophotometer. The spectral measurements were converted to xyY assuming Illuminant C, and then to both CIE Lab and Munsell HVC. The goal of this vignette is to plot the 24 colors as square patches, and then check the calculations from the article.

The spectrophotometer was an ACTA CIII from Beckman Instruments. Relative reflectance data were recorded from 410 to 700nm at 10-nm intervals. The conversion from such data to xyY and Lab is standard, but conversion to Munsell HVC is not explained in the paper. However, in a similar article @OBRIEN1990 by the same authors they state that:

The chromaticity coordinates were converted to Munsell notation by means of graphs ... and the method described by ASTM standard D 1535-80.

so we assume that the same method was used here. Note that the graphical method only applies to Munsell Hue and Chroma; for Munsell Value lookup tables are available in @ASTM-D1535-97. Below we check the published conversions against numerical conversions using munsellinterpol.

Load the required R packages.

library(munsellinterpol)
library(spacesRGB)  # for converting to RGB and plotting the patches
library(spacesXYZ)  # for xyY<->XYZ and Chromatic Adaptation Transform

Featured functions from munsellinterpol are XYZtoMunsell(), MunsellNameFromHVC(), NickersonColorDifference(), and ColorBlockFromMunsell().

Reading and Plotting the Data

Read the published data table.

path = system.file( 'extdata/dental.txt', package='munsellinterpol' )
dental = read.table( path, header=TRUE, sep='\t', stringsAsFactors=FALSE )
dental

Extract xyY, adapt from Illuminant C to D65, convert XYZ to sRGB, and display as a 6x4 grid of patches.

xyY = as.matrix( dental[ c('x','y','Y') ] )
XYZ = XYZfromxyY( xyY ) / 100
# adapt from Illuminant C to the whitepoint of sRGB, which is D65
# make the Chromatic Adaptation Transform
theCAT = spacesXYZ::CAT( 'C', getWhiteXYZ('sRGB',which='display') )
XYZ = adaptXYZ( theCAT, XYZ )
# create data.frame obj for plotting
obj = expand.grid( LEFT=1:6, TOP=1:4 )
obj$WIDTH = 0.9
obj$HEIGHT = 0.9
obj$RGB = RGBfromXYZ( XYZ, space='sRGB' )$RGB     # convert to sRGB
rownames(obj) = rownames(dental)
# plot as square patches
par( omi=c(0,0,0,0), mai=c(0.1,0.1,0.1,0.1) )
plotPatchesRGB( obj, which='signal', labels="bottomleft", adj=c(-0.2,-0.5), cex=0.7 )

This figure is best viewed on a display calibrated for sRGB.

Checking the Conversions

We now recompute Lab and Munsell values, and check against the published values.

Lab  = as.matrix( dental[ c('L','a','b') ] )
XYZ = XYZfromxyY( xyY )
Lab2 = LabfromXYZ( XYZ/100, 'C' )     # recompute Lab
HVC  = HVCfromMunsellName( dental$Munsell )
HVC2 = XYZtoMunsell( XYZ )            # recompute Munsell HVC
comp = data.frame( row.names=rownames(dental) )
comp$Y  = dental$Y
comp$L  = Lab[ ,1]
comp$L2 = round(Lab2[ ,1],4)
comp$Ldiff = round( comp$L - comp$L2, 4 )
comp$DeltaE = round( DeltaE( Lab, Lab2 ), 4 )   # DeltaE is the pairwise color difference
comp$Munsell = dental$Munsell
comp$Munsell2 = MunsellNameFromHVC( HVC2, format='f', digits=2 )
comp$NickersonCD = round( NickersonColorDifference( HVC, HVC2 ), 4 )
comp

The Lab agreement is good, but the published Lightness values are consistently too large and the reason for this is unknown. The exception is B-83 whose published Lightness=71.52 is too small and with the largest DeltaE by far; it appears to be a transcription error.
The Munsell agreement is not bad, but the pubished Munsell Value is too small. This could be due to using magnesium oxide instead of the perfect reflecting diffuser. We can test this by recomputing the value component of HVC2. The newly recomputed Munsell notation is denoted Munsell3.

HVC3 = HVC2
HVC3[ ,2] = VfromY( dental$Y, which='MgO' )
comp$Munsell2 = NULL
comp$NickersonCD = NULL
comp$Munsell3 = MunsellNameFromHVC( HVC3, format='f', digits=2 )
comp$NickersonCD = round( NickersonColorDifference( HVC, HVC3 ), 4 )
comp

The Munsell Value agreement is now much better. The worst Nickerson difference is for B-65 and this is largely because the published Munsell Value is 7.30 instead of 7.40. Note that the Y values for B-65 and B-54 are almost identical and that the Munsell Value for B-54 is correct. So again I think that the problem with B-65 is transcription error. For Munsell Hue and Chroma the agreement is good. It must have been tedious to use the graphical method for all 24 samples.

ISCC-NBS Color Names

I thought it would be interesting to display the ISCC-NBS names for each of the 24 dental shades.

obj = data.frame( row.names=rownames(dental) )
obj$Munsell2 = MunsellNameFromHVC( HVC2, format='f', digits=2 )
block = ColorBlockFromMunsell( HVC2 )
obj[[ "ISCC-NBS Name" ]] = block$Name
obj

All the dental shades are in the same block, except for 2. It would be interesting to turn this into a 3D scatterplot, with the color block boundaries displayed with transparency.