// 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 <cmath>

/*! This utility function evaluates the non-linearity for
	convertXYZToCIELab().*/
inline float nonLinearityXYZToCIELAB(float t) {
	const float delta = 6.0f / 29.0f;
	return (t > delta*delta*delta) ? std::pow(t, 1.0f / 3.0f) : (t / (3.0f*delta*delta) + 4.0f / 29.0f);
}


/*! Maps an XYZ triple to a triple in the CIELAB color space, assuming a D65
	whitepoint. For triples with negative luminance a vector of -1.0f is
	returned.*/
inline std::array<float, 3> convertXYZToCIELAB(const std::array<float, 3>& rXYZ) {
	if (rXYZ[1] < 0.0f) {
		return { -1.0f, -1.0f, -1.0f };
	}
	float nonLinearX = nonLinearityXYZToCIELAB(rXYZ[0] * (1.0f / 0.95047f));
	float nonLinearY = nonLinearityXYZToCIELAB(rXYZ[1]);
	float nonLinearZ = nonLinearityXYZToCIELAB(rXYZ[2] * (1.0f / 1.08883f));
	return {
		116.0f * nonLinearY - 16.0f,
		500.0f * (nonLinearX - nonLinearY),
		200.0f * (nonLinearY - nonLinearZ)
	};
}


/*! This utility function evaluates the non-linearity for
	convertSRGBToXYZ().*/
inline float nonLinearitySRGBToXYZ(float t) {
	return (t > 0.04045f) ? std::pow((t + 0.055f) / 1.055f, 2.4f) : (t / 12.92f);
}


/*! This function converts an sRGB triple to an XYZ triple. The normalization
	is done such that the input triple (1,1,1) corresponds to the XYZ for a
	reflectance spectrum that is constant one, lit by illuminant E or 
	illuminant D65.*/
inline std::array<float, 3> convertSRGBToXYZ(const std::array<float, 3>& rSRGB, bool useD65) {
	std::array<float, 3> linearSRGB = {
		nonLinearitySRGBToXYZ(rSRGB[0]),
		nonLinearitySRGBToXYZ(rSRGB[1]),
		nonLinearitySRGBToXYZ(rSRGB[2])
	};
	if(useD65)
	{
		std::array<float, 3> xyz = {
			0.4124f*linearSRGB[0] + 0.3576f*linearSRGB[1] + 0.1805f*linearSRGB[2],
			0.2126f*linearSRGB[0] + 0.7152f*linearSRGB[1] + 0.0722f*linearSRGB[2],
			0.0193f*linearSRGB[0] + 0.1192f*linearSRGB[1] + 0.9505f*linearSRGB[2]
		};
		const float denormalization = 10567.2949f;
		return { denormalization * xyz[0], denormalization * xyz[1], denormalization * xyz[2] };
	}
	else
	{
		std::array<float, 3> xyz = {
			0.49685929f*linearSRGB[0] + 0.33909388f*linearSRGB[1] + 0.16404638f*linearSRGB[2],
			0.25619309f*linearSRGB[0] + 0.67818808f*linearSRGB[1] + 0.06561853f*linearSRGB[2],
			0.02329024f*linearSRGB[0] + 0.11303116f*linearSRGB[1] + 0.86397868f*linearSRGB[2]
		};
		const float denormalization = 106.85706329345703f;
		return { denormalization * xyz[0], denormalization * xyz[1], denormalization * xyz[2] };
	}
}
