/*
    Copyright (C) 2014-2021 Carl Hetherington <cth@carlh.net>

    This file is part of libdcp.

    libdcp is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    libdcp is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with libdcp.  If not, see <http://www.gnu.org/licenses/>.

    In addition, as a special exception, the copyright holders give
    permission to link the code of portions of this program with the
    OpenSSL library under certain conditions as described in each
    individual source file, and distribute linked combinations
    including the two.

    You must obey the GNU General Public License in all respects
    for all of the code used other than OpenSSL.  If you modify
    file(s) with this exception, you may extend this exception to your
    version of the file(s), but you are not obligated to do so.  If you
    do not wish to do so, delete this exception statement from your
    version.  If you delete this exception statement from all source
    files in the program, then also delete it here.
*/


/** @file  src/colour_conversion.cc
 *  @brief ColourConversion class
 */


#include "colour_conversion.h"
#include "compose.hpp"
#include "dcp_assert.h"
#include "gamma_transfer_function.h"
#include "identity_transfer_function.h"
#include "modified_gamma_transfer_function.h"
#include "openjpeg_image.h"
#include "s_gamut3_transfer_function.h"
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>


using std::make_shared;
using std::max;
using std::min;
using std::shared_ptr;
using boost::optional;
using namespace dcp;


static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;


ColourConversion const &
ColourConversion::srgb_to_dcp()
{
	static auto c = new ColourConversion(
		make_shared<ModifiedGammaTransferFunction>(2.4, 0.04045, 0.055, 12.92),
		YUVToRGB::REC601,
		Chromaticity(0.64, 0.33),
		Chromaticity(0.3, 0.6),
		Chromaticity(0.15, 0.06),
		Chromaticity::D65(),
		optional<Chromaticity>(),
		make_shared<GammaTransferFunction>(2.6),
		make_shared<GammaTransferFunction>(2.4)
		);
	return *c;
}

ColourConversion const &
ColourConversion::rec601_to_dcp()
{
	static auto c = new ColourConversion(
		make_shared<GammaTransferFunction>(2.2),
		YUVToRGB::REC601,
		Chromaticity(0.64, 0.33),
		Chromaticity(0.3, 0.6),
		Chromaticity(0.15, 0.06),
		Chromaticity::D65(),
		optional<Chromaticity>(),
		make_shared<GammaTransferFunction>(2.6),
		make_shared<GammaTransferFunction>(2.4)
		);
	return *c;
}

ColourConversion const &
ColourConversion::rec709_to_dcp()
{
	static auto c = new ColourConversion(
		make_shared<GammaTransferFunction>(2.2),
		YUVToRGB::REC709,
		Chromaticity(0.64, 0.33),
		Chromaticity(0.3, 0.6),
		Chromaticity(0.15, 0.06),
		Chromaticity::D65(),
		optional<Chromaticity>(),
		make_shared<GammaTransferFunction>(2.6),
		make_shared<GammaTransferFunction>(2.4)
		);
	return *c;
}

ColourConversion const &
ColourConversion::p3_to_dcp()
{
	static auto c = new ColourConversion(
		make_shared<GammaTransferFunction>(2.6),
		YUVToRGB::REC709,
		Chromaticity(0.68, 0.32),
		Chromaticity(0.265, 0.69),
		Chromaticity(0.15, 0.06),
		Chromaticity(0.314, 0.351),
		optional<Chromaticity>(),
		make_shared<GammaTransferFunction>(2.6),
		make_shared<GammaTransferFunction>(2.4)
		);
	return *c;
}

ColourConversion const &
ColourConversion::rec1886_to_dcp()
{
	/* According to Olivier on DCP-o-matic bug #832, Rec. 1886 is Rec. 709 with
	   2.4 gamma, so here goes ...
	*/
	static auto c = new ColourConversion(
		make_shared<GammaTransferFunction>(2.4),
		YUVToRGB::REC709,
		Chromaticity(0.64, 0.33),
		Chromaticity(0.3, 0.6),
		Chromaticity(0.15, 0.06),
		Chromaticity::D65(),
		optional<Chromaticity>(),
		make_shared<GammaTransferFunction>(2.6),
		make_shared<GammaTransferFunction>(2.4)
		);
	return *c;
}

ColourConversion const &
ColourConversion::rec2020_to_dcp()
{
	static auto c = new ColourConversion(
		make_shared<GammaTransferFunction>(2.4),
		YUVToRGB::REC2020,
		Chromaticity(0.708, 0.292),
		Chromaticity(0.170, 0.797),
		Chromaticity(0.131, 0.046),
		Chromaticity::D65(),
		optional<Chromaticity>(),
		make_shared<GammaTransferFunction>(2.6),
		make_shared<GammaTransferFunction>(2.4)
		);
	return *c;
}

/** Sony S-Gamut3/S-Log3 */
ColourConversion const &
ColourConversion::s_gamut3_to_dcp()
{
	static auto c = new ColourConversion(
		make_shared<SGamut3TransferFunction>(),
		YUVToRGB::REC709,
		Chromaticity(0.73, 0.280),
		Chromaticity(0.140, 0.855),
		Chromaticity(0.100, -0.050),
		Chromaticity::D65(),
		optional<Chromaticity>(),
		make_shared<IdentityTransferFunction>(),
		make_shared<GammaTransferFunction>(2.4)
		);
	return *c;
}



ColourConversion::ColourConversion(
	shared_ptr<const TransferFunction> in,
	YUVToRGB yuv_to_rgb,
	Chromaticity red,
	Chromaticity green,
	Chromaticity blue,
	Chromaticity white,
	optional<Chromaticity> adjusted_white,
	shared_ptr<const TransferFunction> out_j2k,
	shared_ptr<const TransferFunction> out_mpeg2
	)
	: _in(in)
	, _yuv_to_rgb(yuv_to_rgb)
	, _red(red)
	, _green(green)
	, _blue(blue)
	, _white(white)
	, _adjusted_white(adjusted_white)
	, _out_j2k(out_j2k)
	, _out_mpeg2(out_mpeg2)
{

}


bool
ColourConversion::about_equal(ColourConversion const & other, float epsilon) const
{
	if (!_in->about_equal(other._in, epsilon) ||
	    _yuv_to_rgb != other._yuv_to_rgb ||
	    !_red.about_equal(other._red, epsilon) ||
	    !_green.about_equal(other._green, epsilon) ||
	    !_blue.about_equal(other._blue, epsilon) ||
	    !_white.about_equal(other._white, epsilon) ||
	    (!_out_j2k && other._out_j2k) ||
	    (_out_j2k && !other._out_j2k) ||
	    (_out_j2k && !_out_j2k->about_equal(other._out_j2k, epsilon)) ||
	    (!_out_mpeg2 && other._out_mpeg2) ||
	    (_out_mpeg2 && !other._out_mpeg2) ||
	    (_out_mpeg2 && !_out_mpeg2->about_equal(other._out_mpeg2, epsilon))) {
		return false;
	}

	if (!_adjusted_white && !other._adjusted_white) {
		return true;
	}

	if (
		_adjusted_white && other._adjusted_white &&
		fabs(_adjusted_white.get().x - other._adjusted_white.get().x) < epsilon &&
		fabs(_adjusted_white.get().y - other._adjusted_white.get().y) < epsilon
		) {
		return true;
	}

	/* Otherwise one has an adjusted white and other hasn't, or they both have but different */
	return false;
}


boost::numeric::ublas::matrix<double>
ColourConversion::rgb_to_xyz() const
{
	/* See doc/design/colour.tex */

	double const D = (_red.x - _white.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_red.y - _white.y);
	double const E = (_white.x - _green.x) * (_red.y - _white.y) - (_red.x - _white.x) * (_white.y - _green.y);
	double const F = (_white.x - _green.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_white.y - _green.y);
	double const P = _red.y + _green.y * D / F + _blue.y * E / F;

	boost::numeric::ublas::matrix<double> C(3, 3);
	C(0, 0) = _red.x / P;
	C(0, 1) = _green.x * D / (F * P);
	C(0, 2) = _blue.x * E / (F * P);
	C(1, 0) = _red.y / P;
	C(1, 1) = _green.y * D / (F * P);
	C(1, 2) = _blue.y * E / (F * P);
	C(2, 0) = _red.z() / P;
	C(2, 1) = _green.z() * D / (F * P);
	C(2, 2) = _blue.z() * E / (F * P);
	return C;
}


boost::numeric::ublas::matrix<double>
ColourConversion::xyz_to_rgb() const
{
	boost::numeric::ublas::matrix<double> A(rgb_to_xyz());

	/* permutation matrix for the LU-factorization */
	boost::numeric::ublas::permutation_matrix<std::size_t> pm(A.size1());

	/* perform LU-factorization */
	int const r = lu_factorize(A, pm);
	DCP_ASSERT(r == 0);

	/* create identity matrix of inverse */
	boost::numeric::ublas::matrix<double> xyz_to_rgb(3, 3);
	xyz_to_rgb.assign(boost::numeric::ublas::identity_matrix<double>(A.size1()));

	/* backsubstitute to get the inverse */
	lu_substitute(A, pm, xyz_to_rgb);

	return xyz_to_rgb;
}


boost::numeric::ublas::matrix<double>
ColourConversion::bradford() const
{
	if (!_adjusted_white || fabs(_adjusted_white.get().x) < 1e-6 || fabs(_adjusted_white.get().y) < 1e-6) {
		boost::numeric::ublas::matrix<double> B = boost::numeric::ublas::zero_matrix<double>(3, 3);
		B(0, 0) = 1;
		B(1, 1) = 1;
		B(2, 2) = 1;
		return B;
	}

	/* See doc/design/colour.tex */

	boost::numeric::ublas::matrix<double> M(3, 3);
	M(0, 0) = 0.8951;
	M(0, 1) = 0.2664;
	M(0, 2) = -0.1614;
	M(1, 0) = -0.7502;
	M(1, 1) = 1.7135;
	M(1, 2) = 0.0367;
	M(2, 0) = 0.0389;
	M(2, 1) = -0.0685;
	M(2, 2) = 1.0296;

	boost::numeric::ublas::matrix<double> Mi(3, 3);
	Mi(0, 0) = 0.9869929055;
	Mi(0, 1) = -0.1470542564;
	Mi(0, 2) = 0.1599626517;
	Mi(1, 0) = 0.4323052697;
	Mi(1, 1) = 0.5183602715;
	Mi(1, 2) = 0.0492912282;
	Mi(2, 0) = -0.0085286646;
	Mi(2, 1) = 0.0400428217;
	Mi(2, 2) = 0.9684866958;

	boost::numeric::ublas::matrix<double> Gp(3, 1);
	Gp(0, 0) = _white.x / _white.y;
	Gp(1, 0) = 1;
	Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;

	boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod(M, Gp);

	boost::numeric::ublas::matrix<double> Hp(3, 1);
	Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
	Hp(1, 0) = 1;
	Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;

	boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod(M, Hp);

	boost::numeric::ublas::matrix<double> C = boost::numeric::ublas::zero_matrix<double>(3, 3);
	C(0, 0) = H(0, 0) / G(0, 0);
	C(1, 1) = H(1, 0) / G(1, 0);
	C(2, 2) = H(2, 0) / G(2, 0);

	boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod(C, M);
	return boost::numeric::ublas::prod(Mi, CM);
}




void
dcp::xyz_to_rgba(
	std::shared_ptr<const OpenJPEGImage> xyz_image,
	ColourConversion const & conversion,
	uint8_t* argb,
	int stride
	)
{
	int const max_colour = pow(2, 16) - 1;

	struct {
		double x, y, z;
	} s;

	struct {
		double r, g, b;
	} d;

	int* xyz_x = xyz_image->data(0);
	int* xyz_y = xyz_image->data(1);
	int* xyz_z = xyz_image->data(2);

	auto lut_in = conversion.out_j2k()->double_lut(0, 1, 12, false);
	auto lut_out = conversion.in()->double_lut(0, 1, 16, true);
	boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb();

	double fast_matrix[9] = {
		matrix(0, 0), matrix(0, 1), matrix(0, 2),
		matrix(1, 0), matrix(1, 1), matrix(1, 2),
		matrix(2, 0), matrix(2, 1), matrix(2, 2)
	};

	int const height = xyz_image->size().height;
	int const width = xyz_image->size().width;

	for (int y = 0; y < height; ++y) {
		uint8_t* argb_line = argb;
		for (int x = 0; x < width; ++x) {

			DCP_ASSERT(*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);

			/* In gamma LUT */
			s.x = lut_in[*xyz_x++];
			s.y = lut_in[*xyz_y++];
			s.z = lut_in[*xyz_z++];

			/* DCI companding */
			s.x /= DCI_COEFFICIENT;
			s.y /= DCI_COEFFICIENT;
			s.z /= DCI_COEFFICIENT;

			/* XYZ to RGB */
			d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
			d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
			d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));

			d.r = min(d.r, 1.0);
			d.r = max(d.r, 0.0);

			d.g = min(d.g, 1.0);
			d.g = max(d.g, 0.0);

			d.b = min(d.b, 1.0);
			d.b = max(d.b, 0.0);

			/* Out gamma LUT */
			*argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
			*argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
			*argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
			*argb_line++ = 0xff;
		}

		argb += stride;
	}
}


void
dcp::xyz_to_rgb(
	shared_ptr<const OpenJPEGImage> xyz_image,
	ColourConversion const & conversion,
	uint8_t* rgb,
	int stride,
	optional<NoteHandler> note
	)
{
	struct {
		double x, y, z;
	} s;

	struct {
		double r, g, b;
	} d;

	/* These should be 12-bit values from 0-4095 */
	int* xyz_x = xyz_image->data(0);
	int* xyz_y = xyz_image->data(1);
	int* xyz_z = xyz_image->data(2);

	auto lut_in = conversion.out_j2k()->double_lut(0, 1, 12, false);
	auto lut_out = conversion.in()->double_lut(0, 1, 16, true);
	auto const matrix = conversion.xyz_to_rgb();

	double fast_matrix[9] = {
		matrix(0, 0), matrix(0, 1), matrix(0, 2),
		matrix(1, 0), matrix(1, 1), matrix(1, 2),
		matrix(2, 0), matrix(2, 1), matrix(2, 2)
	};

	int const height = xyz_image->size().height;
	int const width = xyz_image->size().width;

	for (int y = 0; y < height; ++y) {
		auto rgb_line = reinterpret_cast<uint16_t*>(rgb + y * stride);
		for (int x = 0; x < width; ++x) {

			int cx = *xyz_x++;
			int cy = *xyz_y++;
			int cz = *xyz_z++;

			if (cx < 0 || cx > 4095) {
				if (note) {
					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cx));
				}
				cx = max(min(cx, 4095), 0);
			}

			if (cy < 0 || cy > 4095) {
				if (note) {
					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cy));
				}
				cy = max(min(cy, 4095), 0);
			}

			if (cz < 0 || cz > 4095) {
				if (note) {
					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cz));
				}
				cz = max(min(cz, 4095), 0);
			}

			/* In gamma LUT */
			s.x = lut_in[cx];
			s.y = lut_in[cy];
			s.z = lut_in[cz];

			/* DCI companding */
			s.x /= DCI_COEFFICIENT;
			s.y /= DCI_COEFFICIENT;
			s.z /= DCI_COEFFICIENT;

			/* XYZ to RGB */
			d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
			d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
			d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));

			d.r = min(d.r, 1.0);
			d.r = max(d.r, 0.0);

			d.g = min(d.g, 1.0);
			d.g = max(d.g, 0.0);

			d.b = min(d.b, 1.0);
			d.b = max(d.b, 0.0);

			*rgb_line++ = lrint(lut_out[lrint(d.r * 65535)] * 65535);
			*rgb_line++ = lrint(lut_out[lrint(d.g * 65535)] * 65535);
			*rgb_line++ = lrint(lut_out[lrint(d.b * 65535)] * 65535);
		}
	}
}

void
dcp::combined_rgb_to_xyz(ColourConversion const & conversion, double* matrix)
{
	auto const rgb_to_xyz = conversion.rgb_to_xyz();
	auto const bradford = conversion.bradford();

	matrix[0] = (bradford(0, 0) * rgb_to_xyz(0, 0) + bradford(0, 1) * rgb_to_xyz(1, 0) + bradford(0, 2) * rgb_to_xyz(2, 0)) * DCI_COEFFICIENT;
	matrix[1] = (bradford(0, 0) * rgb_to_xyz(0, 1) + bradford(0, 1) * rgb_to_xyz(1, 1) + bradford(0, 2) * rgb_to_xyz(2, 1)) * DCI_COEFFICIENT;
	matrix[2] = (bradford(0, 0) * rgb_to_xyz(0, 2) + bradford(0, 1) * rgb_to_xyz(1, 2) + bradford(0, 2) * rgb_to_xyz(2, 2)) * DCI_COEFFICIENT;
	matrix[3] = (bradford(1, 0) * rgb_to_xyz(0, 0) + bradford(1, 1) * rgb_to_xyz(1, 0) + bradford(1, 2) * rgb_to_xyz(2, 0)) * DCI_COEFFICIENT;
	matrix[4] = (bradford(1, 0) * rgb_to_xyz(0, 1) + bradford(1, 1) * rgb_to_xyz(1, 1) + bradford(1, 2) * rgb_to_xyz(2, 1)) * DCI_COEFFICIENT;
	matrix[5] = (bradford(1, 0) * rgb_to_xyz(0, 2) + bradford(1, 1) * rgb_to_xyz(1, 2) + bradford(1, 2) * rgb_to_xyz(2, 2)) * DCI_COEFFICIENT;
	matrix[6] = (bradford(2, 0) * rgb_to_xyz(0, 0) + bradford(2, 1) * rgb_to_xyz(1, 0) + bradford(2, 2) * rgb_to_xyz(2, 0)) * DCI_COEFFICIENT;
	matrix[7] = (bradford(2, 0) * rgb_to_xyz(0, 1) + bradford(2, 1) * rgb_to_xyz(1, 1) + bradford(2, 2) * rgb_to_xyz(2, 1)) * DCI_COEFFICIENT;
	matrix[8] = (bradford(2, 0) * rgb_to_xyz(0, 2) + bradford(2, 1) * rgb_to_xyz(1, 2) + bradford(2, 2) * rgb_to_xyz(2, 2)) * DCI_COEFFICIENT;
}


PiecewiseLUT2
dcp::make_inverse_gamma_lut(shared_ptr<const TransferFunction> fn, int scale)
{
	/* The parameters here were chosen by trial and error to reduce errors when running
	 * inverse_j2k_lut_test and inverse_mpeg2_lut_test.
	 */
	return PiecewiseLUT2(fn, 0.062, 16, 12, true, scale);
}


template <class T>
void
rgb_to_xyz_internal(
	uint8_t const* rgb,
	T*& xyz_x,
	T*& xyz_y,
	T*& xyz_z,
	dcp::Size size,
	int stride,
	ColourConversion const& conversion
	)
{
	struct {
		double r, g, b;
	} s;

	struct {
		double x, y, z;
	} d;

	auto lut_in = conversion.in()->double_lut(0, 1, 12, false);
	auto lut_out = make_inverse_gamma_lut(conversion.out_j2k(), 4095);

	/* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
	double fast_matrix[9];
	combined_rgb_to_xyz(conversion, fast_matrix);

	for (int y = 0; y < size.height; ++y) {
		auto p = reinterpret_cast<uint16_t const *>(rgb + y * stride);
		for (int x = 0; x < size.width; ++x) {

			/* In gamma LUT (converting 16-bit to 12-bit) */
			s.r = lut_in[*p++ >> 4];
			s.g = lut_in[*p++ >> 4];
			s.b = lut_in[*p++ >> 4];

			/* RGB to XYZ, Bradford transform and DCI companding */
			d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
			d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
			d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];

			/* Clamp */
			d.x = max(0.0, d.x);
			d.y = max(0.0, d.y);
			d.z = max(0.0, d.z);
			d.x = min(1.0, d.x);
			d.y = min(1.0, d.y);
			d.z = min(1.0, d.z);

			/* Out gamma LUT */
			*xyz_x++ = lut_out.lookup(d.x);
			*xyz_y++ = lut_out.lookup(d.y);
			*xyz_z++ = lut_out.lookup(d.z);
		}
	}
}


shared_ptr<dcp::OpenJPEGImage>
dcp::rgb_to_xyz(
	uint8_t const * rgb,
	dcp::Size size,
	int stride,
	ColourConversion const & conversion
	)
{
	auto xyz = make_shared<OpenJPEGImage>(size);

	int* xyz_x = xyz->data(0);
	int* xyz_y = xyz->data(1);
	int* xyz_z = xyz->data(2);

	rgb_to_xyz_internal(rgb, xyz_x, xyz_y, xyz_z, size, stride, conversion);

	return xyz;
}


void
dcp::rgb_to_xyz(
	uint8_t const * rgb,
	uint16_t* dst,
	dcp::Size size,
	int stride,
	ColourConversion const & conversion
	)
{
	rgb_to_xyz_internal(rgb, dst, dst, dst, size, stride, conversion);
}


void
dcp::convert_rgb_gamma(
	uint8_t* rgb,
	dcp::Size size,
	int stride,
	ColourConversion const& conversion
	)
{
	auto lut_in = conversion.in()->double_lut(0, 1, 8, false);
	auto lut_out = make_inverse_gamma_lut(conversion.out_mpeg2(), 255);

	for (int y = 0; y < size.height; ++y) {
		auto p = reinterpret_cast<uint8_t*>(rgb + y * stride);
		for (int x = 0; x < size.width; ++x) {

			/* In gamma LUT */
			auto r = lut_in[p[0]];
			auto g = lut_in[p[1]];
			auto b = lut_in[p[2]];

			/* Clamp */
			r = max(0.0, min(1.0, r));
			g = max(0.0, min(1.0, g));
			b = max(0.0, min(1.0, b));

			/* Out gamma LUT */
			*p++ = lut_out.lookup(r);
			*p++ = lut_out.lookup(g);
			*p++ = lut_out.lookup(b);
		}
	}
}