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]->outcfg.formats)) < 0)
 143         return ret;
 144
 145     if ((ret = ff_formats_ref(ff_make_format_list(out_pix_fmts), &ctx->outputs[0]->incfg.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 = (w - 1) * wup;
1194         wy = (h - 1) - ((int) ((h - 1) * wvp));
1195     } else if (cie == UCS) {
1196         double wu, wv;
1197         xy_to_uv(cs->xWhite, cs->yWhite, &wu, &wv);
1198         wx = (w - 1) * wu;
1199         wy = (h - 1) - ((int) ((h - 1) * wv));
1200     } else if (cie == XYY) {
1201         wx = (w - 1) * cs->xWhite;
1202         wy = (h - 1) - ((int) ((h - 1) * cs->yWhite));
1203     } else {
1204         av_assert0(0);
1205     }
1206
1207     draw_rline(pixels, linesize,
1208                wx + Sz(3), wy, wx + Sz(10), wy,
1209                w, h);
1210     draw_rline(pixels, linesize,
1211                wx - Sz(3), wy, wx - Sz(10), wy,
1212                w, h);
1213     draw_rline(pixels, linesize,
1214                wx, wy + Sz(3), wx, wy + Sz(10),
1215                w, h);
1216     draw_rline(pixels, linesize,
1217                wx, wy - Sz(3), wx, wy - Sz(10),
1218                w, h);
1219 }
1220
1221 static int draw_background(AVFilterContext *ctx)
1222 {
1223     CiescopeContext *s = ctx->priv;
1224     const struct ColorSystem *cs = &color_systems[s->color_system];
1225     AVFilterLink *outlink = ctx->outputs[0];
1226     int w = s->size;
1227     int h = s->size;
1228     uint16_t *pixels;
1229
1230     if ((s->f = ff_get_video_buffer(outlink, outlink->w, outlink->h)) == NULL)
1231         return AVERROR(ENOMEM);
1232     pixels = (uint16_t *)s->f->data[0];
1233
1234     tongue_outline(pixels, s->f->linesize[0] / 2, w, h, 65535, s->cie);
1235
1236     fill_in_tongue(pixels, s->f->linesize[0] / 2, w, h, 65535, cs, (const double (*)[3])s->i, s->cie,
1237                    s->correct_gamma, s->contrast);
1238
1239     return 0;
1240 }
1241
1242 static void filter_rgb48(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1243 {
1244     CiescopeContext *s = ctx->priv;
1245     const uint16_t* src = (const uint16_t*)(in->data[0] + in->linesize[0] * y + x * 6);
1246     double r = src[0] / 65535.;
1247     double g = src[1] / 65535.;
1248     double b = src[2] / 65535.;
1249     double cz;
1250
1251     rgb_to_xy(r, g, b, cx, cy, &cz, (const double (*)[3])s->m);
1252 }
1253
1254 static void filter_rgba64(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1255 {
1256     CiescopeContext *s = ctx->priv;
1257     const uint16_t* src = (const uint16_t*)(in->data[0] + in->linesize[0] * y + x * 8);
1258     double r = src[0] / 65535.;
1259     double g = src[1] / 65535.;
1260     double b = src[2] / 65535.;
1261     double cz;
1262
1263     rgb_to_xy(r, g, b, cx, cy, &cz, (const double (*)[3])s->m);
1264 }
1265
1266 static void filter_rgb24(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1267 {
1268     CiescopeContext *s = ctx->priv;
1269     const uint8_t* src = in->data[0] + in->linesize[0] * y + x * 3;
1270     double r = src[0] / 255.;
1271     double g = src[1] / 255.;
1272     double b = src[2] / 255.;
1273     double cz;
1274
1275     rgb_to_xy(r, g, b, cx, cy, &cz, (const double (*)[3])s->m);
1276 }
1277
1278 static void filter_rgba(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1279 {
1280     CiescopeContext *s = ctx->priv;
1281     const uint8_t* src = in->data[0] + in->linesize[0] * y + x * 4;
1282     double r = src[0] / 255.;
1283     double g = src[1] / 255.;
1284     double b = src[2] / 255.;
1285     double cz;
1286
1287     rgb_to_xy(r, g, b, cx, cy, &cz, (const double (*)[3])s->m);
1288 }
1289
1290 static void filter_xyz(AVFilterContext *ctx, AVFrame *in, double *cx, double *cy, int x, int y)
1291 {
1292     CiescopeContext *s = ctx->priv;
1293     const uint16_t* src = (uint16_t *)(in->data[0] + in->linesize[0] * y + x * 6);
1294     double lx = s->log2lin[src[0]];
1295     double ly = s->log2lin[src[1]];
1296     double lz = s->log2lin[src[2]];
1297     double sum = lx + ly + lz;
1298
1299     if (sum == 0)
1300         sum = 1;
1301     *cx = lx / sum;
1302     *cy = ly / sum;
1303 }
1304
1305 static void plot_gamuts(uint16_t *pixels, int linesize, int w, int h,
1306                         int cie, int gamuts)
1307 {
1308     int i;
1309
1310     for (i = 0; i < NB_CS; i++) {
1311         const struct ColorSystem *cs = &color_systems[i];
1312         int rx, ry, gx, gy, bx, by;
1313
1314         if (!((1 << i) & gamuts))
1315             continue;
1316         if (cie == LUV) {
1317             double wup, wvp;
1318             xy_to_upvp(cs->xRed, cs->yRed, &wup, &wvp);
1319             rx = (w - 1) * wup;
1320             ry = (h - 1) - ((int) ((h - 1) * wvp));
1321             xy_to_upvp(cs->xGreen, cs->yGreen, &wup, &wvp);
1322             gx = (w - 1) * wup;
1323             gy = (h - 1) - ((int) ((h - 1) * wvp));
1324             xy_to_upvp(cs->xBlue, cs->yBlue, &wup, &wvp);
1325             bx = (w - 1) * wup;
1326             by = (h - 1) - ((int) ((h - 1) * wvp));
1327         } else if (cie == UCS) {
1328             double wu, wv;
1329             xy_to_uv(cs->xRed, cs->yRed, &wu, &wv);
1330             rx = (w - 1) * wu;
1331             ry = (h - 1) - ((int) ((h - 1) * wv));
1332             xy_to_uv(cs->xGreen, cs->yGreen, &wu, &wv);
1333             gx = (w - 1) * wu;
1334             gy = (h - 1) - ((int) ((h - 1) * wv));
1335             xy_to_uv(cs->xBlue, cs->yBlue, &wu, &wv);
1336             bx = (w - 1) * wu;
1337             by = (h - 1) - ((int) ((h - 1) * wv));
1338         } else if (cie == XYY) {
1339             rx = (w - 1) * cs->xRed;
1340             ry = (h - 1) - ((int) ((h - 1) * cs->yRed));
1341             gx = (w - 1) * cs->xGreen;
1342             gy = (h - 1) - ((int) ((h - 1) * cs->yGreen));
1343             bx = (w - 1) * cs->xBlue;
1344             by = (h - 1) - ((int) ((h - 1) * cs->yBlue));
1345         } else {
1346             av_assert0(0);
1347         }
1348
1349         draw_rline(pixels, linesize, rx, ry, gx, gy, w, h);
1350         draw_rline(pixels, linesize, gx, gy, bx, by, w, h);
1351         draw_rline(pixels, linesize, bx, by, rx, ry, w, h);
1352     }
1353 }
1354
1355 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
1356 {
1357     AVFilterContext *ctx  = inlink->dst;
1358     CiescopeContext *s = ctx->priv;
1359     AVFilterLink *outlink = ctx->outputs[0];
1360     int i = s->intensity * 65535;
1361     int w = outlink->w;
1362     int h = outlink->h;
1363     AVFrame *out;
1364     int ret, x, y;
1365
1366     out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
1367     if (!out) {
1368         av_frame_free(&in);
1369         return AVERROR(ENOMEM);
1370     }
1371     out->pts = in->pts;
1372
1373     if (!s->background) {
1374         ret = draw_background(ctx);
1375         if (ret < 0) {
1376             av_frame_free(&out);
1377             return ret;
1378         }
1379         s->background = 1;
1380     }
1381     for (y = 0; y < outlink->h; y++) {
1382         memset(out->data[0] + y * out->linesize[0], 0, outlink->w * 8);
1383     }
1384
1385     for (y = 0; y < in->height; y++) {
1386         for (x = 0; x < in->width; x++) {
1387             double cx, cy;
1388             uint16_t *dst;
1389             int wx, wy;
1390
1391             s->filter(ctx, in, &cx, &cy, x, y);
1392
1393             if (s->cie == LUV) {
1394                 double up, vp;
1395                 xy_to_upvp(cx, cy, &up, &vp);
1396                 cx = up;
1397                 cy = vp;
1398             } else if (s->cie == UCS) {
1399                 double u, v;
1400                 xy_to_uv(cx, cy, &u, &v);
1401                 cx = u;
1402                 cy = v;
1403             }
1404
1405             wx = (w - 1) * cx;
1406             wy = (h - 1) - ((h - 1) * cy);
1407
1408             if (wx < 0 || wx >= w ||
1409                 wy < 0 || wy >= h)
1410                 continue;
1411
1412             dst = (uint16_t *)(out->data[0] + wy * out->linesize[0] + wx * 8 + 0);
1413             dst[0] = FFMIN(dst[0] + i, 65535);
1414             dst[1] = FFMIN(dst[1] + i, 65535);
1415             dst[2] = FFMIN(dst[2] + i, 65535);
1416             dst[3] = 65535;
1417         }
1418     }
1419
1420     for (y = 0; y < outlink->h; y++) {
1421         uint16_t *dst = (uint16_t *)(out->data[0] + y * out->linesize[0]);
1422         const uint16_t *src = (const uint16_t *)(s->f->data[0] + y * s->f->linesize[0]);
1423         for (x = 0; x < outlink->w; x++) {
1424             const int xx = x * 4;
1425             if (dst[xx + 3] == 0) {
1426                 dst[xx + 0] = src[xx + 0];
1427                 dst[xx + 1] = src[xx + 1];
1428                 dst[xx + 2] = src[xx + 2];
1429                 dst[xx + 3] = src[xx + 3];
1430             }
1431         }
1432     }
1433
1434     if (s->show_white)
1435         plot_white_point((uint16_t *)out->data[0], out->linesize[0] / 2,
1436                          outlink->w, outlink->h, 65535,
1437                          s->color_system, s->cie);
1438
1439     plot_gamuts((uint16_t *)out->data[0], out->linesize[0] / 2,
1440                 outlink->w, outlink->h,
1441                 s->cie, s->gamuts);
1442
1443     av_frame_free(&in);
1444     return ff_filter_frame(outlink, out);
1445 }
1446
1447 static void av_cold uninit(AVFilterContext *ctx)
1448 {
1449     CiescopeContext *s = ctx->priv;
1450
1451     av_frame_free(&s->f);
1452 }
1453
1454 static int config_input(AVFilterLink *inlink)
1455 {
1456     CiescopeContext *s = inlink->dst->priv;
1457     int i;
1458
1459     get_rgb2xyz_matrix(color_systems[s->color_system], s->m);
1460     invert_matrix3x3(s->m, s->i);
1461
1462     switch (inlink->format) {
1463     case AV_PIX_FMT_RGB24:
1464         s->filter = filter_rgb24;
1465         break;
1466     case AV_PIX_FMT_RGBA:
1467         s->filter = filter_rgba;
1468         break;
1469     case AV_PIX_FMT_RGB48:
1470         s->filter = filter_rgb48;
1471         break;
1472     case AV_PIX_FMT_RGBA64:
1473         s->filter = filter_rgba64;
1474         break;
1475     case AV_PIX_FMT_XYZ12:
1476         s->filter = filter_xyz;
1477         for (i = 0; i < 65536; i++)
1478             s->log2lin[i] = pow(i / 65535., s->igamma) * 65535.;
1479         break;
1480     default:
1481         av_assert0(0);
1482     }
1483
1484     return 0;
1485 }
1486
1487 static const AVFilterPad inputs[] = {
1488     {
1489         .name         = "default",
1490         .type         = AVMEDIA_TYPE_VIDEO,
1491         .filter_frame = filter_frame,
1492         .config_props = config_input,
1493     },
1494     { NULL }
1495 };
1496
1497 static const AVFilterPad outputs[] = {
1498     {
1499         .name         = "default",
1500         .type         = AVMEDIA_TYPE_VIDEO,
1501         .config_props = config_output,
1502     },
1503     { NULL }
1504 };
1505
1506 const AVFilter ff_vf_ciescope = {
1507     .name          = "ciescope",
1508     .description   = NULL_IF_CONFIG_SMALL("Video CIE scope."),
1509     .priv_size     = sizeof(CiescopeContext),
1510     .priv_class    = &ciescope_class,
1511     .query_formats = query_formats,
1512     .uninit        = uninit,
1513     .inputs        = inputs,
1514     .outputs       = outputs,
1515 };