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
 776     out[0][0] =  (m11 * m22 - m21 * m12);
 777     out[0][1] = -(m01 * m22 - m21 * m02);
 778     out[0][2] =  (m01 * m12 - m11 * m02);
 779     out[1][0] = -(m10 * m22 - m20 * m12);
 780     out[1][1] =  (m00 * m22 - m20 * m02);
 781     out[1][2] = -(m00 * m12 - m10 * m02);
 782     out[2][0] =  (m10 * m21 - m20 * m11);
 783     out[2][1] = -(m00 * m21 - m20 * m01);
 784     out[2][2] =  (m00 * m11 - m10 * m01);
 785
 786     double det = m00 * out[0][0] + m10 * out[0][1] + m20 * out[0][2];
 787     det = 1.0 / det;
 788
 789     for (i = 0; i < 3; i++) {
 790         for (j = 0; j < 3; j++)
 791             out[i][j] *= det;
 792     }
 793 }
 794
 795 static void get_rgb2xyz_matrix(struct ColorSystem system, double m[3][3])
 796 {
 797     double S[3], X[4], Z[4];
 798     int i;
 799
 800     X[0] = system.xRed   / system.yRed;
 801     X[1] = system.xGreen / system.yGreen;
 802     X[2] = system.xBlue  / system.yBlue;
 803     X[3] = system.xWhite / system.yWhite;
 804
 805     Z[0] = (1 - system.xRed   - system.yRed)   / system.yRed;
 806     Z[1] = (1 - system.xGreen - system.yGreen) / system.yGreen;
 807     Z[2] = (1 - system.xBlue  - system.yBlue)  / system.yBlue;
 808     Z[3] = (1 - system.xWhite - system.yWhite) / system.yWhite;
 809
 810     for (i = 0; i < 3; i++) {
 811         m[0][i] = X[i];
 812         m[1][i] = 1;
 813         m[2][i] = Z[i];
 814     }
 815
 816     invert_matrix3x3(m, m);
 817
 818     for (i = 0; i < 3; i++)
 819         S[i] = m[i][0] * X[3] + m[i][1] * 1 + m[i][2] * Z[3];
 820
 821     for (i = 0; i < 3; i++) {
 822         m[0][i] = S[i] * X[i];
 823         m[1][i] = S[i] * 1;
 824         m[2][i] = S[i] * Z[i];
 825     }
 826 }
 827
 828 static void
 829 rgb_to_xy(double rc,
 830           double gc,
 831           double bc,
 832           double * const x,
 833           double * const y,
 834           double * const z,
 835           const double m[3][3])
 836 {
 837     double sum;
 838
 839     *x = m[0][0] * rc + m[0][1] * gc + m[0][2] * bc;
 840     *y = m[1][0] * rc + m[1][1] * gc + m[1][2] * bc;
 841     *z = m[2][0] * rc + m[2][1] * gc + m[2][2] * bc;
 842
 843     sum = *x + *y + *z;
 844
 845     *x = *x / sum;
 846     *y = *y / sum;
 847 }
 848
 849 static int
 850 constrain_rgb(double * const r,
 851               double * const g,
 852               double * const b)
 853 {
 854 /*----------------------------------------------------------------------------
 855     If  the  requested RGB shade contains a negative weight for one of
 856     the primaries, it lies outside the color  gamut  accessible  from
 857     the  given  triple  of  primaries.  Desaturate it by adding white,
 858     equal quantities of R, G, and B, enough to make RGB all positive.
 859 -----------------------------------------------------------------------------*/
 860     double w;
 861
 862     /* Amount of white needed is w = - min(0, *r, *g, *b) */
 863     w = (0 < *r) ? 0 : *r;
 864     w = (w < *g) ? w : *g;
 865     w = (w < *b) ? w : *b;
 866     w = - w;
 867
 868     /* Add just enough white to make r, g, b all positive. */
 869     if (w > 0) {
 870         *r += w;  *g += w; *b += w;
 871
 872         return 1;                     /* Color modified to fit RGB gamut */
 873     }
 874
 875     return 0;                         /* Color within RGB gamut */
 876 }
 877
 878 static void
 879 gamma_correct(const struct ColorSystem * const cs,
 880               double *                   const c)
 881 {
 882 /*----------------------------------------------------------------------------
 883     Transform linear RGB values to nonlinear RGB values.
 884
 885     Rec. 709 is ITU-R Recommendation BT. 709 (1990)
 886     ``Basic Parameter Values for the HDTV Standard for the Studio and for
 887     International Programme Exchange'', formerly CCIR Rec. 709.
 888
 889     For details see
 890        http://www.inforamp.net/~poynton/ColorFAQ.html
 891        http://www.inforamp.net/~poynton/GammaFAQ.html
 892 -----------------------------------------------------------------------------*/
 893     double gamma;
 894     double cc;
 895
 896     gamma = cs->gamma;
 897
 898     if (gamma == 0.) {
 899         /* Rec. 709 gamma correction. */
 900         cc = 0.018;
 901         if (*c < cc) {
 902             *c *= (1.099 * pow(cc, 0.45) - 0.099) / cc;
 903         } else {
 904             *c = 1.099 * pow(*c, 0.45) - 0.099;
 905         }
 906     } else {
 907     /* Nonlinear color = (Linear color)^(1/gamma) */
 908         *c = pow(*c, 1./gamma);
 909     }
 910 }
 911
 912
 913
 914 static void
 915 gamma_correct_rgb(const struct ColorSystem * const cs,
 916                   double * const r,
 917                   double * const g,
 918                   double * const b)
 919 {
 920     gamma_correct(cs, r);
 921     gamma_correct(cs, g);
 922     gamma_correct(cs, b);
 923 }
 924
 925 /* Sz(X) is the displacement in pixels of a displacement of X normalized
 926    distance units.  (A normalized distance unit is 1/512 of the smaller
 927    dimension of the canvas)
 928 */
 929 #define Sz(x) (((x) * (int)FFMIN(w, h)) / 512)
 930
 931 static void
 932 monochrome_color_location(double waveLength, int w, int h,
 933                           int cie, int *xP, int *yP)
 934 {
 935     const int ix = waveLength - 360;
 936     const double pX = spectral_chromaticity[ix][0];
 937     const double pY = spectral_chromaticity[ix][1];
 938     const double pZ = spectral_chromaticity[ix][2];
 939     const double px = pX / (pX + pY + pZ);
 940     const double py = pY / (pX + pY + pZ);
 941
 942     if (cie == LUV) {
 943         double up, vp;
 944
 945         xy_to_upvp(px, py, &up, &vp);
 946         *xP = up * (w - 1);
 947         *yP = (h - 1) - vp * (h - 1);
 948     } else if (cie == UCS) {
 949         double u, v;
 950
 951         xy_to_uv(px, py, &u, &v);
 952         *xP = u * (w - 1);
 953         *yP = (h - 1) - v * (h - 1);
 954     } else if (cie == XYY) {
 955         *xP = px * (w - 1);
 956         *yP = (h - 1) - py * (h - 1);
 957     } else {
 958         av_assert0(0);
 959     }
 960 }
 961
 962 static void
 963 find_tongue(uint16_t* const pixels,
 964             int       const w,
 965             int       const linesize,
 966             int       const row,
 967             int *     const presentP,
 968             int *     const leftEdgeP,
 969             int *     const rightEdgeP)
 970 {
 971     int i;
 972
 973     for (i = 0; i < w && pixels[row * linesize + i * 4 + 0] == 0; i++)
 974         ;
 975
 976     if (i >= w) {
 977         *presentP = 0;
 978     } else {
 979         int j;
 980         int const leftEdge = i;
 981
 982         *presentP = 1;
 983
 984         for (j = w - 1; j >= leftEdge && pixels[row * linesize + j * 4 + 0] == 0; j--)
 985             ;
 986
 987         *rightEdgeP = j;
 988         *leftEdgeP = leftEdge;
 989     }
 990 }
 991
 992 static void draw_line(uint16_t *const pixels, int linesize,
 993                       int x0, int y0, int x1, int y1,
 994                       int w, int h,
 995                       const uint16_t *const rgbcolor)
 996 {
 997     int dx  = FFABS(x1 - x0), sx = x0 < x1 ? 1 : -1;
 998     int dy  = FFABS(y1 - y0), sy = y0 < y1 ? 1 : -1;
 999     int err = (dx > dy ? dx : -dy) / 2, e2;
1000
1001     for (;;) {
1002         pixels[y0 * linesize + x0 * 4 + 0] = rgbcolor[0];
1003         pixels[y0 * linesize + x0 * 4 + 1] = rgbcolor[1];
1004         pixels[y0 * linesize + x0 * 4 + 2] = rgbcolor[2];
1005         pixels[y0 * linesize + x0 * 4 + 3] = rgbcolor[3];
1006
1007         if (x0 == x1 && y0 == y1)
1008             break;
1009
1010         e2 = err;
1011
1012         if (e2 >-dx) {
1013             err -= dy;
1014             x0 += sx;
1015         }
1016
1017         if (e2 < dy) {
1018             err += dx;
1019             y0 += sy;
1020         }
1021     }
1022 }
1023
1024 static void draw_rline(uint16_t *const pixels, int linesize,
1025                        int x0, int y0, int x1, int y1,
1026                        int w, int h)
1027 {
1028     int dx  = FFABS(x1 - x0), sx = x0 < x1 ? 1 : -1;
1029     int dy  = FFABS(y1 - y0), sy = y0 < y1 ? 1 : -1;
1030     int err = (dx > dy ? dx : -dy) / 2, e2;
1031
1032     for (;;) {
1033         pixels[y0 * linesize + x0 * 4 + 0] = 65535 - pixels[y0 * linesize + x0 * 4 + 0];
1034         pixels[y0 * linesize + x0 * 4 + 1] = 65535 - pixels[y0 * linesize + x0 * 4 + 1];
1035         pixels[y0 * linesize + x0 * 4 + 2] = 65535 - pixels[y0 * linesize + x0 * 4 + 2];
1036         pixels[y0 * linesize + x0 * 4 + 3] = 65535;
1037
1038         if (x0 == x1 && y0 == y1)
1039             break;
1040
1041         e2 = err;
1042
1043         if (e2 >-dx) {
1044             err -= dy;
1045             x0 += sx;
1046         }
1047
1048         if (e2 < dy) {
1049             err += dx;
1050             y0 += sy;
1051         }
1052     }
1053 }
1054
1055 static void
1056 tongue_outline(uint16_t* const pixels,
1057                int       const linesize,
1058                int       const w,
1059                int       const h,
1060                uint16_t  const maxval,
1061                int       const cie)
1062 {
1063     const uint16_t rgbcolor[4] = { maxval, maxval, maxval, maxval };
1064     int wavelength;
1065     int lx, ly;
1066     int fx, fy;
1067
1068     for (wavelength = 360; wavelength <= 830; wavelength++) {
1069         int icx, icy;
1070
1071         monochrome_color_location(wavelength, w, h, cie,
1072                                   &icx, &icy);
1073
1074         if (wavelength > 360)
1075             draw_line(pixels, linesize, lx, ly, icx, icy, w, h, rgbcolor);
1076         else {
1077             fx = icx;
1078             fy = icy;
1079         }
1080         lx = icx;
1081         ly = icy;
1082     }
1083     draw_line(pixels, linesize, lx, ly, fx, fy, w, h, rgbcolor);
1084 }
1085
1086 static void
1087 fill_in_tongue(uint16_t*                  const pixels,
1088                int                        const linesize,
1089                int                        const w,
1090                int                        const h,
1091                uint16_t                   const maxval,
1092                const struct ColorSystem * const cs,
1093                double                     const m[3][3],
1094                int                        const cie,
1095                int                        const correct_gamma,
1096                float                      const contrast)
1097 {
1098     int y;
1099
1100     /* Scan the image line by line and  fill  the  tongue  outline
1101        with the RGB values determined by the color system for the x-y
1102        co-ordinates within the tongue.
1103     */
1104
1105     for (y = 0; y < h; ++y) {
1106         int  present;  /* There is some tongue on this line */
1107         int leftEdge; /* x position of leftmost pixel in tongue on this line */
1108         int rightEdge; /* same, but rightmost */
1109
1110         find_tongue(pixels, w, linesize, y, &present, &leftEdge, &rightEdge);
1111
1112         if (present) {
1113             int x;
1114
1115             for (x = leftEdge; x <= rightEdge; ++x) {
1116                 double cx, cy, cz, jr, jg, jb, jmax;
1117                 int r, g, b, mx = maxval;
1118
1119                 if (cie == LUV) {
1120                     double up, vp;
1121                     up = ((double) x) / (w - 1);
1122                     vp = 1.0 - ((double) y) / (h - 1);
1123                     upvp_to_xy(up, vp, &cx, &cy);
1124                     cz = 1.0 - (cx + cy);
1125                 } else if (cie == UCS) {
1126                     double u, v;
1127                     u = ((double) x) / (w - 1);
1128                     v = 1.0 - ((double) y) / (h - 1);
1129                     uv_to_xy(u, v, &cx, &cy);
1130                     cz = 1.0 - (cx + cy);
1131                 } else if (cie == XYY) {
1132                     cx = ((double) x) / (w - 1);
1133                     cy = 1.0 - ((double) y) / (h - 1);
1134                     cz = 1.0 - (cx + cy);
1135                 } else {
1136                     av_assert0(0);
1137                 }
1138
1139                 xyz_to_rgb(m, cx, cy, cz, &jr, &jg, &jb);
1140
1141                 /* Check whether the requested color  is  within  the
1142                    gamut  achievable with the given color system.  If
1143                    not, draw it in a reduced  intensity,  interpolated
1144                    by desaturation to the closest within-gamut color. */
1145
1146                 if (constrain_rgb(&jr, &jg, &jb))
1147                     mx *= contrast;
1148
1149                 jmax = FFMAX3(jr, jg, jb);
1150                 if (jmax > 0) {
1151                     jr = jr / jmax;
1152                     jg = jg / jmax;
1153                     jb = jb / jmax;
1154                 }
1155                 /* gamma correct from linear rgb to nonlinear rgb. */
1156                 if (correct_gamma)
1157                     gamma_correct_rgb(cs, &jr, &jg, &jb);
1158                 r = mx * jr;
1159                 g = mx * jg;
1160                 b = mx * jb;
1161                 pixels[y * linesize + x * 4 + 0] = r;
1162                 pixels[y * linesize + x * 4 + 1] = g;
1163                 pixels[y * linesize + x * 4 + 2] = b;
1164                 pixels[y * linesize + x * 4 + 3] = 65535;
1165             }
1166         }
1167     }
1168 }
1169
1170 static void
1171 plot_white_point(uint16_t*      pixels,
1172                  int      const linesize,
1173                  int      const w,
1174                  int      const h,
1175                  int      const maxval,
1176                  int      const color_system,
1177                  int      const cie)
1178 {
1179     const struct ColorSystem *cs = &color_systems[color_system];
1180     int wx, wy;
1181
1182     if (cie == LUV) {
1183         double wup, wvp;
1184         xy_to_upvp(cs->xWhite, cs->yWhite, &wup, &wvp);
1185         wx = wup;
1186         wy = wvp;
1187         wx = (w - 1) * wup;
1188         wy = (h - 1) - ((int) ((h - 1) * wvp));
1189     } else if (cie == UCS) {
1190         double wu, wv;
1191         xy_to_uv(cs->xWhite, cs->yWhite, &wu, &wv);
1192         wx = wu;
1193         wy = wv;
1194         wx = (w - 1) * wu;
1195         wy = (h - 1) - ((int) ((h - 1) * wv));
1196     } else if (cie == XYY) {
1197         wx = (w - 1) * cs->xWhite;
1198         wy = (h - 1) - ((int) ((h - 1) * cs->yWhite));
1199     } else {
1200         av_assert0(0);
1201     }
1202
1203     draw_rline(pixels, linesize,
1204                wx + Sz(3), wy, wx + Sz(10), wy,
1205                w, h);
1206     draw_rline(pixels, linesize,
1207                wx - Sz(3), wy, wx - Sz(10), wy,
1208                w, h);
1209     draw_rline(pixels, linesize,
1210                wx, wy + Sz(3), wx, wy + Sz(10),
1211                w, h);
1212     draw_rline(pixels, linesize,
1213                wx, wy - Sz(3), wx, wy - Sz(10),
1214                w, h);
1215 }
1216
1217 static int draw_background(AVFilterContext *ctx)
1218 {
1219     CiescopeContext *s = ctx->priv;
1220     const struct ColorSystem *cs = &color_systems[s->color_system];
1221     AVFilterLink *outlink = ctx->outputs[0];
1222     int w = s->size;
1223     int h = s->size;
1224     uint16_t *pixels;
1225
1226     if ((s->f = ff_get_video_buffer(outlink, outlink->w, outlink->h)) == NULL)
1227         return AVERROR(ENOMEM);
1228     pixels = (uint16_t *)s->f->data[0];
1229
1230     tongue_outline(pixels, s->f->linesize[0] / 2, w, h, 65535, s->cie);
1231
1232     fill_in_tongue(pixels, s->f->linesize[0] / 2, w, h, 65535, cs, s->i, s->cie,
1233                    s->correct_gamma, s->contrast);
1234
1235     return 0;
1236 }
1237
1238 static void filter_rgb48(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1239 {
1240     CiescopeContext *s = ctx->priv;
1241     const uint16_t* src = (const uint16_t*)(in->data[0] + in->linesize[0] * y + x * 6);
1242     double r = src[0] / 65535.;
1243     double g = src[1] / 65535.;
1244     double b = src[2] / 65535.;
1245     double cz;
1246
1247     rgb_to_xy(r, g, b, cx, cy, &cz, s->m);
1248 }
1249
1250 static void filter_rgba64(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1251 {
1252     CiescopeContext *s = ctx->priv;
1253     const uint16_t* src = (const uint16_t*)(in->data[0] + in->linesize[0] * y + x * 8);
1254     double r = src[0] / 65535.;
1255     double g = src[1] / 65535.;
1256     double b = src[2] / 65535.;
1257     double cz;
1258
1259     rgb_to_xy(r, g, b, cx, cy, &cz, s->m);
1260 }
1261
1262 static void filter_rgb24(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1263 {
1264     CiescopeContext *s = ctx->priv;
1265     const uint8_t* src = in->data[0] + in->linesize[0] * y + x * 3;
1266     double r = src[0] / 255.;
1267     double g = src[1] / 255.;
1268     double b = src[2] / 255.;
1269     double cz;
1270
1271     rgb_to_xy(r, g, b, cx, cy, &cz, s->m);
1272 }
1273
1274 static void filter_rgba(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1275 {
1276     CiescopeContext *s = ctx->priv;
1277     const uint8_t* src = in->data[0] + in->linesize[0] * y + x * 4;
1278     double r = src[0] / 255.;
1279     double g = src[1] / 255.;
1280     double b = src[2] / 255.;
1281     double cz;
1282
1283     rgb_to_xy(r, g, b, cx, cy, &cz, s->m);
1284 }
1285
1286 static void filter_xyz(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1287 {
1288     CiescopeContext *s = ctx->priv;
1289     const uint16_t* src = (uint16_t *)(in->data[0] + in->linesize[0] * y + x * 6);
1290     double lx = s->log2lin[src[0]];
1291     double ly = s->log2lin[src[1]];
1292     double lz = s->log2lin[src[2]];
1293     double sum = lx + ly + lz;
1294
1295     if (sum == 0)
1296         sum = 1;
1297     *cx = lx / sum;
1298     *cy = ly / sum;
1299 }
1300
1301 static void plot_gamuts(uint16_t *pixels, int linesize, int w, int h,
1302                         int cie, int gamuts)
1303 {
1304     int i;
1305
1306     for (i = 0; i < NB_CS; i++) {
1307         const struct ColorSystem *cs = &color_systems[i];
1308         int rx, ry, gx, gy, bx, by;
1309
1310         if (!((1 << i) & gamuts))
1311             continue;
1312         if (cie == LUV) {
1313             double wup, wvp;
1314             xy_to_upvp(cs->xRed, cs->yRed, &wup, &wvp);
1315             rx = (w - 1) * wup;
1316             ry = (h - 1) - ((int) ((h - 1) * wvp));
1317             xy_to_upvp(cs->xGreen, cs->yGreen, &wup, &wvp);
1318             gx = (w - 1) * wup;
1319             gy = (h - 1) - ((int) ((h - 1) * wvp));
1320             xy_to_upvp(cs->xBlue, cs->yBlue, &wup, &wvp);
1321             bx = (w - 1) * wup;
1322             by = (h - 1) - ((int) ((h - 1) * wvp));
1323         } else if (cie == UCS) {
1324             double wu, wv;
1325             xy_to_uv(cs->xRed, cs->yRed, &wu, &wv);
1326             rx = (w - 1) * wu;
1327             ry = (h - 1) - ((int) ((h - 1) * wv));
1328             xy_to_uv(cs->xGreen, cs->yGreen, &wu, &wv);
1329             gx = (w - 1) * wu;
1330             gy = (h - 1) - ((int) ((h - 1) * wv));
1331             xy_to_uv(cs->xBlue, cs->yBlue, &wu, &wv);
1332             bx = (w - 1) * wu;
1333             by = (h - 1) - ((int) ((h - 1) * wv));
1334         } else if (cie == XYY) {
1335             rx = (w - 1) * cs->xRed;
1336             ry = (h - 1) - ((int) ((h - 1) * cs->yRed));
1337             gx = (w - 1) * cs->xGreen;
1338             gy = (h - 1) - ((int) ((h - 1) * cs->yGreen));
1339             bx = (w - 1) * cs->xBlue;
1340             by = (h - 1) - ((int) ((h - 1) * cs->yBlue));
1341         } else {
1342             av_assert0(0);
1343         }
1344
1345         draw_rline(pixels, linesize, rx, ry, gx, gy, w, h);
1346         draw_rline(pixels, linesize, gx, gy, bx, by, w, h);
1347         draw_rline(pixels, linesize, bx, by, rx, ry, w, h);
1348     }
1349 }
1350
1351 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
1352 {
1353     AVFilterContext *ctx  = inlink->dst;
1354     CiescopeContext *s = ctx->priv;
1355     AVFilterLink *outlink = ctx->outputs[0];
1356     int i = s->intensity * 65535;
1357     int w = outlink->w;
1358     int h = outlink->h;
1359     AVFrame *out;
1360     int ret, x, y;
1361
1362     out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
1363     if (!out) {
1364         av_frame_free(&in);
1365         return AVERROR(ENOMEM);
1366     }
1367     out->pts = in->pts;
1368
1369     if (!s->background) {
1370         ret = draw_background(ctx);
1371         if (ret < 0)
1372             return ret;
1373         s->background = 1;
1374     }
1375     for (y = 0; y < outlink->h; y++) {
1376         memset(out->data[0] + y * out->linesize[0], 0, outlink->w * 8);
1377     }
1378
1379     for (y = 0; y < in->height; y++) {
1380         for (x = 0; x < in->width; x++) {
1381             double cx, cy;
1382             uint16_t *dst;
1383             int wx, wy;
1384
1385             s->filter(ctx, in, &cx, &cy, x, y);
1386
1387             if (s->cie == LUV) {
1388                 double up, vp;
1389                 xy_to_upvp(cx, cy, &up, &vp);
1390                 cx = up;
1391                 cy = vp;
1392             } else if (s->cie == UCS) {
1393                 double u, v;
1394                 xy_to_uv(cx, cy, &u, &v);
1395                 cx = u;
1396                 cy = v;
1397             }
1398
1399             wx = (w - 1) * cx;
1400             wy = (h - 1) - ((h - 1) * cy);
1401
1402             if (wx < 0 || wx >= w ||
1403                 wy < 0 || wy >= h)
1404                 continue;
1405
1406             dst = (uint16_t *)(out->data[0] + wy * out->linesize[0] + wx * 8 + 0);
1407             dst[0] = FFMIN(dst[0] + i, 65535);
1408             dst[1] = FFMIN(dst[1] + i, 65535);
1409             dst[2] = FFMIN(dst[2] + i, 65535);
1410             dst[3] = 65535;
1411         }
1412     }
1413
1414     for (y = 0; y < outlink->h; y++) {
1415         uint16_t *dst = (uint16_t *)(out->data[0] + y * out->linesize[0]);
1416         const uint16_t *src = (const uint16_t *)(s->f->data[0] + y * s->f->linesize[0]);
1417         for (x = 0; x < outlink->w; x++) {
1418             const int xx = x * 4;
1419             if (dst[xx + 3] == 0) {
1420                 dst[xx + 0] = src[xx + 0];
1421                 dst[xx + 1] = src[xx + 1];
1422                 dst[xx + 2] = src[xx + 2];
1423                 dst[xx + 3] = src[xx + 3];
1424             }
1425         }
1426     }
1427
1428     if (s->show_white)
1429         plot_white_point((uint16_t *)out->data[0], out->linesize[0] / 2,
1430                          outlink->w, outlink->h, 65535,
1431                          s->color_system, s->cie);
1432
1433     plot_gamuts((uint16_t *)out->data[0], out->linesize[0] / 2,
1434                 outlink->w, outlink->h,
1435                 s->cie, s->gamuts);
1436
1437     av_frame_free(&in);
1438     return ff_filter_frame(outlink, out);
1439 }
1440
1441 static void av_cold uninit(AVFilterContext *ctx)
1442 {
1443     CiescopeContext *s = ctx->priv;
1444
1445     av_frame_free(&s->f);
1446 }
1447
1448 static int config_input(AVFilterLink *inlink)
1449 {
1450     CiescopeContext *s = inlink->dst->priv;
1451     int i;
1452
1453     get_rgb2xyz_matrix(color_systems[s->color_system], s->m);
1454     invert_matrix3x3(s->m, s->i);
1455
1456     switch (inlink->format) {
1457     case AV_PIX_FMT_RGB24:
1458         s->filter = filter_rgb24;
1459         break;
1460     case AV_PIX_FMT_RGBA:
1461         s->filter = filter_rgba;
1462         break;
1463     case AV_PIX_FMT_RGB48:
1464         s->filter = filter_rgb48;
1465         break;
1466     case AV_PIX_FMT_RGBA64:
1467         s->filter = filter_rgba64;
1468         break;
1469     case AV_PIX_FMT_XYZ12:
1470         s->filter = filter_xyz;
1471         for (i = 0; i < 65536; i++)
1472             s->log2lin[i] = pow(i / 65535., s->igamma) * 65535.;
1473         break;
1474     default:
1475         av_assert0(0);
1476     }
1477
1478     return 0;
1479 }
1480
1481 static const AVFilterPad inputs[] = {
1482     {
1483         .name         = "default",
1484         .type         = AVMEDIA_TYPE_VIDEO,
1485         .filter_frame = filter_frame,
1486         .config_props = config_input,
1487     },
1488     { NULL }
1489 };
1490
1491 static const AVFilterPad outputs[] = {
1492     {
1493         .name         = "default",
1494         .type         = AVMEDIA_TYPE_VIDEO,
1495         .config_props = config_output,
1496     },
1497     { NULL }
1498 };
1499
1500 AVFilter ff_vf_ciescope = {
1501     .name          = "ciescope",
1502     .description   = NULL_IF_CONFIG_SMALL("Video CIE scope."),
1503     .priv_size     = sizeof(CiescopeContext),
1504     .priv_class    = &ciescope_class,
1505     .query_formats = query_formats,
1506     .uninit        = uninit,
1507     .inputs        = inputs,
1508     .outputs       = outputs,
1509 };