Extract, test and slightly adjust the piecewise LUT for inverse gamma.

4 files changed, 108 insertions, 12 deletions

diff --git a/src/piecewise_lut.h b/src/piecewise_lut.h
new file mode 100644
index 00000000..ba8f5c13
--- /dev/null
+++ b/src/piecewise_lut.h
@@ -0,0 +1,76 @@
+/*
+    Copyright (C) 2013-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.
+*/
+
+
+#ifndef LIBDCP_PIECEWISE_LUT_H
+#define LIBDCP_PIECEWISE_LUT_H
+
+
+#include "transfer_function.h"
+#include <memory>
+#include <vector>
+
+
+namespace dcp {
+
+
+class PiecewiseLUT2
+{
+public:
+	PiecewiseLUT2(std::shared_ptr<const TransferFunction> fn, double boundary, int low_bits, int high_bits, bool inverse)
+		: _boundary(boundary)
+		, _low(fn->lut(0, boundary, low_bits, inverse))
+		, _high(fn->lut(boundary, 1, high_bits, inverse))
+		, _low_scale(static_cast<int>(std::pow(2.0f, low_bits)) - 1)
+		, _high_scale(static_cast<int>(std::pow(2.0f, high_bits)) - 1)
+	{
+
+	}
+
+	inline double lookup(double x) const {
+		return x < _boundary ?  _low[lrint((x / _boundary) * _low_scale)] : _high[lrint(((x - _boundary) / (1 - _boundary)) * _high_scale)];
+	}
+
+private:
+	double _boundary;
+	std::vector<double> _low;
+	std::vector<double> _high;
+	int _low_scale;
+	int _high_scale;
+};
+
+
+}
+
+#endif
+
diff --git a/src/rgb_xyz.cc b/src/rgb_xyz.cc
index 5c01f1c5..988f9bac 100644
--- a/src/rgb_xyz.cc
+++ b/src/rgb_xyz.cc
@@ -41,6 +41,7 @@
 #include "compose.hpp"
 #include "dcp_assert.h"
 #include "openjpeg_image.h"
+#include "piecewise_lut.h"
 #include "rgb_xyz.h"
 #include "transfer_function.h"
 #include <cmath>
@@ -257,6 +258,14 @@ dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
 }
 
 
+PiecewiseLUT2
+dcp::make_inverse_gamma_lut(shared_ptr<const TransferFunction> fn)
+{
+	/* The parameters here were chosen by trial and error to reduce errors when running rgb_xyz_lut_test */
+	return PiecewiseLUT2(fn, 0.062, 16, 12, true);
+}
+
+
 shared_ptr<dcp::OpenJPEGImage>
 dcp::rgb_to_xyz (
 	uint8_t const * rgb,
@@ -277,13 +286,7 @@ dcp::rgb_to_xyz (
 	} d;
 
 	auto lut_in = conversion.in()->lut(0, 1, 12, false);
-
-	/* Use 2 separate LUTs for the output gamma to keep precision high enough for small values
-	 * where the curve is steep.
-	 */
-	auto constexpr piece = 0.02;
-	auto lut_out_low = conversion.out()->lut(0, piece, 16, true);
-	auto lut_out_high = conversion.out()->lut(piece, 1, 12, true);
+	auto lut_out = make_inverse_gamma_lut(conversion.out());
 
 	/* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 	double fast_matrix[9];
@@ -321,9 +324,9 @@ dcp::rgb_to_xyz (
 			d.z = min (1.0, d.z);
 
 			/* Out gamma LUT */
-			*xyz_x++ = lrint((d.x < piece ? lut_out_low[lrint((d.x / piece) * 65535)] : lut_out_high[lrint(((d.x - piece) / (1 - piece)) * 4095)]) * 4095);
-			*xyz_y++ = lrint((d.y < piece ? lut_out_low[lrint((d.y / piece) * 65535)] : lut_out_high[lrint(((d.y - piece) / (1 - piece)) * 4095)]) * 4095);
-			*xyz_z++ = lrint((d.z < piece ? lut_out_low[lrint((d.z / piece) * 65535)] : lut_out_high[lrint(((d.z - piece) / (1 - piece)) * 4095)]) * 4095);
+			*xyz_x++ = lrint(lut_out.lookup(d.x) * 4095);
+			*xyz_y++ = lrint(lut_out.lookup(d.y) * 4095);
+			*xyz_z++ = lrint(lut_out.lookup(d.z) * 4095);
 		}
 	}
 
diff --git a/src/rgb_xyz.h b/src/rgb_xyz.h
index 4a920efb..f0623940 100644
--- a/src/rgb_xyz.h
+++ b/src/rgb_xyz.h
@@ -37,6 +37,7 @@
  */
 
 
+#include "piecewise_lut.h"
 #include "types.h"
 #include <memory>
 #include <boost/optional.hpp>
@@ -92,6 +93,9 @@ extern void xyz_to_rgb (
 	);
 
 
+extern PiecewiseLUT2 make_inverse_gamma_lut(std::shared_ptr<const TransferFunction> fn);
+
+
 /** @param rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
  *  with the 2-byte value for each R/G/B component stored as
  *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
diff --git a/test/rgb_xyz_test.cc b/test/rgb_xyz_test.cc
index f64debc4..00787938 100644
--- a/test/rgb_xyz_test.cc
+++ b/test/rgb_xyz_test.cc
@@ -32,9 +32,10 @@
 */
 
 
-#include "rgb_xyz.h"
-#include "openjpeg_image.h"
 #include "colour_conversion.h"
+#include "openjpeg_image.h"
+#include "piecewise_lut.h"
+#include "rgb_xyz.h"
 #include "stream_operators.h"
 #include <boost/bind.hpp>
 #include <boost/random.hpp>
@@ -149,6 +150,18 @@ BOOST_AUTO_TEST_CASE (rgb_xyz_test)
 }
 
 
+/** Check the piecewise LUT that is used for inverse gamma calculation */
+BOOST_AUTO_TEST_CASE (rgb_xyz_lut_test)
+{
+	auto conversion = dcp::ColourConversion::rec709_to_xyz();
+	auto lut = dcp::make_inverse_gamma_lut(conversion.out());
+
+	for (double x = 0; x < 1; x += 0.000001) {
+		BOOST_CHECK(std::abs(lrint(lut.lookup(x) * 4095) - lrint(pow(x, 1 / 2.6) * 4095)) < 2);
+	}
+}
+
+
 static list<string> notes;
 
 static void