src/common/ccolor.c

#ifndef lint
static const char RCSid[] = "$Id: ccolor.c,v 3.5 2012/06/26 17:59:16 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}
};

const C_COLOR   c_dfcolor = C_DEFCOLOR;

                                /* 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[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 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] + .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 += 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;
                clr->cx = x / z;
                clr->cy = y / z;
                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] + .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 += c_y31.ssamp[i] * clr->ssamp[i];
                        clr->eff = C_CLPWM * y / clr->ssum;
                } 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 ((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);                  /* 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) + .5;
        }
        clr->flags = C_CDSPEC|C_CSSPEC;
        return(1);
}

#undef  C1
#undef  C2
#undef  bbsp
#undef  bblm
Revision:	3.6
Committed:	Tue Jul 17 23:40:24 2012 UTC (11 years, 9 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R2P2, rad4R2, rad4R2P1
Changes since 3.5:	+7 -4 lines
Log Message:	Additional error-checking
#	User	Rev	Content
1	greg	3.1	#ifndef lint
2	greg	3.6	static const char RCSid[] = "$Id: ccolor.c,v 3.5 2012/06/26 17:59:16 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	greg	3.5	const C_COLOR c_dfcolor = C_DEFCOLOR;
26	greg	3.1
27			/* CIE 1931 Standard Observer curves */
28	greg	3.5	const C_COLOR c_x31 = { 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.5	const C_COLOR c_y31 = { 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.5	const C_COLOR c_z31 = { 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.6	/* assign arbitrary spectrum and return Y value */
116	greg	3.3	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	greg	3.6	if ((nwl <= 1) \| (wlmin >= wlmax))
138			return(0.);
139	greg	3.3	imax = nwl; /* box filter if necessary */
140			boxpos = 0;
141			boxstep = 1;
142			if (wlstep < C_CWLI) {
143			imax = (wlmax - wlmin)/C_CWLI + 1e-7;
144			boxpos = (wlmin - C_CMINWL)/C_CWLI;
145			boxstep = wlstep/C_CWLI;
146			wlstep = C_CWLI;
147			}
148			scale = 0.; /* get values and maximum */
149			yval = 0.;
150			pos = 0;
151			for (i = 0; i < imax; i++) {
152			va[i] = 0.; n = 0;
153			while (boxpos < i+.5 && pos < nwl) {
154			va[i] += spec[pos++];
155			n++;
156			boxpos += boxstep;
157			}
158			if (n > 1)
159			va[i] /= (double)n;
160			if (va[i] > scale)
161			scale = va[i];
162			else if (va[i] < -scale)
163			scale = -va[i];
164	greg	3.5	yval += va[i] * c_y31.ssamp[i];
165	greg	3.3	}
166			if (scale <= 1e-7)
167			return(0.);
168	greg	3.5	yval /= (double)c_y31.ssum;
169	greg	3.3	scale = C_CMAXV / scale;
170			clr->ssum = 0; /* convert to our spacing */
171			wl0 = wlmin;
172			pos = 0;
173			for (i = 0, wl = C_CMINWL; i < C_CNSS; i++, wl += C_CWLI)
174			if ((wl < wlmin) \| (wl > wlmax))
175			clr->ssamp[i] = 0;
176			else {
177			while (wl0 + wlstep < wl+1e-7) {
178			wl0 += wlstep;
179			pos++;
180			}
181			if ((wl+1e-7 >= wl0) & (wl-1e-7 <= wl0))
182			clr->ssamp[i] = scale*va[pos] + .5;
183			else /* interpolate if necessary */
184			clr->ssamp[i] = .5 + scale / wlstep *
185			( va[pos]*(wl0+wlstep - wl) +
186			va[pos+1]*(wl - wl0) );
187			clr->ssum += clr->ssamp[i];
188			}
189			clr->flags = C_CDSPEC\|C_CSSPEC;
190			return(yval);
191			}
192
193	greg	3.1	/* check if color is grey */
194			int
195			c_isgrey(C_COLOR *clr)
196			{
197			if (!(clr->flags & (C_CSXY\|C_CSSPEC)))
198			return(1); /* no settings == grey */
199
200			c_ccvt(clr, C_CSXY);
201
202			return(clr->cx >= .323 && clr->cx <= .343 &&
203			clr->cy >= .323 && clr->cy <= .343);
204			}
205
206			/* convert color representations */
207			void
208			c_ccvt(C_COLOR *clr, int fl)
209			{
210			double x, y, z;
211			int i;
212
213			fl &= ~clr->flags; /* ignore what's done */
214			if (!fl) /* everything's done! */
215			return;
216			if (!(clr->flags & (C_CSXY\|C_CSSPEC))) {
217			clr = c_dfcolor; / nothing was set! */
218			return;
219			}
220			if (fl & C_CSXY) { /* cspec -> cxy */
221			x = y = z = 0.;
222			for (i = 0; i < C_CNSS; i++) {
223	greg	3.5	x += c_x31.ssamp[i] * clr->ssamp[i];
224			y += c_y31.ssamp[i] * clr->ssamp[i];
225			z += c_z31.ssamp[i] * clr->ssamp[i];
226	greg	3.1	}
227	greg	3.5	x /= (double)c_x31.ssum;
228			y /= (double)c_y31.ssum;
229			z /= (double)c_z31.ssum;
230	greg	3.1	z += x + y;
231			clr->cx = x / z;
232			clr->cy = y / z;
233			clr->flags \|= C_CSXY;
234	greg	3.6	}
235			if (fl & C_CSSPEC) { /* cxy -> cspec */
236	greg	3.1	x = clr->cx;
237			y = clr->cy;
238			z = 1. - x - y;
239			clr->ssum = 0;
240			for (i = 0; i < C_CNSS; i++) {
241			clr->ssamp[i] = xcie_xp.ssamp[i] + ycie_yp.ssamp[i]
242			+ z*cie_zp.ssamp[i] + .5;
243			if (clr->ssamp[i] < 0) /* out of gamut! */
244			clr->ssamp[i] = 0;
245			else
246			clr->ssum += clr->ssamp[i];
247			}
248			clr->flags \|= C_CSSPEC;
249			}
250			if (fl & C_CSEFF) { /* compute efficacy */
251			if (clr->flags & C_CSSPEC) { /* from spectrum */
252			y = 0.;
253			for (i = 0; i < C_CNSS; i++)
254	greg	3.5	y += c_y31.ssamp[i] * clr->ssamp[i];
255	greg	3.1	clr->eff = C_CLPWM * y / clr->ssum;
256			} else /* clr->flags & C_CSXY / { / from (x,y) */
257	greg	3.5	clr->eff = clr->cxc_x31.eff + clr->cyc_y31.eff +
258			(1. - clr->cx - clr->cy)*c_z31.eff;
259	greg	3.1	}
260			clr->flags \|= C_CSEFF;
261			}
262			}
263
264	greg	3.2	/* mix two colors according to weights given -- cres may be c1 or c2 */
265	greg	3.1	void
266			c_cmix(C_COLOR cres, double w1, C_COLOR c1, double w2, C_COLOR *c2)
267			{
268			double scale;
269			int i;
270
271			if ((c1->flags\|c2->flags) & C_CDSPEC) { /* spectral mixing */
272	greg	3.2	float cmix[C_CNSS];
273
274	greg	3.1	c_ccvt(c1, C_CSSPEC\|C_CSEFF);
275			c_ccvt(c2, C_CSSPEC\|C_CSEFF);
276			w1 /= c1->eff*c1->ssum;
277			w2 /= c2->eff*c2->ssum;
278			scale = 0.;
279			for (i = 0; i < C_CNSS; i++) {
280			cmix[i] = w1c1->ssamp[i] + w2c2->ssamp[i];
281			if (cmix[i] > scale)
282			scale = cmix[i];
283	greg	3.3	else if (cmix[i] < -scale)
284			scale = -cmix[i];
285	greg	3.1	}
286			scale = C_CMAXV / scale;
287			cres->ssum = 0;
288			for (i = 0; i < C_CNSS; i++)
289			cres->ssum += cres->ssamp[i] = scale*cmix[i] + .5;
290			cres->flags = C_CDSPEC\|C_CSSPEC;
291			} else { /* CIE xy mixing */
292			c_ccvt(c1, C_CSXY);
293			c_ccvt(c2, C_CSXY);
294			w1 *= c2->cy;
295			w2 *= c1->cy;
296			scale = w1 + w2;
297			if (scale == 0.) {
298			*cres = c_dfcolor;
299			return;
300			}
301			scale = 1. / scale;
302			cres->cx = (w1c1->cx + w2c2->cx) * scale;
303			cres->cy = (w1c1->cy + w2c2->cy) * scale;
304			cres->flags = C_CDXY\|C_CSXY;
305			}
306			}
307
308	greg	3.2	/* multiply two colors -- cres may be c1 or c2 */
309			double
310			c_cmult(C_COLOR cres, C_COLOR c1, double y1, C_COLOR *c2, double y2)
311			{
312			double yres;
313			int i;
314
315			if ((c1->flags\|c2->flags) & C_CDSPEC) {
316			long cmix[C_CNSS], cmax = 0; /* spectral multiply */
317
318			c_ccvt(c1, C_CSSPEC\|C_CSEFF);
319			c_ccvt(c2, C_CSSPEC\|C_CSEFF);
320			for (i = 0; i < C_CNSS; i++) {
321			cmix[i] = c1->ssamp[i] * c2->ssamp[i];
322			if (cmix[i] > cmax)
323			cmax = cmix[i];
324	greg	3.3	else if (cmix[i] < -cmax)
325			cmax = -cmix[i];
326	greg	3.2	}
327			cmax /= C_CMAXV;
328			if (!cmax) {
329			*cres = c_dfcolor;
330			return(0.);
331			}
332			cres->ssum = 0;
333			for (i = 0; i < C_CNSS; i++)
334			cres->ssum += cres->ssamp[i] = cmix[i] / cmax;
335			cres->flags = C_CDSPEC\|C_CSSPEC;
336
337	greg	3.6	c_ccvt(cres, C_CSEFF); /* below is correct */
338	greg	3.5	yres = y1 * y2 * c_y31.ssum * C_CLPWM /
339	greg	3.2	(c1->effc1->ssum c2->effc2->ssum)
340			cres->eff( cres->ssum(double)cmax +
341			C_CNSS/2.0*(cmax-1) );
342			} else {
343			float rgb1[3], rgb2[3]; /* CIE xy multiply */
344
345			c_toSharpRGB(c1, y1, rgb1);
346			c_toSharpRGB(c2, y2, rgb2);
347
348			rgb2[0] = rgb1[0]; rgb2[1] = rgb1[1]; rgb2[2] *= rgb1[2];
349
350			yres = c_fromSharpRGB(rgb2, cres);
351			}
352			return(yres);
353			}
354
355	greg	3.1	#define C1 3.741832e-16 /* W-m^2 */
356			#define C2 1.4388e-2 /* m-K */
357
358			#define bbsp(l,t) (C1/((l)(l)(l)(l)(l)(exp(C2/((t)(l)))-1.)))
359	greg	3.2	#define bblm(t) (C2*0.2/(t))
360	greg	3.1
361			/* set black body spectrum */
362			int
363			c_bbtemp(C_COLOR *clr, double tk)
364			{
365			double sf, wl;
366			int i;
367
368			if (tk < 1000)
369			return(0);
370			wl = bblm(tk); /* scalefactor based on peak */
371			if (wl < C_CMINWL*1e-9)
372			wl = C_CMINWL*1e-9;
373			else if (wl > C_CMAXWL*1e-9)
374			wl = C_CMAXWL*1e-9;
375			sf = C_CMAXV/bbsp(wl,tk);
376			clr->ssum = 0;
377			for (i = 0; i < C_CNSS; i++) {
378			wl = (C_CMINWL + iC_CWLI)1e-9;
379			clr->ssum += clr->ssamp[i] = sf*bbsp(wl,tk) + .5;
380			}
381			clr->flags = C_CDSPEC\|C_CSSPEC;
382			return(1);
383			}
384
385			#undef C1
386			#undef C2
387			#undef bbsp
388			#undef bblm