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