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