1use crate::color::Matrix3;
2use crate::error::{LuxError, LuxResult};
3use crate::spectrum::Spectrum;
4
5const ASANO_LMS_ABSORBANCE: &str = include_str!("../data/indvcmf/asano_cie2006_Alms.dat");
6const ASANO_RELATIVE_MACULAR_DENSITY: &str =
7 include_str!("../data/indvcmf/asano_cie2006_RelativeMacularDensity.dat");
8const ASANO_OCULAR_DENSITY: &str = include_str!("../data/indvcmf/asano_cie2006_docul.dat");
9const WAVELENGTH_START: usize = 390;
10const WAVELENGTH_END: usize = 780;
11const WAVELENGTH_STEP: usize = 5;
12const FIELD_SIZE_MIN: f64 = 2.0;
13const FIELD_SIZE_MAX: f64 = 10.0;
14const S_CONE_CUTOFF: f64 = 620.0;
15const LMS_TO_XYZ_2_DEG: Matrix3 = [
16 [0.4151, -0.2424, 0.0425],
17 [0.1355, 0.0833, -0.0043],
18 [-0.0093, 0.0125, 0.2136],
19];
20const LMS_TO_XYZ_10_DEG: Matrix3 = [
21 [0.4499, -0.2630, 0.0460],
22 [0.1617, 0.0726, -0.0011],
23 [-0.0036, 0.0054, 0.2291],
24];
25
26#[derive(Debug, Clone, Copy, PartialEq)]
27pub struct IndividualObserverParameters {
28 pub age: f64,
29 pub field_size: f64,
30 pub lens_density_variation: f64,
31 pub macular_density_variation: f64,
32 pub cone_density_variation: [f64; 3],
33 pub cone_peak_shift: [f64; 3],
34 pub allow_negative_xyz_values: bool,
35}
36
37impl Default for IndividualObserverParameters {
38 fn default() -> Self {
39 Self {
40 age: 32.0,
41 field_size: 10.0,
42 lens_density_variation: 0.0,
43 macular_density_variation: 0.0,
44 cone_density_variation: [0.0, 0.0, 0.0],
45 cone_peak_shift: [0.0, 0.0, 0.0],
46 allow_negative_xyz_values: false,
47 }
48 }
49}
50
51#[derive(Debug, Clone, Copy, PartialEq)]
52pub struct IndividualObserverStdDevs {
53 pub lens_density: f64,
54 pub macular_density: f64,
55 pub cone_density: [f64; 3],
56 pub cone_peak_shift: [f64; 3],
57}
58
59#[derive(Debug, Clone, PartialEq)]
60pub struct IndividualObserverCmf {
61 pub lms: Spectrum,
62 pub xyz: Spectrum,
63 pub lens_transmission: Spectrum,
64 pub macular_transmission: Spectrum,
65 pub photopigment_sensitivity: Spectrum,
66 pub lms_to_xyz_matrix: Matrix3,
67}
68
69pub fn individual_observer_default_std_devs() -> IndividualObserverStdDevs {
70 IndividualObserverStdDevs {
71 lens_density: 19.1,
72 macular_density: 37.2,
73 cone_density: [17.9, 17.9, 14.7],
74 cone_peak_shift: [4.0, 3.0, 2.5],
75 }
76}
77
78pub fn individual_observer_lms_to_xyz_matrix(field_size: f64) -> Matrix3 {
79 let clamped = field_size.clamp(FIELD_SIZE_MIN, FIELD_SIZE_MAX);
80 let a = (FIELD_SIZE_MAX - clamped) / (FIELD_SIZE_MAX - FIELD_SIZE_MIN);
81 interpolate_matrix3(LMS_TO_XYZ_2_DEG, LMS_TO_XYZ_10_DEG, 1.0 - a)
82}
83
84pub fn individual_observer_lms_to_xyz(
85 lms: &Spectrum,
86 field_size: f64,
87 allow_negative_values: bool,
88) -> LuxResult<Spectrum> {
89 if lms.spectrum_count() != 3 {
90 return Err(LuxError::InvalidInput(
91 "individual observer LMS input must contain exactly 3 spectra",
92 ));
93 }
94
95 let matrix = individual_observer_lms_to_xyz_matrix(field_size);
96 let wavelengths = lms.wavelengths().to_vec();
97 let mut xyz = (0..3)
98 .map(|_| Vec::with_capacity(wavelengths.len()))
99 .collect::<Vec<_>>();
100
101 for index in 0..wavelengths.len() {
102 let lms_sample = [
103 lms.spectra()[0][index],
104 lms.spectra()[1][index],
105 lms.spectra()[2][index],
106 ];
107 let mut xyz_sample = multiply_matrix3_vector3(matrix, lms_sample);
108 if !allow_negative_values {
109 for value in &mut xyz_sample {
110 if *value < 0.0 {
111 *value = 0.0;
112 }
113 }
114 }
115 for axis in 0..3 {
116 xyz[axis].push(xyz_sample[axis]);
117 }
118 }
119
120 Spectrum::new(wavelengths, xyz)
121}
122
123pub fn individual_observer_cmf(
124 parameters: IndividualObserverParameters,
125) -> LuxResult<IndividualObserverCmf> {
126 validate_parameters(parameters)?;
127
128 if parameters.cone_peak_shift.iter().any(|shift| *shift != 0.0) {
129 return Err(LuxError::InvalidInput(
130 "non-zero cone peak shifts are not yet supported in the Rust indvcmf slice",
131 ));
132 }
133
134 let wavelengths = base_wavelengths();
135 let lms_absorbance = parse_columns(ASANO_LMS_ABSORBANCE, 3)?;
136 let relative_macular_density = parse_columns(ASANO_RELATIVE_MACULAR_DENSITY, 1)?
137 .into_iter()
138 .next()
139 .ok_or(LuxError::ParseError("missing macular density data"))?;
140 let ocular_density = parse_columns(ASANO_OCULAR_DENSITY, 2)?;
141
142 ensure_len(&wavelengths, &relative_macular_density)?;
143 ensure_len(&wavelengths, &ocular_density[0])?;
144 ensure_len(&wavelengths, &ocular_density[1])?;
145 ensure_len(&wavelengths, &lms_absorbance[0])?;
146 ensure_len(&wavelengths, &lms_absorbance[1])?;
147 ensure_len(&wavelengths, &lms_absorbance[2])?;
148
149 let fs = parameters.field_size;
150 let peak_macular_density =
151 0.485 * (-fs / 6.132).exp() * (1.0 + parameters.macular_density_variation / 100.0);
152 let corrected_macular_density: Vec<f64> = relative_macular_density
153 .iter()
154 .map(|value| value * peak_macular_density)
155 .collect::<Vec<_>>();
156
157 let age_scale = if parameters.age <= 60.0 {
158 1.0 + 0.02 * (parameters.age - 32.0)
159 } else {
160 1.56 + 0.0667 * (parameters.age - 60.0)
161 };
162 let corrected_ocular_density: Vec<f64> = ocular_density[0]
163 .iter()
164 .zip(ocular_density[1].iter())
165 .map(|(first, second)| {
166 (first * age_scale + second) * (1.0 + parameters.lens_density_variation / 100.0)
167 })
168 .collect::<Vec<_>>();
169
170 let cone_peak_density = [
171 (0.38 + 0.54 * (-fs / 1.333).exp()) * (1.0 + parameters.cone_density_variation[0] / 100.0),
172 (0.38 + 0.54 * (-fs / 1.333).exp()) * (1.0 + parameters.cone_density_variation[1] / 100.0),
173 (0.30 + 0.45 * (-fs / 1.333).exp()) * (1.0 + parameters.cone_density_variation[2] / 100.0),
174 ];
175
176 let mut alpha_lms = vec![vec![0.0; wavelengths.len()]; 3];
177 for axis in 0..3 {
178 for (index, wavelength) in wavelengths.iter().enumerate() {
179 alpha_lms[axis][index] = 1.0
180 - 10f64.powf(-cone_peak_density[axis] * 10f64.powf(lms_absorbance[axis][index]));
181 if axis == 2 && *wavelength >= S_CONE_CUTOFF {
182 alpha_lms[axis][index] = 0.0;
183 }
184 }
185 }
186
187 let mut lms = vec![vec![0.0; wavelengths.len()]; 3];
188 let mut photopigment_sensitivity = vec![vec![0.0; wavelengths.len()]; 3];
189 for axis in 0..3 {
190 for (index, wavelength) in wavelengths.iter().enumerate() {
191 let lms_quantal = alpha_lms[axis][index]
192 * 10f64.powf(-corrected_macular_density[index] - corrected_ocular_density[index]);
193 lms[axis][index] = lms_quantal * wavelength;
194 photopigment_sensitivity[axis][index] = alpha_lms[axis][index] * wavelength;
195 }
196 let area: f64 = lms[axis].iter().sum();
197 if area == 0.0 {
198 return Err(LuxError::InvalidInput(
199 "individual observer LMS normalization area must be non-zero",
200 ));
201 }
202 for value in &mut lms[axis] {
203 *value = 100.0 * *value / area;
204 }
205 }
206
207 let lms_matrix = Spectrum::new(wavelengths.clone(), lms)?;
208 let xyz_matrix = individual_observer_lms_to_xyz(
209 &lms_matrix,
210 parameters.field_size,
211 parameters.allow_negative_xyz_values,
212 )?;
213 let lens_transmission = Spectrum::new(
214 wavelengths.clone(),
215 corrected_ocular_density
216 .iter()
217 .map(|value| 10f64.powf(-value))
218 .collect::<Vec<_>>(),
219 )?;
220 let macular_transmission = Spectrum::new(
221 wavelengths.clone(),
222 corrected_macular_density
223 .iter()
224 .map(|value| 10f64.powf(-value))
225 .collect::<Vec<_>>(),
226 )?;
227 let photopigment_sensitivity = Spectrum::new(wavelengths, photopigment_sensitivity)?;
228
229 Ok(IndividualObserverCmf {
230 lms: lms_matrix,
231 xyz: xyz_matrix,
232 lens_transmission,
233 macular_transmission,
234 photopigment_sensitivity,
235 lms_to_xyz_matrix: individual_observer_lms_to_xyz_matrix(parameters.field_size),
236 })
237}
238
239fn validate_parameters(parameters: IndividualObserverParameters) -> LuxResult<()> {
240 if !parameters.age.is_finite() || parameters.age <= 0.0 {
241 return Err(LuxError::InvalidInput("age must be positive and finite"));
242 }
243 if !parameters.field_size.is_finite()
244 || parameters.field_size < FIELD_SIZE_MIN
245 || parameters.field_size > FIELD_SIZE_MAX
246 {
247 return Err(LuxError::InvalidInput(
248 "field size must be finite and within 2..=10 degrees",
249 ));
250 }
251 if !parameters.lens_density_variation.is_finite()
252 || !parameters.macular_density_variation.is_finite()
253 || parameters.lens_density_variation <= -100.0
254 || parameters.macular_density_variation <= -100.0
255 {
256 return Err(LuxError::InvalidInput(
257 "lens and macular density variations must be finite and greater than -100%",
258 ));
259 }
260 if parameters
261 .cone_density_variation
262 .iter()
263 .any(|value| !value.is_finite() || *value <= -100.0)
264 {
265 return Err(LuxError::InvalidInput(
266 "cone density variations must be finite and greater than -100%",
267 ));
268 }
269 if parameters
270 .cone_peak_shift
271 .iter()
272 .any(|value| !value.is_finite())
273 {
274 return Err(LuxError::InvalidInput(
275 "cone peak shifts must be finite when provided",
276 ));
277 }
278 Ok(())
279}
280
281fn base_wavelengths() -> Vec<f64> {
282 (WAVELENGTH_START..=WAVELENGTH_END)
283 .step_by(WAVELENGTH_STEP)
284 .map(|value| value as f64)
285 .collect::<Vec<_>>()
286}
287
288fn parse_columns(data: &str, expected_columns: usize) -> LuxResult<Vec<Vec<f64>>> {
289 let mut columns = vec![Vec::new(); expected_columns];
290
291 for line in data.split(['\n', '\r']) {
292 let trimmed = line.trim();
293 if trimmed.is_empty() {
294 continue;
295 }
296 let values: Vec<f64> = trimmed
297 .split(|char: char| char == ',' || char.is_ascii_whitespace())
298 .filter(|part| !part.is_empty())
299 .map(|part| {
300 part.parse::<f64>()
301 .map_err(|_| LuxError::ParseError("invalid indvcmf numeric value"))
302 })
303 .collect::<LuxResult<Vec<f64>>>()?;
304
305 if values.len() > expected_columns {
306 return Err(LuxError::ParseError("unexpected indvcmf column count"));
307 }
308 let mut padded_values = values;
309 while padded_values.len() < expected_columns {
310 padded_values.push(0.0);
311 }
312 for (column, value) in columns.iter_mut().zip(padded_values.into_iter()) {
313 column.push(value);
314 }
315 }
316
317 Ok(columns)
318}
319
320fn ensure_len(wavelengths: &[f64], values: &[f64]) -> LuxResult<()> {
321 if wavelengths.len() != values.len() {
322 return Err(LuxError::MismatchedLengths {
323 wavelengths: wavelengths.len(),
324 values: values.len(),
325 });
326 }
327 Ok(())
328}
329
330fn interpolate_matrix3(lhs: Matrix3, rhs: Matrix3, lhs_weight: f64) -> Matrix3 {
331 let rhs_weight = 1.0 - lhs_weight;
332 let mut matrix = [[0.0; 3]; 3];
333 for row in 0..3 {
334 for col in 0..3 {
335 matrix[row][col] = lhs[row][col] * lhs_weight + rhs[row][col] * rhs_weight;
336 }
337 }
338 matrix
339}
340
341fn multiply_matrix3_vector3(matrix: Matrix3, vector: [f64; 3]) -> [f64; 3] {
342 [
343 matrix[0][0] * vector[0] + matrix[0][1] * vector[1] + matrix[0][2] * vector[2],
344 matrix[1][0] * vector[0] + matrix[1][1] * vector[1] + matrix[1][2] * vector[2],
345 matrix[2][0] * vector[0] + matrix[2][1] * vector[1] + matrix[2][2] * vector[2],
346 ]
347}