src/common/ccolor.c

#ifndef lint
static const char RCSid[] = "$Id: ccolor.c,v 3.10 2017/04/05 00:54:50 greg Exp $";
#endif
/*
 * Spectral color handling routines
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "ccolor.h"

#undef frand
#define frand() (rand()*(1./(RAND_MAX+.5)))

                                        /* Sharp primary matrix */
float   XYZtoSharp[3][3] = {
        { 1.2694, -0.0988, -0.1706},
        {-0.8364,  1.8006,  0.0357},
        { 0.0297, -0.0315,  1.0018}
};
                                        /* inverse Sharp primary matrix */
float   XYZfromSharp[3][3] = {
        { 0.8156,  0.0472,  0.1372},
        { 0.3791,  0.5769,  0.0440},
        {-0.0123,  0.0167,  0.9955}
};

const C_COLOR   c_dfcolor = C_DEFCOLOR;

const C_CHROMA  c_dfchroma = 49750;     /* c_encodeChroma(&c_dfcolor) */

                                /* CIE 1931 Standard Observer curves */
const C_COLOR   c_x31 = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
                        {14,42,143,435,1344,2839,3483,3362,2908,1954,956,
                        320,49,93,633,1655,2904,4334,5945,7621,9163,10263,
                        10622,10026,8544,6424,4479,2835,1649,874,468,227,
                        114,58,29,14,7,3,2,1,0}, 106836L, .467, .368, 362.230
                        };
const C_COLOR   c_y31 = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
                        {0,1,4,12,40,116,230,380,600,910,1390,2080,3230,
                        5030,7100,8620,9540,9950,9950,9520,8700,7570,6310,
                        5030,3810,2650,1750,1070,610,320,170,82,41,21,10,
                        5,2,1,1,0,0}, 106856L, .398, .542, 493.525
                        };
const C_COLOR   c_z31 = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
                        {65,201,679,2074,6456,13856,17471,17721,16692,
                        12876,8130,4652,2720,1582,782,422,203,87,39,21,17,
                        11,8,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
                        106770L, .147, .077, 54.363
                        };
                                /* Derived CIE 1931 Primaries (imaginary) */
static const C_COLOR    cie_xp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
                        {-174,-198,-195,-197,-202,-213,-235,-272,-333,
                        -444,-688,-1232,-2393,-4497,-6876,-6758,-5256,
                        -3100,-815,1320,3200,4782,5998,6861,7408,7754,
                        7980,8120,8199,8240,8271,8292,8309,8283,8469,
                        8336,8336,8336,8336,8336,8336},
                        127424L, 1., .0,
                        };
static const C_COLOR    cie_yp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
                        {-451,-431,-431,-430,-427,-417,-399,-366,-312,
                        -204,57,691,2142,4990,8810,9871,9122,7321,5145,
                        3023,1123,-473,-1704,-2572,-3127,-3474,-3704,
                        -3846,-3927,-3968,-3999,-4021,-4038,-4012,-4201,
                        -4066,-4066,-4066,-4066,-4066,-4066},
                        -23035L, .0, 1.,
                        };
static const C_COLOR    cie_zp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
                        {4051,4054,4052,4053,4054,4056,4059,4064,4071,
                        4074,4056,3967,3677,2933,1492,313,-440,-795,
                        -904,-918,-898,-884,-869,-863,-855,-855,-851,
                        -848,-847,-846,-846,-846,-845,-846,-843,-845,
                        -845,-845,-845,-845,-845},
                        36057L, .0, .0,
                        };


/* convert to sharpened RGB color for low-error operations */
void
c_toSharpRGB(C_COLOR *cin, double cieY, float cout[3])
{
        double  xyz[3];

        c_ccvt(cin, C_CSXY);

        xyz[0] = cin->cx/cin->cy * cieY;
        xyz[1] = cieY;
        xyz[2] = (1. - cin->cx - cin->cy)/cin->cy * cieY;

        cout[0] = XYZtoSharp[0][0]*xyz[0] + XYZtoSharp[0][1]*xyz[1] +
                                XYZtoSharp[0][2]*xyz[2];
        cout[1] = XYZtoSharp[1][0]*xyz[0] + XYZtoSharp[1][1]*xyz[1] +
                                XYZtoSharp[1][2]*xyz[2];
        cout[2] = XYZtoSharp[2][0]*xyz[0] + XYZtoSharp[2][1]*xyz[1] +
                                XYZtoSharp[2][2]*xyz[2];
}

/* convert back from sharpened RGB color */
double
c_fromSharpRGB(float cin[3], C_COLOR *cout)
{
        double  xyz[3], sf;

        xyz[1] = XYZfromSharp[1][0]*cin[0] + XYZfromSharp[1][1]*cin[1] +
                                XYZfromSharp[1][2]*cin[2];
        if (xyz[1] <= 1e-6) {
                *cout = c_dfcolor;      /* punting, here... */
                return xyz[1];
        }
        xyz[0] = XYZfromSharp[0][0]*cin[0] + XYZfromSharp[0][1]*cin[1] +
                                XYZfromSharp[0][2]*cin[2];
        xyz[2] = XYZfromSharp[2][0]*cin[0] + XYZfromSharp[2][1]*cin[1] +
                                XYZfromSharp[2][2]*cin[2];
        
        sf = 1./(xyz[0] + xyz[1] + xyz[2]);

        cout->cx = xyz[0] * sf;
        cout->cy = xyz[1] * sf;
        cout->flags = C_CDXY|C_CSXY;
        
        return(xyz[1]);
}

/* assign arbitrary spectrum and return Y value */
double
c_sset(C_COLOR *clr, double wlmin, double wlmax, const float spec[], int nwl)
{
        double  yval, scale;
        float   va[C_CNSS];
        int     i, pos, n, imax, wl;
        double  wl0, wlstep;
        double  boxpos, boxstep;
                                        /* check arguments */
        if ((nwl <= 1) | (spec == NULL) | (wlmin >= C_CMAXWL) |
                        (wlmax <= C_CMINWL) | (wlmin >= wlmax))
                return(0.);
        wlstep = (wlmax - wlmin)/(nwl-1);
        while (wlmin < C_CMINWL) {
                wlmin += wlstep;
                --nwl; ++spec;
        }
        while (wlmax > C_CMAXWL) {
                wlmax -= wlstep;
                --nwl;
        }
        if ((nwl <= 1) | (wlmin >= wlmax))
                return(0.);
        imax = nwl;                     /* box filter if necessary */
        boxpos = 0;
        boxstep = 1;
        if (wlstep < C_CWLI) {
                imax = (wlmax - wlmin)/C_CWLI + 1e-7;
                boxpos = (wlmin - C_CMINWL)/C_CWLI;
                boxstep = wlstep/C_CWLI;
                wlstep = C_CWLI;
        }
        scale = 0.;                     /* get values and maximum */
        yval = 0.;
        pos = 0;
        for (i = 0; i < imax; i++) {
                va[i] = 0.; n = 0;
                while (boxpos < i+.5 && pos < nwl) {
                        va[i] += spec[pos++];
                        n++;
                        boxpos += boxstep;
                }
                if (n > 1)
                        va[i] /= (double)n;
                if (va[i] > scale)
                        scale = va[i];
                else if (va[i] < -scale)
                        scale = -va[i];
                yval += va[i] * c_y31.ssamp[i];
        }
        if (scale <= 1e-7)
                return(0.);
        yval /= (double)c_y31.ssum;
        scale = C_CMAXV / scale;
        clr->ssum = 0;                  /* convert to our spacing */
        wl0 = wlmin;
        pos = 0;
        for (i = 0, wl = C_CMINWL; i < C_CNSS; i++, wl += C_CWLI)
                if ((wl < wlmin) | (wl > wlmax))
                        clr->ssamp[i] = 0;
                else {
                        while (wl0 + wlstep < wl+1e-7) {
                                wl0 += wlstep;
                                pos++;
                        }
                        if ((wl+1e-7 >= wl0) & (wl-1e-7 <= wl0))
                                clr->ssamp[i] = scale*va[pos] + frand();
                        else            /* interpolate if necessary */
                                clr->ssamp[i] = frand() + scale / wlstep *
                                                ( va[pos]*(wl0+wlstep - wl) +
                                                        va[pos+1]*(wl - wl0) );
                        clr->ssum += clr->ssamp[i];
                }
        clr->flags = C_CDSPEC|C_CSSPEC;
        return(yval);
}

/* check if color is grey */
int
c_isgrey(C_COLOR *clr)
{
        if (!(clr->flags & (C_CSXY|C_CSSPEC)))
                return(1);              /* no settings == grey */

        c_ccvt(clr, C_CSXY);

        return(clr->cx >= .323 && clr->cx <= .343 &&
                        clr->cy >= .323 && clr->cy <= .343);
}

/* convert color representations */
void
c_ccvt(C_COLOR *clr, int fl)
{
        double  x, y, z;
        int     i;

        fl &= ~clr->flags;                      /* ignore what's done */
        if (!fl)                                /* everything's done! */
                return;
        if (!(clr->flags & (C_CSXY|C_CSSPEC))) {
                *clr = c_dfcolor;               /* nothing was set! */
                return;
        }
        if (fl & C_CSXY) {              /* cspec -> cxy */
                x = y = z = 0.;
                for (i = 0; i < C_CNSS; i++) {
                        x += c_x31.ssamp[i] * clr->ssamp[i];
                        y += c_y31.ssamp[i] * clr->ssamp[i];
                        z += c_z31.ssamp[i] * clr->ssamp[i];
                }
                x /= (double)c_x31.ssum;
                y /= (double)c_y31.ssum;
                z /= (double)c_z31.ssum;
                z += x + y;
                if (z > 1e-6) {
                        clr->cx = x / z;
                        clr->cy = y / z;
                } else
                        clr->cx = clr->cy = 1./3.;
                clr->flags |= C_CSXY;
        }
        if (fl & C_CSSPEC) {            /* cxy -> cspec */
                x = clr->cx;
                y = clr->cy;
                z = 1. - x - y;
                clr->ssum = 0;
                for (i = 0; i < C_CNSS; i++) {
                        clr->ssamp[i] = x*cie_xp.ssamp[i] + y*cie_yp.ssamp[i]
                                        + z*cie_zp.ssamp[i] + frand();
                        if (clr->ssamp[i] < 0)          /* out of gamut! */
                                clr->ssamp[i] = 0;
                        else
                                clr->ssum += clr->ssamp[i];
                }
                clr->flags |= C_CSSPEC;
        }
        if (fl & C_CSEFF) {             /* compute efficacy */
                if (clr->flags & C_CSSPEC) {            /* from spectrum */
                        y = 0.;
                        for (i = 0; i < C_CNSS; i++)
                                y += c_y31.ssamp[i] * clr->ssamp[i];
                        clr->eff = C_CLPWM * y / (clr->ssum + 0.0001);
                } else /* clr->flags & C_CSXY */ {      /* from (x,y) */
                        clr->eff = clr->cx*c_x31.eff + clr->cy*c_y31.eff +
                                        (1. - clr->cx - clr->cy)*c_z31.eff;
                }
                clr->flags |= C_CSEFF;
        }
}

/* mix two colors according to weights given -- cres may be c1 or c2 */
void
c_cmix(C_COLOR *cres, double w1, C_COLOR *c1, double w2, C_COLOR *c2)
{
        double  scale;
        int     i;

        if (w1 == 0) {
                *cres = *c2;
                return;
        }
        if (w2 == 0) {
                *cres = *c1;
                return;
        }
        if ((c1->flags|c2->flags) & C_CDSPEC) {         /* spectral mixing */
                float   cmix[C_CNSS];

                c_ccvt(c1, C_CSSPEC|C_CSEFF);
                c_ccvt(c2, C_CSSPEC|C_CSEFF);
                if (c1->ssum*c2->ssum == 0) {
                        *cres = c_dfcolor;
                        return;
                }
                w1 /= c1->eff*c1->ssum;
                w2 /= c2->eff*c2->ssum;
                scale = 0.;
                for (i = 0; i < C_CNSS; i++) {
                        cmix[i] = w1*c1->ssamp[i] + w2*c2->ssamp[i];
                        if (cmix[i] > scale)
                                scale = cmix[i];
                        else if (cmix[i] < -scale)
                                scale = -cmix[i];
                }
                scale = C_CMAXV / scale;
                cres->ssum = 0;
                for (i = 0; i < C_CNSS; i++)
                        cres->ssum += cres->ssamp[i] = scale*cmix[i] + frand();
                cres->flags = C_CDSPEC|C_CSSPEC;
        } else {                                        /* CIE xy mixing */
                c_ccvt(c1, C_CSXY);
                c_ccvt(c2, C_CSXY);
                w1 *= c2->cy;
                w2 *= c1->cy;
                scale = w1 + w2;
                if (scale == 0.) {
                        *cres = c_dfcolor;
                        return;
                }
                scale = 1. / scale;
                cres->cx = (w1*c1->cx + w2*c2->cx) * scale;
                cres->cy = (w1*c1->cy + w2*c2->cy) * scale;
                cres->flags = C_CDXY|C_CSXY;
        }
}

/* multiply two colors -- cres may be c1 or c2 */
double
c_cmult(C_COLOR *cres, C_COLOR *c1, double y1, C_COLOR *c2, double y2)
{
        double  yres;
        int     i;

        if ((c1->flags|c2->flags) & C_CDSPEC) {
                long    cmix[C_CNSS], cmax = 0;         /* spectral multiply */

                c_ccvt(c1, C_CSSPEC|C_CSEFF);
                c_ccvt(c2, C_CSSPEC|C_CSEFF);
                for (i = 0; i < C_CNSS; i++) {
                        cmix[i] = c1->ssamp[i] * c2->ssamp[i];
                        if (cmix[i] > cmax)
                                cmax = cmix[i];
                        else if (cmix[i] < -cmax)
                                cmax = -cmix[i];
                }
                cmax /= C_CMAXV;
                if (!cmax) {
                        *cres = c_dfcolor;
                        return(0.);
                }
                cres->ssum = 0;
                for (i = 0; i < C_CNSS; i++)
                        cres->ssum += cres->ssamp[i] = cmix[i] / cmax;
                cres->flags = C_CDSPEC|C_CSSPEC;

                c_ccvt(cres, C_CSEFF);                  /* below is correct */
                yres = y1 * y2 * c_y31.ssum * C_CLPWM /
                        (c1->eff*c1->ssum * c2->eff*c2->ssum) *
                        cres->eff*( cres->ssum*(double)cmax +
                                                C_CNSS/2.0*(cmax-1) );
        } else {
                float   rgb1[3], rgb2[3];               /* CIE xy multiply */

                c_toSharpRGB(c1, y1, rgb1);
                c_toSharpRGB(c2, y2, rgb2);
                
                rgb2[0] *= rgb1[0]; rgb2[1] *= rgb1[1]; rgb2[2] *= rgb1[2];

                yres = c_fromSharpRGB(rgb2, cres);
        }
        return(yres);
}

#define C1              3.741832e-16    /* W-m^2 */
#define C2              1.4388e-2       /* m-K */

#define bbsp(l,t)       (C1/((l)*(l)*(l)*(l)*(l)*(exp(C2/((t)*(l)))-1.)))
#define bblm(t)         (C2*0.2/(t))

/* set black body spectrum */
int
c_bbtemp(C_COLOR *clr, double tk)
{
        double  sf, wl;
        int     i;

        if (tk < 1000)
                return(0);
        wl = bblm(tk);                  /* scalefactor based on peak */
        if (wl < C_CMINWL*1e-9)
                wl = C_CMINWL*1e-9;
        else if (wl > C_CMAXWL*1e-9)
                wl = C_CMAXWL*1e-9;
        sf = C_CMAXV/bbsp(wl,tk);
        clr->ssum = 0;
        for (i = 0; i < C_CNSS; i++) {
                wl = (C_CMINWL + i*C_CWLI)*1e-9;
                clr->ssum += clr->ssamp[i] = sf*bbsp(wl,tk) + frand();
        }
        clr->flags = C_CDSPEC|C_CSSPEC;
        return(1);
}

#undef  C1
#undef  C2
#undef  bbsp
#undef  bblm

#define UV_NORMF        410.

/* encode (x,y) chromaticity */
C_CHROMA
c_encodeChroma(C_COLOR *clr)
{
        double  df;
        int     ub, vb;

        c_ccvt(clr, C_CSXY);
        df = UV_NORMF/(-2.*clr->cx + 12.*clr->cy + 3.);
        ub = 4.*clr->cx*df + frand();
        if (ub > 0xff) ub = 0xff;
        else ub *= (ub > 0);
        vb = 9.*clr->cy*df + frand();
        if (vb > 0xff) vb = 0xff;
        else vb *= (vb > 0);

        return(vb<<8 | ub);
}

/* decode (x,y) chromaticity */
void
c_decodeChroma(C_COLOR *cres, C_CHROMA ccode)
{
        double  up = (ccode & 0xff)*(1./UV_NORMF);
        double  vp = (ccode>>8 & 0xff)*(1./UV_NORMF);
        double  df = 1./(6.*up - 16.*vp + 12.);

        cres->cx = 9.*up * df;
        cres->cy = 4.*vp * df;
        cres->flags = C_CDXY|C_CSXY;
}

#undef  UV_NORMF
Revision:	3.11
Committed:	Thu Apr 13 00:42:01 2017 UTC (7 years ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, rad5R2, rad5R1, rad5R3
Changes since 3.10:	+19 -4 lines
Log Message:	Added checks against divide-by-zero
#	User	Rev	Content
1	greg	3.1	#ifndef lint
2	greg	3.11	static const char RCSid[] = "$Id: ccolor.c,v 3.10 2017/04/05 00:54:50 greg Exp $";
3	greg	3.1	#endif
4			/*
5			* Spectral color handling routines
6			*/
7
8			#include <stdio.h>
9	greg	3.9	#include <stdlib.h>
10	greg	3.1	#include <math.h>
11			#include "ccolor.h"
12
13	greg	3.9	#undef frand
14			#define frand() (rand()*(1./(RAND_MAX+.5)))
15
16	greg	3.4	/* Sharp primary matrix */
17			float XYZtoSharp[3][3] = {
18			{ 1.2694, -0.0988, -0.1706},
19			{-0.8364, 1.8006, 0.0357},
20			{ 0.0297, -0.0315, 1.0018}
21			};
22			/* inverse Sharp primary matrix */
23			float XYZfromSharp[3][3] = {
24			{ 0.8156, 0.0472, 0.1372},
25			{ 0.3791, 0.5769, 0.0440},
26			{-0.0123, 0.0167, 0.9955}
27			};
28	greg	3.1
29	greg	3.5	const C_COLOR c_dfcolor = C_DEFCOLOR;
30	greg	3.1
31	greg	3.9	const C_CHROMA c_dfchroma = 49750; /* c_encodeChroma(&c_dfcolor) */
32
33	greg	3.1	/* CIE 1931 Standard Observer curves */
34	greg	3.5	const C_COLOR c_x31 = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY\|C_CSEFF,
35	greg	3.1	{14,42,143,435,1344,2839,3483,3362,2908,1954,956,
36			320,49,93,633,1655,2904,4334,5945,7621,9163,10263,
37			10622,10026,8544,6424,4479,2835,1649,874,468,227,
38			114,58,29,14,7,3,2,1,0}, 106836L, .467, .368, 362.230
39			};
40	greg	3.5	const C_COLOR c_y31 = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY\|C_CSEFF,
41	greg	3.1	{0,1,4,12,40,116,230,380,600,910,1390,2080,3230,
42			5030,7100,8620,9540,9950,9950,9520,8700,7570,6310,
43			5030,3810,2650,1750,1070,610,320,170,82,41,21,10,
44			5,2,1,1,0,0}, 106856L, .398, .542, 493.525
45			};
46	greg	3.5	const C_COLOR c_z31 = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY\|C_CSEFF,
47	greg	3.1	{65,201,679,2074,6456,13856,17471,17721,16692,
48			12876,8130,4652,2720,1582,782,422,203,87,39,21,17,
49			11,8,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
50			106770L, .147, .077, 54.363
51			};
52			/* Derived CIE 1931 Primaries (imaginary) */
53	greg	3.2	static const C_COLOR cie_xp = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY,
54	greg	3.1	{-174,-198,-195,-197,-202,-213,-235,-272,-333,
55			-444,-688,-1232,-2393,-4497,-6876,-6758,-5256,
56			-3100,-815,1320,3200,4782,5998,6861,7408,7754,
57			7980,8120,8199,8240,8271,8292,8309,8283,8469,
58			8336,8336,8336,8336,8336,8336},
59			127424L, 1., .0,
60			};
61	greg	3.2	static const C_COLOR cie_yp = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY,
62	greg	3.1	{-451,-431,-431,-430,-427,-417,-399,-366,-312,
63			-204,57,691,2142,4990,8810,9871,9122,7321,5145,
64			3023,1123,-473,-1704,-2572,-3127,-3474,-3704,
65			-3846,-3927,-3968,-3999,-4021,-4038,-4012,-4201,
66			-4066,-4066,-4066,-4066,-4066,-4066},
67			-23035L, .0, 1.,
68			};
69	greg	3.2	static const C_COLOR cie_zp = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY,
70	greg	3.1	{4051,4054,4052,4053,4054,4056,4059,4064,4071,
71			4074,4056,3967,3677,2933,1492,313,-440,-795,
72			-904,-918,-898,-884,-869,-863,-855,-855,-851,
73			-848,-847,-846,-846,-846,-845,-846,-843,-845,
74			-845,-845,-845,-845,-845},
75			36057L, .0, .0,
76			};
77
78	greg	3.2
79	greg	3.4	/* convert to sharpened RGB color for low-error operations */
80			void
81	greg	3.2	c_toSharpRGB(C_COLOR *cin, double cieY, float cout[3])
82			{
83			double xyz[3];
84
85			c_ccvt(cin, C_CSXY);
86
87			xyz[0] = cin->cx/cin->cy * cieY;
88			xyz[1] = cieY;
89			xyz[2] = (1. - cin->cx - cin->cy)/cin->cy * cieY;
90
91	greg	3.4	cout[0] = XYZtoSharp[0][0]xyz[0] + XYZtoSharp[0][1]xyz[1] +
92			XYZtoSharp[0][2]*xyz[2];
93			cout[1] = XYZtoSharp[1][0]xyz[0] + XYZtoSharp[1][1]xyz[1] +
94			XYZtoSharp[1][2]*xyz[2];
95			cout[2] = XYZtoSharp[2][0]xyz[0] + XYZtoSharp[2][1]xyz[1] +
96			XYZtoSharp[2][2]*xyz[2];
97	greg	3.2	}
98
99	greg	3.4	/* convert back from sharpened RGB color */
100			double
101	greg	3.2	c_fromSharpRGB(float cin[3], C_COLOR *cout)
102			{
103			double xyz[3], sf;
104
105	greg	3.8	xyz[1] = XYZfromSharp[1][0]cin[0] + XYZfromSharp[1][1]cin[1] +
106			XYZfromSharp[1][2]*cin[2];
107			if (xyz[1] <= 1e-6) {
108			cout = c_dfcolor; / punting, here... */
109			return xyz[1];
110			}
111	greg	3.4	xyz[0] = XYZfromSharp[0][0]cin[0] + XYZfromSharp[0][1]cin[1] +
112			XYZfromSharp[0][2]*cin[2];
113			xyz[2] = XYZfromSharp[2][0]cin[0] + XYZfromSharp[2][1]cin[1] +
114			XYZfromSharp[2][2]*cin[2];
115	greg	3.8
116	greg	3.2	sf = 1./(xyz[0] + xyz[1] + xyz[2]);
117
118			cout->cx = xyz[0] * sf;
119			cout->cy = xyz[1] * sf;
120			cout->flags = C_CDXY\|C_CSXY;
121
122			return(xyz[1]);
123			}
124
125	greg	3.6	/* assign arbitrary spectrum and return Y value */
126	greg	3.3	double
127			c_sset(C_COLOR *clr, double wlmin, double wlmax, const float spec[], int nwl)
128			{
129			double yval, scale;
130			float va[C_CNSS];
131			int i, pos, n, imax, wl;
132			double wl0, wlstep;
133			double boxpos, boxstep;
134			/* check arguments */
135			if ((nwl <= 1) \| (spec == NULL) \| (wlmin >= C_CMAXWL) \|
136			(wlmax <= C_CMINWL) \| (wlmin >= wlmax))
137			return(0.);
138			wlstep = (wlmax - wlmin)/(nwl-1);
139			while (wlmin < C_CMINWL) {
140			wlmin += wlstep;
141			--nwl; ++spec;
142			}
143			while (wlmax > C_CMAXWL) {
144			wlmax -= wlstep;
145			--nwl;
146			}
147	greg	3.6	if ((nwl <= 1) \| (wlmin >= wlmax))
148			return(0.);
149	greg	3.3	imax = nwl; /* box filter if necessary */
150			boxpos = 0;
151			boxstep = 1;
152			if (wlstep < C_CWLI) {
153			imax = (wlmax - wlmin)/C_CWLI + 1e-7;
154			boxpos = (wlmin - C_CMINWL)/C_CWLI;
155			boxstep = wlstep/C_CWLI;
156			wlstep = C_CWLI;
157			}
158			scale = 0.; /* get values and maximum */
159			yval = 0.;
160			pos = 0;
161			for (i = 0; i < imax; i++) {
162			va[i] = 0.; n = 0;
163			while (boxpos < i+.5 && pos < nwl) {
164			va[i] += spec[pos++];
165			n++;
166			boxpos += boxstep;
167			}
168			if (n > 1)
169			va[i] /= (double)n;
170			if (va[i] > scale)
171			scale = va[i];
172			else if (va[i] < -scale)
173			scale = -va[i];
174	greg	3.5	yval += va[i] * c_y31.ssamp[i];
175	greg	3.3	}
176			if (scale <= 1e-7)
177			return(0.);
178	greg	3.5	yval /= (double)c_y31.ssum;
179	greg	3.3	scale = C_CMAXV / scale;
180			clr->ssum = 0; /* convert to our spacing */
181			wl0 = wlmin;
182			pos = 0;
183			for (i = 0, wl = C_CMINWL; i < C_CNSS; i++, wl += C_CWLI)
184			if ((wl < wlmin) \| (wl > wlmax))
185			clr->ssamp[i] = 0;
186			else {
187			while (wl0 + wlstep < wl+1e-7) {
188			wl0 += wlstep;
189			pos++;
190			}
191			if ((wl+1e-7 >= wl0) & (wl-1e-7 <= wl0))
192	greg	3.9	clr->ssamp[i] = scale*va[pos] + frand();
193	greg	3.3	else /* interpolate if necessary */
194	greg	3.9	clr->ssamp[i] = frand() + scale / wlstep *
195	greg	3.3	( va[pos]*(wl0+wlstep - wl) +
196			va[pos+1]*(wl - wl0) );
197			clr->ssum += clr->ssamp[i];
198			}
199			clr->flags = C_CDSPEC\|C_CSSPEC;
200			return(yval);
201			}
202
203	greg	3.1	/* check if color is grey */
204			int
205			c_isgrey(C_COLOR *clr)
206			{
207			if (!(clr->flags & (C_CSXY\|C_CSSPEC)))
208			return(1); /* no settings == grey */
209
210			c_ccvt(clr, C_CSXY);
211
212			return(clr->cx >= .323 && clr->cx <= .343 &&
213			clr->cy >= .323 && clr->cy <= .343);
214			}
215
216			/* convert color representations */
217			void
218			c_ccvt(C_COLOR *clr, int fl)
219			{
220			double x, y, z;
221			int i;
222
223			fl &= ~clr->flags; /* ignore what's done */
224			if (!fl) /* everything's done! */
225			return;
226			if (!(clr->flags & (C_CSXY\|C_CSSPEC))) {
227			clr = c_dfcolor; / nothing was set! */
228			return;
229			}
230			if (fl & C_CSXY) { /* cspec -> cxy */
231			x = y = z = 0.;
232			for (i = 0; i < C_CNSS; i++) {
233	greg	3.5	x += c_x31.ssamp[i] * clr->ssamp[i];
234			y += c_y31.ssamp[i] * clr->ssamp[i];
235			z += c_z31.ssamp[i] * clr->ssamp[i];
236	greg	3.1	}
237	greg	3.5	x /= (double)c_x31.ssum;
238			y /= (double)c_y31.ssum;
239			z /= (double)c_z31.ssum;
240	greg	3.1	z += x + y;
241	greg	3.11	if (z > 1e-6) {
242			clr->cx = x / z;
243			clr->cy = y / z;
244			} else
245			clr->cx = clr->cy = 1./3.;
246	greg	3.1	clr->flags \|= C_CSXY;
247	greg	3.6	}
248			if (fl & C_CSSPEC) { /* cxy -> cspec */
249	greg	3.1	x = clr->cx;
250			y = clr->cy;
251			z = 1. - x - y;
252			clr->ssum = 0;
253			for (i = 0; i < C_CNSS; i++) {
254			clr->ssamp[i] = xcie_xp.ssamp[i] + ycie_yp.ssamp[i]
255	greg	3.9	+ z*cie_zp.ssamp[i] + frand();
256	greg	3.1	if (clr->ssamp[i] < 0) /* out of gamut! */
257			clr->ssamp[i] = 0;
258			else
259			clr->ssum += clr->ssamp[i];
260			}
261			clr->flags \|= C_CSSPEC;
262			}
263			if (fl & C_CSEFF) { /* compute efficacy */
264			if (clr->flags & C_CSSPEC) { /* from spectrum */
265			y = 0.;
266			for (i = 0; i < C_CNSS; i++)
267	greg	3.5	y += c_y31.ssamp[i] * clr->ssamp[i];
268	greg	3.11	clr->eff = C_CLPWM * y / (clr->ssum + 0.0001);
269	greg	3.1	} else /* clr->flags & C_CSXY / { / from (x,y) */
270	greg	3.5	clr->eff = clr->cxc_x31.eff + clr->cyc_y31.eff +
271			(1. - clr->cx - clr->cy)*c_z31.eff;
272	greg	3.1	}
273			clr->flags \|= C_CSEFF;
274			}
275			}
276
277	greg	3.2	/* mix two colors according to weights given -- cres may be c1 or c2 */
278	greg	3.1	void
279			c_cmix(C_COLOR cres, double w1, C_COLOR c1, double w2, C_COLOR *c2)
280			{
281			double scale;
282			int i;
283
284	greg	3.11	if (w1 == 0) {
285			cres = c2;
286			return;
287			}
288			if (w2 == 0) {
289			cres = c1;
290			return;
291			}
292	greg	3.1	if ((c1->flags\|c2->flags) & C_CDSPEC) { /* spectral mixing */
293	greg	3.2	float cmix[C_CNSS];
294
295	greg	3.1	c_ccvt(c1, C_CSSPEC\|C_CSEFF);
296			c_ccvt(c2, C_CSSPEC\|C_CSEFF);
297	greg	3.11	if (c1->ssum*c2->ssum == 0) {
298			*cres = c_dfcolor;
299			return;
300			}
301	greg	3.1	w1 /= c1->eff*c1->ssum;
302			w2 /= c2->eff*c2->ssum;
303			scale = 0.;
304			for (i = 0; i < C_CNSS; i++) {
305			cmix[i] = w1c1->ssamp[i] + w2c2->ssamp[i];
306			if (cmix[i] > scale)
307			scale = cmix[i];
308	greg	3.3	else if (cmix[i] < -scale)
309			scale = -cmix[i];
310	greg	3.1	}
311			scale = C_CMAXV / scale;
312			cres->ssum = 0;
313			for (i = 0; i < C_CNSS; i++)
314	greg	3.9	cres->ssum += cres->ssamp[i] = scale*cmix[i] + frand();
315	greg	3.1	cres->flags = C_CDSPEC\|C_CSSPEC;
316			} else { /* CIE xy mixing */
317			c_ccvt(c1, C_CSXY);
318			c_ccvt(c2, C_CSXY);
319			w1 *= c2->cy;
320			w2 *= c1->cy;
321			scale = w1 + w2;
322			if (scale == 0.) {
323			*cres = c_dfcolor;
324			return;
325			}
326			scale = 1. / scale;
327			cres->cx = (w1c1->cx + w2c2->cx) * scale;
328			cres->cy = (w1c1->cy + w2c2->cy) * scale;
329			cres->flags = C_CDXY\|C_CSXY;
330			}
331			}
332
333	greg	3.2	/* multiply two colors -- cres may be c1 or c2 */
334			double
335			c_cmult(C_COLOR cres, C_COLOR c1, double y1, C_COLOR *c2, double y2)
336			{
337			double yres;
338			int i;
339
340			if ((c1->flags\|c2->flags) & C_CDSPEC) {
341			long cmix[C_CNSS], cmax = 0; /* spectral multiply */
342
343			c_ccvt(c1, C_CSSPEC\|C_CSEFF);
344			c_ccvt(c2, C_CSSPEC\|C_CSEFF);
345			for (i = 0; i < C_CNSS; i++) {
346			cmix[i] = c1->ssamp[i] * c2->ssamp[i];
347			if (cmix[i] > cmax)
348			cmax = cmix[i];
349	greg	3.3	else if (cmix[i] < -cmax)
350			cmax = -cmix[i];
351	greg	3.2	}
352			cmax /= C_CMAXV;
353			if (!cmax) {
354			*cres = c_dfcolor;
355			return(0.);
356			}
357			cres->ssum = 0;
358			for (i = 0; i < C_CNSS; i++)
359			cres->ssum += cres->ssamp[i] = cmix[i] / cmax;
360			cres->flags = C_CDSPEC\|C_CSSPEC;
361
362	greg	3.6	c_ccvt(cres, C_CSEFF); /* below is correct */
363	greg	3.5	yres = y1 * y2 * c_y31.ssum * C_CLPWM /
364	greg	3.2	(c1->effc1->ssum c2->effc2->ssum)
365			cres->eff( cres->ssum(double)cmax +
366			C_CNSS/2.0*(cmax-1) );
367			} else {
368			float rgb1[3], rgb2[3]; /* CIE xy multiply */
369
370			c_toSharpRGB(c1, y1, rgb1);
371			c_toSharpRGB(c2, y2, rgb2);
372
373			rgb2[0] = rgb1[0]; rgb2[1] = rgb1[1]; rgb2[2] *= rgb1[2];
374
375			yres = c_fromSharpRGB(rgb2, cres);
376			}
377			return(yres);
378			}
379
380	greg	3.1	#define C1 3.741832e-16 /* W-m^2 */
381			#define C2 1.4388e-2 /* m-K */
382
383			#define bbsp(l,t) (C1/((l)(l)(l)(l)(l)(exp(C2/((t)(l)))-1.)))
384	greg	3.2	#define bblm(t) (C2*0.2/(t))
385	greg	3.1
386			/* set black body spectrum */
387			int
388			c_bbtemp(C_COLOR *clr, double tk)
389			{
390			double sf, wl;
391			int i;
392
393			if (tk < 1000)
394			return(0);
395			wl = bblm(tk); /* scalefactor based on peak */
396			if (wl < C_CMINWL*1e-9)
397			wl = C_CMINWL*1e-9;
398			else if (wl > C_CMAXWL*1e-9)
399			wl = C_CMAXWL*1e-9;
400			sf = C_CMAXV/bbsp(wl,tk);
401			clr->ssum = 0;
402			for (i = 0; i < C_CNSS; i++) {
403			wl = (C_CMINWL + iC_CWLI)1e-9;
404	greg	3.9	clr->ssum += clr->ssamp[i] = sf*bbsp(wl,tk) + frand();
405	greg	3.1	}
406			clr->flags = C_CDSPEC\|C_CSSPEC;
407			return(1);
408			}
409
410			#undef C1
411			#undef C2
412			#undef bbsp
413			#undef bblm
414	greg	3.7
415			#define UV_NORMF 410.
416
417			/* encode (x,y) chromaticity */
418			C_CHROMA
419			c_encodeChroma(C_COLOR *clr)
420			{
421			double df;
422			int ub, vb;
423
424			c_ccvt(clr, C_CSXY);
425			df = UV_NORMF/(-2.clr->cx + 12.clr->cy + 3.);
426	greg	3.9	ub = 4.clr->cxdf + frand();
427			if (ub > 0xff) ub = 0xff;
428	greg	3.10	else ub *= (ub > 0);
429	greg	3.9	vb = 9.clr->cydf + frand();
430			if (vb > 0xff) vb = 0xff;
431	greg	3.10	else vb *= (vb > 0);
432	greg	3.7
433			return(vb<<8 \| ub);
434			}
435
436			/* decode (x,y) chromaticity */
437			void
438			c_decodeChroma(C_COLOR *cres, C_CHROMA ccode)
439			{
440	greg	3.9	double up = (ccode & 0xff)*(1./UV_NORMF);
441			double vp = (ccode>>8 & 0xff)*(1./UV_NORMF);
442	greg	3.7	double df = 1./(6.up - 16.vp + 12.);
443
444			cres->cx = 9.up df;
445			cres->cy = 4.vp df;
446			cres->flags = C_CDXY\|C_CSXY;
447			}
448
449			#undef UV_NORMF