Matching spectra for a-opic irradiance with STLAB

When designing stimuli with STLAB, you may need to find the settings that most closely match a spectrum you measured elsewhere, which requires some linear algebra. In this example we are aiming to match the spectral output of a NeurOptics PLR-3000 automated pupillometer, which administers light stimuli with 4 white LEDs. We measured the spectral output of the PLR-3000 with an OceanOptics STS-VIS spectrometer at the usual eye position and calibrated using our standard pipeline.

import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
sns.set_context('paper', font_scale=1.5)

from pyplr.CIE import get_CIES026

# Load data
plr3000 = pd.read_csv(
    '../data/PLR-3000_oo_calibrated_spectra.csv', index_col='uW')
plr3000.columns = plr3000.columns.astype('int')

# Plot spectra
specs = (plr3000.reset_index()
                .melt(id_vars='uW',
                      var_name='Wavelength',
                      value_name='w/m2/nm')
                .sort_values(['uW','Wavelength'])
                .reset_index(drop=True))
ax = sns.lineplot(
    data=specs,
    x='Wavelength',
    y='w/m2/nm',
    units='uW',
    hue='uW',
    estimator=None)
ax.set_ylabel('W/m$^2$/nm')
ax.set_title('PLR-3000 spectral measurements')

# Plot a-opic irradiances
sss = get_CIES026(asdf=True)
sss = sss.fillna(0)
plr_3k_ao = plr3000.dot(sss)
data = (plr_3k_ao.reset_index()
                 .melt(id_vars=['uW'],
                       var_name=['aopic'],
                       value_name='irradiance'))
f2, ax2 = plt.subplots()
sns.barplot(
    data=data, x='uW', y='irradiance', hue='aopic', ax=ax2)
ax2.set_ylabel('W/m$^2$')
ax2.set_title('PLR-3000 $a$-opic irradiances');

In order to find the right STLAB settings, we need to make a CalibrationContext from our calibrated OceanOptics data.

from pyplr.calibrate import CalibrationContext

cc = CalibrationContext(
    '../data/S2_corrected_oo_spectra.csv', binwidth=1)
_ = cc.plot_calibrated_spectra()

Now we can start searching for the STLAB settings that match the PLR-3000 settings for a-opic irradiance. We will focus on the 180 uW setting to keep things simple. For all possible combinations in 10 choose 5 (one LED for each photoreceptor class), we use linear algebra to work out the input fraction of the chosen LEDs that is required for matching the spectrum. Only those solutions where the input fractions are between 0 and 1 are valid, because eventually we will need to convert these to 12-bit integers for STLAB.

import itertools

import numpy as np

# The PLR-3000 setting we care about
uW_setting = 180

# Photoreceptors classes we are aiming to match
opsins = ['L','M','S','Mel', 'Rods']

# An LED for each photoreceptors
num_leds = 5

# List to store valid settings
keep = []

# Loop through all possible combinations in 10 choose 5
for choose in itertools.combinations(range(10), num_leds):

    # Get the irradiances for each LED at maximum
    settings_to_irradiances = cc.aopic.loc[[(led, 4095) for led in choose], opsins]

    # Take the inverse
    irradiances_to_settings = np.linalg.inv(settings_to_irradiances)

    # Calculate the required input fraction for the chosen LEDs
    plr_irradiances = plr_3k_ao.loc[uW_setting, opsins].T
    settings = plr_irradiances.dot(irradiances_to_settings)

    # Keep the settings where all values are greater than 0 and less then 1
    if all(settings < 1) and all(settings > 0):
        keep.append((uW_setting, choose, settings))

print('\n> ' + str(len(keep)) + ' settings found')
keep

> 18 settings found

[(180,
  (0, 1, 2, 6, 7),
  array([0.00517166, 0.32173031, 0.02026771, 0.19441607, 0.0164306 ])),
 (180,
  (0, 1, 2, 6, 8),
  array([0.31550785, 0.05958875, 0.10219134, 0.20210005, 0.06132067])),
 (180,
  (0, 1, 2, 6, 9),
  array([0.3583914 , 0.02262412, 0.11388467, 0.20229605, 0.07338902])),
 (180,
  (0, 1, 3, 6, 7),
  array([0.01684065, 0.31693341, 0.01842986, 0.19397591, 0.01666438])),
 (180,
  (0, 1, 3, 6, 8),
  array([0.40155686, 0.01397228, 0.09859531, 0.20033016, 0.06598839])),
 (180,
  (0, 2, 3, 6, 9),
  array([0.40011401, 0.06656192, 0.04605868, 0.20147622, 0.07599868])),
 (180,
  (0, 2, 4, 6, 7),
  array([0.41946373, 0.08731701, 0.04332805, 0.17244175, 0.02821919])),
 (180,
  (0, 2, 4, 6, 8),
  array([0.38995726, 0.11689858, 0.00506435, 0.200176  , 0.06646312])),
 (180,
  (0, 2, 4, 6, 9),
  array([0.38655782, 0.11954637, 0.00182771, 0.2016076 , 0.07561016])),
 (180,
  (0, 2, 5, 6, 7),
  array([0.40016009, 0.10775283, 0.05472756, 0.15665175, 0.03607427])),
 (180,
  (0, 2, 5, 6, 8),
  array([0.3873756 , 0.11958843, 0.00513457, 0.19941887, 0.06819887])),
 (180,
  (0, 2, 5, 6, 9),
  array([0.38560989, 0.12052573, 0.00182263, 0.20135211, 0.0763111 ])),
 (180,
  (0, 3, 4, 6, 7),
  array([0.44201214, 0.07566702, 0.04067575, 0.17197971, 0.0284574 ])),
 (180,
  (0, 3, 4, 6, 8),
  array([0.41989687, 0.1015923 , 0.00106964, 0.19986998, 0.06721641])),
 (180,
  (0, 3, 5, 6, 7),
  array([0.42903038, 0.09205737, 0.05065177, 0.15726557, 0.03577909])),
 (180,
  (0, 3, 5, 6, 8),
  array([0.41950423, 0.10207723, 0.00106514, 0.19971146, 0.06758008])),
 (180,
  (0, 3, 6, 7, 9),
  array([0.41847277, 0.10241046, 0.19974393, 0.00152053, 0.07239993])),
 (180,
  (0, 3, 6, 8, 9),
  array([0.419022  , 0.10233208, 0.200328  , 0.0458547 , 0.02430788]))]

So there are 18 solutions to the problem using 5 LEDs. Next, we need to convert the input fractions to 12-bit integers for STLAB.

# Lists to store settings and predicted spectra
settings = []
predicted = []

# Remove the index
keep = [val[1:] for val in keep if val[0]]

# Loop over settings
for k in keep:
    leds = k[0]

    # Convert to 12-bit integer and save LED settings
    intensities = np.round(k[1] * 4095).astype('int')
    spec = [0]*10
    for led, i in zip(leds, intensities):
        spec[led] = i

    # Get predicted spectrum
    pred = cc.predict_spd(spec)

    # Add to lists
    settings.append(spec)
    predicted.append(pred)

# Make DFs
settings = pd.DataFrame(settings)
predicted = pd.concat(predicted)
predicted.reset_index(inplace=True, drop=True)

# In theory it doesn't matter which one we use, but let's define
# the 'best' solution as the one with the least squared error
best = predicted.sub(
    plr3000.loc[uW_setting].to_numpy()).pow(2).sum(axis=1).idxmin()

optimal_predicted = predicted.loc[best]
optimal_settings = settings.loc[best]

Time to visualise the solutions.

fig, ax = plt.subplots(nrows=1, sharex=True, figsize=(12,4))
plr3000.loc[uW_setting].plot(
    label='PLR-3000: {} uW'.format(
        uW_setting), ax=ax, color='k', lw=2, linestyle='--')
ax.plot(plr3000.columns,
        optimal_predicted,
        label='Optimal STLAB: ' + str(optimal_settings.to_list()),
        linestyle='--')
ax.legend()
for idx, p in predicted.iterrows():
    p.plot(lw=.2)
ax.set_xlabel('Wavelength (nm)')
ax.set_ylabel('SPD (W/m$^2$/nm)')
fig.savefig('../img/PLR-3000-STLAB-stimuli.tiff', dpi=300, bbox_inches='tight')

All that remains is to make a video file for use with STLAB.

from pyplr.stlab import pulse_protocol, video_file_to_dict
pulse_protocol(pulse_spec=optimal_settings.to_list(),
               pulse_duration=1000,
               fname='PLR-3000-{}-mw'.format(uW_setting))

"PLR-3000-180-mw.dsf" saved in the current working directory.

vf = video_file_to_dict('PLR-3000-180-mw.dsf')
vf

{'header': {'version': 1,
  'model': 'VEGA10',
  'channels': 10,
  'spectracount': 4,
  'transitionsCount': 4,
  'fluxReference': 0,
  'repeats': 1},
 'metadata': {'creation_time': '2021-05-27 08:33:12.766758',
  'creator': 'jtm',
  'protocol': 'pulse',
  'pulse_spec': '[21, 1317, 83, 0, 0, 0, 796, 67, 0, 0]',
  'pulse_duration': '1000'},
 'spectra': [[21, 1317, 83, 0, 0, 0, 796, 67, 0, 0],
  [21, 1317, 83, 0, 0, 0, 796, 67, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]],
 'transitions': [{'spectrum': 0, 'power': 100, 'time': 0, 'flags': 0},
  {'spectrum': 1, 'power': 100, 'time': 1000, 'flags': 0},
  {'spectrum': 2, 'power': 100, 'time': 1000, 'flags': 0},
  {'spectrum': 3, 'power': 100, 'time': 1100, 'flags': 0}]}