git.sesse.net Git - ffmpeg/blob - libavfilter/vf_ciescope.c

   1 /*
   2  * Copyright (c) 2000 John Walker
   3  * Copyright (c) 2016 Paul B Mahol
   4  *
   5  * This file is part of FFmpeg.
   6  *
   7  * FFmpeg is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU Lesser General Public
   9  * License as published by the Free Software Foundation; either
  10  * version 2.1 of the License, or (at your option) any later version.
  11  *
  12  * FFmpeg is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15  * Lesser General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public
  18  * License along with FFmpeg; if not, write to the Free Software
  19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20  */
  21
  22 #include "libavutil/avassert.h"
  23 #include "libavutil/intreadwrite.h"
  24 #include "libavutil/opt.h"
  25 #include "libavutil/parseutils.h"
  26 #include "libavutil/pixdesc.h"
  27 #include "avfilter.h"
  28 #include "formats.h"
  29 #include "internal.h"
  30 #include "video.h"
  31
  32 enum CieSystem {
  33     XYY,
  34     UCS,
  35     LUV,
  36     NB_CIE
  37 };
  38
  39 enum ColorsSystems {
  40     NTSCsystem,
  41     EBUsystem,
  42     SMPTEsystem,
  43     SMPTE240Msystem,
  44     APPLEsystem,
  45     wRGBsystem,
  46     CIE1931system,
  47     Rec709system,
  48     Rec2020system,
  49     NB_CS
  50 };
  51
  52 typedef struct CiescopeContext {
  53     const AVClass *class;
  54     int color_system;
  55     unsigned gamuts;
  56     int size;
  57     int show_white;
  58     int correct_gamma;
  59     int cie;
  60     float intensity;
  61     float contrast;
  62     int background;
  63
  64     double log2lin[65536];
  65     double igamma;
  66     double i[3][3];
  67     double m[3][3];
  68     AVFrame *f;
  69     void (*filter)(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y);
  70 } CiescopeContext;
  71
  72 #define OFFSET(x) offsetof(CiescopeContext, x)
  73 #define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
  74
  75 static const AVOption ciescope_options[] = {
  76     { "system",     "set color system", OFFSET(color_system), AV_OPT_TYPE_INT, {.i64=Rec709system}, 0, NB_CS-1, FLAGS, "system" },
  77     {   "ntsc",       "NTSC 1953 Y'I'O' (ITU-R BT.470 System M)", 0, AV_OPT_TYPE_CONST, {.i64=NTSCsystem},     0, 0, FLAGS, "system" },
  78     {   "470m",       "NTSC 1953 Y'I'O' (ITU-R BT.470 System M)", 0, AV_OPT_TYPE_CONST, {.i64=NTSCsystem},     0, 0, FLAGS, "system" },
  79     {   "ebu",        "EBU Y'U'V' (PAL/SECAM) (ITU-R BT.470 System B, G)", 0, AV_OPT_TYPE_CONST, {.i64=EBUsystem},      0, 0, FLAGS, "system" },
  80     {   "470bg",      "EBU Y'U'V' (PAL/SECAM) (ITU-R BT.470 System B, G)", 0, AV_OPT_TYPE_CONST, {.i64=EBUsystem},      0, 0, FLAGS, "system" },
  81     {   "smpte",      "SMPTE-C RGB",            0, AV_OPT_TYPE_CONST, {.i64=SMPTEsystem},    0, 0, FLAGS, "system" },
  82     {   "240m",       "SMPTE-240M Y'PbPr",      0, AV_OPT_TYPE_CONST, {.i64=SMPTE240Msystem},0, 0, FLAGS, "system" },
  83     {   "apple",      "Apple RGB",              0, AV_OPT_TYPE_CONST, {.i64=APPLEsystem},    0, 0, FLAGS, "system" },
  84     {   "widergb",    "Adobe Wide Gamut RGB",   0, AV_OPT_TYPE_CONST, {.i64=wRGBsystem},     0, 0, FLAGS, "system" },
  85     {   "cie1931",    "CIE 1931 RGB",           0, AV_OPT_TYPE_CONST, {.i64=CIE1931system},  0, 0, FLAGS, "system" },
  86     {   "hdtv",       "ITU.BT-709 Y'CbCr",      0, AV_OPT_TYPE_CONST, {.i64=Rec709system},   0, 0, FLAGS, "system" },
  87     {   "rec709",     "ITU.BT-709 Y'CbCr",      0, AV_OPT_TYPE_CONST, {.i64=Rec709system},   0, 0, FLAGS, "system" },
  88     {   "uhdtv",      "ITU-R.BT-2020",          0, AV_OPT_TYPE_CONST, {.i64=Rec2020system},  0, 0, FLAGS, "system" },
  89     {   "rec2020",    "ITU-R.BT-2020",          0, AV_OPT_TYPE_CONST, {.i64=Rec2020system},  0, 0, FLAGS, "system" },
  90     { "cie",        "set cie system", OFFSET(cie), AV_OPT_TYPE_INT,   {.i64=XYY}, 0, NB_CIE-1, FLAGS, "cie" },
  91     {   "xyy",      "CIE 1931 xyY", 0, AV_OPT_TYPE_CONST, {.i64=XYY}, 0, 0, FLAGS, "cie" },
  92     {   "ucs",      "CIE 1960 UCS", 0, AV_OPT_TYPE_CONST, {.i64=UCS}, 0, 0, FLAGS, "cie" },
  93     {   "luv",      "CIE 1976 Luv", 0, AV_OPT_TYPE_CONST, {.i64=LUV}, 0, 0, FLAGS, "cie" },
  94     { "gamuts",     "set what gamuts to draw", OFFSET(gamuts), AV_OPT_TYPE_FLAGS, {.i64=0}, 0, 0xFFF, FLAGS, "gamuts" },
  95     {   "ntsc",     NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<NTSCsystem},      0, 0, FLAGS, "gamuts" },
  96     {   "470m",     NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<NTSCsystem},      0, 0, FLAGS, "gamuts" },
  97     {   "ebu",      NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<EBUsystem},       0, 0, FLAGS, "gamuts" },
  98     {   "470bg",    NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<EBUsystem},       0, 0, FLAGS, "gamuts" },
  99     {   "smpte",    NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<SMPTEsystem},     0, 0, FLAGS, "gamuts" },
 100     {   "240m",     NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<SMPTE240Msystem}, 0, 0, FLAGS, "gamuts" },
 101     {   "apple",    NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<APPLEsystem},     0, 0, FLAGS, "gamuts" },
 102     {   "widergb",  NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<wRGBsystem},      0, 0, FLAGS, "gamuts" },
 103     {   "cie1931",  NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<CIE1931system},   0, 0, FLAGS, "gamuts" },
 104     {   "hdtv",     NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<Rec709system},    0, 0, FLAGS, "gamuts" },
 105     {   "rec709",   NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<Rec709system},    0, 0, FLAGS, "gamuts" },
 106     {   "uhdtv",    NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<Rec2020system},   0, 0, FLAGS, "gamuts" },
 107     {   "rec2020",  NULL, 0, AV_OPT_TYPE_CONST, {.i64=1<<Rec2020system},   0, 0, FLAGS, "gamuts" },
 108     { "size",       "set ciescope size", OFFSET(size), AV_OPT_TYPE_INT, {.i64=512}, 256, 8192, FLAGS },
 109     { "s",          "set ciescope size", OFFSET(size), AV_OPT_TYPE_INT, {.i64=512}, 256, 8192, FLAGS },
 110     { "intensity",  "set ciescope intensity", OFFSET(intensity), AV_OPT_TYPE_FLOAT, {.dbl=0.001}, 0, 1, FLAGS },
 111     { "i",          "set ciescope intensity", OFFSET(intensity), AV_OPT_TYPE_FLOAT, {.dbl=0.001}, 0, 1, FLAGS },
 112     { "contrast",   NULL, OFFSET(contrast), AV_OPT_TYPE_FLOAT, {.dbl=0.75},  0, 1, FLAGS },
 113     { "corrgamma",  NULL, OFFSET(correct_gamma), AV_OPT_TYPE_BOOL, {.i64=1}, 0, 1, FLAGS },
 114     { "showwhite",  NULL, OFFSET(show_white), AV_OPT_TYPE_BOOL, {.i64=0}, 0, 1, FLAGS },
 115     { "gamma",      NULL, OFFSET(igamma), AV_OPT_TYPE_DOUBLE, {.dbl=2.6}, 0.1, 6, FLAGS },
 116     { NULL }
 117 };
 118
 119 AVFILTER_DEFINE_CLASS(ciescope);
 120
 121 static const enum AVPixelFormat in_pix_fmts[] = {
 122     AV_PIX_FMT_RGB24,
 123     AV_PIX_FMT_RGBA,
 124     AV_PIX_FMT_RGB48,
 125     AV_PIX_FMT_RGBA64,
 126     AV_PIX_FMT_XYZ12,
 127     AV_PIX_FMT_NONE
 128 };
 129
 130 static const enum AVPixelFormat out_pix_fmts[] = {
 131     AV_PIX_FMT_RGBA64,
 132     AV_PIX_FMT_NONE
 133 };
 134
 135 static int query_formats(AVFilterContext *ctx)
 136 {
 137     int ret;
 138
 139     if ((ret = ff_formats_ref(ff_make_format_list(in_pix_fmts), &ctx->inputs[0]->out_formats)) < 0)
 140         return ret;
 141
 142     if ((ret = ff_formats_ref(ff_make_format_list(out_pix_fmts), &ctx->outputs[0]->in_formats)) < 0)
 143         return ret;
 144
 145     return 0;
 146 }
 147
 148 static int config_output(AVFilterLink *outlink)
 149 {
 150     CiescopeContext *s = outlink->src->priv;
 151
 152     outlink->h = outlink->w = s->size;
 153     outlink->sample_aspect_ratio = (AVRational){1,1};
 154
 155     return 0;
 156 }
 157
 158 /* A  color  system is defined by the CIE x and y  coordinates of its
 159    three primary illuminants and the x and y coordinates of the  white
 160    point. */
 161
 162 struct ColorSystem {
 163     double xRed, yRed,                /* Red primary illuminant */
 164            xGreen, yGreen,            /* Green primary illuminant */
 165            xBlue, yBlue,              /* Blue primary illuminant */
 166            xWhite, yWhite,            /* White point */
 167            gamma;             /* gamma of nonlinear correction */
 168 };
 169
 170 static float const spectral_chromaticity[][3] = {
 171     { 0.175560, 0.005294, 0.819146 },
 172     { 0.175483, 0.005286, 0.819231 },
 173     { 0.175400, 0.005279, 0.819321 },
 174     { 0.175317, 0.005271, 0.819412 },
 175     { 0.175237, 0.005263, 0.819500 },
 176     { 0.175161, 0.005256, 0.819582 },
 177     { 0.175088, 0.005247, 0.819665 },
 178     { 0.175015, 0.005236, 0.819749 },
 179     { 0.174945, 0.005226, 0.819829 },
 180     { 0.174880, 0.005221, 0.819899 },
 181     { 0.174821, 0.005221, 0.819959 },
 182     { 0.174770, 0.005229, 0.820001 },
 183     { 0.174722, 0.005238, 0.820040 },
 184     { 0.174665, 0.005236, 0.820098 },
 185     { 0.174595, 0.005218, 0.820187 },
 186     { 0.174510, 0.005182, 0.820309 },
 187     { 0.174409, 0.005127, 0.820464 },
 188     { 0.174308, 0.005068, 0.820624 },
 189     { 0.174222, 0.005017, 0.820761 },
 190     { 0.174156, 0.004981, 0.820863 },
 191     { 0.174112, 0.004964, 0.820924 },
 192     { 0.174088, 0.004964, 0.820948 },
 193     { 0.174073, 0.004973, 0.820955 },
 194     { 0.174057, 0.004982, 0.820961 },
 195     { 0.174036, 0.004986, 0.820978 },
 196     { 0.174008, 0.004981, 0.821012 },
 197     { 0.173972, 0.004964, 0.821064 },
 198     { 0.173932, 0.004943, 0.821125 },
 199     { 0.173889, 0.004926, 0.821185 },
 200     { 0.173845, 0.004916, 0.821239 },
 201     { 0.173801, 0.004915, 0.821284 },
 202     { 0.173754, 0.004925, 0.821321 },
 203     { 0.173705, 0.004937, 0.821358 },
 204     { 0.173655, 0.004944, 0.821401 },
 205     { 0.173606, 0.004940, 0.821454 },
 206     { 0.173560, 0.004923, 0.821517 },
 207     { 0.173514, 0.004895, 0.821590 },
 208     { 0.173468, 0.004865, 0.821667 },
 209     { 0.173424, 0.004836, 0.821740 },
 210     { 0.173380, 0.004813, 0.821807 },
 211     { 0.173337, 0.004797, 0.821866 },
 212     { 0.173291, 0.004786, 0.821923 },
 213     { 0.173238, 0.004779, 0.821983 },
 214     { 0.173174, 0.004775, 0.822051 },
 215     { 0.173101, 0.004774, 0.822125 },
 216     { 0.173021, 0.004775, 0.822204 },
 217     { 0.172934, 0.004781, 0.822285 },
 218     { 0.172843, 0.004791, 0.822366 },
 219     { 0.172751, 0.004799, 0.822450 },
 220     { 0.172662, 0.004802, 0.822536 },
 221     { 0.172577, 0.004799, 0.822624 },
 222     { 0.172489, 0.004795, 0.822715 },
 223     { 0.172396, 0.004796, 0.822808 },
 224     { 0.172296, 0.004803, 0.822901 },
 225     { 0.172192, 0.004815, 0.822993 },
 226     { 0.172087, 0.004833, 0.823081 },
 227     { 0.171982, 0.004855, 0.823163 },
 228     { 0.171871, 0.004889, 0.823240 },
 229     { 0.171741, 0.004939, 0.823319 },
 230     { 0.171587, 0.005010, 0.823402 },
 231     { 0.171407, 0.005102, 0.823490 },
 232     { 0.171206, 0.005211, 0.823583 },
 233     { 0.170993, 0.005334, 0.823674 },
 234     { 0.170771, 0.005470, 0.823759 },
 235     { 0.170541, 0.005621, 0.823838 },
 236     { 0.170301, 0.005789, 0.823911 },
 237     { 0.170050, 0.005974, 0.823976 },
 238     { 0.169786, 0.006177, 0.824037 },
 239     { 0.169505, 0.006398, 0.824097 },
 240     { 0.169203, 0.006639, 0.824158 },
 241     { 0.168878, 0.006900, 0.824222 },
 242     { 0.168525, 0.007184, 0.824291 },
 243     { 0.168146, 0.007491, 0.824363 },
 244     { 0.167746, 0.007821, 0.824433 },
 245     { 0.167328, 0.008175, 0.824496 },
 246     { 0.166895, 0.008556, 0.824549 },
 247     { 0.166446, 0.008964, 0.824589 },
 248     { 0.165977, 0.009402, 0.824622 },
 249     { 0.165483, 0.009865, 0.824652 },
 250     { 0.164963, 0.010351, 0.824687 },
 251     { 0.164412, 0.010858, 0.824731 },
 252     { 0.163828, 0.011385, 0.824787 },
 253     { 0.163210, 0.011937, 0.824853 },
 254     { 0.162552, 0.012520, 0.824928 },
 255     { 0.161851, 0.013137, 0.825011 },
 256     { 0.161105, 0.013793, 0.825102 },
 257     { 0.160310, 0.014491, 0.825199 },
 258     { 0.159466, 0.015232, 0.825302 },
 259     { 0.158573, 0.016015, 0.825412 },
 260     { 0.157631, 0.016840, 0.825529 },
 261     { 0.156641, 0.017705, 0.825654 },
 262     { 0.155605, 0.018609, 0.825786 },
 263     { 0.154525, 0.019556, 0.825920 },
 264     { 0.153397, 0.020554, 0.826049 },
 265     { 0.152219, 0.021612, 0.826169 },
 266     { 0.150985, 0.022740, 0.826274 },
 267     { 0.149691, 0.023950, 0.826359 },
 268     { 0.148337, 0.025247, 0.826416 },
 269     { 0.146928, 0.026635, 0.826437 },
 270     { 0.145468, 0.028118, 0.826413 },
 271     { 0.143960, 0.029703, 0.826337 },
 272     { 0.142405, 0.031394, 0.826201 },
 273     { 0.140796, 0.033213, 0.825991 },
 274     { 0.139121, 0.035201, 0.825679 },
 275     { 0.137364, 0.037403, 0.825233 },
 276     { 0.135503, 0.039879, 0.824618 },
 277     { 0.133509, 0.042692, 0.823798 },
 278     { 0.131371, 0.045876, 0.822753 },
 279     { 0.129086, 0.049450, 0.821464 },
 280     { 0.126662, 0.053426, 0.819912 },
 281     { 0.124118, 0.057803, 0.818079 },
 282     { 0.121469, 0.062588, 0.815944 },
 283     { 0.118701, 0.067830, 0.813468 },
 284     { 0.115807, 0.073581, 0.810612 },
 285     { 0.112776, 0.079896, 0.807328 },
 286     { 0.109594, 0.086843, 0.803563 },
 287     { 0.106261, 0.094486, 0.799253 },
 288     { 0.102776, 0.102864, 0.794360 },
 289     { 0.099128, 0.112007, 0.788865 },
 290     { 0.095304, 0.121945, 0.782751 },
 291     { 0.091294, 0.132702, 0.776004 },
 292     { 0.087082, 0.144317, 0.768601 },
 293     { 0.082680, 0.156866, 0.760455 },
 294     { 0.078116, 0.170420, 0.751464 },
 295     { 0.073437, 0.185032, 0.741531 },
 296     { 0.068706, 0.200723, 0.730571 },
 297     { 0.063993, 0.217468, 0.718539 },
 298     { 0.059316, 0.235254, 0.705430 },
 299     { 0.054667, 0.254096, 0.691238 },
 300     { 0.050031, 0.274002, 0.675967 },
 301     { 0.045391, 0.294976, 0.659633 },
 302     { 0.040757, 0.316981, 0.642262 },
 303     { 0.036195, 0.339900, 0.623905 },
 304     { 0.031756, 0.363598, 0.604646 },
 305     { 0.027494, 0.387921, 0.584584 },
 306     { 0.023460, 0.412703, 0.563837 },
 307     { 0.019705, 0.437756, 0.542539 },
 308     { 0.016268, 0.462955, 0.520777 },
 309     { 0.013183, 0.488207, 0.498610 },
 310     { 0.010476, 0.513404, 0.476120 },
 311     { 0.008168, 0.538423, 0.453409 },
 312     { 0.006285, 0.563068, 0.430647 },
 313     { 0.004875, 0.587116, 0.408008 },
 314     { 0.003982, 0.610447, 0.385570 },
 315     { 0.003636, 0.633011, 0.363352 },
 316     { 0.003859, 0.654823, 0.341318 },
 317     { 0.004646, 0.675898, 0.319456 },
 318     { 0.006011, 0.696120, 0.297869 },
 319     { 0.007988, 0.715342, 0.276670 },
 320     { 0.010603, 0.733413, 0.255984 },
 321     { 0.013870, 0.750186, 0.235943 },
 322     { 0.017766, 0.765612, 0.216622 },
 323     { 0.022244, 0.779630, 0.198126 },
 324     { 0.027273, 0.792104, 0.180623 },
 325     { 0.032820, 0.802926, 0.164254 },
 326     { 0.038852, 0.812016, 0.149132 },
 327     { 0.045328, 0.819391, 0.135281 },
 328     { 0.052177, 0.825164, 0.122660 },
 329     { 0.059326, 0.829426, 0.111249 },
 330     { 0.066716, 0.832274, 0.101010 },
 331     { 0.074302, 0.833803, 0.091894 },
 332     { 0.082053, 0.834090, 0.083856 },
 333     { 0.089942, 0.833289, 0.076769 },
 334     { 0.097940, 0.831593, 0.070468 },
 335     { 0.106021, 0.829178, 0.064801 },
 336     { 0.114161, 0.826207, 0.059632 },
 337     { 0.122347, 0.822770, 0.054882 },
 338     { 0.130546, 0.818928, 0.050526 },
 339     { 0.138702, 0.814774, 0.046523 },
 340     { 0.146773, 0.810395, 0.042832 },
 341     { 0.154722, 0.805864, 0.039414 },
 342     { 0.162535, 0.801238, 0.036226 },
 343     { 0.170237, 0.796519, 0.033244 },
 344     { 0.177850, 0.791687, 0.030464 },
 345     { 0.185391, 0.786728, 0.027881 },
 346     { 0.192876, 0.781629, 0.025495 },
 347     { 0.200309, 0.776399, 0.023292 },
 348     { 0.207690, 0.771055, 0.021255 },
 349     { 0.215030, 0.765595, 0.019375 },
 350     { 0.222337, 0.760020, 0.017643 },
 351     { 0.229620, 0.754329, 0.016051 },
 352     { 0.236885, 0.748524, 0.014591 },
 353     { 0.244133, 0.742614, 0.013253 },
 354     { 0.251363, 0.736606, 0.012031 },
 355     { 0.258578, 0.730507, 0.010916 },
 356     { 0.265775, 0.724324, 0.009901 },
 357     { 0.272958, 0.718062, 0.008980 },
 358     { 0.280129, 0.711725, 0.008146 },
 359     { 0.287292, 0.705316, 0.007391 },
 360     { 0.294450, 0.698842, 0.006708 },
 361     { 0.301604, 0.692308, 0.006088 },
 362     { 0.308760, 0.685712, 0.005528 },
 363     { 0.315914, 0.679063, 0.005022 },
 364     { 0.323066, 0.672367, 0.004566 },
 365     { 0.330216, 0.665628, 0.004156 },
 366     { 0.337363, 0.658848, 0.003788 },
 367     { 0.344513, 0.652028, 0.003459 },
 368     { 0.351664, 0.645172, 0.003163 },
 369     { 0.358814, 0.638287, 0.002899 },
 370     { 0.365959, 0.631379, 0.002662 },
 371     { 0.373102, 0.624451, 0.002448 },
 372     { 0.380244, 0.617502, 0.002254 },
 373     { 0.387379, 0.610542, 0.002079 },
 374     { 0.394507, 0.603571, 0.001922 },
 375     { 0.401626, 0.596592, 0.001782 },
 376     { 0.408736, 0.589607, 0.001657 },
 377     { 0.415836, 0.582618, 0.001546 },
 378     { 0.422921, 0.575631, 0.001448 },
 379     { 0.429989, 0.568649, 0.001362 },
 380     { 0.437036, 0.561676, 0.001288 },
 381     { 0.444062, 0.554714, 0.001224 },
 382     { 0.451065, 0.547766, 0.001169 },
 383     { 0.458041, 0.540837, 0.001123 },
 384     { 0.464986, 0.533930, 0.001084 },
 385     { 0.471899, 0.527051, 0.001051 },
 386     { 0.478775, 0.520202, 0.001023 },
 387     { 0.485612, 0.513389, 0.001000 },
 388     { 0.492405, 0.506615, 0.000980 },
 389     { 0.499151, 0.499887, 0.000962 },
 390     { 0.505845, 0.493211, 0.000944 },
 391     { 0.512486, 0.486591, 0.000923 },
 392     { 0.519073, 0.480029, 0.000899 },
 393     { 0.525600, 0.473527, 0.000872 },
 394     { 0.532066, 0.467091, 0.000843 },
 395     { 0.538463, 0.460725, 0.000812 },
 396     { 0.544787, 0.454434, 0.000779 },
 397     { 0.551031, 0.448225, 0.000744 },
 398     { 0.557193, 0.442099, 0.000708 },
 399     { 0.563269, 0.436058, 0.000673 },
 400     { 0.569257, 0.430102, 0.000641 },
 401     { 0.575151, 0.424232, 0.000616 },
 402     { 0.580953, 0.418447, 0.000601 },
 403     { 0.586650, 0.412758, 0.000591 },
 404     { 0.592225, 0.407190, 0.000586 },
 405     { 0.597658, 0.401762, 0.000580 },
 406     { 0.602933, 0.396497, 0.000571 },
 407     { 0.608035, 0.391409, 0.000556 },
 408     { 0.612977, 0.386486, 0.000537 },
 409     { 0.617779, 0.381706, 0.000516 },
 410     { 0.622459, 0.377047, 0.000493 },
 411     { 0.627037, 0.372491, 0.000472 },
 412     { 0.631521, 0.368026, 0.000453 },
 413     { 0.635900, 0.363665, 0.000435 },
 414     { 0.640156, 0.359428, 0.000416 },
 415     { 0.644273, 0.355331, 0.000396 },
 416     { 0.648233, 0.351395, 0.000372 },
 417     { 0.652028, 0.347628, 0.000344 },
 418     { 0.655669, 0.344018, 0.000313 },
 419     { 0.659166, 0.340553, 0.000281 },
 420     { 0.662528, 0.337221, 0.000251 },
 421     { 0.665764, 0.334011, 0.000226 },
 422     { 0.668874, 0.330919, 0.000207 },
 423     { 0.671859, 0.327947, 0.000194 },
 424     { 0.674720, 0.325095, 0.000185 },
 425     { 0.677459, 0.322362, 0.000179 },
 426     { 0.680079, 0.319747, 0.000174 },
 427     { 0.682582, 0.317249, 0.000170 },
 428     { 0.684971, 0.314863, 0.000167 },
 429     { 0.687250, 0.312586, 0.000164 },
 430     { 0.689426, 0.310414, 0.000160 },
 431     { 0.691504, 0.308342, 0.000154 },
 432     { 0.693490, 0.306366, 0.000145 },
 433     { 0.695389, 0.304479, 0.000133 },
 434     { 0.697206, 0.302675, 0.000119 },
 435     { 0.698944, 0.300950, 0.000106 },
 436     { 0.700606, 0.299301, 0.000093 },
 437     { 0.702193, 0.297725, 0.000083 },
 438     { 0.703709, 0.296217, 0.000074 },
 439     { 0.705163, 0.294770, 0.000067 },
 440     { 0.706563, 0.293376, 0.000061 },
 441     { 0.707918, 0.292027, 0.000055 },
 442     { 0.709231, 0.290719, 0.000050 },
 443     { 0.710500, 0.289453, 0.000047 },
 444     { 0.711724, 0.288232, 0.000044 },
 445     { 0.712901, 0.287057, 0.000041 },
 446     { 0.714032, 0.285929, 0.000040 },
 447     { 0.715117, 0.284845, 0.000038 },
 448     { 0.716159, 0.283804, 0.000036 },
 449     { 0.717159, 0.282806, 0.000035 },
 450     { 0.718116, 0.281850, 0.000034 },
 451     { 0.719033, 0.280935, 0.000032 },
 452     { 0.719912, 0.280058, 0.000030 },
 453     { 0.720753, 0.279219, 0.000028 },
 454     { 0.721555, 0.278420, 0.000026 },
 455     { 0.722315, 0.277662, 0.000023 },
 456     { 0.723032, 0.276948, 0.000020 },
 457     { 0.723702, 0.276282, 0.000016 },
 458     { 0.724328, 0.275660, 0.000012 },
 459     { 0.724914, 0.275078, 0.000007 },
 460     { 0.725467, 0.274530, 0.000003 },
 461     { 0.725992, 0.274008, 0.000000 },
 462     { 0.726495, 0.273505, 0.000000 },
 463     { 0.726975, 0.273025, 0.000000 },
 464     { 0.727432, 0.272568, 0.000000 },
 465     { 0.727864, 0.272136, 0.000000 },
 466     { 0.728272, 0.271728, 0.000000 },
 467     { 0.728656, 0.271344, 0.000000 },
 468     { 0.729020, 0.270980, 0.000000 },
 469     { 0.729361, 0.270639, 0.000000 },
 470     { 0.729678, 0.270322, 0.000000 },
 471     { 0.729969, 0.270031, 0.000000 },
 472     { 0.730234, 0.269766, 0.000000 },
 473     { 0.730474, 0.269526, 0.000000 },
 474     { 0.730693, 0.269307, 0.000000 },
 475     { 0.730896, 0.269104, 0.000000 },
 476     { 0.731089, 0.268911, 0.000000 },
 477     { 0.731280, 0.268720, 0.000000 },
 478     { 0.731467, 0.268533, 0.000000 },
 479     { 0.731650, 0.268350, 0.000000 },
 480     { 0.731826, 0.268174, 0.000000 },
 481     { 0.731993, 0.268007, 0.000000 },
 482     { 0.732150, 0.267850, 0.000000 },
 483     { 0.732300, 0.267700, 0.000000 },
 484     { 0.732443, 0.267557, 0.000000 },
 485     { 0.732581, 0.267419, 0.000000 },
 486     { 0.732719, 0.267281, 0.000000 },
 487     { 0.732859, 0.267141, 0.000000 },
 488     { 0.733000, 0.267000, 0.000000 },
 489     { 0.733142, 0.266858, 0.000000 },
 490     { 0.733281, 0.266719, 0.000000 },
 491     { 0.733417, 0.266583, 0.000000 },
 492     { 0.733551, 0.266449, 0.000000 },
 493     { 0.733683, 0.266317, 0.000000 },
 494     { 0.733813, 0.266187, 0.000000 },
 495     { 0.733936, 0.266064, 0.000000 },
 496     { 0.734047, 0.265953, 0.000000 },
 497     { 0.734143, 0.265857, 0.000000 },
 498     { 0.734221, 0.265779, 0.000000 },
 499     { 0.734286, 0.265714, 0.000000 },
 500     { 0.734341, 0.265659, 0.000000 },
 501     { 0.734390, 0.265610, 0.000000 },
 502     { 0.734438, 0.265562, 0.000000 },
 503     { 0.734482, 0.265518, 0.000000 },
 504     { 0.734523, 0.265477, 0.000000 },
 505     { 0.734560, 0.265440, 0.000000 },
 506     { 0.734592, 0.265408, 0.000000 },
 507     { 0.734621, 0.265379, 0.000000 },
 508     { 0.734649, 0.265351, 0.000000 },
 509     { 0.734673, 0.265327, 0.000000 },
 510     { 0.734690, 0.265310, 0.000000 },
 511     { 0.734690, 0.265310, 0.000000 },
 512     { 0.734690, 0.265310, 0.000000 },
 513     { 0.734690, 0.265310, 0.000000 },
 514     { 0.734690, 0.265310, 0.000000 },
 515     { 0.734690, 0.265310, 0.000000 },
 516     { 0.734690, 0.265310, 0.000000 },
 517     { 0.734690, 0.265310, 0.000000 },
 518     { 0.734690, 0.265310, 0.000000 },
 519     { 0.734690, 0.265310, 0.000000 },
 520     { 0.734690, 0.265310, 0.000000 },
 521     { 0.734690, 0.265310, 0.000000 },
 522     { 0.734690, 0.265310, 0.000000 },
 523     { 0.734690, 0.265310, 0.000000 },
 524     { 0.734690, 0.265310, 0.000000 },
 525     { 0.734690, 0.265310, 0.000000 },
 526     { 0.734690, 0.265310, 0.000000 },
 527     { 0.734690, 0.265310, 0.000000 },
 528     { 0.734690, 0.265310, 0.000000 },
 529     { 0.734690, 0.265310, 0.000000 },
 530     { 0.734690, 0.265310, 0.000000 },
 531     { 0.734690, 0.265310, 0.000000 },
 532     { 0.734690, 0.265310, 0.000000 },
 533     { 0.734690, 0.265310, 0.000000 },
 534     { 0.734690, 0.265310, 0.000000 },
 535     { 0.734690, 0.265310, 0.000000 },
 536     { 0.734690, 0.265310, 0.000000 },
 537     { 0.734690, 0.265310, 0.000000 },
 538     { 0.734690, 0.265310, 0.000000 },
 539     { 0.734690, 0.265310, 0.000000 },
 540     { 0.734690, 0.265310, 0.000000 },
 541     { 0.734690, 0.265310, 0.000000 },
 542     { 0.734690, 0.265310, 0.000000 },
 543     { 0.734690, 0.265310, 0.000000 },
 544     { 0.734690, 0.265310, 0.000000 },
 545     { 0.734690, 0.265310, 0.000000 },
 546     { 0.734690, 0.265310, 0.000000 },
 547     { 0.734690, 0.265310, 0.000000 },
 548     { 0.734690, 0.265310, 0.000000 },
 549     { 0.734690, 0.265310, 0.000000 },
 550     { 0.734690, 0.265310, 0.000000 },
 551     { 0.734690, 0.265310, 0.000000 },
 552     { 0.734690, 0.265310, 0.000000 },
 553     { 0.734690, 0.265310, 0.000000 },
 554     { 0.734690, 0.265310, 0.000000 },
 555     { 0.734690, 0.265310, 0.000000 },
 556     { 0.734690, 0.265310, 0.000000 },
 557     { 0.734690, 0.265310, 0.000000 },
 558     { 0.734690, 0.265310, 0.000000 },
 559     { 0.734690, 0.265310, 0.000000 },
 560     { 0.734690, 0.265310, 0.000000 },
 561     { 0.734690, 0.265310, 0.000000 },
 562     { 0.734690, 0.265310, 0.000000 },
 563     { 0.734690, 0.265310, 0.000000 },
 564     { 0.734690, 0.265310, 0.000000 },
 565     { 0.734690, 0.265310, 0.000000 },
 566     { 0.734690, 0.265310, 0.000000 },
 567     { 0.734690, 0.265310, 0.000000 },
 568     { 0.734690, 0.265310, 0.000000 },
 569     { 0.734690, 0.265310, 0.000000 },
 570     { 0.734690, 0.265310, 0.000000 },
 571     { 0.734690, 0.265310, 0.000000 },
 572     { 0.734690, 0.265310, 0.000000 },
 573     { 0.734690, 0.265310, 0.000000 },
 574     { 0.734690, 0.265310, 0.000000 },
 575     { 0.734690, 0.265310, 0.000000 },
 576     { 0.734690, 0.265310, 0.000000 },
 577     { 0.734690, 0.265310, 0.000000 },
 578     { 0.734690, 0.265310, 0.000000 },
 579     { 0.734690, 0.265310, 0.000000 },
 580     { 0.734690, 0.265310, 0.000000 },
 581     { 0.734690, 0.265310, 0.000000 },
 582     { 0.734690, 0.265310, 0.000000 },
 583     { 0.734690, 0.265310, 0.000000 },
 584     { 0.734690, 0.265310, 0.000000 },
 585     { 0.734690, 0.265310, 0.000000 },
 586     { 0.734690, 0.265310, 0.000000 },
 587     { 0.734690, 0.265310, 0.000000 },
 588     { 0.734690, 0.265310, 0.000000 },
 589     { 0.734690, 0.265310, 0.000000 },
 590     { 0.734690, 0.265310, 0.000000 },
 591     { 0.734690, 0.265310, 0.000000 },
 592     { 0.734690, 0.265310, 0.000000 },
 593     { 0.734690, 0.265310, 0.000000 },
 594     { 0.734690, 0.265310, 0.000000 },
 595     { 0.734690, 0.265310, 0.000000 },
 596     { 0.734690, 0.265310, 0.000000 },
 597     { 0.734690, 0.265310, 0.000000 },
 598     { 0.734690, 0.265310, 0.000000 },
 599     { 0.734690, 0.265310, 0.000000 },
 600     { 0.734690, 0.265310, 0.000000 },
 601     { 0.734690, 0.265310, 0.000000 },
 602     { 0.734690, 0.265310, 0.000000 },
 603     { 0.734690, 0.265310, 0.000000 },
 604     { 0.734690, 0.265310, 0.000000 },
 605     { 0.734690, 0.265310, 0.000000 },
 606     { 0.734690, 0.265310, 0.000000 },
 607     { 0.734690, 0.265310, 0.000000 },
 608     { 0.734690, 0.265310, 0.000000 },
 609     { 0.734690, 0.265310, 0.000000 },
 610     { 0.734690, 0.265310, 0.000000 },
 611     { 0.734690, 0.265310, 0.000000 },
 612     { 0.734690, 0.265310, 0.000000 },
 613     { 0.734690, 0.265310, 0.000000 },
 614     { 0.734690, 0.265310, 0.000000 },
 615     { 0.734690, 0.265310, 0.000000 },
 616     { 0.734690, 0.265310, 0.000000 },
 617     { 0.734690, 0.265310, 0.000000 },
 618     { 0.734690, 0.265310, 0.000000 },
 619     { 0.734690, 0.265310, 0.000000 },
 620     { 0.734690, 0.265310, 0.000000 },
 621     { 0.734690, 0.265310, 0.000000 },
 622     { 0.734690, 0.265310, 0.000000 },
 623     { 0.734690, 0.265310, 0.000000 },
 624     { 0.734690, 0.265310, 0.000000 },
 625     { 0.734690, 0.265310, 0.000000 },
 626     { 0.734690, 0.265310, 0.000000 },
 627     { 0.734690, 0.265310, 0.000000 },
 628     { 0.734690, 0.265310, 0.000000 },
 629     { 0.734690, 0.265310, 0.000000 },
 630     { 0.734690, 0.265310, 0.000000 },
 631     { 0.734690, 0.265310, 0.000000 },
 632     { 0.734690, 0.265310, 0.000000 },
 633     { 0.734690, 0.265310, 0.000000 },
 634     { 0.734690, 0.265310, 0.000000 },
 635     { 0.734690, 0.265310, 0.000000 },
 636     { 0.734690, 0.265310, 0.000000 },
 637     { 0.734690, 0.265310, 0.000000 },
 638     { 0.734690, 0.265310, 0.000000 },
 639     { 0.734690, 0.265310, 0.000000 },
 640     { 0.734690, 0.265310, 0.000000 },
 641     { 0.734690, 0.265310, 0.000000 },
 642 };
 643
 644
 645 /* Standard white point chromaticities. */
 646
 647 #define C     0.310063, 0.316158
 648 #define E     1.0/3.0, 1.0/3.0
 649 #define D50   0.34570, 0.3585
 650 #define D65   0.312713, 0.329016
 651
 652 /* Gamma of nonlinear correction.
 653    See Charles Poynton's ColorFAQ Item 45 and GammaFAQ Item 6 at
 654    http://www.inforamp.net/~poynton/ColorFAQ.html
 655    http://www.inforamp.net/~poynton/GammaFAQ.html
 656 */
 657
 658 #define GAMMA_REC709    0.      /* Rec. 709 */
 659
 660 static const struct ColorSystem color_systems[] = {
 661     [NTSCsystem] = {
 662         0.67,  0.33,  0.21,  0.71,  0.14,  0.08,
 663         C, GAMMA_REC709
 664     },
 665     [EBUsystem] = {
 666         0.64,  0.33,  0.29,  0.60,  0.15,  0.06,
 667         D65, GAMMA_REC709
 668     },
 669     [SMPTEsystem] = {
 670         0.630, 0.340, 0.310, 0.595, 0.155, 0.070,
 671         D65, GAMMA_REC709
 672     },
 673     [SMPTE240Msystem] = {
 674         0.670, 0.330, 0.210, 0.710, 0.150, 0.060,
 675         D65, GAMMA_REC709
 676     },
 677     [APPLEsystem] = {
 678         0.625, 0.340, 0.280, 0.595, 0.115, 0.070,
 679         D65, GAMMA_REC709
 680     },
 681     [wRGBsystem] = {
 682         0.7347, 0.2653, 0.1152, 0.8264, 0.1566, 0.0177,
 683         D50, GAMMA_REC709
 684     },
 685     [CIE1931system] = {
 686         0.7347, 0.2653, 0.2738, 0.7174, 0.1666, 0.0089,
 687         E, GAMMA_REC709
 688     },
 689     [Rec709system] = {
 690         0.64,  0.33,  0.30,  0.60,  0.15,  0.06,
 691         D65, GAMMA_REC709
 692     },
 693     [Rec2020system] = {
 694         0.708,  0.292,  0.170,  0.797,  0.131,  0.046,
 695         D65, GAMMA_REC709
 696     },
 697 };
 698
 699 /*
 700 static struct ColorSystem CustomSystem = {
 701     "Custom",
 702     0.64,  0.33,  0.30,  0.60,  0.15,  0.06,
 703     D65, GAMMA_REC709
 704 };
 705 */
 706
 707 static void
 708 uv_to_xy(double   const u,
 709          double   const v,
 710          double * const xc,
 711          double * const yc)
 712 {
 713 /*
 714     Given 1970 coordinates u, v, determine 1931 chromaticities x, y
 715 */
 716     *xc = 3*u / (2*u - 8*v + 4);
 717     *yc = 2*v / (2*u - 8*v + 4);
 718 }
 719
 720 static void
 721 upvp_to_xy(double   const up,
 722            double   const vp,
 723            double * const xc,
 724            double * const yc)
 725 {
 726 /*
 727     Given 1976 coordinates u', v', determine 1931 chromaticities x, y
 728 */
 729     *xc = 9*up / (6*up - 16*vp + 12);
 730     *yc = 4*vp / (6*up - 16*vp + 12);
 731 }
 732
 733 static void
 734 xy_to_upvp(double xc,
 735            double yc,
 736            double * const up,
 737            double * const vp)
 738 {
 739 /*
 740     Given 1931 chromaticities x, y, determine 1976 coordinates u', v'
 741 */
 742     *up = 4*xc / (- 2*xc + 12*yc + 3);
 743     *vp = 9*yc / (- 2*xc + 12*yc + 3);
 744 }
 745
 746 static void
 747 xy_to_uv(double xc,
 748          double yc,
 749          double * const u,
 750          double * const v)
 751 {
 752 /*
 753     Given 1931 chromaticities x, y, determine 1960 coordinates u, v
 754 */
 755     *u = 4*xc / (- 2*xc + 12*yc + 3);
 756     *v = 6*yc / (- 2*xc + 12*yc + 3);
 757 }
 758
 759 static void
 760 xyz_to_rgb(const double m[3][3],
 761            double xc, double yc, double zc,
 762            double * const r, double * const g, double * const b)
 763 {
 764     *r = m[0][0]*xc + m[0][1]*yc + m[0][2]*zc;
 765     *g = m[1][0]*xc + m[1][1]*yc + m[1][2]*zc;
 766     *b = m[2][0]*xc + m[2][1]*yc + m[2][2]*zc;
 767 }
 768
 769 static void invert_matrix3x3(double in[3][3], double out[3][3])
 770 {
 771     double m00 = in[0][0], m01 = in[0][1], m02 = in[0][2],
 772            m10 = in[1][0], m11 = in[1][1], m12 = in[1][2],
 773            m20 = in[2][0], m21 = in[2][1], m22 = in[2][2];
 774     int i, j;
 775     double det;
 776
 777     out[0][0] =  (m11 * m22 - m21 * m12);
 778     out[0][1] = -(m01 * m22 - m21 * m02);
 779     out[0][2] =  (m01 * m12 - m11 * m02);
 780     out[1][0] = -(m10 * m22 - m20 * m12);
 781     out[1][1] =  (m00 * m22 - m20 * m02);
 782     out[1][2] = -(m00 * m12 - m10 * m02);
 783     out[2][0] =  (m10 * m21 - m20 * m11);
 784     out[2][1] = -(m00 * m21 - m20 * m01);
 785     out[2][2] =  (m00 * m11 - m10 * m01);
 786
 787     det = m00 * out[0][0] + m10 * out[0][1] + m20 * out[0][2];
 788     det = 1.0 / det;
 789
 790     for (i = 0; i < 3; i++) {
 791         for (j = 0; j < 3; j++)
 792             out[i][j] *= det;
 793     }
 794 }
 795
 796 static void get_rgb2xyz_matrix(struct ColorSystem system, double m[3][3])
 797 {
 798     double S[3], X[4], Z[4];
 799     int i;
 800
 801     X[0] = system.xRed   / system.yRed;
 802     X[1] = system.xGreen / system.yGreen;
 803     X[2] = system.xBlue  / system.yBlue;
 804     X[3] = system.xWhite / system.yWhite;
 805
 806     Z[0] = (1 - system.xRed   - system.yRed)   / system.yRed;
 807     Z[1] = (1 - system.xGreen - system.yGreen) / system.yGreen;
 808     Z[2] = (1 - system.xBlue  - system.yBlue)  / system.yBlue;
 809     Z[3] = (1 - system.xWhite - system.yWhite) / system.yWhite;
 810
 811     for (i = 0; i < 3; i++) {
 812         m[0][i] = X[i];
 813         m[1][i] = 1;
 814         m[2][i] = Z[i];
 815     }
 816
 817     invert_matrix3x3(m, m);
 818
 819     for (i = 0; i < 3; i++)
 820         S[i] = m[i][0] * X[3] + m[i][1] * 1 + m[i][2] * Z[3];
 821
 822     for (i = 0; i < 3; i++) {
 823         m[0][i] = S[i] * X[i];
 824         m[1][i] = S[i] * 1;
 825         m[2][i] = S[i] * Z[i];
 826     }
 827 }
 828
 829 static void
 830 rgb_to_xy(double rc,
 831           double gc,
 832           double bc,
 833           double * const x,
 834           double * const y,
 835           double * const z,
 836           const double m[3][3])
 837 {
 838     double sum;
 839
 840     *x = m[0][0] * rc + m[0][1] * gc + m[0][2] * bc;
 841     *y = m[1][0] * rc + m[1][1] * gc + m[1][2] * bc;
 842     *z = m[2][0] * rc + m[2][1] * gc + m[2][2] * bc;
 843
 844     sum = *x + *y + *z;
 845
 846     *x = *x / sum;
 847     *y = *y / sum;
 848 }
 849
 850 static int
 851 constrain_rgb(double * const r,
 852               double * const g,
 853               double * const b)
 854 {
 855 /*----------------------------------------------------------------------------
 856     If  the  requested RGB shade contains a negative weight for one of
 857     the primaries, it lies outside the color  gamut  accessible  from
 858     the  given  triple  of  primaries.  Desaturate it by adding white,
 859     equal quantities of R, G, and B, enough to make RGB all positive.
 860 -----------------------------------------------------------------------------*/
 861     double w;
 862
 863     /* Amount of white needed is w = - min(0, *r, *g, *b) */
 864     w = (0 < *r) ? 0 : *r;
 865     w = (w < *g) ? w : *g;
 866     w = (w < *b) ? w : *b;
 867     w = - w;
 868
 869     /* Add just enough white to make r, g, b all positive. */
 870     if (w > 0) {
 871         *r += w;  *g += w; *b += w;
 872
 873         return 1;                     /* Color modified to fit RGB gamut */
 874     }
 875
 876     return 0;                         /* Color within RGB gamut */
 877 }
 878
 879 static void
 880 gamma_correct(const struct ColorSystem * const cs,
 881               double *                   const c)
 882 {
 883 /*----------------------------------------------------------------------------
 884     Transform linear RGB values to nonlinear RGB values.
 885
 886     Rec. 709 is ITU-R Recommendation BT. 709 (1990)
 887     ``Basic Parameter Values for the HDTV Standard for the Studio and for
 888     International Programme Exchange'', formerly CCIR Rec. 709.
 889
 890     For details see
 891        http://www.inforamp.net/~poynton/ColorFAQ.html
 892        http://www.inforamp.net/~poynton/GammaFAQ.html
 893 -----------------------------------------------------------------------------*/
 894     double gamma;
 895     double cc;
 896
 897     gamma = cs->gamma;
 898
 899     if (gamma == 0.) {
 900         /* Rec. 709 gamma correction. */
 901         cc = 0.018;
 902         if (*c < cc) {
 903             *c *= (1.099 * pow(cc, 0.45) - 0.099) / cc;
 904         } else {
 905             *c = 1.099 * pow(*c, 0.45) - 0.099;
 906         }
 907     } else {
 908     /* Nonlinear color = (Linear color)^(1/gamma) */
 909         *c = pow(*c, 1./gamma);
 910     }
 911 }
 912
 913
 914
 915 static void
 916 gamma_correct_rgb(const struct ColorSystem * const cs,
 917                   double * const r,
 918                   double * const g,
 919                   double * const b)
 920 {
 921     gamma_correct(cs, r);
 922     gamma_correct(cs, g);
 923     gamma_correct(cs, b);
 924 }
 925
 926 /* Sz(X) is the displacement in pixels of a displacement of X normalized
 927    distance units.  (A normalized distance unit is 1/512 of the smaller
 928    dimension of the canvas)
 929 */
 930 #define Sz(x) (((x) * (int)FFMIN(w, h)) / 512)
 931
 932 static void
 933 monochrome_color_location(double waveLength, int w, int h,
 934                           int cie, int *xP, int *yP)
 935 {
 936     const int ix = waveLength - 360;
 937     const double pX = spectral_chromaticity[ix][0];
 938     const double pY = spectral_chromaticity[ix][1];
 939     const double pZ = spectral_chromaticity[ix][2];
 940     const double px = pX / (pX + pY + pZ);
 941     const double py = pY / (pX + pY + pZ);
 942
 943     if (cie == LUV) {
 944         double up, vp;
 945
 946         xy_to_upvp(px, py, &up, &vp);
 947         *xP = up * (w - 1);
 948         *yP = (h - 1) - vp * (h - 1);
 949     } else if (cie == UCS) {
 950         double u, v;
 951
 952         xy_to_uv(px, py, &u, &v);
 953         *xP = u * (w - 1);
 954         *yP = (h - 1) - v * (h - 1);
 955     } else if (cie == XYY) {
 956         *xP = px * (w - 1);
 957         *yP = (h - 1) - py * (h - 1);
 958     } else {
 959         av_assert0(0);
 960     }
 961 }
 962
 963 static void
 964 find_tongue(uint16_t* const pixels,
 965             int       const w,
 966             int       const linesize,
 967             int       const row,
 968             int *     const presentP,
 969             int *     const leftEdgeP,
 970             int *     const rightEdgeP)
 971 {
 972     int i;
 973
 974     for (i = 0; i < w && pixels[row * linesize + i * 4 + 0] == 0; i++)
 975         ;
 976
 977     if (i >= w) {
 978         *presentP = 0;
 979     } else {
 980         int j;
 981         int const leftEdge = i;
 982
 983         *presentP = 1;
 984
 985         for (j = w - 1; j >= leftEdge && pixels[row * linesize + j * 4 + 0] == 0; j--)
 986             ;
 987
 988         *rightEdgeP = j;
 989         *leftEdgeP = leftEdge;
 990     }
 991 }
 992
 993 static void draw_line(uint16_t *const pixels, int linesize,
 994                       int x0, int y0, int x1, int y1,
 995                       int w, int h,
 996                       const uint16_t *const rgbcolor)
 997 {
 998     int dx  = FFABS(x1 - x0), sx = x0 < x1 ? 1 : -1;
 999     int dy  = FFABS(y1 - y0), sy = y0 < y1 ? 1 : -1;
1000     int err = (dx > dy ? dx : -dy) / 2, e2;
1001
1002     for (;;) {
1003         pixels[y0 * linesize + x0 * 4 + 0] = rgbcolor[0];
1004         pixels[y0 * linesize + x0 * 4 + 1] = rgbcolor[1];
1005         pixels[y0 * linesize + x0 * 4 + 2] = rgbcolor[2];
1006         pixels[y0 * linesize + x0 * 4 + 3] = rgbcolor[3];
1007
1008         if (x0 == x1 && y0 == y1)
1009             break;
1010
1011         e2 = err;
1012
1013         if (e2 >-dx) {
1014             err -= dy;
1015             x0 += sx;
1016         }
1017
1018         if (e2 < dy) {
1019             err += dx;
1020             y0 += sy;
1021         }
1022     }
1023 }
1024
1025 static void draw_rline(uint16_t *const pixels, int linesize,
1026                        int x0, int y0, int x1, int y1,
1027                        int w, int h)
1028 {
1029     int dx  = FFABS(x1 - x0), sx = x0 < x1 ? 1 : -1;
1030     int dy  = FFABS(y1 - y0), sy = y0 < y1 ? 1 : -1;
1031     int err = (dx > dy ? dx : -dy) / 2, e2;
1032
1033     for (;;) {
1034         pixels[y0 * linesize + x0 * 4 + 0] = 65535 - pixels[y0 * linesize + x0 * 4 + 0];
1035         pixels[y0 * linesize + x0 * 4 + 1] = 65535 - pixels[y0 * linesize + x0 * 4 + 1];
1036         pixels[y0 * linesize + x0 * 4 + 2] = 65535 - pixels[y0 * linesize + x0 * 4 + 2];
1037         pixels[y0 * linesize + x0 * 4 + 3] = 65535;
1038
1039         if (x0 == x1 && y0 == y1)
1040             break;
1041
1042         e2 = err;
1043
1044         if (e2 >-dx) {
1045             err -= dy;
1046             x0 += sx;
1047         }
1048
1049         if (e2 < dy) {
1050             err += dx;
1051             y0 += sy;
1052         }
1053     }
1054 }
1055
1056 static void
1057 tongue_outline(uint16_t* const pixels,
1058                int       const linesize,
1059                int       const w,
1060                int       const h,
1061                uint16_t  const maxval,
1062                int       const cie)
1063 {
1064     const uint16_t rgbcolor[4] = { maxval, maxval, maxval, maxval };
1065     int wavelength;
1066     int lx, ly;
1067     int fx, fy;
1068
1069     for (wavelength = 360; wavelength <= 830; wavelength++) {
1070         int icx, icy;
1071
1072         monochrome_color_location(wavelength, w, h, cie,
1073                                   &icx, &icy);
1074
1075         if (wavelength > 360)
1076             draw_line(pixels, linesize, lx, ly, icx, icy, w, h, rgbcolor);
1077         else {
1078             fx = icx;
1079             fy = icy;
1080         }
1081         lx = icx;
1082         ly = icy;
1083     }
1084     draw_line(pixels, linesize, lx, ly, fx, fy, w, h, rgbcolor);
1085 }
1086
1087 static void
1088 fill_in_tongue(uint16_t*                  const pixels,
1089                int                        const linesize,
1090                int                        const w,
1091                int                        const h,
1092                uint16_t                   const maxval,
1093                const struct ColorSystem * const cs,
1094                double                     const m[3][3],
1095                int                        const cie,
1096                int                        const correct_gamma,
1097                float                      const contrast)
1098 {
1099     int y;
1100
1101     /* Scan the image line by line and  fill  the  tongue  outline
1102        with the RGB values determined by the color system for the x-y
1103        co-ordinates within the tongue.
1104     */
1105
1106     for (y = 0; y < h; ++y) {
1107         int  present;  /* There is some tongue on this line */
1108         int leftEdge; /* x position of leftmost pixel in tongue on this line */
1109         int rightEdge; /* same, but rightmost */
1110
1111         find_tongue(pixels, w, linesize, y, &present, &leftEdge, &rightEdge);
1112
1113         if (present) {
1114             int x;
1115
1116             for (x = leftEdge; x <= rightEdge; ++x) {
1117                 double cx, cy, cz, jr, jg, jb, jmax;
1118                 int r, g, b, mx = maxval;
1119
1120                 if (cie == LUV) {
1121                     double up, vp;
1122                     up = ((double) x) / (w - 1);
1123                     vp = 1.0 - ((double) y) / (h - 1);
1124                     upvp_to_xy(up, vp, &cx, &cy);
1125                     cz = 1.0 - (cx + cy);
1126                 } else if (cie == UCS) {
1127                     double u, v;
1128                     u = ((double) x) / (w - 1);
1129                     v = 1.0 - ((double) y) / (h - 1);
1130                     uv_to_xy(u, v, &cx, &cy);
1131                     cz = 1.0 - (cx + cy);
1132                 } else if (cie == XYY) {
1133                     cx = ((double) x) / (w - 1);
1134                     cy = 1.0 - ((double) y) / (h - 1);
1135                     cz = 1.0 - (cx + cy);
1136                 } else {
1137                     av_assert0(0);
1138                 }
1139
1140                 xyz_to_rgb(m, cx, cy, cz, &jr, &jg, &jb);
1141
1142                 /* Check whether the requested color  is  within  the
1143                    gamut  achievable with the given color system.  If
1144                    not, draw it in a reduced  intensity,  interpolated
1145                    by desaturation to the closest within-gamut color. */
1146
1147                 if (constrain_rgb(&jr, &jg, &jb))
1148                     mx *= contrast;
1149
1150                 jmax = FFMAX3(jr, jg, jb);
1151                 if (jmax > 0) {
1152                     jr = jr / jmax;
1153                     jg = jg / jmax;
1154                     jb = jb / jmax;
1155                 }
1156                 /* gamma correct from linear rgb to nonlinear rgb. */
1157                 if (correct_gamma)
1158                     gamma_correct_rgb(cs, &jr, &jg, &jb);
1159                 r = mx * jr;
1160                 g = mx * jg;
1161                 b = mx * jb;
1162                 pixels[y * linesize + x * 4 + 0] = r;
1163                 pixels[y * linesize + x * 4 + 1] = g;
1164                 pixels[y * linesize + x * 4 + 2] = b;
1165                 pixels[y * linesize + x * 4 + 3] = 65535;
1166             }
1167         }
1168     }
1169 }
1170
1171 static void
1172 plot_white_point(uint16_t*      pixels,
1173                  int      const linesize,
1174                  int      const w,
1175                  int      const h,
1176                  int      const maxval,
1177                  int      const color_system,
1178                  int      const cie)
1179 {
1180     const struct ColorSystem *cs = &color_systems[color_system];
1181     int wx, wy;
1182
1183     if (cie == LUV) {
1184         double wup, wvp;
1185         xy_to_upvp(cs->xWhite, cs->yWhite, &wup, &wvp);
1186         wx = wup;
1187         wy = wvp;
1188         wx = (w - 1) * wup;
1189         wy = (h - 1) - ((int) ((h - 1) * wvp));
1190     } else if (cie == UCS) {
1191         double wu, wv;
1192         xy_to_uv(cs->xWhite, cs->yWhite, &wu, &wv);
1193         wx = wu;
1194         wy = wv;
1195         wx = (w - 1) * wu;
1196         wy = (h - 1) - ((int) ((h - 1) * wv));
1197     } else if (cie == XYY) {
1198         wx = (w - 1) * cs->xWhite;
1199         wy = (h - 1) - ((int) ((h - 1) * cs->yWhite));
1200     } else {
1201         av_assert0(0);
1202     }
1203
1204     draw_rline(pixels, linesize,
1205                wx + Sz(3), wy, wx + Sz(10), wy,
1206                w, h);
1207     draw_rline(pixels, linesize,
1208                wx - Sz(3), wy, wx - Sz(10), wy,
1209                w, h);
1210     draw_rline(pixels, linesize,
1211                wx, wy + Sz(3), wx, wy + Sz(10),
1212                w, h);
1213     draw_rline(pixels, linesize,
1214                wx, wy - Sz(3), wx, wy - Sz(10),
1215                w, h);
1216 }
1217
1218 static int draw_background(AVFilterContext *ctx)
1219 {
1220     CiescopeContext *s = ctx->priv;
1221     const struct ColorSystem *cs = &color_systems[s->color_system];
1222     AVFilterLink *outlink = ctx->outputs[0];
1223     int w = s->size;
1224     int h = s->size;
1225     uint16_t *pixels;
1226
1227     if ((s->f = ff_get_video_buffer(outlink, outlink->w, outlink->h)) == NULL)
1228         return AVERROR(ENOMEM);
1229     pixels = (uint16_t *)s->f->data[0];
1230
1231     tongue_outline(pixels, s->f->linesize[0] / 2, w, h, 65535, s->cie);
1232
1233     fill_in_tongue(pixels, s->f->linesize[0] / 2, w, h, 65535, cs, (const double (*)[3])s->i, s->cie,
1234                    s->correct_gamma, s->contrast);
1235
1236     return 0;
1237 }
1238
1239 static void filter_rgb48(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1240 {
1241     CiescopeContext *s = ctx->priv;
1242     const uint16_t* src = (const uint16_t*)(in->data[0] + in->linesize[0] * y + x * 6);
1243     double r = src[0] / 65535.;
1244     double g = src[1] / 65535.;
1245     double b = src[2] / 65535.;
1246     double cz;
1247
1248     rgb_to_xy(r, g, b, cx, cy, &cz, (const double (*)[3])s->m);
1249 }
1250
1251 static void filter_rgba64(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1252 {
1253     CiescopeContext *s = ctx->priv;
1254     const uint16_t* src = (const uint16_t*)(in->data[0] + in->linesize[0] * y + x * 8);
1255     double r = src[0] / 65535.;
1256     double g = src[1] / 65535.;
1257     double b = src[2] / 65535.;
1258     double cz;
1259
1260     rgb_to_xy(r, g, b, cx, cy, &cz, (const double (*)[3])s->m);
1261 }
1262
1263 static void filter_rgb24(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1264 {
1265     CiescopeContext *s = ctx->priv;
1266     const uint8_t* src = in->data[0] + in->linesize[0] * y + x * 3;
1267     double r = src[0] / 255.;
1268     double g = src[1] / 255.;
1269     double b = src[2] / 255.;
1270     double cz;
1271
1272     rgb_to_xy(r, g, b, cx, cy, &cz, (const double (*)[3])s->m);
1273 }
1274
1275 static void filter_rgba(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1276 {
1277     CiescopeContext *s = ctx->priv;
1278     const uint8_t* src = in->data[0] + in->linesize[0] * y + x * 4;
1279     double r = src[0] / 255.;
1280     double g = src[1] / 255.;
1281     double b = src[2] / 255.;
1282     double cz;
1283
1284     rgb_to_xy(r, g, b, cx, cy, &cz, (const double (*)[3])s->m);
1285 }
1286
1287 static void filter_xyz(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1288 {
1289     CiescopeContext *s = ctx->priv;
1290     const uint16_t* src = (uint16_t *)(in->data[0] + in->linesize[0] * y + x * 6);
1291     double lx = s->log2lin[src[0]];
1292     double ly = s->log2lin[src[1]];
1293     double lz = s->log2lin[src[2]];
1294     double sum = lx + ly + lz;
1295
1296     if (sum == 0)
1297         sum = 1;
1298     *cx = lx / sum;
1299     *cy = ly / sum;
1300 }
1301
1302 static void plot_gamuts(uint16_t *pixels, int linesize, int w, int h,
1303                         int cie, int gamuts)
1304 {
1305     int i;
1306
1307     for (i = 0; i < NB_CS; i++) {
1308         const struct ColorSystem *cs = &color_systems[i];
1309         int rx, ry, gx, gy, bx, by;
1310
1311         if (!((1 << i) & gamuts))
1312             continue;
1313         if (cie == LUV) {
1314             double wup, wvp;
1315             xy_to_upvp(cs->xRed, cs->yRed, &wup, &wvp);
1316             rx = (w - 1) * wup;
1317             ry = (h - 1) - ((int) ((h - 1) * wvp));
1318             xy_to_upvp(cs->xGreen, cs->yGreen, &wup, &wvp);
1319             gx = (w - 1) * wup;
1320             gy = (h - 1) - ((int) ((h - 1) * wvp));
1321             xy_to_upvp(cs->xBlue, cs->yBlue, &wup, &wvp);
1322             bx = (w - 1) * wup;
1323             by = (h - 1) - ((int) ((h - 1) * wvp));
1324         } else if (cie == UCS) {
1325             double wu, wv;
1326             xy_to_uv(cs->xRed, cs->yRed, &wu, &wv);
1327             rx = (w - 1) * wu;
1328             ry = (h - 1) - ((int) ((h - 1) * wv));
1329             xy_to_uv(cs->xGreen, cs->yGreen, &wu, &wv);
1330             gx = (w - 1) * wu;
1331             gy = (h - 1) - ((int) ((h - 1) * wv));
1332             xy_to_uv(cs->xBlue, cs->yBlue, &wu, &wv);
1333             bx = (w - 1) * wu;
1334             by = (h - 1) - ((int) ((h - 1) * wv));
1335         } else if (cie == XYY) {
1336             rx = (w - 1) * cs->xRed;
1337             ry = (h - 1) - ((int) ((h - 1) * cs->yRed));
1338             gx = (w - 1) * cs->xGreen;
1339             gy = (h - 1) - ((int) ((h - 1) * cs->yGreen));
1340             bx = (w - 1) * cs->xBlue;
1341             by = (h - 1) - ((int) ((h - 1) * cs->yBlue));
1342         } else {
1343             av_assert0(0);
1344         }
1345
1346         draw_rline(pixels, linesize, rx, ry, gx, gy, w, h);
1347         draw_rline(pixels, linesize, gx, gy, bx, by, w, h);
1348         draw_rline(pixels, linesize, bx, by, rx, ry, w, h);
1349     }
1350 }
1351
1352 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
1353 {
1354     AVFilterContext *ctx  = inlink->dst;
1355     CiescopeContext *s = ctx->priv;
1356     AVFilterLink *outlink = ctx->outputs[0];
1357     int i = s->intensity * 65535;
1358     int w = outlink->w;
1359     int h = outlink->h;
1360     AVFrame *out;
1361     int ret, x, y;
1362
1363     out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
1364     if (!out) {
1365         av_frame_free(&in);
1366         return AVERROR(ENOMEM);
1367     }
1368     out->pts = in->pts;
1369
1370     if (!s->background) {
1371         ret = draw_background(ctx);
1372         if (ret < 0) {
1373             av_frame_free(&out);
1374             return ret;
1375         }
1376         s->background = 1;
1377     }
1378     for (y = 0; y < outlink->h; y++) {
1379         memset(out->data[0] + y * out->linesize[0], 0, outlink->w * 8);
1380     }
1381
1382     for (y = 0; y < in->height; y++) {
1383         for (x = 0; x < in->width; x++) {
1384             double cx, cy;
1385             uint16_t *dst;
1386             int wx, wy;
1387
1388             s->filter(ctx, in, &cx, &cy, x, y);
1389
1390             if (s->cie == LUV) {
1391                 double up, vp;
1392                 xy_to_upvp(cx, cy, &up, &vp);
1393                 cx = up;
1394                 cy = vp;
1395             } else if (s->cie == UCS) {
1396                 double u, v;
1397                 xy_to_uv(cx, cy, &u, &v);
1398                 cx = u;
1399                 cy = v;
1400             }
1401
1402             wx = (w - 1) * cx;
1403             wy = (h - 1) - ((h - 1) * cy);
1404
1405             if (wx < 0 || wx >= w ||
1406                 wy < 0 || wy >= h)
1407                 continue;
1408
1409             dst = (uint16_t *)(out->data[0] + wy * out->linesize[0] + wx * 8 + 0);
1410             dst[0] = FFMIN(dst[0] + i, 65535);
1411             dst[1] = FFMIN(dst[1] + i, 65535);
1412             dst[2] = FFMIN(dst[2] + i, 65535);
1413             dst[3] = 65535;
1414         }
1415     }
1416
1417     for (y = 0; y < outlink->h; y++) {
1418         uint16_t *dst = (uint16_t *)(out->data[0] + y * out->linesize[0]);
1419         const uint16_t *src = (const uint16_t *)(s->f->data[0] + y * s->f->linesize[0]);
1420         for (x = 0; x < outlink->w; x++) {
1421             const int xx = x * 4;
1422             if (dst[xx + 3] == 0) {
1423                 dst[xx + 0] = src[xx + 0];
1424                 dst[xx + 1] = src[xx + 1];
1425                 dst[xx + 2] = src[xx + 2];
1426                 dst[xx + 3] = src[xx + 3];
1427             }
1428         }
1429     }
1430
1431     if (s->show_white)
1432         plot_white_point((uint16_t *)out->data[0], out->linesize[0] / 2,
1433                          outlink->w, outlink->h, 65535,
1434                          s->color_system, s->cie);
1435
1436     plot_gamuts((uint16_t *)out->data[0], out->linesize[0] / 2,
1437                 outlink->w, outlink->h,
1438                 s->cie, s->gamuts);
1439
1440     av_frame_free(&in);
1441     return ff_filter_frame(outlink, out);
1442 }
1443
1444 static void av_cold uninit(AVFilterContext *ctx)
1445 {
1446     CiescopeContext *s = ctx->priv;
1447
1448     av_frame_free(&s->f);
1449 }
1450
1451 static int config_input(AVFilterLink *inlink)
1452 {
1453     CiescopeContext *s = inlink->dst->priv;
1454     int i;
1455
1456     get_rgb2xyz_matrix(color_systems[s->color_system], s->m);
1457     invert_matrix3x3(s->m, s->i);
1458
1459     switch (inlink->format) {
1460     case AV_PIX_FMT_RGB24:
1461         s->filter = filter_rgb24;
1462         break;
1463     case AV_PIX_FMT_RGBA:
1464         s->filter = filter_rgba;
1465         break;
1466     case AV_PIX_FMT_RGB48:
1467         s->filter = filter_rgb48;
1468         break;
1469     case AV_PIX_FMT_RGBA64:
1470         s->filter = filter_rgba64;
1471         break;
1472     case AV_PIX_FMT_XYZ12:
1473         s->filter = filter_xyz;
1474         for (i = 0; i < 65536; i++)
1475             s->log2lin[i] = pow(i / 65535., s->igamma) * 65535.;
1476         break;
1477     default:
1478         av_assert0(0);
1479     }
1480
1481     return 0;
1482 }
1483
1484 static const AVFilterPad inputs[] = {
1485     {
1486         .name         = "default",
1487         .type         = AVMEDIA_TYPE_VIDEO,
1488         .filter_frame = filter_frame,
1489         .config_props = config_input,
1490     },
1491     { NULL }
1492 };
1493
1494 static const AVFilterPad outputs[] = {
1495     {
1496         .name         = "default",
1497         .type         = AVMEDIA_TYPE_VIDEO,
1498         .config_props = config_output,
1499     },
1500     { NULL }
1501 };
1502
1503 AVFilter ff_vf_ciescope = {
1504     .name          = "ciescope",
1505     .description   = NULL_IF_CONFIG_SMALL("Video CIE scope."),
1506     .priv_size     = sizeof(CiescopeContext),
1507     .priv_class    = &ciescope_class,
1508     .query_formats = query_formats,
1509     .uninit        = uninit,
1510     .inputs        = inputs,
1511     .outputs       = outputs,
1512 };