Bug 1799258 - Add prereq Colorspaces stuff, including generic gamma->linear LUT inversion approximation. r=gfx-reviewers,bradwerth

Differential Revision: https://phabricator.services.mozilla.com/D163664
This commit is contained in:
Kelsey Gilbert 2023-03-13 21:04:10 +00:00
Родитель 3bd71468e2
Коммит 000ff9b4e5
5 изменённых файлов: 841 добавлений и 9 удалений

Просмотреть файл

@ -498,7 +498,8 @@ enum class YUVRangedColorSpace : uint8_t {
// one.
// Some times Worse Is Better.
enum class ColorSpace2 : uint8_t {
UNKNOWN, // Really "DISPLAY". Eventually we will remove this.
Display,
UNKNOWN = Display, // We feel sufficiently bad about this TODO.
SRGB,
DISPLAY_P3,
BT601_525, // aka smpte170m NTSC
@ -506,7 +507,7 @@ enum class ColorSpace2 : uint8_t {
BT601_625 =
BT709, // aka bt470bg PAL. Basically BT709, just Xg is 0.290 not 0.300.
BT2020,
_First = UNKNOWN,
_First = Display,
_Last = BT2020,
};

Просмотреть файл

@ -7,6 +7,9 @@
#include "Colorspaces.h"
#include "nsDebug.h"
#include "qcms.h"
namespace mozilla::color {
// tf = { k * linear | linear < b
@ -76,6 +79,12 @@ mat4 YuvFromYcbcr(const YcbcrDesc& d) {
return yuvFromYcbcr;
}
inline vec3 CIEXYZ_from_CIExyY(const vec2 xy, const float Y = 1) {
const auto xyz = vec3(xy, 1 - xy.x() - xy.y());
const auto XYZ = xyz * (Y / xy.y());
return XYZ;
}
mat3 XyzFromLinearRgb(const Chromaticities& c) {
// http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
@ -115,7 +124,7 @@ mat3 XyzFromLinearRgb(const Chromaticities& c) {
const auto XYZrgb = mat3({Xrgb, Yrgb, Zrgb});
const auto XYZrgb_inv = inverse(XYZrgb);
const auto XYZwhitepoint = vec3({c.wx, c.wy, 1 - c.wx - c.wy});
const auto XYZwhitepoint = vec3({c.wx, c.wy, 1 - c.wx - c.wy}) / c.wy;
const auto Srgb = XYZrgb_inv * XYZwhitepoint;
const auto M = mat3({Srgb * Xrgb, Srgb * Yrgb, Srgb * Zrgb});
@ -183,6 +192,29 @@ vec3 ColorspaceTransform::DstFromSrc(const vec3 src) const {
// -
mat3 XyzAFromXyzB_BradfordLinear(const vec2 xyA, const vec2 xyB) {
// This is what ICC profiles use to do whitepoint transforms,
// because ICC also requires D50 for the Profile Connection Space.
// From https://www.color.org/specification/ICC.1-2022-05.pdf
// E.3 "Linearized Bradford transformation":
constexpr auto M_BFD = mat3{{
vec3{{0.8951, 0.2664f, -0.1614f}},
vec3{{-0.7502f, 1.7135f, 0.0367f}},
vec3{{0.0389f, -0.0685f, 1.0296f}},
}};
// NB: They use rho/gamma/beta, but we'll use R/G/B here.
const auto XYZDst = CIEXYZ_from_CIExyY(xyA); // "XYZ_W", WP of PCS
const auto XYZSrc = CIEXYZ_from_CIExyY(xyB); // "XYZ_NAW", WP of src
const auto rgbSrc = M_BFD * XYZSrc; // "RGB_SRC"
const auto rgbDst = M_BFD * XYZDst; // "RGB_PCS"
const auto rgbDstOverSrc = rgbDst / rgbSrc;
const auto M_dstOverSrc = mat3::Scale(rgbDstOverSrc);
const auto M_adapt = inverse(M_BFD) * M_dstOverSrc * M_BFD;
return M_adapt;
}
std::optional<mat4> ColorspaceTransform::ToMat4() const {
mat4 fromSrc = srcRgbTfFromSrc;
if (srcTf) return {};
@ -226,4 +258,169 @@ vec3 Lut3::Sample(const vec3 in01) const {
return fxyz;
}
// -
ColorProfileDesc ColorProfileDesc::From(const ColorspaceDesc& cspace) {
auto ret = ColorProfileDesc{};
if (cspace.yuv) {
const auto yuvFromYcbcr = YuvFromYcbcr(cspace.yuv->ycbcr);
const auto yuvFromRgb = YuvFromRgb(cspace.yuv->yCoeffs);
const auto rgbFromYuv = inverse(yuvFromRgb);
ret.rgbFromYcbcr = mat4(rgbFromYuv) * yuvFromYcbcr;
}
if (cspace.tf) {
const size_t tableSize = 256;
auto& tableR = ret.linearFromTf.r;
tableR.resize(tableSize);
for (size_t i = 0; i < tableR.size(); i++) {
const float tfVal = i / float(tableR.size() - 1);
const float linearVal = LinearFromTf(*cspace.tf, tfVal);
tableR[i] = linearVal;
}
ret.linearFromTf.g = tableR;
ret.linearFromTf.b = tableR;
}
ret.xyzd65FromLinearRgb = XyzFromLinearRgb(cspace.chrom);
return ret;
}
// -
template <class T>
constexpr inline T NewtonEstimateX(const T x1, const T y1, const T dydx,
const T y2 = 0) {
// Estimate x s.t. y=0
// y = y0 + x*dydx;
// y0 = y - x*dydx;
// y1 - x1*dydx = y2 - x2*dydx
// x2*dydx = y2 - y1 + x1*dydx
// x2 = (y2 - y1)/dydx + x1
return (y2 - y1) / dydx + x1;
}
float GuessGamma(const std::vector<float>& vals, float exp_guess) {
// Approximate (signed) error = 0.0.
constexpr float d_exp = 0.001;
constexpr float error_tolerance = 0.001;
struct Samples {
float y1, y2;
};
const auto Sample = [&](const float exp) {
int i = -1;
auto samples = Samples{};
for (const auto& expected : vals) {
i += 1;
const auto in = i / float(vals.size() - 1);
samples.y1 += powf(in, exp) - expected;
samples.y2 += powf(in, exp + d_exp) - expected;
}
samples.y1 /= vals.size(); // Normalize by val count.
samples.y2 /= vals.size();
return samples;
};
constexpr int MAX_ITERS = 10;
for (int i = 1;; i++) {
const auto err = Sample(exp_guess);
const auto derr = err.y2 - err.y1;
exp_guess = NewtonEstimateX(exp_guess, err.y1, derr / d_exp);
// Check if we were close before, because then this last round of estimation
// should get us pretty much right on it.
if (std::abs(err.y1) < error_tolerance) {
return exp_guess;
}
if (i >= MAX_ITERS) {
printf_stderr("GuessGamma() -> %f after %i iterations (avg err %f)\n",
exp_guess, i, err.y1);
MOZ_ASSERT(false, "GuessGamma failed.");
return exp_guess;
}
}
}
// -
ColorProfileDesc ColorProfileDesc::From(const qcms_profile& qcmsProfile) {
ColorProfileDesc ret;
qcms_profile_data data = {};
qcms_profile_get_data(&qcmsProfile, &data);
auto xyzd50FromLinearRgb = mat3{};
// X contributions from [R,G,B]
xyzd50FromLinearRgb.at(0, 0) = data.red_colorant_xyzd50[0];
xyzd50FromLinearRgb.at(1, 0) = data.green_colorant_xyzd50[0];
xyzd50FromLinearRgb.at(2, 0) = data.blue_colorant_xyzd50[0];
// Y contributions from [R,G,B]
xyzd50FromLinearRgb.at(0, 1) = data.red_colorant_xyzd50[1];
xyzd50FromLinearRgb.at(1, 1) = data.green_colorant_xyzd50[1];
xyzd50FromLinearRgb.at(2, 1) = data.blue_colorant_xyzd50[1];
// Z contributions from [R,G,B]
xyzd50FromLinearRgb.at(0, 2) = data.red_colorant_xyzd50[2];
xyzd50FromLinearRgb.at(1, 2) = data.green_colorant_xyzd50[2];
xyzd50FromLinearRgb.at(2, 2) = data.blue_colorant_xyzd50[2];
const auto d65FromD50 = XyzAFromXyzB_BradfordLinear(D65, D50);
ret.xyzd65FromLinearRgb = d65FromD50 * xyzd50FromLinearRgb;
// -
const auto Fn = [&](std::vector<float>* const linearFromTf,
int32_t claimed_samples,
const qcms_color_channel channel) {
if (claimed_samples == 0) return; // No tf.
if (claimed_samples == -1) {
claimed_samples = 4096; // Ask it to generate a bunch.
claimed_samples = 256; // Ask it to generate a bunch.
}
linearFromTf->resize(AssertedCast<size_t>(claimed_samples));
const auto begin = linearFromTf->data();
qcms_profile_get_lut(&qcmsProfile, channel, begin,
begin + linearFromTf->size());
};
Fn(&ret.linearFromTf.r, data.linear_from_trc_red_samples,
qcms_color_channel::Red);
Fn(&ret.linearFromTf.b, data.linear_from_trc_blue_samples,
qcms_color_channel::Blue);
Fn(&ret.linearFromTf.g, data.linear_from_trc_green_samples,
qcms_color_channel::Green);
// -
return ret;
}
// -
ColorProfileConversionDesc ColorProfileConversionDesc::From(
const FromDesc& desc) {
const auto dstLinearRgbFromXyzd65 = inverse(desc.dst.xyzd65FromLinearRgb);
auto ret = ColorProfileConversionDesc{
.srcRgbFromSrcYuv = desc.src.rgbFromYcbcr,
.srcLinearFromSrcTf = desc.src.linearFromTf,
.dstLinearFromSrcLinear =
dstLinearRgbFromXyzd65 * desc.src.xyzd65FromLinearRgb,
.dstTfFromDstLinear = {},
};
const auto Invert = [](const std::vector<float>& linearFromTf,
std::vector<float>* const tfFromLinear) {
const auto size = linearFromTf.size();
MOZ_ASSERT(size != 1); // Less than two is uninvertable.
if (size < 2) return;
(*tfFromLinear).resize(size);
InvertLut(linearFromTf, &*tfFromLinear);
};
Invert(desc.dst.linearFromTf.r, &ret.dstTfFromDstLinear.r);
Invert(desc.dst.linearFromTf.g, &ret.dstTfFromDstLinear.g);
Invert(desc.dst.linearFromTf.b, &ret.dstTfFromDstLinear.b);
return ret;
}
} // namespace mozilla::color

Просмотреть файл

@ -15,12 +15,16 @@
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <functional>
#include <optional>
#include <vector>
#include "AutoMappable.h"
#include "mozilla/Assertions.h"
#include "mozilla/Attributes.h"
#include "mozilla/Span.h"
#ifdef DEBUG
# define ASSERT(EXPR) \
@ -33,6 +37,9 @@
# define ASSERT(EXPR) (void)(EXPR)
#endif
struct _qcms_profile;
typedef struct _qcms_profile qcms_profile;
namespace mozilla::color {
struct YuvLumaCoeffs final {
@ -74,8 +81,7 @@ struct PiecewiseGammaDesc final {
4.5,
};
}
static constexpr auto Rec2020_10bit() { return Rec709(); }
// static constexpr auto Rec2020_10bit() { return Rec709(); }
static constexpr auto Rec2020_12bit() {
return PiecewiseGammaDesc{
1.0993,
@ -290,6 +296,7 @@ struct avec final {
return eq;
}
};
using vec2 = avec<float, 2>;
using vec3 = avec<float, 3>;
using vec4 = avec<float, 4>;
using ivec3 = avec<int32_t, 3>;
@ -364,10 +371,15 @@ struct mat final {
static constexpr auto Identity() {
auto ret = mat{};
for (int x = 0; x < x_cols; x++) {
for (int y = 0; y < y_rows; y++) {
ret.at(x, y) = (x == y ? 1 : 0);
}
for (int i = 0; i < std::min(x_cols, y_rows); i++) {
ret.at(i, i) = 1;
}
return ret;
}
static constexpr auto Scale(const avec<float, std::min(x_cols, y_rows)>& v) {
auto ret = mat{};
for (int i = 0; i < v.N; i++) {
ret.at(i, i) = v[i];
}
return ret;
}
@ -427,7 +439,30 @@ struct mat final {
}
return c;
}
// For e.g. similarity evaluation
friend auto operator-(const mat& a, const mat& b) {
mat c;
for (int y = 0; y < y_rows; y++) {
c.rows[y] = a.rows[y] - b.rows[y];
}
return c;
}
};
template <class M>
inline float dotDifference(const M& a, const M& b) {
const auto c = a - b;
const auto d = c * avec<float, M::x_cols>(1);
const auto d2 = dot(d, d);
return d2;
}
template <class M>
inline bool approx(const M& a, const M& b, const float eps = 0.0001) {
const auto errSquared = dotDifference(a, b);
return errSquared <= (eps * eps);
}
using mat3 = mat<3, 3>;
using mat4 = mat<4, 4>;
@ -651,6 +686,291 @@ struct ColorspaceTransform final {
}
};
// -
struct RgbTransferTables {
std::vector<float> r;
std::vector<float> g;
std::vector<float> b;
};
float GuessGamma(const std::vector<float>& vals, float exp_guess = 1.0);
static constexpr auto D65 = vec2{{0.3127, 0.3290}};
static constexpr auto D50 = vec2{{0.34567, 0.35850}};
mat3 XyzAFromXyzB_BradfordLinear(const vec2 xyA, const vec2 xyB);
// -
struct ColorProfileDesc {
// ICC profiles are phrased as PCS-from-encoded (PCS is CIEXYZ-D50)
// However, all of our colorspaces are D65, so let's normalize to that,
// even though it's a reverseable transform.
color::mat4 rgbFromYcbcr = color::mat4::Identity();
RgbTransferTables linearFromTf;
color::mat3 xyzd65FromLinearRgb = color::mat3::Identity();
static ColorProfileDesc From(const ColorspaceDesc&);
static ColorProfileDesc From(const qcms_profile&);
};
template <class C>
inline float SampleOutByIn(const C& outByIn, const float in) {
MOZ_ASSERT(outByIn.size() >= 2);
const auto begin = outByIn.begin();
const auto in0i = size_t(floorf(in * (outByIn.size() - 1)));
const auto out0_itr = begin + std::min(in0i, outByIn.size() - 2);
const auto in0 = float(out0_itr - begin) / (outByIn.size() - 1);
const auto out0 = *out0_itr;
const auto d_in = float(1) / (outByIn.size() - 1);
const auto d_out = *(out0_itr + 1) - *out0_itr;
const auto out = out0 + (d_out / d_in) * (in - in0);
// printf("SampleOutByIn(%f)->%f\n", in, out);
return out;
}
template <class C>
inline float SampleInByOut(const C& outByIn, const float out) {
MOZ_ASSERT(outByIn.size() >= 2);
const auto begin = outByIn.begin();
const auto out0_itr = std::lower_bound(begin + 1, outByIn.end() - 1, out) - 1;
const auto in0 = float(out0_itr - begin) / (outByIn.size() - 1);
const auto out0 = *out0_itr;
const auto d_in = float(1) / (outByIn.size() - 1);
const auto d_out = *(out0_itr + 1) - *out0_itr;
// printf("%f - (%f / %f) * (%f - %f)\n", in0, d_in, d_out, out, out0);
const auto in = in0 + (d_in / d_out) * (out - out0);
// printf("SampleInByOut(%f)->%f\n", out, in);
return in;
}
template <class C, class FnLessEqualT = std::less_equal<typename C::value_type>>
inline bool IsMonotonic(const C& vals, const FnLessEqualT& LessEqual = {}) {
bool ok = true;
const auto begin = vals.begin();
for (size_t i = 1; i < vals.size(); i++) {
const auto itr = begin + i;
ok &= LessEqual(*(itr - 1), *itr);
// Assert(true, [&]() {
// return prints("[%zu]->%f <= [%zu]->%f", i-1, *(itr-1), i, *itr);
// });
}
return ok;
}
template <class T, class I>
inline std::optional<I> SeekNeq(const T& ref, const I first, const I last) {
const auto inc = (last - first) > 0 ? 1 : -1;
auto itr = first;
while (true) {
if (*itr != ref) return itr;
if (itr == last) return {};
itr += inc;
}
}
template <class T>
struct TwoPoints {
struct {
T x;
T y;
} p0;
struct {
T x;
T y;
} p1;
T y(const T x) const {
return p0.y + (x - p0.x) / (p1.x - p0.x) * (p1.y - p0.y);
}
};
template <class T>
static void LinearFill(T& vals, const TwoPoints<float>& line) {
float x = -1;
for (auto& val : vals) {
x += 1;
val = line.y(x);
}
}
// -
inline void DequantizeMonotonic(const Span<float> vals) {
MOZ_ASSERT(IsMonotonic(vals));
const auto first = vals.begin();
const auto end = vals.end();
if (first == end) return;
const auto last = end - 1;
if (first == last) return;
// Three monotonic cases:
// 1. [0,0,0,0]
// 2. [0,0,1,1]
// 3. [0,1,1,2]
const auto body_first = SeekNeq(*first, first, last);
if (!body_first) {
// E.g. [0,0,0,0]
return;
}
const auto body_last = SeekNeq(*last, last, *body_first);
if (!body_last) {
// E.g. [0,0,1,1]
// This isn't the most accurate, but close enough.
// print("#2: %s", to_str(vals).c_str());
LinearFill(vals, {
{0, *first},
{float(vals.size() - 1), *last},
});
// print(" -> %s\n", to_str(vals).c_str());
return;
}
// E.g. [0,1,1,2]
// ^^^ body
// => f(0.5)->0.5, f(2.5)->1.5
// => f(x) = f(x0) + (x-x0) * (f(x1) - f(x0)) / (x1-x1)
// => f(x) = f(x0) + (x-x0) * dfdx
const auto head_end = *body_first;
const auto head = vals.subspan(0, head_end - vals.begin());
const auto tail_begin = *body_last + 1;
const auto tail = vals.subspan(tail_begin - vals.begin());
// print("head tail: %s %s\n",
// to_str(head).c_str(),
// to_str(tail).c_str());
// const auto body = vals->subspan(head.size(), vals->size()-tail.size());
auto next_part_first = head_end;
while (next_part_first != tail_begin) {
const auto part_first = next_part_first;
// print("part_first: %f\n", *part_first);
next_part_first = *SeekNeq(*part_first, part_first, tail_begin);
// print("next_part_first: %f\n", *next_part_first);
const auto part =
Span<float>{part_first, size_t(next_part_first - part_first)};
// print("part: %s\n", to_str(part).c_str());
const auto prev_part_last = part_first - 1;
const auto part_last = next_part_first - 1;
const auto line = TwoPoints<float>{
{-0.5, (*prev_part_last + *part_first) / 2},
{part.size() - 0.5f, (*part_last + *next_part_first) / 2},
};
LinearFill(part, line);
}
static constexpr bool INFER_HEAD_TAIL_FROM_BODY_EDGE = false;
// Basically ignore contents of head and tail, and infer from edges of body.
// print("3: %s\n", to_str(vals).c_str());
if (!IsMonotonic(head, std::less<float>{})) {
if (!INFER_HEAD_TAIL_FROM_BODY_EDGE) {
LinearFill(head,
{
{0, *head.begin()},
{head.size() - 0.5f, (*(head.end() - 1) + *head_end) / 2},
});
} else {
LinearFill(head, {
{head.size() + 0.0f, *head_end},
{head.size() + 1.0f, *(head_end + 1)},
});
}
}
if (!IsMonotonic(tail, std::less<float>{})) {
if (!INFER_HEAD_TAIL_FROM_BODY_EDGE) {
LinearFill(tail, {
{-0.5, (*(tail_begin - 1) + *tail.begin()) / 2},
{tail.size() - 1.0f, *(tail.end() - 1)},
});
} else {
LinearFill(tail, {
{-2.0f, *(tail_begin - 2)},
{-1.0f, *(tail_begin - 1)},
});
}
}
// print("3: %s\n", to_str(vals).c_str());
MOZ_ASSERT(IsMonotonic(vals, std::less<float>{}));
// Rescale, because we tend to lose range.
static constexpr bool RESCALE = false;
if (RESCALE) {
const auto firstv = *first;
const auto lastv = *last;
for (auto& val : vals) {
val = (val - firstv) / (lastv - firstv);
}
}
// print("4: %s\n", to_str(vals).c_str());
}
template <class In, class Out>
static void InvertLut(const In& lut, Out* const out_invertedLut) {
MOZ_ASSERT(IsMonotonic(lut));
auto plut = &lut;
auto vec = std::vector<float>{};
if (!IsMonotonic(lut, std::less<float>{})) {
// print("Not strictly monotonic...\n");
vec.assign(lut.begin(), lut.end());
DequantizeMonotonic(vec);
plut = &vec;
// print(" Now strictly monotonic: %i: %s\n",
// int(IsMonotonic(*plut, std::less<float>{})), to_str(*plut).c_str());
MOZ_ASSERT(IsMonotonic(*plut, std::less<float>{}));
}
MOZ_ASSERT(plut->size() >= 2);
auto& ret = *out_invertedLut;
for (size_t i_out = 0; i_out < ret.size(); i_out++) {
const auto f_out = i_out / float(ret.size() - 1);
const auto f_in = SampleInByOut(*plut, f_out);
ret[i_out] = f_in;
}
MOZ_ASSERT(IsMonotonic(ret));
MOZ_ASSERT(IsMonotonic(ret, std::less<float>{}));
}
// -
struct ColorProfileConversionDesc {
// ICC profiles are phrased as PCS-from-encoded (PCS is CIEXYZ-D50)
color::mat4 srcRgbFromSrcYuv = color::mat4::Identity();
RgbTransferTables srcLinearFromSrcTf;
color::mat3 dstLinearFromSrcLinear = color::mat3::Identity();
RgbTransferTables dstTfFromDstLinear;
struct FromDesc {
ColorProfileDesc src;
ColorProfileDesc dst;
};
static ColorProfileConversionDesc From(const FromDesc&);
vec3 Apply(const vec3 src) const {
const auto srcRgb = vec3(srcRgbFromSrcYuv * vec4(src, 1));
const auto srcLinear = vec3{{
SampleOutByIn(srcLinearFromSrcTf.r, srcRgb.x()),
SampleOutByIn(srcLinearFromSrcTf.g, srcRgb.y()),
SampleOutByIn(srcLinearFromSrcTf.b, srcRgb.z()),
}};
const auto dstLinear = dstLinearFromSrcLinear * srcLinear;
const auto dstRgb = vec3{{
SampleOutByIn(dstTfFromDstLinear.r, dstLinear.x()),
SampleOutByIn(dstTfFromDstLinear.g, dstLinear.y()),
SampleOutByIn(dstTfFromDstLinear.b, dstLinear.z()),
}};
return dstRgb;
}
};
} // namespace mozilla::color
#undef ASSERT

Просмотреть файл

@ -6,12 +6,14 @@
#include "gtest/gtest.h"
#include "Colorspaces.h"
#include <array>
#include <limits>
namespace mozilla::color {
mat4 YuvFromYcbcr(const YcbcrDesc&);
float TfFromLinear(const PiecewiseGammaDesc&, float linear);
float LinearFromTf(const PiecewiseGammaDesc&, float tf);
mat3 XyzFromLinearRgb(const Chromaticities&);
} // namespace mozilla::color
using namespace mozilla::color;
@ -314,6 +316,24 @@ TEST(Colorspaces, SrgbFromDisplayP3)
// -
template <class Fn, class Tuple, size_t... I>
constexpr auto map_tups_seq(const Tuple& a, const Tuple& b, const Fn& fn,
std::index_sequence<I...>) {
return std::tuple{fn(std::get<I>(a), std::get<I>(b))...};
}
template <class Fn, class Tuple>
constexpr auto map_tups(const Tuple& a, const Tuple& b, const Fn& fn) {
return map_tups_seq(a, b, fn,
std::make_index_sequence<std::tuple_size_v<Tuple>>{});
}
template <class Fn, class Tuple>
constexpr auto cmp_tups_all(const Tuple& a, const Tuple& b, const Fn& fn) {
bool all = true;
map_tups(a, b, [&](const auto& a, const auto& b) { return all &= fn(a, b); });
return all;
}
struct Stats {
double mean = 0;
double variance = 0;
@ -329,6 +349,7 @@ struct Stats {
ret.max = std::max(ret.max, cur);
}
ret.mean /= iterable.size();
// Gather mean first before we can calc variance.
for (const auto& cur : iterable) {
ret.variance += pow(cur - ret.mean, 2);
}
@ -336,8 +357,61 @@ struct Stats {
return ret;
}
template <class T, class U>
static Stats Diff(const T& a, const U& b) {
MOZ_ASSERT(a.size() == b.size());
std::vector<double> diff;
diff.reserve(a.size());
for (size_t i = 0; i < diff.capacity(); i++) {
diff.push_back(a[i] - b[i]);
}
return Stats::For(diff);
}
double standardDeviation() const { return sqrt(variance); }
friend std::ostream& operator<<(std::ostream& s, const Stats& a) {
return s << "Stats"
<< "{ mean:" << a.mean << ", stddev:" << a.standardDeviation()
<< ", min:" << a.min << ", max:" << a.max << " }";
}
struct Error {
double absmean = std::numeric_limits<double>::infinity();
double stddev = std::numeric_limits<double>::infinity();
double absmax = std::numeric_limits<double>::infinity();
constexpr auto Fields() const { return std::tie(absmean, stddev, absmax); }
template <class Fn>
friend constexpr bool cmp_all(const Error& a, const Error& b,
const Fn& fn) {
return cmp_tups_all(a.Fields(), b.Fields(), fn);
}
friend constexpr bool operator<(const Error& a, const Error& b) {
return cmp_all(a, b, [](const auto& a, const auto& b) { return a < b; });
}
friend constexpr bool operator<=(const Error& a, const Error& b) {
return cmp_all(a, b, [](const auto& a, const auto& b) { return a <= b; });
}
friend std::ostream& operator<<(std::ostream& s, const Error& a) {
return s << "Stats::Error"
<< "{ absmean:" << a.absmean << ", stddev:" << a.stddev
<< ", absmax:" << a.absmax << " }";
}
};
operator Error() const {
return {abs(mean), standardDeviation(), std::max(abs(min), abs(max))};
}
};
static_assert(Stats::Error{0, 0, 0} < Stats::Error{1, 1, 1});
static_assert(!(Stats::Error{0, 1, 0} < Stats::Error{1, 1, 1}));
static_assert(Stats::Error{0, 1, 0} <= Stats::Error{1, 1, 1});
static_assert(!(Stats::Error{0, 2, 0} <= Stats::Error{1, 1, 1}));
// -
static Stats StatsForLutError(const ColorspaceTransform& ct,
const ivec3 srcQuants, const ivec3 dstQuants) {
@ -403,3 +477,222 @@ TEST(Colorspaces, LutError_Rec709Full_Srgb)
EXPECT_NEAR(stats.min, 0, 0.001);
EXPECT_NEAR(stats.max, 17, 0.001);
}
// -
// https://www.reedbeta.com/blog/python-like-enumerate-in-cpp17/
template <typename T, typename TIter = decltype(std::begin(std::declval<T>())),
typename = decltype(std::end(std::declval<T>()))>
constexpr auto enumerate(T&& iterable) {
struct iterator {
size_t i;
TIter iter;
bool operator!=(const iterator& other) const { return iter != other.iter; }
void operator++() {
++i;
++iter;
}
auto operator*() const { return std::tie(i, *iter); }
};
struct iterable_wrapper {
T iterable;
auto begin() { return iterator{0, std::begin(iterable)}; }
auto end() { return iterator{0, std::end(iterable)}; }
};
return iterable_wrapper{std::forward<T>(iterable)};
}
inline auto MakeLinear(const float from, const float to, const int n) {
std::vector<float> ret;
ret.resize(n);
for (auto [i, val] : enumerate(ret)) {
const auto t = i / float(ret.size() - 1);
val = from + (to - from) * t;
}
return ret;
};
inline auto MakeGamma(const float exp, const int n) {
std::vector<float> ret;
ret.resize(n);
for (auto [i, val] : enumerate(ret)) {
const auto t = i / float(ret.size() - 1);
val = powf(t, exp);
}
return ret;
};
// -
TEST(Colorspaces, GuessGamma)
{
EXPECT_NEAR(GuessGamma(MakeGamma(1, 11)), 1.0, 0);
EXPECT_NEAR(GuessGamma(MakeGamma(2.2, 11)), 2.2, 4.8e-8);
EXPECT_NEAR(GuessGamma(MakeGamma(1 / 2.2, 11)), 1 / 2.2, 1.1e-7);
}
// -
template <class T, class U>
float StdDev(const T& test, const U& ref) {
float sum = 0;
for (size_t i = 0; i < test.size(); i++) {
const auto diff = test[i] - ref[i];
sum += diff * diff;
}
const auto variance = sum / test.size();
return sqrt(variance);
}
template <class T>
inline void AutoLinearFill(T& vals) {
LinearFill(vals, {
{0, 0},
{vals.size() - 1.0f, 1},
});
}
template <class T, class... More>
auto MakeArray(const T& a0, const More&... args) {
return std::array<T, 1 + sizeof...(More)>{a0, static_cast<float>(args)...};
}
TEST(Colorspaces, LinearFill)
{
EXPECT_NEAR(StdDev(MakeLinear(0, 1, 3), MakeArray<float>(0, 0.5, 1)), 0,
0.001);
auto vals = std::vector<float>(3);
LinearFill(vals, {
{0, 0},
{vals.size() - 1.0f, 1},
});
EXPECT_NEAR(StdDev(vals, MakeArray<float>(0, 0.5, 1)), 0, 0.001);
LinearFill(vals, {
{0, 1},
{vals.size() - 1.0f, 0},
});
EXPECT_NEAR(StdDev(vals, MakeArray<float>(1, 0.5, 0)), 0, 0.001);
}
TEST(Colorspaces, DequantizeMonotonic)
{
auto orig = std::vector<float>{0, 0, 0, 1, 1, 2};
auto vals = orig;
EXPECT_TRUE(IsMonotonic(vals));
EXPECT_TRUE(!IsMonotonic(vals, std::less<float>{}));
DequantizeMonotonic(vals);
EXPECT_TRUE(IsMonotonic(vals, std::less<float>{}));
EXPECT_LT(StdDev(vals, orig),
StdDev(MakeLinear(orig.front(), orig.back(), vals.size()), orig));
}
TEST(Colorspaces, InvertLut)
{
const auto linear = MakeLinear(0, 1, 256);
auto linearFromSrgb = linear;
for (auto& val : linearFromSrgb) {
val = powf(val, 2.2);
}
auto srgbFromLinearExpected = linear;
for (auto& val : srgbFromLinearExpected) {
val = powf(val, 1 / 2.2);
}
auto srgbFromLinearViaInvert = linearFromSrgb;
InvertLut(linearFromSrgb, &srgbFromLinearViaInvert);
// I just want to appreciate that InvertLut is a non-analytical approximation,
// and yet it's extraordinarily close to the analytical inverse.
EXPECT_LE(Stats::Diff(srgbFromLinearViaInvert, srgbFromLinearExpected),
(Stats::Error{3e-6, 3e-6, 3e-5}));
const auto srcSrgb = MakeLinear(0, 1, 256);
auto roundtripSrgb = srcSrgb;
for (auto& srgb : roundtripSrgb) {
const auto linear = SampleOutByIn(linearFromSrgb, srgb);
const auto srgb2 = SampleOutByIn(srgbFromLinearViaInvert, linear);
// printf("[%f] %f -> %f -> %f\n", srgb2-srgb, srgb, linear, srgb2);
srgb = srgb2;
}
EXPECT_LE(Stats::Diff(roundtripSrgb, srcSrgb),
(Stats::Error{0.0013, 0.0046, 0.023}));
}
TEST(Colorspaces, XyzFromLinearRgb)
{
const auto xyzd65FromLinearRgb = XyzFromLinearRgb(Chromaticities::Srgb());
// http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
const auto XYZD65_FROM_LINEAR_RGB = mat3({
vec3{{0.4124564, 0.3575761, 0.1804375}},
vec3{{0.2126729, 0.7151522, 0.0721750}},
vec3{{0.0193339, 0.1191920, 0.9503041}},
});
EXPECT_NEAR(sqrt(dotDifference(xyzd65FromLinearRgb, XYZD65_FROM_LINEAR_RGB)),
0, 0.001);
}
TEST(Colorspaces, ColorProfileConversionDesc_SrgbFromRec709)
{
const auto srgb = ColorProfileDesc::From({
Chromaticities::Srgb(),
PiecewiseGammaDesc::Srgb(),
});
const auto rec709 = ColorProfileDesc::From({
Chromaticities::Rec709(),
PiecewiseGammaDesc::Rec709(),
});
{
const auto conv = ColorProfileConversionDesc::From({
.src = srgb,
.dst = srgb,
});
auto src = vec3(16.0);
auto dst = conv.Apply(src / 255) * 255;
const auto tfa = PiecewiseGammaDesc::Srgb();
const auto tfb = PiecewiseGammaDesc::Srgb();
const auto expected =
TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
printf("%f %f %f\n", src.x(), src.y(), src.z());
printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{0.42}));
}
{
const auto conv = ColorProfileConversionDesc::From({
.src = rec709,
.dst = rec709,
});
auto src = vec3(16.0);
auto dst = conv.Apply(src / 255) * 255;
const auto tfa = PiecewiseGammaDesc::Rec709();
const auto tfb = PiecewiseGammaDesc::Rec709();
const auto expected =
TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
printf("%f %f %f\n", src.x(), src.y(), src.z());
printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{1e-6}));
}
{
const auto conv = ColorProfileConversionDesc::From({
.src = rec709,
.dst = srgb,
});
auto src = vec3(16.0);
auto dst = conv.Apply(src / 255) * 255;
const auto tfa = PiecewiseGammaDesc::Rec709();
const auto tfb = PiecewiseGammaDesc::Srgb();
const auto expected =
TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
printf("expected: %f\n", expected);
printf("%f %f %f\n", src.x(), src.y(), src.z());
printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{0.12}));
}
}

Просмотреть файл

@ -366,6 +366,7 @@ class Span {
public:
// constants and types
using element_type = ElementType;
using value_type = std::remove_cv_t<element_type>;
using index_type = size_t;
using pointer = element_type*;
using reference = element_type&;
@ -420,6 +421,16 @@ class Span {
aEnd)
: storage_(aBegin == aEnd ? nullptr : &*aBegin, aEnd - aBegin) {}
/**
* Constructor for {iterator,size_t}
*/
template <typename OtherElementType, size_t OtherExtent, bool IsConst>
constexpr Span(
span_details::span_iterator<Span<OtherElementType, OtherExtent>, IsConst>
aBegin,
index_type aLength)
: storage_(!aLength ? nullptr : &*aBegin, aLength) {}
/**
* Constructor for C array.
*/
@ -664,6 +675,16 @@ class Span {
return Subspan(0, aEnd);
}
/// std::span-compatible method name
constexpr auto subspan(index_type aStart,
index_type aLength = dynamic_extent) const {
return Subspan(aStart, aLength);
}
/// std::span-compatible method name
constexpr auto from(index_type aStart) const { return From(aStart); }
/// std::span-compatible method name
constexpr auto to(index_type aEnd) const { return To(aEnd); }
/**
* Subspan with run-time start index and exclusive end index.
* (Rust's &foo[start..end])