% *************************************************************************
% varsa_2023.m
% *************************************************************************
%
% Supplemental data accompanying "Rendering the Bluish Appearance of Snow:
% When Light Transmission Matters" by Varsa and Baranoski 2023. Published
% in IEEE Computer Graphics and Applications. Vol. XXX, No. YYY.
%
% This dataset provides simulated spectral responses resulting in the
% process of light transmission through snow. For detailed information
% regarding the simulation methodology, please see "SPLITSnow: A spectral
% light transport model for snow" by Varsa, Baranoski and Kimmel published
% in Remote Sensing of Environment, Vol. 255, pp. 112272:1-20, March, 2021.
% For more information regarding how this data may be employed in rendering
% applications please see the accompanying article mentioned above.
% Inline doccumentation has been provided to explain the data format and
% layout.
%
% Varsa et al. 2021:        https://doi.org/10.1016/j.rse.2020.112272
% Varsa and Baranoski 2023: https://doi.org/10.1109/mcg.2023.3307517
%
% Prepared by Petri Varsa, September 2023.
clear;

waves = getLambda();
M = getSpectra();

% Define of the snow sample propeties
GRAIN_SIZE = 400;
FACETNESS = 0.3;
DENSITY = 300;
INCIDENT_ANGLE = 20;
THICKNESS = '15';

% Compute the Index of the spectrum
i = getIndex(GRAIN_SIZE, DENSITY, FACETNESS, INCIDENT_ANGLE, THICKNESS);

% Retrieve the spectrum
Mi = M(i,:);

% Given a sample characterization, compute the corresponding index for that
% sample in the matrix returned by calling getSpectra().
%
%   grainsize - The grain size of the sample.
%   density   - The density of the sample.
%   facetness - The facetness of the sample.
%   angle     - The angle of incidence at the surface.
%   depth     - The sample depth.
%
function index = getIndex(grain_size, density, facetness, incident_angle, thickness)
    c = getSampleConstants();

    % Determine the individual indexes for each paramter from the list of
    % possible values for each sample characterization parameter
    f_index = find(c.facetnesses == facetness) - 1;
    gs_index = find(c.grain_sizes == grain_size) - 1;
    den_index = find(c.densities == density) - 1;
    a_index = find(c.incident_angles == incident_angle) - 1;
    dep_index = find(strcmp(c.thicknesses, thickness), 1) - 1;

    % The following line is a five-dimensional array index computation
    % One must be added to the result since Matlab uses one-based arrays
    index = f_index * c.f_step + gs_index * c.gs_step + ...
        den_index * c.den_step + a_index * c.a_step + dep_index + 1;
end

% Return a structure that defines the various sample characterization
% constants.
function constants = getSampleConstants()

    % Valid parameter values
    constants.facetnesses = [ 0.1 0.3 0.5 0.7 0.9 ];
    constants.grain_sizes = [ 100 150 200 250 300 350 400 450 500 ];
    constants.densities = [ 200 250 300 350 400 ];
    constants.incident_angles = [ 0 10 20 30 40 50 60 70 80 89 ];
    constants.thicknesses = { '1' '1.5' '2' '3' '5' '7' '10' '15' '25' '35' '50' };
    
    % Size of each valid parameter array
    constants.f_size = size(constants.facetnesses,2);
    constants.gs_size = size(constants.grain_sizes,2);
    constants.den_size = size(constants.densities,2);
    constants.a_size = size(constants.incident_angles,2);
    constants.dep_size = size(constants.thicknesses,2);

    % Compute some helpful values, such as total size and step sizes
    constants.num_rows = constants.f_size * constants.gs_size * constants.den_size ...
        * constants.a_size * constants.dep_size;
    constants.f_step = constants.gs_size * constants.den_size * constants.a_size * constants.dep_size;
    constants.gs_step = constants.den_size * constants.a_size * constants.dep_size;
    constants.den_step = constants.a_size * constants.dep_size;
    constants.a_step = constants.dep_size;
    constants.dep_step = 1;
end


% Return the wavelengths of light associated with each column of the matrix
% returned by the getSpectra() function.
function L = getLambda()
	L = 400:10:700;
end

% Return the raw data of simulation.
% A single matrix containing 24750 rows and 31 columns is returned. Each row
% represents a single spectrum of 31 wavelengths. The wavelengths are provided
% by the getLambda() function.
function M = getSpectra()
    FILE='varsa_2023.csv';
    M = readmatrix(FILE);
end
