src/common/ccolor.c

#ifndef lint
static const char RCSid[] = "$Id: ccolor.c,v 3.3 2012/05/17 17:10:23 greg Exp $";
#endif
/*
 * Spectral color handling routines
 */

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

                                        /* 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}
};

C_COLOR         c_dfcolor = C_DEFCOLOR;

                                /* CIE 1931 Standard Observer curves */
static const C_COLOR    cie_xf = { 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
                        };
static const C_COLOR    cie_yf = { 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
                        };
static const C_COLOR    cie_zf = { 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[0] = XYZfromSharp[0][0]*cin[0] + XYZfromSharp[0][1]*cin[1] +
                                XYZfromSharp[0][2]*cin[2];
        xyz[1] = XYZfromSharp[1][0]*cin[0] + XYZfromSharp[1][1]*cin[1] +
                                XYZfromSharp[1][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 */
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;
        }
        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] * cie_yf.ssamp[i];
        }
        if (scale <= 1e-7)
                return(0.);
        yval /= (double)cie_yf.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] + .5;
                        else            /* interpolate if necessary */
                                clr->ssamp[i] = .5 + 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 += cie_xf.ssamp[i] * clr->ssamp[i];
                        y += cie_yf.ssamp[i] * clr->ssamp[i];
                        z += cie_zf.ssamp[i] * clr->ssamp[i];
                }
                x /= (double)cie_xf.ssum;
                y /= (double)cie_yf.ssum;
                z /= (double)cie_zf.ssum;
                z += x + y;
                clr->cx = x / z;
                clr->cy = y / z;
                clr->flags |= C_CSXY;
        } else 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] + .5;
                        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 += cie_yf.ssamp[i] * clr->ssamp[i];
                        clr->eff = C_CLPWM * y / clr->ssum;
                } else /* clr->flags & C_CSXY */ {      /* from (x,y) */
                        clr->eff = clr->cx*cie_xf.eff + clr->cy*cie_yf.eff +
                                        (1. - clr->cx - clr->cy)*cie_zf.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 ((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);
                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] + .5;
                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);                  /* nasty, but true */
                yres = y1 * y2 * cie_yf.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) + .5;
        }
        clr->flags = C_CDSPEC|C_CSSPEC;
        return(1);
}

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