// 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 "momentBasedSpectra.real.hpp"

/*! \return true iff the Toeplitz matrix associated with the given moments is
			positive definite.*/
inline bool is2x2ToeplitzMatrixPositiveDefinite(const float_complex pMoment[2]) {
	// The approach is based on Levinson's algorithm to ensure that the result 
	// is equivalent to the result of the reconstruction code even in presence of 
	// rounding errors
	float_complex pOutSolution[1];
	if(pMoment[0].x <= 0.0f) {
		return false;
	}
	pOutSolution[0] = float_complex(1.0f / pMoment[0].x, 0.0f);
	float_complex dotProduct;
	float_complex flippedSolution1;
	float factor;
	dotProduct = pOutSolution[0].x * pMoment[1];
	factor = 1.0f / (1.0f - absSqr(dotProduct));
	if(factor <= 0.0f) {
		return false;
	}
	return true;
}

/*! \return true iff the Toeplitz matrix associated with the given moments is
			positive definite.*/
inline bool is3x3ToeplitzMatrixPositiveDefinite(const float_complex pMoment[3]) {
	// The approach is based on Levinson's algorithm to ensure that the result 
	// is equivalent to the result of the reconstruction code even in presence of 
	// rounding errors
	float_complex pOutSolution[2];
	if(pMoment[0].x <= 0.0f) {
		return false;
	}
	pOutSolution[0] = float_complex(1.0f / pMoment[0].x, 0.0f);
	float_complex dotProduct;
	float_complex flippedSolution1;
	float_complex flippedSolution2;
	float factor;
	dotProduct = pOutSolution[0].x * pMoment[1];
	factor = 1.0f / (1.0f - absSqr(dotProduct));
	if(factor <= 0.0f) {
		return false;
	}
	flippedSolution1 = float_complex(pOutSolution[0].x, 0.0f);
	pOutSolution[0] = float_complex(factor*pOutSolution[0].x, 0.0f);
	pOutSolution[1] = factor * (-flippedSolution1.x * dotProduct);
	dotProduct = pOutSolution[0].x * pMoment[2] + pOutSolution[1] * pMoment[1];
	factor = 1.0f / (1.0f - absSqr(dotProduct));
	if(factor <= 0.0f) {
		return false;
	}
	return true;
}


/*! This class template implements evaluation of a reflectance spectrum
	described by three Fourier coefficients. The instance is specific to one
	set of Fourier coefficients. Once the instance is constructed, evaluation
	for many wavelengths will be fast.*/
template<bool mirrorSignal, class WavelengthToPhaseType>
class FourierSpectrum {
public:
	/*! Constructs an object that can be used to evaluate the reflectance
		spectrum for the given Fourier coefficients. If mirrorSignal is false,
		fourierCoefficients[2] provides the imaginary part of Fourier
		coefficient 1.*/
	inline FourierSpectrum(const WavelengthToPhaseType& wavelengthToPhase, std::array<float, 3> fourierCoefficients) :
		m_wavelengthToPhase(wavelengthToPhase)
	{
		if(mirrorSignal) {
			float_complex pExponentialMoment[3];
			float_complex pEvaluationPolynomial[3];
			prepareReflectanceSpectrumReal3(pExponentialMoment, pEvaluationPolynomial, fourierCoefficients.data());
			m_valid = is3x3ToeplitzMatrixPositiveDefinite(pExponentialMoment);
			if(m_valid) {
				float_complex pAutocorrelation[3];
				computeAutocorrelation3(pAutocorrelation, pEvaluationPolynomial);
				pExponentialMoment[0] = 0.5f * pExponentialMoment[0];
				computeImaginaryCorrelation3(m_pLagrangeMultiplier, pAutocorrelation, pExponentialMoment);
				float normalizationFactor = 1.0f / (3.14159265f * pEvaluationPolynomial[0].x);
				for(std::size_t i = 0; i != 3; ++i) {
					m_pLagrangeMultiplier[i] = normalizationFactor * m_pLagrangeMultiplier[i];
				}
			}
		}
		else {
			float_complex pComplexFourierCoefficient[2] = { float_complex(fourierCoefficients[0],0.0f),float_complex(fourierCoefficients[1],fourierCoefficients[2]) };
			float_complex pExponentialMoment[2];
			float_complex pEvaluationPolynomial[2];
			prepareReflectanceSpectrum2(pExponentialMoment, pEvaluationPolynomial, pComplexFourierCoefficient);
			m_valid = is2x2ToeplitzMatrixPositiveDefinite(pExponentialMoment);
			if(m_valid) {
				float_complex pAutocorrelation[2];
				computeAutocorrelation2(pAutocorrelation, pEvaluationPolynomial);
				pExponentialMoment[0] = 0.5f * pExponentialMoment[0];
				float_complex pLagrangeMultiplier[2];
				computeCorrelation2(pLagrangeMultiplier, pAutocorrelation, pExponentialMoment);
				float normalizationFactor = 1.0f / (3.14159265f*pEvaluationPolynomial[0].x);
				m_pLagrangeMultiplier[0] = normalizationFactor * pLagrangeMultiplier[0].y;
				m_pLagrangeMultiplier[1] = normalizationFactor * pLagrangeMultiplier[1].y;
				m_pLagrangeMultiplier[2] = -normalizationFactor * pLagrangeMultiplier[1].x;
			}
		}
	}

	/*! \return true iff the provided Fourier coefficients admit a matching
		reflectance spectrum.*/
	inline bool isValid() const {
		return m_valid;
	}

	/*! Given a wavelength in nanometers, this operator returns the
		corresponding reflectance between zero and one.*/
	inline float operator()(float wavelength) const {
		float phase;
		if(m_wavelengthToPhase(wavelength, phase)) {
			if(mirrorSignal) {
				return evaluateReflectanceSpectrumLagrangeReal3(phase, m_pLagrangeMultiplier);
			}
			else {
				float_complex pLagrangeMultiplier[2] = { float_complex(m_pLagrangeMultiplier[0],0.0f),float_complex(m_pLagrangeMultiplier[1],m_pLagrangeMultiplier[2]) };
				return evaluateReflectanceSpectrumLagrange2(phase, pLagrangeMultiplier);
			}
		}
		return 0.0f;
	}

private:
	/*! true iff the provided Fourier coefficients admit a matching reflectance
		spectrum.*/
	bool m_valid;
	/*! The Lagrange multipliers that describe the reconstruction, packed like
		the input Fourier coefficients.*/
	float m_pLagrangeMultiplier[3];
	/*! A bool(float,float&) functor mapping wavelength in nanometers to a
		boolean indicating whether the reflectance at this wavelength may be
		non-zero (i.e. it is not out of range) and a phase.*/
	WavelengthToPhaseType m_wavelengthToPhase;
};
