Load data

# For calculation, I use wavelengths between 390 and 780nm in steps of 2nm
wavelengths <- seq(390, 780, 2)

# Read data from Manchester
manchester <- data.table(read.xlsx(here("other_data", "all manchester spectra.xlsx")))
names(manchester)[1] <- "Wavelength"
fundamentals_manchester <- manchester[Wavelength >= 390, .(Wavelength, R, S, M, L)]
names(fundamentals_manchester) <- c("Wavelength", "rod", "scone", "mcone", "lcone")

led_spectra <- manchester[Wavelength >= 390, .(Wavelength, Red, Green, Blue, Amber)]

# Read fundamentals from Erlangen
load(here("data", "photoreceptor_sensitivities.rda"))
fundamentals_cord <- data.table(Wavelength = wavelengths, 
                                photoreceptor_spectral_sensitivities)
fundamentals_cord$melanopsin <- NULL

Compare fundamentals graphically

The red lines are the fundamentals from Erlangen, the black lines from Manchester.

ggplot(melt(fundamentals_manchester, id.vars = "Wavelength"), 
       aes(x = Wavelength, y = value, group = variable)) +
  geom_line() +
  geom_line(data = melt(fundamentals_cord, id.vars = "Wavelength"), color = "red")

Comparison with the Stockman fundamentals

There is also a Stockman & Sharpe (2000) table (“2-deg fundamentals based on the Stiles & Burch 10-deg CMFs (adjusted to 2-deg)”) on cvrl.org. The cone fundamentals from Jan’s Excel sheet are apparently identical to the ones from Manchester and similar but not identical to these Stockman&Sharpe 2° fundamentals (green). The differ at lower wavelengths.

data.table(Wavelength = wavelengths, analyzeTCS::ConeFundamentals2) %>%
  melt(id.vars = "Wavelength") %>%
  ggplot(aes(x = Wavelength, y = value, group = variable)) +
  geom_line(color = "green") +
  geom_line(data = melt(fundamentals_manchester, id.vars = "Wavelength"))

LED spectra from Manchester

ggplot(melt(led_spectra, id.vars = "Wavelength"), 
       aes(x = Wavelength, y = value, group = variable)) +
  geom_line() 

I’m not sure what luminances to use for calculating the spectra. However, this should not matter, so I will just pick 10cd/m^2 for each LED.

led_luminances <- c(Red = 10, Amber = 10, Green = 10, Blue = 10)

p_led <- 
  led_spectra %>%
  get_normalized_spectrum_matrix() %>%
  calculate_power_spectra(luminances = led_luminances, spectra = .)

This results in the following LED powers [W/m^2].

apply(p_led, 2, sum)
##         Red       Amber       Green        Blue 
## 0.030810957 0.006616563 0.015643687 0.364513300

Calculate A-Matrices

Fundamentals from Erlangen

print(create_A_matrix)
## function (P_led, s_receptors) 
## {
##     stopifnot(is_valid_spectra_matrix(P_led), is.matrix(s_receptors))
##     a_matrix <- crossprod(P_led, s_receptors)
##     return(t(t(a_matrix)/apply(a_matrix, 2, sum)))
## }

The main part of this function is: a_matrix <- crossprod(P_led, s_receptors)

s_receptors <- 
  fundamentals_cord %>%
  get_normalized_spectrum_matrix() 

p_led <- 
  led_spectra %>%
  get_normalized_spectrum_matrix() %>%
  calculate_power_spectra(luminances = led_luminances, spectra = .)

a_matrix_cord <- create_A_matrix(P_led = p_led, s_receptors = s_receptors)

print(a_matrix_cord)
##            lcone      mcone          rod        scone
## Red   0.18371781 0.02747917 0.0008675821 2.948406e-05
## Amber 0.09375166 0.03681379 0.0023429799 6.087550e-07
## Green 0.16939207 0.17126117 0.0724423315 2.606042e-03
## Blue  0.55313846 0.76444587 0.9243471065 9.973639e-01

In Jan’s spreadsheet, the A-Matrix is normalized, so that all colums add up to one. As a consequence, calculations with photoreceptor contrasts should be possible.

# Calculate column sums
apply(a_matrix_cord, 2, sum)
## lcone mcone   rod scone 
##     1     1     1     1

Fundamentals from Manchester

s_receptors <- 
  fundamentals_manchester %>%
  get_normalized_spectrum_matrix() 

p_led <- 
  led_spectra %>%
  get_normalized_spectrum_matrix() %>%
  calculate_power_spectra(luminances = led_luminances, spectra = .)

a_matrix_manchester <- create_A_matrix(P_led = p_led, s_receptors = s_receptors)

print(a_matrix_manchester)
##           lcone      mcone         rod        scone
## Red   0.2484116 0.04439978 0.000867162 5.394000e-05
## Amber 0.1181298 0.05421146 0.002342095 1.634294e-05
## Green 0.1944796 0.21780200 0.072440387 4.355952e-03
## Blue  0.4389790 0.68358676 0.924350356 9.955738e-01

Compare resulting LED contrasts

LED contrasts are calculated by calculating the crossproduct of the photoreceptor contrasts and the inverse of the A-Matrix.

print(find_LED_contrasts)
## function (photoreceptor_contrasts, A_matrix) 
## {
##     stopifnot(is.vector(photoreceptor_contrasts), is.matrix(A_matrix), 
##         !is.null(names(photoreceptor_contrasts)), colnames(A_matrix) == 
##             names(photoreceptor_contrasts), all(photoreceptor_contrasts <= 
##             1))
##     contrast_matrix <- crossprod(photoreceptor_contrasts, solve(A_matrix))
##     contrasts <- as.vector(contrast_matrix)
##     names(contrasts) <- colnames(contrast_matrix)
##     return(contrasts)
## }

The main part of this function is: contrast_matrix <- crossprod(photoreceptor_contrasts, solve(A_matrix))

L-cone contrast of 10%

Comparison of the Michelson contrast in %. A negative signs means modulation in counterphase.

find_LED_contrasts(c(lcone = .1, mcone = 0, rod = 0, scone = 0), a_matrix_cord) * 100
##           Red         Amber         Green          Blue 
##  91.114205494 -74.366745795   1.394296684  -0.006291326
find_LED_contrasts(c(lcone = .1, mcone = 0, rod = 0, scone = 0), a_matrix_manchester) * 100
##           Red         Amber         Green          Blue 
##  68.194255904 -60.783901546   1.253276218  -0.008180431

M-cone contrast of 10%

Apparently, there is a considerable difference between calculations. Note that amber modulation is > 100% here, so that 10% M-cone contrast is above instrument gamut, but this should not matter here.

find_LED_contrasts(c(lcone = 0, mcone = 0.1, rod = 0, scone = 0), a_matrix_cord) * 100
##           Red         Amber         Green          Blue 
## -257.77986489  531.41057832  -14.68265279    0.04566088
find_LED_contrasts(c(lcone = 0, mcone = 0.1, rod = 0, scone = 0), a_matrix_manchester) * 100
##           Red         Amber         Green          Blue 
## -162.25620441  357.88015275  -10.23717817    0.04770711

Rod contrast of 10%

find_LED_contrasts(c(lcone = 0, mcone = 0, rod = 0.1, scone = 0), a_matrix_cord) * 100
##           Red         Amber         Green          Blue 
##   404.5641391 -1106.1275420   174.9476544    -0.4684106
find_LED_contrasts(c(lcone = 0, mcone = 0, rod = 0.1, scone = 0), a_matrix_manchester) * 100
##          Red        Amber        Green         Blue 
##  317.6062436 -952.8639350  174.8306824   -0.7665059

S-cone contrast of 10%

find_LED_contrasts(c(lcone = 0, mcone = 0, rod = 0, scone = 0.1), a_matrix_cord) * 100
##        Red      Amber      Green       Blue 
## -227.89848  659.08371 -151.65930   10.42904
find_LED_contrasts(c(lcone = 0, mcone = 0, rod = 0, scone = 0.1), a_matrix_manchester) * 100
##        Red      Amber      Green       Blue 
## -213.54430  665.76768 -155.84678   10.72698