/* Generate a PPM file representing a CIE colour gamut chart Author: John Walker 92 Caro Mio 1 Sark (via Guernsey) Channel Islands GY9 0SE GREAT BRITAIN Internet: kelvint@netcom.com Fax: ++44 481 832515 Voice: ++44 481 832613 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. */ #include /* To avoid atof() clash in math.h */ #include #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 { double xRed, yRed, /* Red primary illuminant */ xGreen, yGreen, /* Green primary illuminant */ xBlue, yBlue, /* Blue primary illuminant */ xWhite, yWhite; /* White point */ }; /* Prototypes */ static void xyz_to_rgb ARGS((struct colourSystem *cs, double xc, double yc, double zc, double *r, double *g, double *b)); static int inside_gamut ARGS((double r, double g, double b)); static int constrain_rgb ARGS((struct colourSystem *cs, double *x, double *y, double *z, double *r, double *g, double *b)); /* 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 */ static int interpwp = FALSE; /* Interpolate to white point ? */ /* Standard white point chromaticities. */ #define IlluminantC 0.3101, 0.3162 /* For NTSC television */ #define IlluminantD65 0.3127, 0.3291 /* For EBU and SMPTE */ static struct colourSystem NTSCsystem = { 0.67, 0.33, 0.21, 0.71, 0.14, 0.08, IlluminantC }, EBUsystem = { 0.64, 0.33, 0.29, 0.60, 0.15, 0.06, IlluminantD65 }, SMPTEsystem = { 0.630, 0.340, 0.310, 0.595, 0.155, 0.070, IlluminantD65 }, HDTVsystem = { 0.670, 0.330, 0.210, 0.710, 0.150, 0.060, IlluminantD65 }, CIEsystem = { 0.7355,0.2645,0.2658,0.7243,0.1669,0.0085, 0.3894,0.3324}, Customsystem = { 0.64, 0.33, 0.29, 0.60, 0.15, 0.06, IlluminantD65 }; static struct colourSystem *cs = &EBUsystem; /* 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) struct colourSystem *cs; double xc, yc, zc; double *r, *g, *b; { double xr, yr, zr, xg, yg, zg, xb, yb, zb; 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); *r = (-xg*yc*zb + xc*yg*zb + xg*yb*zc - xb*yg*zc - xc*yb*zg + xb*yc*zg) / (xr*yg*zb - xg*yr*zb - xr*yb*zg + xb*yr*zg + xg*yb*zr - xb*yg*zr); *g = (xr*yc*zb - xc*yr*zb - xr*yb*zc + xb*yr*zc + xc*yb*zr - xb*yc*zr) / (xr*yg*zb - xg*yr*zb - xr*yb*zg + xb*yr*zg + xg*yb*zr - xb*yg*zr); *b = (xr*yg*zc - xg*yr*zc - xr*yc*zg + xc*yr*zg + xg*yc*zr - xc*yg*zr) / (xr*yg*zb - xg*yr*zb - xr*yb*zg + xb*yr*zg + xg*yb*zr - xb*yg*zr); } /* 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 mixing with the white point of the colour system so as to reduce the primary with the negative weight to zero. This is equivalent to finding the intersection on the CIE diagram of a line drawn between the white point and and the requested color and the edge of the Maxwell triangle formed by the three primaries. If INTERPWP is nonzero, the white point defined by the sum of the three primary illuminants is used instead of the colour system's actual white point. While indefensible from the standpoint of colour theory, this produces much better looking charts on computer monitors, since the white points of most colour systems are well to the blue of the white defined by the R = G = B = 1. */ static int constrain_rgb(cs, x, y, z, r, g, b) struct colourSystem *cs; double *x, *y, *z, *r, *g, *b; { /* Is the contribution of one of the primaries negative ? */ if (!inside_gamut(*r, *g, *b)) { double par, wr, wg, wb, xw, yw; /* Determine the white point used in interpolating out of gamut colours. If INTERPWP is set, the colour system's white point is used. Otherwise, the white defined by an equal mix of the three illuminants is taken as the origin of the interpolation line drawn to the out of gamut colour. */ if (interpwp) { xw = cs->xWhite; yw = cs->yWhite; } else { xw = (cs->xRed + cs->xGreen + cs->xBlue) / 3; yw = (cs->yRed + cs->yGreen + cs->yBlue) / 3; } /* Yes. Find the RGB mixing weights of the white point (we assume the white point is in the gamut!). */ xyz_to_rgb(cs, xw, yw, 1 - (xw + yw), &wr, &wg, &wb); /* Find the primary with negative weight and calculate the parameter of the point on the vector from the white point to the original requested colour in RGB space. */ if (*r < *g && *r < *b) { par = wr / (wr - *r); } else if (*g < *r && *g < *b) { par = wg / (wg - *g); } else { par = wb / (wb - *b); } /* Since XYZ space is a linear transformation of RGB space, we can find the XYZ space coordinates of the point where the edge of the gamut intersects the vector from the white point to the original colour by multiplying the parameter in RGB space by the difference vector in XYZ space. */ *x = CLAMP(xw + par * (*x - xw), 0, 1); *y = CLAMP(yw + par * (*y - yw), 0, 1); *z = CLAMP(1 - (*x + *y), 0, 1); /* Now finally calculate the gamut-constrained RGB weights. */ *r = CLAMP(wr + par * (*r - wr), 0, 1); *g = CLAMP(wg + par * (*g - wg), 0, 1); *b = CLAMP(wb + par * (*b - wb), 0, 1); return 1; /* Colour modified to fit RGB gamut */ } return 0; /* Colour within RGB gamut */ } /* Main program. */ int main(argc, argv) int argc; char *argv[]; { int argn, x, y; char *usage = "[-[no]black] [-[no]wpoint] [-full] [-interpwp]\n\ [-ntsc|-ebu|-smpte|-hdtv|-cie]\n\ [-red -green -blue -white ]\n\ [-size ] [-xsize|-width ] [-ysize|-height ]"; int widspec = FALSE, hgtspec = FALSE, ix, icx, icy, fx, fy, lx, ly; pixel rgbcolour; /* Pixel used to clear pixmap */ int showWhite = TRUE; /* Show white point ? */ int showBlack = TRUE; /* Show black body curve ? */ int fullChart = FALSE; /* Fill entire tongue ? */ ppm_init(&argc, argv); argn = 1; while (argn < argc && argv[argn][0] == '-' && argv[argn][1] != '\0') { 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], "-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], "-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", 3)) { 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", 2)) { 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 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; /* Show black body curve */ } else if (pm_keymatch(argv[argn], "-nowpoint", 3)) { showWhite = FALSE; /* Show white point of colour system */ } else if (pm_keymatch(argv[argn], "-full", 1)) { fullChart = TRUE; /* Fill whole tongue full-intensity */ } else if (pm_keymatch(argv[argn], "-interpwp", 1)) { interpwp = TRUE; /* Interpolate colours to white point */ } else { pm_usage(usage); } argn++; } if (argn != argc) { /* Extra bogus arguments ? */ pm_usage(usage); } /* 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); /* 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]; icx = px * (pixcols - 1); icy = (pixrows - 1) - py * (pixrows - 1); if (x > 380) { ppmd_line(pixels, pixcols, pixrows, Maxval, lx, ly, icx, icy, PPMD_NULLDRAWPROC, (char *) &rgbcolour); } else { fx = icx; fy = icy; } lx = icx; ly = icy; } ppmd_line(pixels, pixcols, pixrows, Maxval, lx, ly, 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 < pixrows; y++) { int xe; /* Find horizontal extents of tongue on this line. */ for (x = 0; x < pixcols; x++) { if (PPM_GETR(pixels[y][x]) != 0) { for (xe = pixcols - 1; xe >= x; xe--) { if (PPM_GETR(pixels[y][xe]) != 0) { break; } } break; } } if (x < pixcols) { for ( ; x <= xe; x++) { double cx, cy, cz, jr, jg, jb, jmax; int r, g, b, mx; cx = ((double) x) / (pixcols - 1); cy = 1.0 - ((double) y) / (pixrows - 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(cs, &cx, &cy, &cz, &jr, &jg, &jb)) { mx = fullChart ? Maxval : ((Maxval + 1) * 3) / 4; } if (jr <= 1 || jg <= 1 || jb <= 1) { jmax = max(jr, max(jg, jb)); if (!(jr < 0 || jg < 0 || jb < 0)) { r = mx * (jr / jmax); g = mx * (jg / jmax); b = mx * (jb / jmax); PPM_ASSIGN(pixels[y][x], (pixval) r, (pixval) g, (pixval) b); } } } } } /* Draw the axes and bounding line. */ PPM_ASSIGN(rgbcolour, Maxval, Maxval, Maxval); ppmd_line(pixels, pixcols, pixrows, Maxval, 0, 0, 0, pixrows - 1, PPMD_NULLDRAWPROC, (char *) &rgbcolour); ppmd_line(pixels, pixcols, pixrows, Maxval, 0, pixrows - 1, pixcols - 1, pixrows - 1, PPMD_NULLDRAWPROC, (char *) &rgbcolour); ppmd_line(pixels, pixcols, pixrows, Maxval, 0, 0, pixcols - 1, pixrows - 1, PPMD_NULLDRAWPROC, (char *) &rgbcolour); /* Draw tick marks on X and Y axes every 0.2 units. */ for (y = 2; y <= 8; y += 2) { ppmd_line(pixels, pixcols, pixrows, Maxval, (y * (pixcols - 1)) / 10, pixrows - 1, (y * (pixcols - 1)) / 10, pixrows - 4, PPMD_NULLDRAWPROC, (char *) &rgbcolour); ppmd_line(pixels, pixcols, pixrows, Maxval, 0, (y * (pixrows - 1)) / 10, 3, (y * (pixrows - 1)) / 10, PPMD_NULLDRAWPROC, (char *) &rgbcolour); } /* Plot the white point of the chosen colour system. */ if (showWhite) { int wx = (pixcols - 1) * cs->xWhite, wy = (pixrows - 1) - ((int) ((pixrows - 1) * cs->yWhite)); PPM_ASSIGN(rgbcolour, 0, 0, 0); ppmd_line(pixels, pixcols, pixrows, Maxval, wx + 1, wy, wx + 3, wy, PPMD_NULLDRAWPROC, (char *) &rgbcolour); ppmd_line(pixels, pixcols, pixrows, Maxval, wx - 1, wy, wx - 3, wy, PPMD_NULLDRAWPROC, (char *) &rgbcolour); ppmd_line(pixels, pixcols, pixrows, Maxval, wx, wy + 1, wx, wy + 3, PPMD_NULLDRAWPROC, (char *) &rgbcolour); ppmd_line(pixels, pixcols, pixrows, Maxval, wx, wy - 1, wx, wy - 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]; } xb = (pixcols - 1) * X / (X + Y + Z); yb = (pixrows - 1) - ((pixrows - 1) * Y / (X + Y + Z)); if (t > 1000) { ppmd_line(pixels, pixcols, pixrows, Maxval, lx, ly, xb, yb, PPMD_NULLDRAWPROC, (char *) &rgbcolour); /* Draw tick mark every 1000 degrees Kelvin */ if ((((int) t) % 1000) == 0) { ppmd_line(pixels, pixcols, pixrows, Maxval, lx, ly - 2, lx, ly + 2, PPMD_NULLDRAWPROC, (char *) &rgbcolour); } } lx = xb; ly = yb; } } ppm_writeppm(stdout, pixels, pixcols, pixrows, Maxval, FALSE); return 0; }