Perceptual Distances with coldist()

pavo implements the noise-weighted colour distance calculations in the function coldist(), assuming that raw receptor quantum catches are provided (via relative = FALSE in vismodel()). For the achromatic contrast, coldist() uses n4 to calculate \(w\) for the achromatic contrast. Note that even if \(Q_i\) is chosen, values are still log-transformed. This option is available in case the user wants to specify a data frame of quantum catches that was not generated by vismodel as an input. In this case, the argument qcatch should be used to inform the function if \(Q_i\) or \(f_i\) values are being used (note that if the input to coldist() is an object generated using the vismodel() function, this argument is ignored.) The type of noise to be calculated can be selected from the coldist() argument noise (which accepts either "neural" or "quantum").

coldist(vismod1,
  noise = "neural", achromatic = TRUE, n = c(1, 2, 2, 4),
  weber = 0.1, weber.achro = 0.1
)

##      patch1   patch2        dS        dL
## 1  cardinal   jacana  8.423965 3.3449627
## 2  cardinal   oriole  7.905876 3.5110820
## 3  cardinal parakeet  7.999730 0.5177495
## 4  cardinal    robin  5.805562 2.5606167
## 5  cardinal  tanager 10.100311 1.9068807
## 6    jacana   oriole 10.151253 0.1661193
## 7    jacana parakeet  4.231570 2.8272132
## 8    jacana    robin  3.503769 5.9055794
## 9    jacana  tanager  5.126266 1.4380820
## 10   oriole parakeet  8.135793 2.9933325
## 11   oriole    robin  7.220505 6.0716987
## 12   oriole  tanager 13.608253 1.6042013
## 13 parakeet    robin  4.606847 3.0783661
## 14 parakeet  tanager  5.969760 1.3891312
## 15    robin  tanager  7.688707 4.4674974

coldist(vismod.idi, n = c(1, 2), weber = 0.1)

##      patch1   patch2        dS
## 1  cardinal   jacana 2.0993283
## 2  cardinal   oriole 4.6567462
## 3  cardinal parakeet 2.2575457
## 4  cardinal    robin 3.7236780
## 5  cardinal  tanager 3.5417700
## 6    jacana   oriole 2.5574178
## 7    jacana parakeet 4.3568740
## 8    jacana    robin 1.6243496
## 9    jacana  tanager 5.6410983
## 10   oriole parakeet 6.9142918
## 11   oriole    robin 0.9330682
## 12   oriole  tanager 8.1985162
## 13 parakeet    robin 5.9812236
## 14 parakeet  tanager 1.2842244
## 15    robin  tanager 7.2654480

Where dS is the chromatic contrast (\(\Delta S\)) and dL is the achromatic contrast (\(\Delta L\)). Note that, by default, achromatic = FALSE, so dL isn’t included in the second result (this is to maintain consistency since, in the vismodel() function, the achromatic argument defaults to none). As expected, values are really high under the avian colour vision, since the colours of these species are quite different and because of the enhanced discriminatory ability with four compared to two cones.

coldist() also has a subset argument, which is useful if only certain comparisons are of interest (for example, of colour patches against a background, or only comparisons among a species or body patch). subset can be a vector of length one or two. If only one subsetting option is passed, all comparisons against the matching argument are returned (useful in the case of comparing to a background, for example). If two values are passed, comparisons will only be made between samples that match that rule (partial string matching and regular expressions are accepted). For example, compare:

coldist(vismod1, subset = 'cardinal')

to:

coldist(vismod1, subset = c('cardinal', 'jacana'))

Converting receptor noise-corrected distances to Cartesian coordinates

You can convert distances in JND back to Cartesian position coordinates with the function jnd2xyz(). Note that the actual position of these points in the XYZ space is arbitrary; therefore the function allows you to rotate the data so that, for example, the vector leading to the long-wavelength cone aligns with the X axis:

fakedata1 <- vapply(
  seq(100, 500, by = 20),
  function(x) rowSums(cbind(
      dnorm(300:700, x, 30),
      dnorm(300:700, x + 400, 30)
    )),
    numeric(401)
)

# Creating idealized specs with varying saturation
fakedata2 <- vapply(
  c(500, 300, 150, 105, 75, 55, 40, 30),
  function(x) dnorm(300:700, 550, x),
  numeric(401)
)

fakedata1 <- as.rspec(data.frame(wl = 300:700, fakedata1))

## wavelengths found in column 1

fakedata1 <- procspec(fakedata1, "max")

## processing options applied:
## Scaling spectra to a maximum value of 1

fakedata2 <- as.rspec(data.frame(wl = 300:700, fakedata2))

## wavelengths found in column 1

fakedata2 <- procspec(fakedata2, "sum")

## processing options applied:
## Scaling spectra to a total area of 1

fakedata2 <- procspec(fakedata2, "min")

## processing options applied:
## Scaling spectra to a minimum value of zero

# Converting reflectance to percentage
fakedata1[, -1] <- fakedata1[, -1] * 100
fakedata2[, -1] <- fakedata2[, -1] / max(fakedata2[, -1]) * 100

# Combining and converting to rspec
fakedata.c <- data.frame(wl = 300:700, fakedata1[, -1], fakedata2[, -1])
fakedata.c <- as.rspec(fakedata.c)

## wavelengths found in column 1

# Visual model and colour distances
fakedata.vm <- vismodel(fakedata.c, relative = FALSE, achromatic = TRUE)
fakedata.cd <- coldist(fakedata.vm,
  noise = "neural", n = c(1, 2, 2, 4),
  weber = 0.1, achromatic = TRUE
)

# Converting to Cartesian coordinates
fakedata.cc <- jnd2xyz(fakedata.cd, ref1 = "l", axis1 = c(1, 0, 0), ref2 = NULL)
head(fakedata.cc)

##             x          y         z      lum
## X1 -26.165302   3.692783 -17.71718 1.533258
## X2 -11.212399  -2.085407 -27.97937 1.533258
## X3   3.076291  -7.441594 -37.76761 1.533258
## X4  16.523287 -12.475488 -46.33606 1.533255
## X5  25.812150 -17.848152 -40.71621 1.533213
## X6  31.912643 -23.368560 -24.79660 1.532744

which you can then plot as well:

plot(fakedata.cc, theta = 55, phi = 25, col = spec2rgb(fakedata.c))

Spectral data in a receptor noise-corrected colourspace

The axes in this colourspace are in JND units. For more information on these functions, see ?jnd2xyz, ?jndrot and ?jndplot.

Di-, Tri-, and Tetrachromatic Colourspaces

pavo has extensive modelling and visualisation capabilities for generic di-, tri-, and tetra-chromatic spaces, uniting these approaches in a cohesive workflow. As with most colourspace models, we first estimate relative quantum catches with various assumptions by using the vismodel() function, before converting each set of values to a location in colourspace by using the space argument in colspace() (the function can also be set to try to detect the dimensionality of the colourspace automatically). For di- tri- and tetrachromatic spaces, colspace() calculates the coordinates of stimuli as:

Dichromats: \[x = \frac{1}{\sqrt{2}}(Q_l - Q_s)\]

Trichromats: \[x = \frac{1}{\sqrt{2}}(Q_l - Q_m)\] \[y = \frac{\sqrt{2}}{\sqrt{3}}(Q_s - \frac{Q_l + Q_m}{2})\]

Tetrachromats: \[x = \frac{1}{\sqrt{2}}(Q_l - Q_m)\] \[y = \frac{\sqrt{2}}{\sqrt{3}}(Q_s - \frac{Q_l + Q_m}{2})\] \[z = \frac{\sqrt{3}}{2}(Q_u - \frac{Q_l + Q_m + Q_s}{3})\]

Where Q_u, Q_s, Q_m, and Q_l refer to quantum catch estimates for UV-, short, medium-, and long-wavelength photoreceptors, respectively.

For a dichromatic example, we can model our floral reflectance data using the visual system of the domestic dog Canis familiaris, which has two cones with maximal sensitivity near 440 and 560 nm.

vis.flowers <- vismodel(flowers, visual = 'canis')

di.flowers <- colspace(vis.flowers, space = 'di')

head(di.flowers)

##                                s         l           x      r.vec
## Goodenia_heterophylla 0.52954393 0.4704561 -0.04178143 0.04178143
## Goodenia_geniculata   0.25430150 0.7456985  0.34747015 0.34747015
## Goodenia_gracilis     0.01747832 0.9825217  0.68238870 0.68238870
## Xyris_operculata      0.39433933 0.6056607  0.14942675 0.14942675
## Eucalyptus_sp         0.40628552 0.5937145  0.13253229 0.13253229
## Faradaya_splendida    0.23166580 0.7683342  0.37948187 0.37948187

The output contains values for the relative stimulation of short- and long-wavelength sensitive photoreceptors associated with each flower, along with its single coordinate in dichromatic space and its r.vector (distance from the origin). To visualise where these points lie, we can simply plot them on a segment.

plot(di.flowers, col = spec2rgb(flowers))

Flowers in a dichromatic colourspace, as modelled according to a canid visual system.

For our trichromatic viewer we can use the honeybee Apis mellifera, one of the most significant and widespread pollinators. We’ll also transform our quantum catches according to Fechner’s law by specifying qcatch = 'fi', and will model photoreceptor stimulation under bright conditions by scaling our illuminant with the scale argument.

vis.flowers <- vismodel(flowers, visual = 'apis', qcatch = 'fi', scale = 10000)

tri.flowers <- colspace(vis.flowers, space = 'tri')

head(tri.flowers)

##                               s         m         l            x
## Goodenia_heterophylla 0.2882383 0.3587254 0.3530363 -0.004022828
## Goodenia_geniculata   0.2475326 0.3549832 0.3974842  0.030052720
## Goodenia_gracilis     0.1373117 0.3125304 0.5501579  0.168028001
## Xyris_operculata      0.2395035 0.3729820 0.3875145  0.010275990
## Eucalyptus_sp         0.2760091 0.3548096 0.3691813  0.010162377
## Faradaya_splendida    0.2654951 0.3431266 0.3913783  0.034119147
##                                 y    h.theta      r.vec
## Goodenia_heterophylla -0.05522995 -1.6435057 0.05537626
## Goodenia_geniculata   -0.10508396 -1.2922439 0.10929687
## Goodenia_gracilis     -0.24007647 -0.9601417 0.29303604
## Xyris_operculata      -0.11491764 -1.4816131 0.11537616
## Eucalyptus_sp         -0.07020755 -1.4270471 0.07093922
## Faradaya_splendida    -0.08308452 -1.1811377 0.08981733

As in the case of our dichromat, the output contains relative photoreceptor stimulations, coordinates in the Maxwell triangle, as well as the ‘hue angle’ h.theta and distance from the origin (r.vec).

plot(tri.flowers, pch = 21, bg = spec2rgb(flowers))

Floral reflectance in a Maxwell triangle, considering a honeybee visual system.

Finally, we’ll draw on the blue tit’s visual system to model our floral reflectance spectra in a tetrahedral space, again using log-transformed quantum catches and assuming bright viewing conditions.

vis.flowers <- vismodel(flowers, visual = "bluetit", qcatch = "fi", scale = 10000)

tetra.flowers <- colspace(vis.flowers, space = "tcs")

head(tetra.flowers)

##                               u         s         m         l         u.r
## Goodenia_heterophylla 0.2274310 0.2659669 0.2524170 0.2541851 -0.02256897
## Goodenia_geniculata   0.1625182 0.2721621 0.2828016 0.2825182 -0.08748184
## Goodenia_gracilis     0.0428552 0.2443091 0.3504414 0.3623943 -0.20714480
## Xyris_operculata      0.1766234 0.2709785 0.2642637 0.2881344 -0.07337659
## Eucalyptus_sp         0.2118481 0.2609468 0.2637023 0.2635028 -0.03815190
## Faradaya_splendida    0.1844183 0.2557660 0.2786393 0.2811764 -0.06558169
##                                s.r         m.r         l.r            x
## Goodenia_heterophylla  0.015966893 0.002417012 0.004185066 -0.007214866
## Goodenia_geniculata    0.022162064 0.032801599 0.032518179  0.006341799
## Goodenia_gracilis     -0.005690920 0.100441373 0.112394348  0.072312163
## Xyris_operculata       0.020978454 0.014263689 0.038134446  0.010505857
## Eucalyptus_sp          0.010946788 0.013702317 0.013502794  0.001565228
## Faradaya_splendida     0.005765955 0.028639328 0.031176408  0.015560661
##                                  y           z    h.theta     h.phi
## Goodenia_heterophylla -0.005415708 -0.02256897 -2.4976873 -1.190529
## Goodenia_geniculata    0.003861848 -0.08748184  0.5469755 -1.486123
## Goodenia_gracilis      0.033297418 -0.20714480  0.4315247 -1.203879
## Xyris_operculata      -0.010813615 -0.07337659 -0.7998327 -1.368146
## Eucalyptus_sp          0.001044768 -0.03815190  0.5885699 -1.521510
## Faradaya_splendida     0.007189965 -0.06558169  0.4328380 -1.315140
##                            r.vec     r.max r.achieved
## Goodenia_heterophylla 0.02430520 0.2692325 0.09027589
## Goodenia_geniculata   0.08779638 0.2508989 0.34992737
## Goodenia_gracilis     0.22191605 0.2678272 0.82857920
## Xyris_operculata      0.07490949 0.2552227 0.29350636
## Eucalyptus_sp         0.03819828 0.2503039 0.15260760
## Faradaya_splendida    0.06778487 0.2583986 0.26232676

Tetrahedral data (via colspace(space = 'tcs')) may be plotted in a standard static tetrahedron by using plot(), or can be visualised as part of an interactive tetrahedron by using tcsplot() with the accessory functions tcspoints() and tcsvol() for adding points and convex hulls, respectively. As with other colourspace plots there are a number of associated graphical options, though the theta and phi arguments are particularly useful in this case, as they control the orientation (in degrees) of the tetrahedron in the xy and yz planes, respectively.

When plotting the tetrahedral colourspace, one can also force perspective by changing the size of points relative to their distance from the plane of observation, using the arguments perspective = TRUE and controlling the size range with the range argument. Several other options control the appearance of this plot, you can check these using ?plot.colspace or ?tetraplot

plot(tetra.flowers, pch = 21, bg = spec2rgb(flowers), perspective = TRUE, range = c(1, 2), cex = 0.5)

Flowers in a tetrahedral colourspace modelled using the visual phenotype of the blue tit. Point size is used to force perspective

Two additional functions may help with tetrahedral colourspace plotting:

• axistetra() function can be used to draw arrows showing the direction and magnitude of distortion of x, y and z in the tetrahedral plot.
• legendtetra() allows you to add a legend to a plot.

par(mfrow = c(1, 2), pty = "s")
plot(tetra.flowers, pch = 21, bg = spec2rgb(flowers))
axistetra(x = 0, y = 1.8)
plot(tetra.flowers, theta = 110, phi = 10, pch = 21, bg = spec2rgb(flowers))
axistetra(x = 0, y = 1.8)

Flowers in a tetrahedral colourspace, with varied orientations and perspectives, modelled using the visual phenotype of the blue tit.

For tetrahedral models, another plotting option available is projplot, which projects colour points in the surface of a sphere encompassing the tetrahedron. This plot is particularly useful to see differences in hue.

projplot(tetra.flowers, pch = 20, col = spec2rgb(flowers))

## Warning in seq.default(lim[1], lim[2], len = 100): partial argument match
## of 'len' to 'length.out'

## Warning in seq.default(lim[3], lim[4], len = 100): partial argument match
## of 'len' to 'length.out'

Projection plot from a tetrahedral colour space.

Finally, a useful function for tetrahedral models is voloverlap(), which calculates the overlap in tetrahedral colour volume between two sets of points. This can be useful to explore whether different species occupy similar (overlapping) or different (non- overlapping) “sensory niches”, or to test for mimetism, dichromatism, etc. [11]. To show this function, we will use the sicalis dataset, which includes measurements from the crown (C), throat (T) and breast (B) of seven stripe-tailed yellow finches (Sicalis citrina).

data(sicalis)

We will use this dataset to test for the overlap between the volume determined by the measurements of those body parts from multiple individuals in the tetrahedral colourspace (note the option plot for plotting of the volumes):

par(mfrow = c(1, 2), pty = "s")
tcs.sicalis.C <- subset(colspace(vismodel(sicalis)), "C")
tcs.sicalis.T <- subset(colspace(vismodel(sicalis)), "T")
tcs.sicalis.B <- subset(colspace(vismodel(sicalis)), "B")
voloverlap(tcs.sicalis.T, tcs.sicalis.B, plot = TRUE)

##           vol1         vol2   overlapvol vsmallest      vboth
## 1 5.183721e-06 6.281511e-06 6.904074e-07 0.1331876 0.06407598

voloverlap(tcs.sicalis.T, tcs.sicalis.C, plot = TRUE)

##           vol1         vol2 overlapvol vsmallest vboth
## 1 5.183721e-06 4.739152e-06          0         0     0

The function voloverlap() gives the volume (\(V\)) of the convex hull delimited by the overlap between the two original volumes, and two proportions are calculated from that: \[\mathtt{vsmallest} = V_{overlap} / min(V_A, V_B)\] and \[\mathtt{vboth} = V_{overlap} / (V_A + V_B)\]. Thus, if one of the volumes is entirely contained in the other, vsmallest will equal 1.

So we can clearly see that there is overlap between the throat and breast colours (of about 6%), but not between the throat and the crown colours (Figures above).

Summary variables for groups of points

Another advantage of colourspace models is that they allow for the calculation of useful summary statistics of groups of points, such as the centroid of the points, the total volume occupied, the mean and variance of hue span and the mean saturation. In pavo, the result of a colspace() call is an object of class colspace, and thus these summary statistics can be calculated simply by calling summary. Note that the summary call can also take a by vector of group identities, so that the variables are calculated for each group separately:

summary(tetra.flowers)

## Colorspace & visual model options:
##  * Colorspace: tcs 
##  * Quantal catch: fi 
##  * Visual system, chromatic: bluetit 
##  * Visual system, achromatic: none 
##  * Illuminant: ideal, scale = 10000 (von Kries color correction not applied) 
##  * Background: ideal 
##  * Relative: TRUE 
## 

##            centroid.u centroid.s centroid.m centroid.l        c.vol
## all.points  0.1948478  0.2554421  0.2704411   0.279269 0.0002510915
##              rel.c.vol  colspan.m   colspan.v huedisp.m huedisp.v
## all.points 0.001159742 0.05706825 0.001682224 0.6413439 0.3157949
##              mean.ra    max.ra
## all.points 0.2428508 0.8285792

The Colour Hexagon

The hexagon colour space of [4] is a generalised colour-opponent model of hymenopteran vision that has found extremely broad use, particularly in studies of bee-flower interactions. It’s also often broadly applied across hymenopteran species, because the photopigments underlying trichromatic vision in Hymenoptera appear to be quite conserved ???. What’s particularly useful is that colour distances within the hexagon have been extensively validated against behaviour, and thus offer a relatively reliable measure of perceptual distance.

In the hexagon, photoreceptor quantum catches are typically hyperbolically transformed (and pavo will return a warning if the transform is not selected), and vonkries correction is often used used to model photoreceptor adaptation to a vegetation background. This can all now be specified in vismodel(). including the optional use of a ‘green’ vegetation background. Note that although this is a colourspace model, we specific relative = FALSE to return raw (albeit transformed) quantum catches, as required for the model.

vis.flowers <- vismodel(flowers, visual = 'apis', qcatch = 'Ei', relative = FALSE, vonkries = TRUE, achromatic = 'l', bkg = 'green')

We can then apply the hexagon model in colspace, which will convert our photoreceptor ‘excitation values’ to coordinates in the hexagon according to:

\[x = \frac{\sqrt{3}}{2(E_g + E_{uv})}\]

\[y = E_b - 0.5(E_{uv} + E_g)\]

hex.flowers <- colspace(vis.flowers, space = 'hexagon')

head(hex.flowers)

##                                s         m         l         x          y
## Goodenia_heterophylla 0.56572097 0.8237141 0.7053057 0.1208839  0.1882007
## Goodenia_geniculata   0.35761236 0.8176153 0.8670134 0.4411542  0.2053024
## Goodenia_gracilis     0.01888788 0.1622766 0.7810589 0.6600594 -0.2376968
## Xyris_operculata      0.21080752 0.7345122 0.6796464 0.4060263  0.2892853
## Eucalyptus_sp         0.55622758 0.8515289 0.8208038 0.2291297  0.1630132
## Faradaya_splendida    0.45855056 0.7828905 0.8565895 0.3447118  0.1253205
##                         h.theta     r.vec sec.fine sec.coarse       lum
## Goodenia_heterophylla  32.71320 0.2236793       30  bluegreen 0.1680153
## Goodenia_geniculata    65.04388 0.4865862       60  bluegreen 0.3548825
## Goodenia_gracilis     109.80467 0.7015541      100      green 0.2313678
## Xyris_operculata       54.53092 0.4985412       50  bluegreen 0.1518319
## Eucalyptus_sp          54.57024 0.2812005       50  bluegreen 0.2787543
## Faradaya_splendida     70.02119 0.3667853       70  bluegreen 0.3351008

Again, the output includes the photoreceptor excitation values for short- medium- and long-wave sensitive photoreceptors, x and y coordinates, and measures of hue and saturation for each stimulus. The hexagon model also outputs two additional measures of of subjective ‘bee-hue’; sec.fine and sec.coarse. sec.fine describes the location of stimuli within one of 36 ‘hue sectors’ that are specified by radially dissecting the hexagon in 10-degree increments. sec.coarse follows a similar principle, though here the hexagon is divided into only five ‘bee-hue’ sectors: UV, UV-blue, blue, blue-green, green, and UV-green. These can easily be visualised by specifying sectors = 'coarse or sectors = 'fine' in a call to plot after modelling.

plot(hex.flowers, sectors = 'coarse', pch = 21, bg = spec2rgb(flowers))

Flowers as modelled in the hymenopteran colour hexagon of Chittka (1992), overlain with coarse bee-hue sectors.

The Colour Opponent Coding (COC) Space

The colour opponent coding (coc) space is an earlier hymenopteran visual model [5] which, although now seldom used in favour of the hexagon, may prove useful for comparative work. While the initial estimation of photoreceptor excitation is similar to that in the hexagon, the coc subsequently specifies A and B coordinates based on empirically-derived weights for the output from each photoreceptor:

\[A = -9.86E_u + 7.70E_b + 2.16E_g\] \[B = -5.17E_u + 20.25E_b - 15.08E_g\]

Where \(E_i\) is the excitation value (quantum catch) in photoreceptor \(i\).

vis.flowers <- vismodel(flowers, visual = "apis", qcatch = "Ei", relative = FALSE, vonkries = TRUE, bkg = "green")

coc.flowers <- colspace(vis.flowers, space = "coc")

head(coc.flowers)

##                                s         m         l        x          y
## Goodenia_heterophylla 0.56572097 0.8237141 0.7053057 2.288050  3.1194224
## Goodenia_geniculata   0.35761236 0.8176153 0.8670134 4.642329  1.6332926
## Goodenia_gracilis     0.01888788 0.1622766 0.7810589 2.750382 -8.5899171
## Xyris_operculata      0.21080752 0.7345122 0.6796464 5.045218  3.5349309
## Eucalyptus_sp         0.55622758 0.8515289 0.8208038 2.845304  1.9900420
## Faradaya_splendida    0.45855056 0.7828905 0.8565895 3.357182  0.5654568
##                           r.vec
## Goodenia_heterophylla  5.407473
## Goodenia_geniculata    6.275622
## Goodenia_gracilis     11.340299
## Xyris_operculata       8.580149
## Eucalyptus_sp          4.835346
## Faradaya_splendida     3.922638

The A and B coordinates are designated x and y in the output of coc for consistency, and while the model includes a measure of saturation in r.vec, it contains no associated measure of hue.

plot(coc.flowers, pch = 21, bg = spec2rgb(flowers), yaxt = "n")

Flowers in the colour-opponent-coding space of Backhaus (1991), as modelling according to the honeybee.

CIE Spaces

The CIE (International Commission on Illumination) colourspaces are a suite of models of human colour vision and perception. pavo 1.0 now includes two of the most commonly used: the foundational 1931 CIE XYZ space, and the more modern, perceptually calibrated CIE LAB space and its cylindrical CIE LCh transformation. Tristimulus values in XYZ space are calculated as:

\[X = k\int_{300}^{700}{R(\lambda)I(\lambda)\bar{x}(\lambda)\,d\lambda}\] \[Y = k\int_{300}^{700}{R(\lambda)I(\lambda)\bar{y}(\lambda)\,d\lambda}\] \[Z = k\int_{300}^{700}{R(\lambda)I(\lambda)\bar{z}(\lambda)\,d\lambda}\]

where x, y, and z are the trichromatic colour matching functions for a ‘standard colourimetric viewer’. These functions are designed to describe an average human’s chromatic response within a specified viewing arc in the fovea (to account for the uneven distribution of cones across eye). pavo includes both the CIE 2-degree and the modern 10-degree standard observer, which can be selected in the visual option in the vismodel function. In these equations, k is the normalising factor

\[k = \frac{100}{\int_{300}^{700}{I(\lambda)\bar{y}(\lambda)\,d\lambda}}\]

and the chromaticity coordinates of stimuli are calculated as

\[x = \frac{X}{X + Y + Z}\] \[y = \frac{Y}{X + Y + Z}\] \[z = \frac{Z}{X + Y + Z} = 1 - x - y\]

For modelling in both XYZ and LAB spaces, here we’ll use the CIE 10-degree standard observer, and assume a D65 ‘standard daylight’ illuminant. Again, although they are colourspace models, we need to set relative = FALSE to return raw quantum catch estimates, vonkries = TRUE to account for the required normalising factor (as above), and achromatic = 'none' since there is no applicable way to estimate luminance in the CIE models. As with all models that have particular requirements, vismodel() will output warnings and/or errors if unusual or non-standard arguments are specified.

vis.flowers <- vismodel(flowers, visual = 'cie10', illum = 'D65', vonkries = TRUE, relative = FALSE, achromatic = 'none')

ciexyz.flowers <- colspace(vis.flowers, space = 'ciexyz')
head(ciexyz.flowers)

##                               X         Y         Z         x         y
## Goodenia_heterophylla 0.2074305 0.2067084 0.3013614 0.2899097 0.2889005
## Goodenia_geniculata   0.6472803 0.6752456 0.4240380 0.3706021 0.3866137
## Goodenia_gracilis     0.5153064 0.4587957 0.0153714 0.5207885 0.4636766
## Xyris_operculata      0.2884446 0.2353124 0.2224271 0.3865595 0.3153544
## Eucalyptus_sp         0.4431842 0.4470468 0.4019893 0.3429634 0.3459524
## Faradaya_splendida    0.6636975 0.6553481 0.2869075 0.4132733 0.4080743
##                                z
## Goodenia_heterophylla 0.42118974
## Goodenia_geniculata   0.24278414
## Goodenia_gracilis     0.01553492
## Xyris_operculata      0.29808613
## Eucalyptus_sp         0.31108420
## Faradaya_splendida    0.17865248

The output is simply the tristimulus values and chromaticity coordinates of stimuli, and we can visualise our results (along with a line connecting monochromatic loci, by default) by calling

plot(ciexyz.flowers, pch = 21, bg = spec2rgb(flowers))

Floral reflectance in the CIEXYZ human visual model. Note that this space is not perceptually calibrated, so we cannot make inferences about the similarity or differences of colours based on their relative location.

The Lab space is a more recent development, and is a colour-opponent model that attempts to mimic the nonlinear responses of the human eye. The Lab space is also calibrated (unlike the XYZ space) such that Euclidean distances between points represent relative perceptual distances. This means that two stimuli that are farther apart in Lab space should be perceived as more different than two closer points. As the name suggests, the dimensions of the Lab space are Lightness (i.e. subjective brightness), along with two colour-opponent dimensions designated a and b. The colspace() function, when space = cielab, simply converts points from the XYZ model according to:

\[ L=\begin{cases} 116\left(\frac{Y}{Y_n}\right)^\frac{1}{3} & \text{if } \frac{Y}{Y_n} > 0.008856 \\ 903.3\left(\frac{Y}{Y_n}\right) & \text{if } \frac{Y}{Y_n} \leq 0.008856 \end{cases} \]

\[a = 500\left(f\left(\frac{X}{X_n}\right) - f\left(\frac{Y}{Y_n}\right)\right)\]

\[b = 500\left(f\left(\frac{Y}{Y_n}\right) - f\left(\frac{Z}{Z_n}\right)\right)\]

where

\[ f(x)=\begin{cases} x^\frac{1}{3} & \text{if } x > 0.008856 \\ 7.787x + \frac{4}{29} & \text{if } x \leq 0.008856 \end{cases} \]

Here, \(X_n, Y_n, Z_n\) are neutral point values to model visual adaptation, calculated as:

\[X_n = \int_{300}^{700}{R_n(\lambda)I(\lambda)\bar{x}(\lambda)\,d\lambda}\] \[Y_n = \int_{300}^{700}{R_n(\lambda)I(\lambda)\bar{y}(\lambda)\,d\lambda}\] \[Z_n = \int_{300}^{700}{R_n(\lambda)I(\lambda)\bar{z}(\lambda)\,d\lambda}\]

when \(R_n(\lambda)\) is a perfect diffuse reflector (i.e. 1).

cielab.flowers <- colspace(vis.flowers, space = 'cielab')
head(cielab.flowers)

##                               X         Y         Z        L           a
## Goodenia_heterophylla 0.2074305 0.2067084 0.3013614 16.91253  0.91908831
## Goodenia_geniculata   0.6472803 0.6752456 0.4240380 55.24744 -7.78898966
## Goodenia_gracilis     0.5153064 0.4587957 0.0153714 37.53788 21.57990123
## Xyris_operculata      0.2884446 0.2353124 0.2224271 19.25286 19.66170677
## Eucalyptus_sp         0.4431842 0.4470468 0.4019893 36.57660  0.05744138
## Faradaya_splendida    0.6636975 0.6553481 0.2869075 53.61946  5.07048199
##                                b
## Goodenia_heterophylla -24.336189
## Goodenia_geniculata    19.981497
## Goodenia_gracilis      61.991586
## Xyris_operculata       -6.289287
## Eucalyptus_sp          -8.295538
## Faradaya_splendida     41.517115

Our output now contains the tristimulus XYZ values, as well as their Lab counterparts, which are coordinates in the Lab space. These can also be visualised in three-dimensional Lab space by calling plot:

plot(cielab.flowers, pch = 21, bg = spec2rgb(flowers))

Floral reflectance spectra represented in the CIELab model of human colour sensation.

Categorical Fly Colourspace

The categorical colour vision model of ??? is a model of dipteran vision, based on behavioural data from the blowfly Lucilia sp. It assumes the involvement of all four dipteran photoreceptor classes (R7 & R8 ‘pale’ and ‘yellow’ subtypes), and further posits that colour vision is based on two specific opponent mechanisms (R7p - R8p, and R7y - R8y). The model assumes that all colours are perceptually grouped into one of four colour categories, and that flies are unable to distinguish between colours that fall within the same category.

We’ll use the visual systems of the muscoid fly Musca domestica, and will begin by estimating linear (i.e. untransformed) quantum catches for each of the four photoreceptors.

vis.flowers <- vismodel(flowers, qcatch = 'Qi', visual = 'musca', achromatic = 'none', relative = TRUE)

Our call to colspace() will then simply estimate the location of stimuli in the categorical space as the difference in relative stimulation between ‘pale’ (R7p - R8p) and ‘yellow’ (R7y - R8y) photoreceptor pairs:

\[x = R7_p - R8_p\]

\[y = R7_y - R8_y\]

cat.flowers <- colspace(vis.flowers, space = 'categorical')

head(cat.flowers)

##                              R7p        R7y       R8p       R8y          x
## Goodenia_heterophylla 0.05387620 0.18710675 0.4335989 0.3254182 -0.3797227
## Goodenia_geniculata   0.02080100 0.08276037 0.3737596 0.5226791 -0.3529586
## Goodenia_gracilis     0.01153133 0.02501422 0.1420052 0.8214492 -0.1304739
## Xyris_operculata      0.01270355 0.12523436 0.4488383 0.4132238 -0.4361348
## Eucalyptus_sp         0.03461951 0.14207857 0.3977211 0.4255808 -0.3631016
## Faradaya_splendida    0.03599956 0.09188569 0.3129199 0.5591949 -0.2769203
##                                y category     r.vec
## Goodenia_heterophylla -0.1383114     p-y- 0.4041279
## Goodenia_geniculata   -0.4399187     p-y- 0.5640108
## Goodenia_gracilis     -0.7964350     p-y- 0.8070515
## Xyris_operculata      -0.2879894     p-y- 0.5226389
## Eucalyptus_sp         -0.2835023     p-y- 0.4606694
## Faradaya_splendida    -0.4673092     p-y- 0.5431968

And it is simply the signs of these differences that define the four possible fly-colour categories (p+y+, p-y+, p+y-, p-y-), which we can see in the associated plot.

plot(cat.flowers, pch = 21, bg = spec2rgb(flowers))

Flowers in the categorical colourspace of Troje (1993).

Segment Classification

The segment classification analysis [3] does not assume any particular visual system, but instead tries to classify colours in a manner that captures common properties of many vertebrate (and some invertebrate) visual systems. In essence, it breaks down the reflectance spectrum region of interest into four equally-spaced regions, measuring the relative signal along those regions. This approximates a tetrachromatic system with ultraviolet, short, medium, and long-wavelength sensitive photoreceptors.

Though somewhat simplistic, this model captures many of the properties of other, more complex visual models, but without many of the additional assumptions these make. It also provides results in a fairly intuitive colour space, in which the angle corresponds to hue and the distance from the centre corresponds to chroma (Figure below; in fact, variables S5 and H4 from summary.rspec() are calculated from these relative segments). Note that, while a segment analysis ranging from 300 or 400 nm to 700 nm corresponds quite closely to the human visual system colour wheel, any wavelength range can be analysed in this way, returning a 360° hue space delimited by the range used.

The segment differences or “opponents” are calculated as:

\[LM = \frac{ R_\lambda \sum_{\lambda={Q4}} R_\lambda - \sum_{\lambda={Q2}} R_\lambda }{\sum_{\lambda={min}}^{max}R_\lambda}\] \[MS = \frac{ R_\lambda \sum_{\lambda={Q3}} R_\lambda - \sum_{\lambda={Q1}} R_\lambda }{\sum_{\lambda={min}}^{max}R_\lambda}\]

Where \(Qi\) represent the interquantile distances (e.g. \(Q1=ultraviolet\), \(Q2=blue\), \(Q3=green\) and \(Q4=red\))

The segment classification model is obtained through the colspace() argument space = 'segment', following initial modelling with vismodel(). The example below uses idealized reflectance spectra to illustrate how the avian colour space defined from the segment classification maps to the human colour wheel:

fakedata1 <- vapply(
  seq(100, 500, by = 20),
  function(x) rowSums(cbind(
      dnorm(300:700, x, 30),
      dnorm(300:700, x + 400, 30)
    )), numeric(401)
)

# creating idealized specs with varying saturation
fakedata2 <- vapply(
  c(500, 300, 150, 105, 75, 55, 40, 30),
  function(x) dnorm(300:700, 550, x),
  numeric(401)
)

fakedata1 <- as.rspec(data.frame(wl = 300:700, fakedata1))

## wavelengths found in column 1

fakedata1 <- procspec(fakedata1, "max")

## processing options applied:
## Scaling spectra to a maximum value of 1

fakedata2 <- as.rspec(data.frame(wl = 300:700, fakedata2))

## wavelengths found in column 1

fakedata2 <- procspec(fakedata2, "sum")

## processing options applied:
## Scaling spectra to a total area of 1

fakedata2 <- procspec(fakedata2, "min")

## processing options applied:
## Scaling spectra to a minimum value of zero

# converting reflectance to percentage
fakedata1[, -1] <- fakedata1[, -1] * 100
fakedata2[, -1] <- fakedata2[, -1] / max(fakedata2[, -1]) * 100

# combining and converting to rspec
fakedata.c <- data.frame(wl = 300:700, fakedata1[, -1], fakedata2[, -1])
fakedata.c <- as.rspec(fakedata.c)

## wavelengths found in column 1

# segment classification analysis
seg.vis <- vismodel(fakedata.c, visual = "segment", achromatic = "all")

## Warning: segment analysis chosen, overriding incompatible parameters.

seg.fdc <- colspace(seg.vis, space = "segment")

# plot results
plot(seg.fdc, col = spec2rgb(fakedata.c))

Idealized reflectance spectra and their projection on the axes of segment classification

Colourspace Distances with coldist

Under the colour space framework, colour distances can be calculated simply as Euclidean distances of the relative cone stimulation data, either log-transformed or not, depending on how it was defined. However, these distances cannot be interpreted in terms of JNDs, since no receptor noise is incorporated in the model. Euclidean distances can be computed in R using the coldist() function on the colspace() output:

head(coldist(tetra.flowers))

## Warning: Quantum catch are relative, distances may not be meaningful

##                  patch1              patch2        dS
## 1 Goodenia_heterophylla Goodenia_geniculata 0.4389066
## 2 Goodenia_heterophylla   Goodenia_gracilis 1.4421739
## 3 Goodenia_heterophylla    Xyris_operculata 0.3859043
## 4 Goodenia_heterophylla       Eucalyptus_sp 0.1383047
## 5 Goodenia_heterophylla  Faradaya_splendida 0.3640816
## 6 Goodenia_heterophylla  Gaultheria_hispida 0.2588980

Specialised colourspace models such as the colour-hexagon, coc space, CIELAB and CIELCh models have their own distance measures, which are returned be default when run through coldist(). Be sure to read the ?coldist documentation, and the original publications, to understand what is being returned. As an example, if we run the results of colour-hexagon modelling through coldist():

# Model flower colours according to a honeybee
vis.flowers <- vismodel(flowers, visual = "apis", qcatch = "Ei", relative = FALSE, vonkries = TRUE, achromatic = "l", bkg = "green")
hex.flowers <- colspace(vis.flowers, space = "hexagon")

# Estimate colour distances. No need to specify relative receptor densities, noise etc.,
# which only apply in the case of receptor-noise modelling
dist.flowers <- coldist(hex.flowers)
head(dist.flowers)

##                  patch1              patch2        dS
## 1 Goodenia_heterophylla Goodenia_geniculata 0.3207265
## 2 Goodenia_heterophylla   Goodenia_gracilis 0.6870945
## 3 Goodenia_heterophylla    Xyris_operculata 0.3025298
## 4 Goodenia_heterophylla       Eucalyptus_sp 0.1111375
## 5 Goodenia_heterophylla  Faradaya_splendida 0.2324927
## 6 Goodenia_heterophylla  Gaultheria_hispida 0.2280836

the chromatic contrasts dS and achromatic contrasts dL are expressed as Euclidean distances in the hexagon, known as ‘hexagon units’ [4]. If we had instead used the colour-opponent coding space the units would have been city-bloc distances, while in the CIELAB model the distances would be derived from the CIE colour-distance formula (2000).

Distances in N-Dimensions

The coldist() function is no longer limited to di-, tri- or tetrachromatic visual systems. The function has been generalised, and can now calculate colour distances for n-dimensional visual phenotypes. This means there is no limit on the dimensionality of visual systems that may be input, which may prove useful for modelling nature’s extremes (e.g. butterflies, mantis shrimp) or for simulation-based work. Naturally, since these calculations aren’t largely implemented elsewhere, we recommend caution and validation of results prior to publication.

# Create an arbitrary visual phenotype with 10 photoreceptors
fakemantisshrimp <- sensmodel(c(325, 350, 400, 425, 450, 500, 550, 600, 650, 700), beta = FALSE, integrate = FALSE)

# Convert to percentages, just to colour the plot
fakemantisshrimp.colours <- fakemantisshrimp * 100
fakemantisshrimp.colours[, "wl"] <- fakemantisshrimp[, "wl"]

plot(fakemantisshrimp, col = spec2rgb(fakemantisshrimp.colours), lwd = 2, ylab = "Absorbance")

Visual system of a pretend mantis shrimp with 10 cones

# Run visual model and calculate colour distances
vm.fms <- vismodel(flowers, visual = fakemantisshrimp, relative = FALSE, achromatic = FALSE)

JND.fms <- coldist(vm.fms, n = c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5))

head(JND.fms)

##                  patch1              patch2        dS
## 1 Goodenia_heterophylla Goodenia_geniculata 16.297557
## 2 Goodenia_heterophylla   Goodenia_gracilis 48.390309
## 3 Goodenia_heterophylla    Xyris_operculata 30.760482
## 4 Goodenia_heterophylla       Eucalyptus_sp  6.044184
## 5 Goodenia_heterophylla  Faradaya_splendida 15.901435
## 6 Goodenia_heterophylla  Gaultheria_hispida 11.188246