#include <iostream>
#include <Eigen/Dense>
#include <vector>
using namespace Eigen;
struct Display {
std::string name;
Matrix3d rgb2xyz;
Matrix3d xyz2rgb;
double gain;
Vector3d offset;
Vector3d gamma;
VectorXd gammatest;
};
Display o2Ics() {
Display disp;
Vector2d xyW(0.3127, 0.3290);
Vector2d xyR(0.708, 0.292);
Vector2d xyG(0.170, 0.797);
Vector2d xyB(0.131, 0.046);
Matrix3d rgb2xyz;
rgb2xyz << xyR, xyG, xyB, Vector2d::Ones() - xyR - xyG - xyB;
Vector3d rgbWei = rgb2xyz.inverse() * Vector3d(xyW, Vector3d::Ones() - Vector3d(xyW.sum(), xyW.sum(), xyW.sum()));
rgb2xyz *= rgbWei.asDiagonal();
rgb2xyz /= rgb2xyz.row(1).sum();
double beta = 0.0181;
double alpha = 1.0993;
VectorXd in = VectorXd::LinSpaced(65536, 0, 1);
VectorXd invGamma = (4.5 * in.array()).min(beta).eval() +
((alpha
* in.
array().
pow(0.45) - (alpha
- 1)).
max(0)).
eval(); VectorXd gamma = (in.array() < beta).select(4.5 * in.array(), invGamma);
disp.name = "Rec.2020";
disp.rgb2xyz = rgb2xyz;
disp.xyz2rgb = rgb2xyz.inverse();
disp.gain = 225;
disp.offset.setZero();
disp.gamma = gamma.replicate(3, 1);
disp.gammatest = ((1 / 4.5) * in.array()).min(beta).eval() +
(((1 / alpha
) * in.
array() + beta
/ 4.5).
pow(1 / 0.45)).
max(0).
eval();
return disp;
}
void ConvertSonyPrimaries2ccm(const std::vector<double>& Rx, const std::vector<double>& Ry,
double Luminance, Matrix3d& ccm_colorOptimized,
Matrix3d& ccm_brightnessOptimized) {
double Xw = x(x.size() - 1) * Luminance / y(y.size() - 1);
double Zw = (1 - x(x.size() - 1) - y(y.size() - 1)) * Luminance / y(y.size() - 1);
double Yw = Luminance;
Vector3d z = Vector3d::Ones() - x.head(3) - y.head(3);
Matrix3d A;
A << x.head(3), y.head(3), z;
Vector3d B(Xw, Yw, Zw);
Vector3d S = A.colPivHouseholderQr().solve(B);
Vector3d X = S.array() * x.head(3).array();
Vector3d Y = S.array() * y.head(3).array();
Vector3d Z = S.array() * z.array();
Matrix3d XYZ;
XYZ << X, Y, Z;
double dispGain = Luminance;
Matrix3d dispXYZnorm = XYZ.array() / dispGain;
Vector3d checkWhite = dispXYZnorm.row(1).array().sum() * Vector3d::Ones();
if ((checkWhite - Vector3d::Ones()).cwiseAbs().maxCoeff() > 0.03) {
std::cerr << "*** WARNING: Linearity seems off for white. MUST CHECK. ***" << std::endl;
}
Matrix3d disp2xyz = dispXYZnorm.array().colwise() / checkWhite.array();
// Assuming you have dispSony.rgb2xyz defined elsewhere
Matrix3d dispIcs; // Initialize dispIcs
ccm_brightnessOptimized = dispSony.rgb2xyz.inverse() * dispIcs.rgb2xyz;
Matrix3d mat;
mat << 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1;
ccm_colorOptimized = ccm_brightnessOptimized.array() / (ccm_brightnessOptimized * mat).maxCoeff();
if (ccm_colorOptimized.minCoeff() < -std::numeric_limits<double>::epsilon()) {
std::cerr << "*** WARNING: Some points in source can't be simulated on the destination. ***" << std::endl;
}
}
int main() {
// Given input vectors
std::vector<double> Rx = {0.660, 0.720};
std::vector<double> Ry = {0.277, 0.337};
std::vector<double> Gx = {0.208, 0.268};
std::vector<double> Gy = {0.663, 0.723};
std::vector<double> Bx = {0.117, 0.177};
std::vector<double> By = {0.026, 0.086};
std::vector<double> Wx = {0.293, 0.333};
std::vector<double> Wy = {0.309, 0.349};
// Calculate intermediate values
double xr = Rx[0] + ((Rx[1] - Rx[0]) / 2);
double yr = Ry[0] + ((Ry[1] - Ry[0]) / 2);
double xg = Gx[0] + ((Gx[1] - Gx[0]) / 2);
double yg = Gy[0] + ((Gy[1] - Gy[0]) / 2);
double xb = Bx[0] + ((Bx[1] - Bx[0]) / 2);
double yb = By[0] + ((By[1] - By[0]) / 2);
double xw = Wx[0] + ((Wx[1] - Wx[0]) / 2);
double yw = Wy[0] + ((Wy[1] - Wy[0]) / 2);
// Create vectors x and y
VectorXd x(4), y(4);
x << xr, xg, xb, xw;
x << yr, yg, yb, yw;
double Luminance = 1.0;
Matrix3d ccm_colorOptimized, ccm_brightnessOptimized;
ConvertSonyPrimaries2ccm(Rx, Ry, Luminance, ccm_colorOptimized, ccm_brightnessOptimized);
// Verify the result
std::cout << "ccm_colorOptimized:\n" << ccm_colorOptimized << std::endl;
return 0;
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8RWlnZW4vRGVuc2U+CiNpbmNsdWRlIDx2ZWN0b3I+Cgp1c2luZyBuYW1lc3BhY2UgRWlnZW47CgpzdHJ1Y3QgRGlzcGxheSB7CiAgICBzdGQ6OnN0cmluZyBuYW1lOwogICAgTWF0cml4M2QgcmdiMnh5ejsKICAgIE1hdHJpeDNkIHh5ejJyZ2I7CiAgICBkb3VibGUgZ2FpbjsKICAgIFZlY3RvcjNkIG9mZnNldDsKICAgIFZlY3RvcjNkIGdhbW1hOwogICAgVmVjdG9yWGQgZ2FtbWF0ZXN0Owp9OwoKRGlzcGxheSBvMkljcygpIHsKICAgIERpc3BsYXkgZGlzcDsKICAgIFZlY3RvcjJkIHh5VygwLjMxMjcsIDAuMzI5MCk7CiAgICBWZWN0b3IyZCB4eVIoMC43MDgsIDAuMjkyKTsKICAgIFZlY3RvcjJkIHh5RygwLjE3MCwgMC43OTcpOwogICAgVmVjdG9yMmQgeHlCKDAuMTMxLCAwLjA0Nik7CgogICAgTWF0cml4M2QgcmdiMnh5ejsKICAgIHJnYjJ4eXogPDwgeHlSLCB4eUcsIHh5QiwgVmVjdG9yMmQ6Ok9uZXMoKSAtIHh5UiAtIHh5RyAtIHh5QjsKICAgIFZlY3RvcjNkIHJnYldlaSA9IHJnYjJ4eXouaW52ZXJzZSgpICogVmVjdG9yM2QoeHlXLCBWZWN0b3IzZDo6T25lcygpIC0gVmVjdG9yM2QoeHlXLnN1bSgpLCB4eVcuc3VtKCksIHh5Vy5zdW0oKSkpOwogICAgcmdiMnh5eiAqPSByZ2JXZWkuYXNEaWFnb25hbCgpOwogICAgcmdiMnh5eiAvPSByZ2IyeHl6LnJvdygxKS5zdW0oKTsKCiAgICBkb3VibGUgYmV0YSA9IDAuMDE4MTsKICAgIGRvdWJsZSBhbHBoYSA9IDEuMDk5MzsKICAgIFZlY3RvclhkIGluID0gVmVjdG9yWGQ6OkxpblNwYWNlZCg2NTUzNiwgMCwgMSk7CiAgICBWZWN0b3JYZCBpbnZHYW1tYSA9ICg0LjUgKiBpbi5hcnJheSgpKS5taW4oYmV0YSkuZXZhbCgpICsKICAgICAgICAgICAgICAgICAgICAgICAgKChhbHBoYSAqIGluLmFycmF5KCkucG93KDAuNDUpIC0gKGFscGhhIC0gMSkpLm1heCgwKSkuZXZhbCgpOwogICAgVmVjdG9yWGQgZ2FtbWEgPSAoaW4uYXJyYXkoKSA8IGJldGEpLnNlbGVjdCg0LjUgKiBpbi5hcnJheSgpLCBpbnZHYW1tYSk7CgogICAgZGlzcC5uYW1lID0gIlJlYy4yMDIwIjsKICAgIGRpc3AucmdiMnh5eiA9IHJnYjJ4eXo7CiAgICBkaXNwLnh5ejJyZ2IgPSByZ2IyeHl6LmludmVyc2UoKTsKICAgIGRpc3AuZ2FpbiA9IDIyNTsKICAgIGRpc3Aub2Zmc2V0LnNldFplcm8oKTsKICAgIGRpc3AuZ2FtbWEgPSBnYW1tYS5yZXBsaWNhdGUoMywgMSk7CiAgICBkaXNwLmdhbW1hdGVzdCA9ICgoMSAvIDQuNSkgKiBpbi5hcnJheSgpKS5taW4oYmV0YSkuZXZhbCgpICsKICAgICAgICAgICAgICAgICAgICAgKCgoMSAvIGFscGhhKSAqIGluLmFycmF5KCkgKyBiZXRhIC8gNC41KS5wb3coMSAvIDAuNDUpKS5tYXgoMCkuZXZhbCgpOwoKICAgIHJldHVybiBkaXNwOwp9Cgp2b2lkIENvbnZlcnRTb255UHJpbWFyaWVzMmNjbShjb25zdCBzdGQ6OnZlY3Rvcjxkb3VibGU+JiBSeCwgY29uc3Qgc3RkOjp2ZWN0b3I8ZG91YmxlPiYgUnksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGRvdWJsZSBMdW1pbmFuY2UsIE1hdHJpeDNkJiBjY21fY29sb3JPcHRpbWl6ZWQsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIE1hdHJpeDNkJiBjY21fYnJpZ2h0bmVzc09wdGltaXplZCkgewogICAgZG91YmxlIFh3ID0geCh4LnNpemUoKSAtIDEpICogTHVtaW5hbmNlIC8geSh5LnNpemUoKSAtIDEpOwogICAgZG91YmxlIFp3ID0gKDEgLSB4KHguc2l6ZSgpIC0gMSkgLSB5KHkuc2l6ZSgpIC0gMSkpICogTHVtaW5hbmNlIC8geSh5LnNpemUoKSAtIDEpOwogICAgZG91YmxlIFl3ID0gTHVtaW5hbmNlOwoKICAgIFZlY3RvcjNkIHogPSBWZWN0b3IzZDo6T25lcygpIC0geC5oZWFkKDMpIC0geS5oZWFkKDMpOwoKICAgIE1hdHJpeDNkIEE7CiAgICBBIDw8IHguaGVhZCgzKSwgeS5oZWFkKDMpLCB6OwogICAgVmVjdG9yM2QgQihYdywgWXcsIFp3KTsKICAgIFZlY3RvcjNkIFMgPSBBLmNvbFBpdkhvdXNlaG9sZGVyUXIoKS5zb2x2ZShCKTsKCiAgICBWZWN0b3IzZCBYID0gUy5hcnJheSgpICogeC5oZWFkKDMpLmFycmF5KCk7CiAgICBWZWN0b3IzZCBZID0gUy5hcnJheSgpICogeS5oZWFkKDMpLmFycmF5KCk7CiAgICBWZWN0b3IzZCBaID0gUy5hcnJheSgpICogei5hcnJheSgpOwoKICAgIE1hdHJpeDNkIFhZWjsKICAgIFhZWiA8PCBYLCBZLCBaOwoKICAgIGRvdWJsZSBkaXNwR2FpbiA9IEx1bWluYW5jZTsKICAgIE1hdHJpeDNkIGRpc3BYWVpub3JtID0gWFlaLmFycmF5KCkgLyBkaXNwR2FpbjsKCiAgICBWZWN0b3IzZCBjaGVja1doaXRlID0gZGlzcFhZWm5vcm0ucm93KDEpLmFycmF5KCkuc3VtKCkgKiBWZWN0b3IzZDo6T25lcygpOwogICAgaWYgKChjaGVja1doaXRlIC0gVmVjdG9yM2Q6Ok9uZXMoKSkuY3dpc2VBYnMoKS5tYXhDb2VmZigpID4gMC4wMykgewogICAgICAgIHN0ZDo6Y2VyciA8PCAiKioqIFdBUk5JTkc6IExpbmVhcml0eSBzZWVtcyBvZmYgZm9yIHdoaXRlLiBNVVNUIENIRUNLLiAqKioiIDw8IHN0ZDo6ZW5kbDsKICAgIH0KCiAgICBNYXRyaXgzZCBkaXNwMnh5eiA9IGRpc3BYWVpub3JtLmFycmF5KCkuY29sd2lzZSgpIC8gY2hlY2tXaGl0ZS5hcnJheSgpOwoKICAgIC8vIEFzc3VtaW5nIHlvdSBoYXZlIGRpc3BTb255LnJnYjJ4eXogZGVmaW5lZCBlbHNld2hlcmUKICAgIE1hdHJpeDNkIGRpc3BJY3M7IC8vIEluaXRpYWxpemUgZGlzcEljcwoKICAgIGNjbV9icmlnaHRuZXNzT3B0aW1pemVkID0gZGlzcFNvbnkucmdiMnh5ei5pbnZlcnNlKCkgKiBkaXNwSWNzLnJnYjJ4eXo7CgogICAgTWF0cml4M2QgbWF0OwogICAgbWF0IDw8IDAsIDEsIDAsIDEsIDAsIDEsIDAsIDEsIDAsIDAsIDEsIDEsIDAsIDAsIDEsIDEsIDAsIDAsIDAsIDAsIDEsIDEsIDEsIDE7CiAgICBjY21fY29sb3JPcHRpbWl6ZWQgPSBjY21fYnJpZ2h0bmVzc09wdGltaXplZC5hcnJheSgpIC8gKGNjbV9icmlnaHRuZXNzT3B0aW1pemVkICogbWF0KS5tYXhDb2VmZigpOwoKICAgIGlmIChjY21fY29sb3JPcHRpbWl6ZWQubWluQ29lZmYoKSA8IC1zdGQ6Om51bWVyaWNfbGltaXRzPGRvdWJsZT46OmVwc2lsb24oKSkgewogICAgICAgIHN0ZDo6Y2VyciA8PCAiKioqIFdBUk5JTkc6IFNvbWUgcG9pbnRzIGluIHNvdXJjZSBjYW4ndCBiZSBzaW11bGF0ZWQgb24gdGhlIGRlc3RpbmF0aW9uLiAqKioiIDw8IHN0ZDo6ZW5kbDsKICAgIH0KfQoKaW50IG1haW4oKSB7CiAgICAvLyBHaXZlbiBpbnB1dCB2ZWN0b3JzCiAgICBzdGQ6OnZlY3Rvcjxkb3VibGU+IFJ4ID0gezAuNjYwLCAwLjcyMH07CiAgICBzdGQ6OnZlY3Rvcjxkb3VibGU+IFJ5ID0gezAuMjc3LCAwLjMzN307CiAgICBzdGQ6OnZlY3Rvcjxkb3VibGU+IEd4ID0gezAuMjA4LCAwLjI2OH07CiAgICBzdGQ6OnZlY3Rvcjxkb3VibGU+IEd5ID0gezAuNjYzLCAwLjcyM307CiAgICBzdGQ6OnZlY3Rvcjxkb3VibGU+IEJ4ID0gezAuMTE3LCAwLjE3N307CiAgICBzdGQ6OnZlY3Rvcjxkb3VibGU+IEJ5ID0gezAuMDI2LCAwLjA4Nn07CiAgICBzdGQ6OnZlY3Rvcjxkb3VibGU+IFd4ID0gezAuMjkzLCAwLjMzM307CiAgICBzdGQ6OnZlY3Rvcjxkb3VibGU+IFd5ID0gezAuMzA5LCAwLjM0OX07CgogICAgLy8gQ2FsY3VsYXRlIGludGVybWVkaWF0ZSB2YWx1ZXMKICAgIGRvdWJsZSB4ciA9IFJ4WzBdICsgKChSeFsxXSAtIFJ4WzBdKSAvIDIpOwogICAgZG91YmxlIHlyID0gUnlbMF0gKyAoKFJ5WzFdIC0gUnlbMF0pIC8gMik7CiAgICBkb3VibGUgeGcgPSBHeFswXSArICgoR3hbMV0gLSBHeFswXSkgLyAyKTsKICAgIGRvdWJsZSB5ZyA9IEd5WzBdICsgKChHeVsxXSAtIEd5WzBdKSAvIDIpOwogICAgZG91YmxlIHhiID0gQnhbMF0gKyAoKEJ4WzFdIC0gQnhbMF0pIC8gMik7CiAgICBkb3VibGUgeWIgPSBCeVswXSArICgoQnlbMV0gLSBCeVswXSkgLyAyKTsKICAgIGRvdWJsZSB4dyA9IFd4WzBdICsgKChXeFsxXSAtIFd4WzBdKSAvIDIpOwogICAgZG91YmxlIHl3ID0gV3lbMF0gKyAoKFd5WzFdIC0gV3lbMF0pIC8gMik7CgogICAgLy8gQ3JlYXRlIHZlY3RvcnMgeCBhbmQgeQogICAgVmVjdG9yWGQgeCg0KSwgeSg0KTsKICAgIHggPDwgeHIsIHhnLCB4YiwgeHc7CiAgICB4IDw8IHlyLCB5ZywgeWIsIHl3OwogICAgZG91YmxlIEx1bWluYW5jZSA9IDEuMDsKCiAgICBNYXRyaXgzZCBjY21fY29sb3JPcHRpbWl6ZWQsIGNjbV9icmlnaHRuZXNzT3B0aW1pemVkOwogICAgQ29udmVydFNvbnlQcmltYXJpZXMyY2NtKFJ4LCBSeSwgTHVtaW5hbmNlLCBjY21fY29sb3JPcHRpbWl6ZWQsIGNjbV9icmlnaHRuZXNzT3B0aW1pemVkKTsKCiAgICAvLyBWZXJpZnkgdGhlIHJlc3VsdAogICAgc3RkOjpjb3V0IDw8ICJjY21fY29sb3JPcHRpbWl6ZWQ6XG4iIDw8IGNjbV9jb2xvck9wdGltaXplZCA8PCBzdGQ6OmVuZGw7CgogICAgcmV0dXJuIDA7Cn0K