/*

      Generate a PPM file representing a CIE colour gamut chart

               by John Walker  --  kelvin@@fourmilab.ch
               WWW home page: http://www.fourmilab.ch/

    Permission  to  use, copy, modify, and distribute this software and
    its documentation  for  any  purpose  and  without  fee  is  hereby
    granted,  without any conditions or restrictions.  This software is
    provided "as is" without express or implied warranty.

    This program was called cietoppm in Walker's original work.  
    Because "cie" is not a graphics format, Bryan changed the name
    when he integrated it into the Netpbm package in March 2000.
*/

/*
  Modified by
  Andrew J. S. Hamilton 21 May 1999
  Andrew.Hamilton@Colorado.EDU
  http://casa.colorado.edu/~ajsh/

  Corrected XYZ -> RGB transform.
  Introduced gamma correction.
  Introduced option to plot 1976 u' v' chromaticities.
*/

#include <math.h>
#include "ppm.h"
#include "ppmdraw.h"

#define CLAMP(v, l, h)  ((v) < (l) ? (l) : (v) > (h) ? (h) : (v))
#define TRUE    1
#define FALSE   0

#define Maxval  255                   /* Maxval to use in generated pixmaps */

/* A  colour  system is defined by the CIE x and y  coordinates of its
   three primary illuminants and the x and y coordinates of the  white
   point. */

struct colourSystem {
    char *name;                       /* Colour system name */
    double xRed, yRed,                /* Red primary illuminant */
           xGreen, yGreen,            /* Green primary illuminant */
           xBlue, yBlue,              /* Blue primary illuminant */
           xWhite, yWhite,            /* White point */
           gamma;		      /* gamma of nonlinear correction */
};

/* The  following  table  gives  the  CIE  colour  matching  functions
   \bar{x}(\lambda),  \bar{y}(\lambda),  and   \bar{z}(\lambda),   for
   wavelengths  \lambda  at 5 nanometre increments from 380 nm through
   780 nm.  This table is used in conjunction with  Planck's  law  for
   the  energy spectrum of a black body at a given temperature to plot
   the black body curve on the CIE chart. */

static double cie_colour_match[][3] = {
    { 0.0014, 0.0000, 0.0065 },       /* 380 nm */
    { 0.0022, 0.0001, 0.0105 },
    { 0.0042, 0.0001, 0.0201 },
    { 0.0076, 0.0002, 0.0362 },
    { 0.0143, 0.0004, 0.0679 },
    { 0.0232, 0.0006, 0.1102 },
    { 0.0435, 0.0012, 0.2074 },
    { 0.0776, 0.0022, 0.3713 },
    { 0.1344, 0.0040, 0.6456 },
    { 0.2148, 0.0073, 1.0391 },
    { 0.2839, 0.0116, 1.3856 },
    { 0.3285, 0.0168, 1.6230 },
    { 0.3483, 0.0230, 1.7471 },
    { 0.3481, 0.0298, 1.7826 },
    { 0.3362, 0.0380, 1.7721 },
    { 0.3187, 0.0480, 1.7441 },
    { 0.2908, 0.0600, 1.6692 },
    { 0.2511, 0.0739, 1.5281 },
    { 0.1954, 0.0910, 1.2876 },
    { 0.1421, 0.1126, 1.0419 },
    { 0.0956, 0.1390, 0.8130 },
    { 0.0580, 0.1693, 0.6162 },
    { 0.0320, 0.2080, 0.4652 },
    { 0.0147, 0.2586, 0.3533 },
    { 0.0049, 0.3230, 0.2720 },
    { 0.0024, 0.4073, 0.2123 },
    { 0.0093, 0.5030, 0.1582 },
    { 0.0291, 0.6082, 0.1117 },
    { 0.0633, 0.7100, 0.0782 },
    { 0.1096, 0.7932, 0.0573 },
    { 0.1655, 0.8620, 0.0422 },
    { 0.2257, 0.9149, 0.0298 },
    { 0.2904, 0.9540, 0.0203 },
    { 0.3597, 0.9803, 0.0134 },
    { 0.4334, 0.9950, 0.0087 },
    { 0.5121, 1.0000, 0.0057 },
    { 0.5945, 0.9950, 0.0039 },
    { 0.6784, 0.9786, 0.0027 },
    { 0.7621, 0.9520, 0.0021 },
    { 0.8425, 0.9154, 0.0018 },
    { 0.9163, 0.8700, 0.0017 },
    { 0.9786, 0.8163, 0.0014 },
    { 1.0263, 0.7570, 0.0011 },
    { 1.0567, 0.6949, 0.0010 },
    { 1.0622, 0.6310, 0.0008 },
    { 1.0456, 0.5668, 0.0006 },
    { 1.0026, 0.5030, 0.0003 },
    { 0.9384, 0.4412, 0.0002 },
    { 0.8544, 0.3810, 0.0002 },
    { 0.7514, 0.3210, 0.0001 },
    { 0.6424, 0.2650, 0.0000 },
    { 0.5419, 0.2170, 0.0000 },
    { 0.4479, 0.1750, 0.0000 },
    { 0.3608, 0.1382, 0.0000 },
    { 0.2835, 0.1070, 0.0000 },
    { 0.2187, 0.0816, 0.0000 },
    { 0.1649, 0.0610, 0.0000 },
    { 0.1212, 0.0446, 0.0000 },
    { 0.0874, 0.0320, 0.0000 },
    { 0.0636, 0.0232, 0.0000 },
    { 0.0468, 0.0170, 0.0000 },
    { 0.0329, 0.0119, 0.0000 },
    { 0.0227, 0.0082, 0.0000 },
    { 0.0158, 0.0057, 0.0000 },
    { 0.0114, 0.0041, 0.0000 },
    { 0.0081, 0.0029, 0.0000 },
    { 0.0058, 0.0021, 0.0000 },
    { 0.0041, 0.0015, 0.0000 },
    { 0.0029, 0.0010, 0.0000 },
    { 0.0020, 0.0007, 0.0000 },
    { 0.0014, 0.0005, 0.0000 },
    { 0.0010, 0.0004, 0.0000 },
    { 0.0007, 0.0002, 0.0000 },
    { 0.0005, 0.0002, 0.0000 },
    { 0.0003, 0.0001, 0.0000 },
    { 0.0002, 0.0001, 0.0000 },
    { 0.0002, 0.0001, 0.0000 },
    { 0.0001, 0.0000, 0.0000 },
    { 0.0001, 0.0000, 0.0000 },
    { 0.0001, 0.0000, 0.0000 },
    { 0.0000, 0.0000, 0.0000 }        /* 780 nm */
};

/* The following table gives the  spectral  chromaticity  co-ordinates
   x(\lambda) and y(\lambda) for wavelengths in 5 nanometre increments
   from 380 nm through  780  nm.   These  co-ordinates  represent  the
   position in the CIE x-y space of pure spectral colours of the given
   wavelength, and  thus  define  the  outline  of  the  CIE  "tongue"
   diagram. */

static double spectral_chromaticity[81][3] = {
    { 0.1741, 0.0050 },               /* 380 nm */
    { 0.1740, 0.0050 },
    { 0.1738, 0.0049 },
    { 0.1736, 0.0049 },
    { 0.1733, 0.0048 },
    { 0.1730, 0.0048 },
    { 0.1726, 0.0048 },
    { 0.1721, 0.0048 },
    { 0.1714, 0.0051 },
    { 0.1703, 0.0058 },
    { 0.1689, 0.0069 },
    { 0.1669, 0.0086 },
    { 0.1644, 0.0109 },
    { 0.1611, 0.0138 },
    { 0.1566, 0.0177 },
    { 0.1510, 0.0227 },
    { 0.1440, 0.0297 },
    { 0.1355, 0.0399 },
    { 0.1241, 0.0578 },
    { 0.1096, 0.0868 },
    { 0.0913, 0.1327 },
    { 0.0687, 0.2007 },
    { 0.0454, 0.2950 },
    { 0.0235, 0.4127 },
    { 0.0082, 0.5384 },
    { 0.0039, 0.6548 },
    { 0.0139, 0.7502 },
    { 0.0389, 0.8120 },
    { 0.0743, 0.8338 },
    { 0.1142, 0.8262 },
    { 0.1547, 0.8059 },
    { 0.1929, 0.7816 },
    { 0.2296, 0.7543 },
    { 0.2658, 0.7243 },
    { 0.3016, 0.6923 },
    { 0.3373, 0.6589 },
    { 0.3731, 0.6245 },
    { 0.4087, 0.5896 },
    { 0.4441, 0.5547 },
    { 0.4788, 0.5202 },
    { 0.5125, 0.4866 },
    { 0.5448, 0.4544 },
    { 0.5752, 0.4242 },
    { 0.6029, 0.3965 },
    { 0.6270, 0.3725 },
    { 0.6482, 0.3514 },
    { 0.6658, 0.3340 },
    { 0.6801, 0.3197 },
    { 0.6915, 0.3083 },
    { 0.7006, 0.2993 },
    { 0.7079, 0.2920 },
    { 0.7140, 0.2859 },
    { 0.7190, 0.2809 },
    { 0.7230, 0.2770 },
    { 0.7260, 0.2740 },
    { 0.7283, 0.2717 },
    { 0.7300, 0.2700 },
    { 0.7311, 0.2689 },
    { 0.7320, 0.2680 },
    { 0.7327, 0.2673 },
    { 0.7334, 0.2666 },
    { 0.7340, 0.2660 },
    { 0.7344, 0.2656 },
    { 0.7346, 0.2654 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 },
    { 0.7347, 0.2653 }                /* 780 nm */
};

static pixel **pixels;                /* Pixel map */
static int pixcols, pixrows;          /* Pixel map size */
static int sxsize = 512, sysize = 512; /* X, Y size */

/* Standard white point chromaticities. */

#define IlluminantC     0.3101, 0.3162  /* For NTSC television */
#define IlluminantD65   0.3127, 0.3291  /* For EBU and SMPTE */

/* Gamma of nonlinear correction.
   See Charles Poynton's ColorFAQ Item 45 and GammaFAQ Item 6 at
   http://www.inforamp.net/~poynton/ColorFAQ.html
   http://www.inforamp.net/~poynton/GammaFAQ.html
*/

#define GAMMA_REC709	0.		/* Rec. 709 */

static struct colourSystem const
    NTSCsystem  =  { "NTSC",               0.67,  0.33,  0.21,  0.71,  0.14,  0.08,  IlluminantC,   GAMMA_REC709},
    EBUsystem   =  { "EBU (PAL/SECAM)",    0.64,  0.33,  0.29,  0.60,  0.15,  0.06,  IlluminantD65, GAMMA_REC709},
    SMPTEsystem =  { "SMPTE",             0.630, 0.340, 0.310, 0.595, 0.155, 0.070, IlluminantD65,  GAMMA_REC709},
    HDTVsystem  =  { "HDTV",              0.670, 0.330, 0.210, 0.710, 0.150, 0.060, IlluminantD65,  GAMMA_REC709},
/* Huh? -ajsh
    CIEsystem   =  { "CIE",              0.7355,0.2645,0.2658,0.7243,0.1669,0.0085, 0.3894,0.3324,  GAMMA_REC709},
*/
    CIEsystem   =  { "CIE",              0.7355,0.2645,0.2658,0.7243,0.1669,0.0085, 0.33333333,0.33333333,  GAMMA_REC709},
    Rec709system = { "CIE REC 709",        0.64,  0.33,  0.30,  0.60,  0.15,  0.06,    IlluminantD65, GAMMA_REC709};

/* Customsystem  is a variable that is initialized to CIE Rec 709, but
   we modify it with information specified by the user's options.
*/
static struct colourSystem Customsystem = 
{ "Custom",             0.64,  0.33,  0.30,  0.60,  0.15,  0.06,  
  IlluminantD65, GAMMA_REC709};
    


/*
    Given 1976 coordinates u', v', determine 1931 chromaticities x, y
*/
static void upvp_to_xy(up, vp, xc, yc)
  double up, vp;
  double *xc, *yc;
{
    *xc = 9*up / (6*up - 16*vp + 12);
    *yc = 4*vp / (6*up - 16*vp + 12);
}

/*
    Given 1931 chromaticities x, y, determine 1976 coordinates u', v'
*/
static void xy_to_upvp(xc, yc, up, vp)
  double xc, yc;
  double *up, *vp;
{
    *up = 4*xc / (- 2*xc + 12*yc + 3);
    *vp = 9*yc / (- 2*xc + 12*yc + 3);
}

/*                            XYZ_TO_RGB

    Given  an additive tricolour system CS, defined by the CIE x and y
    chromaticities of its three primaries (z is derived  trivially  as
    1-(x+y)),  and  a  desired chromaticity (XC, YC, ZC) in CIE space,
    determine the contribution of each primary in a linear combination
    which   sums  to  the  desired  chromaticity.   If  the  requested
    chromaticity falls outside the  Maxwell  triangle  (colour  gamut)
    formed  by the three primaries, one of the r, g, or b weights will
    be negative.  Use inside_gamut() to test for a  valid  colour  and
    constrain_rgb()  to  desaturate  an  outside-gamut  colour  to the
    closest representation within the available gamut. */

static void xyz_to_rgb(cs, xc, yc, zc, r, g, b)
  const struct colourSystem *cs;
  double xc, yc, zc;
  double *r, *g, *b;
{
    double xr, yr, zr, xg, yg, zg, xb, yb, zb;
    double xw, yw, zw;
    double rx, ry, rz, gx, gy, gz, bx, by, bz;
    double rw, gw, bw;

    xr = cs->xRed;    yr = cs->yRed;    zr = 1 - (xr + yr);
    xg = cs->xGreen;  yg = cs->yGreen;  zg = 1 - (xg + yg);
    xb = cs->xBlue;   yb = cs->yBlue;   zb = 1 - (xb + yb);

    xw = cs->xWhite;  yw = cs->yWhite;  zw = 1 - (xw + yw);

    /* xyz -> rgb matrix, before scaling to white. */
    rx = yg*zb - yb*zg;  ry = xb*zg - xg*zb;  rz = xg*yb - xb*yg;
    gx = yb*zr - yr*zb;  gy = xr*zb - xb*zr;  gz = xb*yr - xr*yb;
    bx = yr*zg - yg*zr;  by = xg*zr - xr*zg;  bz = xr*yg - xg*yr;

    /* White scaling factors.
       Dividing by yw scales the white luminance to unity, as conventional. */
    rw = (rx*xw + ry*yw + rz*zw) / yw;
    gw = (gx*xw + gy*yw + gz*zw) / yw;
    bw = (bx*xw + by*yw + bz*zw) / yw;

    /* xyz -> rgb matrix, correctly scaled to white. */
    rx = rx / rw;  ry = ry / rw;  rz = rz / rw;
    gx = gx / gw;  gy = gy / gw;  gz = gz / gw;
    bx = bx / bw;  by = by / bw;  bz = bz / bw;

    /* rgb of the desired point */
    *r = rx*xc + ry*yc + rz*zc;
    *g = gx*xc + gy*yc + gz*zc;
    *b = bx*xc + by*yc + bz*zc;
}

/*                           INSIDE_GAMUT

    Test whether a requested colour is  within  the  gamut  achievable
    with  the  primaries  of  the current colour system.  This amounts
    simply  to  testing  whether   all   the   primary   weights   are
    non-negative. */

static int inside_gamut(r, g, b)
  double r, g, b;
{
    return r >= 0 && g >= 0 && b >= 0;
}

/*                          CONSTRAIN-RGB

    If  the  requested RGB shade contains a negative weight for one of
    the primaries, it lies outside the colour  gamut  accessible  from
    the  given  triple  of  primaries.  Desaturate it by adding white,
    equal quantities of R, G, and B, enough to make RGB all positive.
*/

static int constrain_rgb(r, g, b)
  double *r, *g, *b;
{
    double w;

    /* Amount of white needed is w = - min(0, *r, *g, *b) */
    w = (0 < *r) ? 0 : *r;
    w = (w < *g) ? w : *g;
    w = (w < *b) ? w : *b;
    w = - w;

    /* Add just enough white to make r, g, b all positive. */
    if (w > 0) {
	*r += w;  *g += w; *b += w;

        return 1;                     /* Colour modified to fit RGB gamut */
    }

    return 0;                         /* Colour within RGB gamut */
}

/*                          GAMMA-CORRECT-RGB

    Transform linear RGB values to nonlinear RGB values.
    Rec. 709 is
    ITU-R Recommendation BT. 709 (1990)
    ``Basic Parameter Values for the HDTV Standard for the Studio and for
    International Programme Exchange'', formerly CCIR Rec. 709.
    For details see
       http://www.inforamp.net/~poynton/ColorFAQ.html
       http://www.inforamp.net/~poynton/GammaFAQ.html
*/

static void gamma_correct(cs, c)
  const struct colourSystem *cs;
  double *c;
{
    double gamma;

    gamma = cs->gamma;

    if (gamma == 0.) {
	/* Rec. 709 gamma correction. */
	double cc = 0.018;
	if (*c < cc) {
	    *c *= (1.099 * pow(cc, 0.45) - 0.099) / cc;
	} else {
	    *c = 1.099 * pow(*c, 0.45) - 0.099;
	}
    } else {
	/* Nonlinear colour = (Linear colour)^(1/gamma) */
	*c = pow(*c, 1./gamma);
    }
}

static void gamma_correct_rgb(cs, r, g, b)
  const struct colourSystem *cs;
  double *r, *g, *b;
{
    gamma_correct(cs, r);
    gamma_correct(cs, g);
    gamma_correct(cs, b);
}

/*  Main program. */

int main(argc, argv)
  int argc;
  char *argv[];
{
    int argn, x, y;
    char *usage = "[-[no]black] [-[no]wpoint] [-[no]label] [-no[axes]] [-full]\n\
                 [-xy|-upvp] [-rec709|-ntsc|-ebu|-smpte|-hdtv|-cie]\n\
                 [-red <x> <y>] [-green <x> <y>] [-blue <x> <y>]\n\
                 [-white <x> <y>] [-gamma <g>]\n\
                 [-size <s>] [-xsize|-width <x>] [-ysize|-height <y>]";
    const struct colourSystem *cs;

    int widspec = FALSE, hgtspec = FALSE,
        ix, icx, icy, fx, fy, lx, ly, xBias, yBias, pxcols, pxrows;
    pixel rgbcolour;                  /* Pixel used to clear pixmap */
    int upvp = FALSE;		      /* xy or u'v' color coordinates? */
    int showWhite = TRUE;             /* Show white point ? */
    int showBlack = TRUE;             /* Show black body curve ? */
    int fullChart = FALSE;            /* Fill entire tongue ? */
    int showLabel = TRUE;             /* Show labels ? */
    int showAxes = TRUE;              /* Plot axes ? */
    char sysdesc[256];                /* Colour system description */

    ppm_init(&argc, argv);
    argn = 1;

    cs = &Rec709system;  /* default */
    while (argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0') {
        if (pm_keymatch(argv[argn], "-xy", 2)) {
            upvp = FALSE;
        } else if (pm_keymatch(argv[argn], "-upvp", 1)) {
            upvp = TRUE;
        } else if (pm_keymatch(argv[argn], "-xsize", 1) ||
                   pm_keymatch(argv[argn], "-width", 2)) {
            if (widspec) {
                pm_error("already specified a size/width/xsize");
            }
            argn++;
            if ((argn == argc) || (sscanf(argv[argn], "%d", &sxsize) != 1))
                pm_usage(usage);
            widspec = TRUE;
        } else if (pm_keymatch(argv[argn], "-ysize", 1) ||
                   pm_keymatch(argv[argn], "-height", 2)) {
            if (hgtspec) {
                pm_error("already specified a size/height/ysize");
            }
            argn++;
            if ((argn == argc) || (sscanf(argv[argn], "%d", &sysize) != 1))
                pm_usage(usage);
            hgtspec = TRUE;
        } else if (pm_keymatch(argv[argn], "-size", 2)) {
            if (hgtspec || widspec) {
                pm_error("already specified a size/height/ysize");
            }
            argn++;
            if ((argn == argc) || (sscanf(argv[argn], "%d", &sysize) != 1))
                pm_usage(usage);
            sxsize = sysize;
            hgtspec = widspec = TRUE;
        } else if (pm_keymatch(argv[argn], "-rec709", 1)) {
            cs = &Rec709system;
        } else if (pm_keymatch(argv[argn], "-ntsc", 1)) {
            cs = &NTSCsystem;
        } else if (pm_keymatch(argv[argn], "-ebu", 1)) {
            cs = &EBUsystem;
        } else if (pm_keymatch(argv[argn], "-smpte", 2)) {
            cs = &SMPTEsystem;
        } else if (pm_keymatch(argv[argn], "-hdtv", 2)) {
            cs = &HDTVsystem;                 
        } else if (pm_keymatch(argv[argn], "-cie", 1)) {
            cs = &CIEsystem;                 
        } else if (pm_keymatch(argv[argn], "-black", 3)) {
            showBlack = TRUE;         /* Show black body curve */
        } else if (pm_keymatch(argv[argn], "-wpoint", 2)) {
            showWhite = TRUE;         /* Show white point of colour system */
        } else if (pm_keymatch(argv[argn], "-noblack", 3)) {
            showBlack = FALSE;        /* Don't show black body curve */
        } else if (pm_keymatch(argv[argn], "-nowpoint", 3)) {
            showWhite = FALSE;        /* Don't show white point of system */
        } else if (pm_keymatch(argv[argn], "-label", 1)) {
            showLabel = TRUE;	      /* Show labels. */
        } else if (pm_keymatch(argv[argn], "-nolabel", 3)) {
            showLabel = FALSE;	      /* Don't show labels */
        } else if (pm_keymatch(argv[argn], "-axes", 1)) {
            showAxes = TRUE;	      /* Show axes. */
        } else if (pm_keymatch(argv[argn], "-noaxes", 3)) {
            showAxes = FALSE;	      /* Don't show axes */
        } else if (pm_keymatch(argv[argn], "-full", 1)) {
            fullChart = TRUE;         /* Fill whole tongue full-intensity */
        } else if (pm_keymatch(argv[argn], "-gamma", 2)) {
            cs = &Customsystem;
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.gamma) != 1))
                pm_usage(usage);
        } else if (pm_keymatch(argv[argn], "-red", 1)) {
            cs = &Customsystem;
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.xRed) != 1))
                pm_usage(usage);
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.yRed) != 1))
                pm_usage(usage);
        } else if (pm_keymatch(argv[argn], "-green", 1)) {
            cs = &Customsystem;
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.xGreen) != 1))
                pm_usage(usage);
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.yGreen) != 1))
                pm_usage(usage);
        } else if (pm_keymatch(argv[argn], "-blue", 1)) {
            cs = &Customsystem;
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.xBlue) != 1))
                pm_usage(usage);
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.yBlue) != 1))
                pm_usage(usage);
        } else if (pm_keymatch(argv[argn], "-white", 1)) {
            cs = &Customsystem;
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.xWhite) != 1))
                pm_usage(usage);
            argn++;
            if ((argn == argc) ||
                (sscanf(argv[argn], "%lf", &Customsystem.yWhite) != 1))
                pm_usage(usage);
        } else {
            pm_usage(usage);
        }
        argn++;
    }

    if (argn != argc) {               /* Extra bogus arguments ? */
        pm_usage(usage);
    }

#define Sz(x) (((x) * MIN(pixcols, pixrows)) / 512)
#define B(x, y) ((x) + xBias), (y)
#define Bixels(y, x) pixels[y][x + xBias]

    /* Allocate image buffer and clear it to black. */

    pixels = ppm_allocarray(pixcols = sxsize, pixrows = sysize);
    PPM_ASSIGN(rgbcolour, 0, 0, 0);
    ppmd_filledrectangle(pixels, pixcols, pixrows, Maxval, 0, 0,
         pixcols, pixrows, PPMD_NULLDRAWPROC,
         (char *) &rgbcolour);

    /* Partition into plot area and axes and establish subwindow. */

    xBias = Sz(32);
    yBias = Sz(20);

    pxcols = pixcols - xBias;
    pxrows = pixrows - yBias;

    /* Draw the CIE tongue outline. */

    PPM_ASSIGN(rgbcolour, Maxval, Maxval, Maxval);
    for (x = 380; x <= 700; x += 5) {
        double px, py;

        ix = (x - 380) / 5;
        px = spectral_chromaticity[ix][0];
        py = spectral_chromaticity[ix][1];
	if (upvp) {
	    double up, vp;
	    xy_to_upvp(px, py, &up, &vp);
	    icx = up * (pxcols - 1);
	    icy = (pxrows - 1) - vp * (pxrows - 1);
	} else {
	    icx = px * (pxcols - 1);
	    icy = (pxrows - 1) - py * (pxrows - 1);
	}
        if (x > 380) {
            ppmd_line(pixels, pixcols, pixrows, Maxval,
                  B(lx, ly), B(icx, icy),
                  PPMD_NULLDRAWPROC, (char *) &rgbcolour);
        } else {
            fx = icx;
            fy = icy;
        }
        lx = icx;
        ly = icy;
    }
    ppmd_line(pixels, pixcols, pixrows, Maxval,
          B(lx, ly), B(fx, fy),
          PPMD_NULLDRAWPROC, (char *) &rgbcolour);

    /* Now scan the image line by line and  fill  the  tongue  outline
       with the RGB values determined by the colour system for the x-y
       co-ordinates within the tongue. */

    for (y = 0; y < pxrows; y++) {
        int xe;

        /* Find horizontal extents of tongue on this line. */

        for (x = 0; x < pxcols; x++) {
            if (PPM_GETR(Bixels(y, x)) != 0) {

                for (xe = pxcols - 1; xe >= x; xe--) {
                    if (PPM_GETR(Bixels(y, xe)) != 0) {
                        break;
                    }
                }
                break;
            }
        }
        if (x < pxcols) {
            for ( ; x <= xe; x++) {
                double cx, cy, cz, jr, jg, jb, jmax;
                int r, g, b, mx;

		if (upvp) {
		    double up, vp;
		    up = ((double) x) / (pxcols - 1);
		    vp = 1.0 - ((double) y) / (pxrows - 1);
		    upvp_to_xy(up, vp, &cx, &cy);
		    cz = 1.0 - (cx + cy);
                } else {
		    cx = ((double) x) / (pxcols - 1);
		    cy = 1.0 - ((double) y) / (pxrows - 1);
		    cz = 1.0 - (cx + cy);
		}

                xyz_to_rgb(cs, cx, cy, cz, &jr, &jg, &jb);

                mx = Maxval;

                /* Check whether the requested colour  is  within  the
                   gamut  achievable with the given colour system.  If
                   not, draw it in a reduced  intensity,  interpolated
                   by desaturation to the closest within-gamut colour. */

                if (constrain_rgb(&jr, &jg, &jb)) {
                    mx = fullChart ? Maxval :
                                     ((Maxval + 1) * 3) / 4;
                }
		/* Scale to max(rgb) = 1. */
		jmax = MAX(jr, MAX(jg, jb));
                if (jmax > 0) {
		    jr = jr / jmax;
		    jg = jg / jmax;
		    jb = jb / jmax;
		}
		/* gamma correct from linear rgb to nonlinear rgb. */
                gamma_correct_rgb(cs, &jr, &jg, &jb);
		r = mx * jr;
		g = mx * jg;
		b = mx * jb;
		PPM_ASSIGN(Bixels(y, x), (pixval) r, (pixval) g, (pixval) b);
            }
        }
    }

    /* Draw the axes and bounding line. */

    if (showAxes) {
	PPM_ASSIGN(rgbcolour, Maxval, Maxval, Maxval);
	ppmd_line(pixels, pixcols, pixrows, Maxval,
	      B(0, 0), B(0, pxrows - 1), PPMD_NULLDRAWPROC, (char *) &rgbcolour);
	ppmd_line(pixels, pixcols, pixrows, Maxval,
              B(0, pxrows - 1), B(pxcols - 1, pxrows - 1),
              PPMD_NULLDRAWPROC, (char *) &rgbcolour);

	/* Draw tick marks on X and Y axes every 0.1 units.  Also label axes. */

	for (y = 1; y <= 9; y += 1) {
	    char s[20];

	    /* X axis ticks */

	    sprintf(s, "0.%d", y);
	    ppmd_line(pixels, pixcols, pixrows, Maxval,
                  B((y * (pxcols - 1)) / 10, pxrows - Sz(1)),
                  B((y * (pxcols - 1)) / 10, pxrows - Sz(4)),
                  PPMD_NULLDRAWPROC, (char *) &rgbcolour);
	    ppmd_text(pixels, pixcols, pixrows, Maxval,
                  B((y * (pxcols - 1)) / 10 - Sz(11), pxrows + Sz(12)),
                  Sz(10), 0, s,
                  PPMD_NULLDRAWPROC, (char *) &rgbcolour);

	    /* Y axis ticks */

	    sprintf(s, "0.%d", 10 - y);
	    ppmd_line(pixels, pixcols, pixrows, Maxval,
                  B(0, (y * (pxrows - 1)) / 10),
                  B(Sz(3), (y * (pxrows - 1)) / 10),
                  PPMD_NULLDRAWPROC, (char *) &rgbcolour);
	    ppmd_text(pixels, pixcols, pixrows, Maxval,
                  B(Sz(-30), (y * (pxrows - 1)) / 10 + Sz(5)),
                  Sz(10), 0, s,
                  PPMD_NULLDRAWPROC, (char *) &rgbcolour);
	}
	ppmd_text(pixels, pixcols, pixrows, Maxval,
              B((98 * (pxcols - 1)) / 100 - Sz(11), pxrows + Sz(12)),
              Sz(10), 0, (upvp ? "u'" : "x"),
              PPMD_NULLDRAWPROC, (char *) &rgbcolour);
	ppmd_text(pixels, pixcols, pixrows, Maxval,
              B(Sz(-22), (2 * (pxrows - 1)) / 100 + Sz(5)),
              Sz(10), 0, (upvp ? "v'" : "y"),
              PPMD_NULLDRAWPROC, (char *) &rgbcolour);
    }

    /* Plot the white point of the chosen colour system. */

    if (showWhite) {
        int wx, wy;

	if (upvp) {
	    double wup, wvp;
	    xy_to_upvp(cs->xWhite, cs->yWhite, &wup, &wvp);
	    wx = wup;
	    wy = wvp;
	    wx = (pxcols - 1) * wup;
            wy = (pxrows - 1) - ((int) ((pxrows - 1) * wvp));
	} else {
	    wx = (pxcols - 1) * cs->xWhite;
            wy = (pxrows - 1) - ((int) ((pxrows - 1) * cs->yWhite));
	}

        PPM_ASSIGN(rgbcolour, 0, 0, 0);
        /* We draw the four arms of the cross separately so as to
           leave the pixel representing the precise white point
           undisturbed. */
        ppmd_line(pixels, pixcols, pixrows, Maxval,
              B(wx + Sz(1), wy), B(wx + Sz(3), wy),
              PPMD_NULLDRAWPROC, (char *) &rgbcolour);
        ppmd_line(pixels, pixcols, pixrows, Maxval,
              B(wx - Sz(1), wy), B(wx - Sz(3), wy),
              PPMD_NULLDRAWPROC, (char *) &rgbcolour);
        ppmd_line(pixels, pixcols, pixrows, Maxval,
              B(wx, wy + Sz(1)), B(wx, wy + Sz(3)),
              PPMD_NULLDRAWPROC, (char *) &rgbcolour);
        ppmd_line(pixels, pixcols, pixrows, Maxval,
              B(wx, wy - Sz(1)), B(wx, wy - Sz(3)),
              PPMD_NULLDRAWPROC, (char *) &rgbcolour);
    }

    /* Plot the black body curve. */

    if (showBlack) {
        double t;

        PPM_ASSIGN(rgbcolour, 0, 0, 0);

        /* Plot black body curve from 1000 to 30000 degrees Kelvin. */

        for (t = 1000.0; t < 30000.0; t += 50.0) {
            double lambda, X = 0, Y = 0, Z = 0;
            int xb, yb;

            /* Determine X, Y, and Z for blackbody by summing colour
               match functions over the visual range. */

            for (ix = 0, lambda = 380; lambda <= 780.0; ix++, lambda += 5) {
                double Me;

                /* Evaluate Planck's black body equation for the
                   power at this wavelength. */

                Me = 3.74183e-16 * pow(lambda * 1e-9, -5.0) /
                        (exp(1.4388e-2/(lambda * 1e-9 * t)) - 1.0);
                X += Me * cie_colour_match[ix][0];
                Y += Me * cie_colour_match[ix][1];
                Z += Me * cie_colour_match[ix][2];
            }
	    if (upvp) {
		double up, vp;
		xy_to_upvp(X / (X + Y + Z), Y / (X + Y + Z), &up, &vp);
		xb = (pxcols - 1) * up;
		yb = (pxrows - 1) - ((pxrows - 1) * vp);
	    } else {
		xb = (pxcols - 1) * X / (X + Y + Z);
		yb = (pxrows - 1) - ((pxrows - 1) * Y / (X + Y + Z));
	    }

            if (t > 1000) {
                ppmd_line(pixels, pixcols, pixrows, Maxval,
                    B(lx, ly), B(xb, yb), PPMD_NULLDRAWPROC, (char *) &rgbcolour);

                /* Draw tick mark every 1000 degrees Kelvin */

                if ((((int) t) % 1000) == 0) {
                    ppmd_line(pixels, pixcols, pixrows, Maxval,
                        B(lx, ly - Sz(2)), B(lx, ly + Sz(2)),
                        PPMD_NULLDRAWPROC, (char *) &rgbcolour);

                    /* Label selected tick marks with decreasing density. */

                    if (t <= 8000.1 || (t > 8000.0 && ((((int) t) % 5000) == 0) && t != 20000.0)) {
                        char bb[20];

                        sprintf(bb, "%g", t);
                        ppmd_text(pixels, pixcols, pixrows, Maxval,
                            B(lx - Sz(12), ly - Sz(4)), Sz(6), 0, bb,
                            PPMD_NULLDRAWPROC, (char *) &rgbcolour);
                    }
  
                }
            }
            lx = xb;
            ly = yb;
        }
    }

    /* Plot wavelengths around periphery of the tongue. */

    if (showAxes) {
	/*for (x = 450; x <= 650; x += (x > 470 && x < 600) ? 5 : 10) {*/
	for (x = (upvp? 420 : 450); x <= 650;
	    x += (upvp? 10 : (x > 470 && x < 600) ? 5 : 10)) {
	    double px, py, cx, cy, cz, jr, jg, jb, jmax;
	    char wl[20];
	    int bx = 0, by = 0, tx, ty, r, g, b;

	    /* Ick.  Drop legends that overlap and twiddle position
	       so they appear at reasonable positions with respect to
	       the tongue. */

	    if ((upvp && (x == 430 || x == 640))
	        || (!upvp && (x == 460 || x == 630 || x == 640))) {
		continue;
	    }
	    if (x < 520) {
		bx = Sz(-22);
		by = Sz(2);
	    } else if (x < 535) {
		bx = Sz(-8);
		by = Sz(-6);
	    } else {
		bx = Sz(4);
	    }

	    ix = (x - 380) / 5;
	    px = spectral_chromaticity[ix][0];
	    py = spectral_chromaticity[ix][1];
	    if (upvp) {
		double up, vp;
		xy_to_upvp(px, py, &up, &vp);
		icx = up * (pxcols - 1);
		icy = (pxrows - 1) - vp * (pxrows - 1);
	    } else {
		icx = px * (pxcols - 1);
		icy = (pxrows - 1) - py * (pxrows - 1);
	    }

	    PPM_ASSIGN(rgbcolour, Maxval, Maxval, Maxval);
	    tx = icx + ((x < 520) ? Sz(-2) : ((x >= 535) ? Sz(2) : 0));
	    ty = icy + ((x < 520) ? 0 : ((x >= 535) ? Sz(-1) : Sz(-2))); 
	    ppmd_line(pixels, pixcols, pixrows, Maxval,
		 B(icx, icy), B(tx, ty),
		 PPMD_NULLDRAWPROC, (char *) &rgbcolour);

	    /* The following flailing about sets the drawing colour to
	       the hue corresponding to the pure wavelength (constrained
	       to the display gamut). */

	    if (upvp) {
		double up, vp;
		up = ((double) icx) / (pxcols - 1);
		vp = 1.0 - ((double) icy) / (pxrows - 1);
		upvp_to_xy(up, vp, &cx, &cy);
		cz = 1.0 - (cx + cy);
	    } else {
		cx = ((double) icx) / (pxcols - 1);
		cy = 1.0 - ((double) icy) / (pxrows - 1);
		cz = 1.0 - (cx + cy);
	    }

	    xyz_to_rgb(cs, cx, cy, cz, &jr, &jg, &jb);
	    (void) constrain_rgb(&jr, &jg, &jb);

	    /* Scale to max(rgb) = 1 */
	    jmax = MAX(jr, MAX(jg, jb));
	    if (jmax > 0) {
		jr = jr / jmax;
		jg = jg / jmax;
		jb = jb / jmax;
	    }
	    /* gamma correct from linear rgb to nonlinear rgb. */
	    gamma_correct_rgb(cs, &jr, &jg, &jb);
	    r = Maxval * jr;
	    g = Maxval * jg;
	    b = Maxval * jb;
	    PPM_ASSIGN(rgbcolour, (pixval) r, (pixval) g, (pixval) b);

	    sprintf(wl, "%d", x);
	    ppmd_text(pixels, pixcols, pixrows, Maxval,
		B(icx + bx, icy + by), Sz(6), 0, wl,
		PPMD_NULLDRAWPROC, (char *) &rgbcolour);
	}
    }

    /* Place description on chart. */

    if (showLabel) {
	PPM_ASSIGN(rgbcolour, Maxval, Maxval, Maxval);
	sprintf(sysdesc, "System: %s\nPrimary illuminants (X, Y)\n\
          Red:  %0.4f, %0.4f\n\
          Green: %0.4f, %0.4f\n\
          Blue:  %0.4f, %0.4f\n\
White point (X, Y): %0.4f, %0.4f",
	    cs->name, cs->xRed, cs->yRed, cs->xGreen, cs->yGreen,
	    cs->xBlue, cs->yBlue, cs->xWhite, cs->yWhite);
	ppmd_text(pixels, pixcols, pixrows, Maxval,
	    pixcols / 3, Sz(24), Sz(12), 0, sysdesc,
	    PPMD_NULLDRAWPROC, (char *) &rgbcolour);
    }

    ppm_writeppm(stdout, pixels, pixcols, pixrows, Maxval, FALSE);
    return 0;
}
