// Copyright (c) 2019, Christoph Peters
// All rights reserved.
// 
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//     * Redistributions of source code must retain the above copyright
//       notice, this list of conditions and the following disclaimer.
//     * Redistributions in binary form must reproduce the above copyright
//       notice, this list of conditions and the following disclaimer in the
//       documentation and/or other materials provided with the distribution.
//     * Neither the name of the Karlsruhe Institute of Technology nor the
//       names of its contributors may be used to endorse or promote products
//       derived from this software without specific prior written permission.
// 
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY
// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


#pragma once
#include <array>
#include <vector>
#include <cmath>
#include "FourierSpectrum.h"

/*! This class implements a quadrature rule to compute tristimulus values
	representing a given continuously defined reflectance spectrum.*/
class TristimulusIntegrator {
public:
	/*! A vector defining a quadrature rule. Each entry consists of a
		wavelength in nanometers and three weights (which may be negative).
		The spectrum is evaluated at each of the wavelengths, multiplied by
		each of the weights and the three results are added to three totals.*/
	typedef std::vector<std::array<float, 4>> QuadratureRuleType;

	/*! This constructor defines a quadrature rule that uses the midpoint rule
		to integrate the product of the reflectance spectrum and three
		piecewise constant functions.
		\param wavelengthMin The wavelength in nanometers for the first sample
				of the piecewise-constant functions in valueVector.
		\param wavelengthDelta The distance in nanometers between two
				subsequent samples. The maximal wavelength then follows from
				the size of valueVector.
		\param valueVector Samples of the three piecewise constant functions.
				On the interval from wavelengthMin+(i-0.5f)*wavelengthDelta to
				wavelengthMin+(i+0.5f)*wavelengthDelta function j takes value
				valueVector[i][j].
		\param nSamplePerInterval The number of equidistant samples that will
				be taken from the reflectance spectrum for each interval where
				the piecewise constant functions are constant.*/
	TristimulusIntegrator(float wavelengthMin, float wavelengthDelta, std::vector<std::array<float, 3>> valueVector, std::size_t nSamplePerInterval = 1) :
		m_quadratureRule(valueVector.size()*nSamplePerInterval)
	{
		float sampleDelta = wavelengthDelta / nSamplePerInterval;
		float wavelengthSample0 = wavelengthMin - 0.5f*wavelengthDelta + 0.5f*sampleDelta;
		for (std::size_t i = 0; i != valueVector.size(); ++i) {
			for (std::size_t j = 0; j != nSamplePerInterval; ++j) {
				float wavelength = wavelengthSample0 + wavelengthDelta * i + sampleDelta * j;
				m_quadratureRule[i*nSamplePerInterval + j][0] = wavelength;
				for (std::size_t k = 0; k != 3; ++k) {
					m_quadratureRule[nSamplePerInterval * i + j][k + 1] = valueVector[i][k] * sampleDelta;
				}
			}
		}
	}

	/*! This constructor enables definition of arbitrary quadrature rules.*/
	TristimulusIntegrator(const QuadratureRuleType& rQuadratureRule) :
		m_quadratureRule(rQuadratureRule) {}

	/*! Applies this quadrature rule to the given Fourier spectrum.*/
	template<bool mirrorSignal, class WavelengthToPhaseType>
	inline std::array<float, 3> integrate(const FourierSpectrum<mirrorSignal, WavelengthToPhaseType>& rFourierSpectrum) const {
		std::array<float, 3> result = { -1.0f, -1.0f, -1.0f };
		if (rFourierSpectrum.isValid()) {
			bool valid = true;
			result = { 0.0f, 0.0f, 0.0f };
			for (std::array<float, 4> sample : m_quadratureRule) {
				float reflectance = rFourierSpectrum(sample[0]);
				for (std::size_t i = 0; i != 3; ++i) {
					result[i] += sample[i + 1] * reflectance;
				}
				// FourierSpectrum::isValid() is slightly unstable for points close to the 
				// boundary. Thus, we do a second check for invalid results.
				valid &= ((reflectance>=0.0f) && (reflectance<=1.0f));
			}
			if(!valid || !std::isfinite(result[0]) || !std::isfinite(result[1]) || !std::isfinite(result[2])){
				return { -1.0f, -1.0f, -1.0f };
			}
		}
		return result;
	}

	/*! \return The integrated values for a spectrum that is constant one.*/
	inline std::array<float, 3> integrateConstantOne() const {
		std::array<float, 3> result = { 0.0f,0.0f,0.0f };
		for (std::array<float, 4> sample : m_quadratureRule) {
			for (std::size_t i = 0; i != 3; ++i) {
				result[i] += sample[i + 1];
			}
		}
		return result;
	}

private:
	/*! The used quadrature rule.*/
	QuadratureRuleType m_quadratureRule;
};


/*! Returns an integrator that computes XYZ values of a reflectance spectrum
	lit by illuminant D65 or illuminant E.
	\sa TristimulusIntegrator(float,float,std::vector<std::array<float,3>>,int)*/
TristimulusIntegrator getXYZIntegrator(std::size_t nSamplePerInterval = 5, bool useD65 = false) {
	static const std::vector<float> X = { 0.0001299f, 0.0002321f, 0.0004149f, 0.0007416f, 0.001368f, 0.002236f, 0.004243f, 0.00765f, 0.01431f, 0.02319f, 0.04351f, 0.07763f, 0.13438f, 0.21477f, 0.2839f, 0.3285f, 0.34828f, 0.34806f, 0.3362f, 0.3187f, 0.2908f, 0.2511f, 0.19536f, 0.1421f, 0.09564f, 0.05795001f, 0.03201f, 0.0147f, 0.0049f, 0.0024f, 0.0093f, 0.0291f, 0.06327f, 0.1096f, 0.1655f, 0.2257499f, 0.2904f, 0.3597f, 0.4334499f, 0.5120501f, 0.5945f, 0.6784f, 0.7621f, 0.8425f, 0.9163f, 0.9786f, 1.0263f, 1.0567f, 1.0622f, 1.0456f, 1.0026f, 0.9384f, 0.8544499f, 0.7514f, 0.6424f, 0.5419f, 0.4479f, 0.3608f, 0.2835f, 0.2187f, 0.1649f, 0.1212f, 0.0874f, 0.0636f, 0.04677f, 0.0329f, 0.0227f, 0.01584f, 0.01135916f, 0.008110916f, 0.005790346f, 0.004109457f, 0.002899327f, 0.00204919f, 0.001439971f, 0.0009999493f, 0.0006900786f, 0.0004760213f, 0.0003323011f, 0.0002348261f, 0.0001661505f, 0.000117413f, 0.00008307527f, 0.00005870652f, 0.00004150994f, 0.00002935326f, 0.00002067383f, 0.00001455977f, 0.00001025398f, 0.000007221456f, 0.000005085868f, 0.000003581652f, 0.000002522525f, 0.000001776509f, 0.000001251141f };
	static const std::vector<float> Y = { 0.000003917f, 0.000006965f, 0.00001239f, 0.00002202f, 0.000039f, 0.000064f, 0.00012f, 0.000217f, 0.000396f, 0.00064f, 0.00121f, 0.00218f, 0.004f, 0.0073f, 0.0116f, 0.01684f, 0.023f, 0.0298f, 0.038f, 0.048f, 0.06f, 0.0739f, 0.09098f, 0.1126f, 0.13902f, 0.1693f, 0.20802f, 0.2586f, 0.323f, 0.4073f, 0.503f, 0.6082f, 0.71f, 0.7932f, 0.862f, 0.9148501f, 0.954f, 0.9803f, 0.9949501f, 1.f, 0.995f, 0.9786f, 0.952f, 0.9154f, 0.87f, 0.8163f, 0.757f, 0.6949f, 0.631f, 0.5668f, 0.503f, 0.4412f, 0.381f, 0.321f, 0.265f, 0.217f, 0.175f, 0.1382f, 0.107f, 0.0816f, 0.061f, 0.04458f, 0.032f, 0.0232f, 0.017f, 0.01192f, 0.00821f, 0.005723f, 0.004102f, 0.002929f, 0.002091f, 0.001484f, 0.001047f, 0.00074f, 0.00052f, 0.0003611f, 0.0002492f, 0.0001719f, 0.00012f, 0.0000848f, 0.00006f, 0.0000424f, 0.00003f, 0.0000212f, 0.00001499f, 0.0000106f, 0.0000074657f, 0.0000052578f, 0.0000037029f, 0.0000026078f, 0.0000018366f, 0.0000012934f, 0.00000091093f, 0.00000064153f, 0.000000451810f };
	static const std::vector<float> Z = { 0.0006061f, 0.001086f, 0.001946f, 0.003486f, 0.006450001f, 0.01054999f, 0.02005001f, 0.03621f, 0.06785001f, 0.1102f, 0.2074f, 0.3713f, 0.6456f, 1.0390501f, 1.3856f, 1.62296f, 1.74706f, 1.7826f, 1.77211f, 1.7441f, 1.6692f, 1.5281f, 1.28764f, 1.0419f, 0.8129501f, 0.6162f, 0.46518f, 0.3533f, 0.272f, 0.2123f, 0.1582f, 0.1117f, 0.07824999f, 0.05725001f, 0.04216f, 0.02984f, 0.0203f, 0.0134f, 0.008749999f, 0.005749999f, 0.0039f, 0.002749999f, 0.0021f, 0.0018f, 0.001650001f, 0.0014f, 0.0011f, 0.001f, 0.0008f, 0.0006f, 0.00034f, 0.00024f, 0.00019f, 0.0001f, 0.00004999999f, 0.00003f, 0.00002f, 0.00001f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f };
	static const std::vector<float> D65 = { 46.6383f, 49.3637f, 52.0891f, 51.0323f, 49.9755f, 52.3118f, 54.6482f, 68.7015f, 82.7549f, 87.1204f, 91.486f, 92.4589f, 93.4318f, 90.057f, 86.6823f, 95.7736f, 104.865f, 110.936f, 117.008f, 117.41f, 117.812f, 116.336f, 114.861f, 115.392f, 115.923f, 112.367f, 108.811f, 109.082f, 109.354f, 108.578f, 107.802f, 106.296f, 104.79f, 106.239f, 107.689f, 106.047f, 104.405f, 104.225f, 104.046f, 102.023f, 100.0f, 98.1671f, 96.3342f, 96.0611f, 95.788f, 92.2368f, 88.6856f, 89.3459f, 90.0062f, 89.8026f, 89.5991f, 88.6489f, 87.6987f, 85.4936f, 83.2886f, 83.4939f, 83.6992f, 81.863f, 80.0268f, 80.1207f, 80.2146f, 81.2462f, 82.2778f, 80.281f, 78.2842f, 74.0027f, 69.7213f, 70.6652f, 71.6091f, 72.979f, 74.349f, 67.9765f, 61.604f, 65.7448f, 69.8856f, 72.4863f, 75.087f, 69.3398f, 63.5927f, 55.0054f, 46.4182f, 56.6118f, 66.8054f, 65.0941f, 63.3828f, 63.8434f, 64.304f, 61.8779f, 59.4519f, 55.7054f, 51.959f, 54.6998f, 57.4406f, 58.8765f, 60.312500 };
	std::vector<std::array<float, 3>> valueVector(X.size());
	for (std::size_t i = 0; i != X.size(); ++i) {
		valueVector[i][0] = X[i] * (useD65 ? D65[i] : 1.0f);
		valueVector[i][1] = Y[i] * (useD65 ? D65[i] : 1.0f);
		valueVector[i][2] = Z[i] * (useD65 ? D65[i] : 1.0f);
	}
	return TristimulusIntegrator(360.0f, 5.0f, valueVector, nSamplePerInterval);
}
