src/common/ccolor.c

#ifndef lint
static const char RCSid[] = "$Id: ccolor.c,v 3.11 2017/04/13 00:42:01 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+.01) {
                                wl0 += wlstep;
                                pos++;
                        }
                        if ((wl+.01 >= wl0) & (wl-.01 <= 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 colors are equivalent */
int
c_equiv(C_COLOR *c1, C_COLOR *c2)
{
        long    thresh;
        int     i;

        if (c1 == c2)
                return(1);
        c_ccvt(c1, C_CSXY);             /* first check chromaticities */
        c_ccvt(c2, C_CSXY);
        if (fabs(c1->cx - c2->cx) + fabs(c1->cy - c2->cy) > .015)
                return(0);              /* mismatch means definitely different */
        if (c1->flags & c2->flags & C_CDXY)
                return(1);              /* done if both defined as (x,y) */
        c_ccvt(c1, C_CSSPEC);           /* else compare spectra */
        c_ccvt(c2, C_CSSPEC);
        thresh = C_CMAXV/200*(c1->ssum + c2->ssum);
        for (i = 0; i < C_CNSS; i++)
                if (labs(c1->ssamp[i]*c2->ssum - c2->ssamp[i]*c1->ssum) > thresh)
                        return(0);
        return(1);
}


/* 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.12
Committed:	Thu Jan 4 01:55:42 2024 UTC (4 months, 1 week ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.11:	+33 -7 lines
Log Message:	feat(mgf2rad): Now handles spectral color translation
#	User	Rev	Content
1	greg	3.1	#ifndef lint
2	greg	3.12	static const char RCSid[] = "$Id: ccolor.c,v 3.11 2017/04/13 00:42:01 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	greg	3.12	wlstep = (wlmax - wlmin)/(nwl-1.);
139	greg	3.3	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	greg	3.12	while (wl0 + wlstep < wl+.01) {
188	greg	3.3	wl0 += wlstep;
189			pos++;
190			}
191	greg	3.12	if ((wl+.01 >= wl0) & (wl-.01 <= 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.12
204			/* check if colors are equivalent */
205			int
206			c_equiv(C_COLOR c1, C_COLOR c2)
207			{
208			long thresh;
209			int i;
210
211			if (c1 == c2)
212			return(1);
213			c_ccvt(c1, C_CSXY); /* first check chromaticities */
214			c_ccvt(c2, C_CSXY);
215			if (fabs(c1->cx - c2->cx) + fabs(c1->cy - c2->cy) > .015)
216			return(0); /* mismatch means definitely different */
217			if (c1->flags & c2->flags & C_CDXY)
218			return(1); /* done if both defined as (x,y) */
219			c_ccvt(c1, C_CSSPEC); /* else compare spectra */
220			c_ccvt(c2, C_CSSPEC);
221			thresh = C_CMAXV/200*(c1->ssum + c2->ssum);
222			for (i = 0; i < C_CNSS; i++)
223			if (labs(c1->ssamp[i]c2->ssum - c2->ssamp[i]c1->ssum) > thresh)
224			return(0);
225			return(1);
226			}
227
228
229	greg	3.1	/* check if color is grey */
230			int
231			c_isgrey(C_COLOR *clr)
232			{
233			if (!(clr->flags & (C_CSXY\|C_CSSPEC)))
234			return(1); /* no settings == grey */
235
236			c_ccvt(clr, C_CSXY);
237
238			return(clr->cx >= .323 && clr->cx <= .343 &&
239			clr->cy >= .323 && clr->cy <= .343);
240			}
241
242			/* convert color representations */
243			void
244			c_ccvt(C_COLOR *clr, int fl)
245			{
246			double x, y, z;
247			int i;
248
249	greg	3.12	fl &= ~clr->flags; /* ignore what's done */
250			if (!fl) /* everything's done! */
251	greg	3.1	return;
252			if (!(clr->flags & (C_CSXY\|C_CSSPEC))) {
253	greg	3.12	clr = c_dfcolor; / nothing was set! */
254	greg	3.1	return;
255			}
256			if (fl & C_CSXY) { /* cspec -> cxy */
257			x = y = z = 0.;
258			for (i = 0; i < C_CNSS; i++) {
259	greg	3.5	x += c_x31.ssamp[i] * clr->ssamp[i];
260			y += c_y31.ssamp[i] * clr->ssamp[i];
261			z += c_z31.ssamp[i] * clr->ssamp[i];
262	greg	3.1	}
263	greg	3.5	x /= (double)c_x31.ssum;
264			y /= (double)c_y31.ssum;
265			z /= (double)c_z31.ssum;
266	greg	3.1	z += x + y;
267	greg	3.11	if (z > 1e-6) {
268			clr->cx = x / z;
269			clr->cy = y / z;
270			} else
271			clr->cx = clr->cy = 1./3.;
272	greg	3.1	clr->flags \|= C_CSXY;
273	greg	3.6	}
274			if (fl & C_CSSPEC) { /* cxy -> cspec */
275	greg	3.1	x = clr->cx;
276			y = clr->cy;
277			z = 1. - x - y;
278			clr->ssum = 0;
279			for (i = 0; i < C_CNSS; i++) {
280			clr->ssamp[i] = xcie_xp.ssamp[i] + ycie_yp.ssamp[i]
281	greg	3.9	+ z*cie_zp.ssamp[i] + frand();
282	greg	3.1	if (clr->ssamp[i] < 0) /* out of gamut! */
283			clr->ssamp[i] = 0;
284			else
285			clr->ssum += clr->ssamp[i];
286			}
287			clr->flags \|= C_CSSPEC;
288			}
289			if (fl & C_CSEFF) { /* compute efficacy */
290			if (clr->flags & C_CSSPEC) { /* from spectrum */
291			y = 0.;
292			for (i = 0; i < C_CNSS; i++)
293	greg	3.5	y += c_y31.ssamp[i] * clr->ssamp[i];
294	greg	3.11	clr->eff = C_CLPWM * y / (clr->ssum + 0.0001);
295	greg	3.1	} else /* clr->flags & C_CSXY / { / from (x,y) */
296	greg	3.5	clr->eff = clr->cxc_x31.eff + clr->cyc_y31.eff +
297			(1. - clr->cx - clr->cy)*c_z31.eff;
298	greg	3.1	}
299			clr->flags \|= C_CSEFF;
300			}
301			}
302
303	greg	3.2	/* mix two colors according to weights given -- cres may be c1 or c2 */
304	greg	3.1	void
305			c_cmix(C_COLOR cres, double w1, C_COLOR c1, double w2, C_COLOR *c2)
306			{
307			double scale;
308			int i;
309
310	greg	3.11	if (w1 == 0) {
311			cres = c2;
312			return;
313			}
314			if (w2 == 0) {
315			cres = c1;
316			return;
317			}
318	greg	3.1	if ((c1->flags\|c2->flags) & C_CDSPEC) { /* spectral mixing */
319	greg	3.2	float cmix[C_CNSS];
320
321	greg	3.1	c_ccvt(c1, C_CSSPEC\|C_CSEFF);
322			c_ccvt(c2, C_CSSPEC\|C_CSEFF);
323	greg	3.11	if (c1->ssum*c2->ssum == 0) {
324			*cres = c_dfcolor;
325			return;
326			}
327	greg	3.1	w1 /= c1->eff*c1->ssum;
328			w2 /= c2->eff*c2->ssum;
329			scale = 0.;
330			for (i = 0; i < C_CNSS; i++) {
331			cmix[i] = w1c1->ssamp[i] + w2c2->ssamp[i];
332			if (cmix[i] > scale)
333			scale = cmix[i];
334	greg	3.3	else if (cmix[i] < -scale)
335			scale = -cmix[i];
336	greg	3.1	}
337			scale = C_CMAXV / scale;
338			cres->ssum = 0;
339			for (i = 0; i < C_CNSS; i++)
340	greg	3.9	cres->ssum += cres->ssamp[i] = scale*cmix[i] + frand();
341	greg	3.1	cres->flags = C_CDSPEC\|C_CSSPEC;
342			} else { /* CIE xy mixing */
343			c_ccvt(c1, C_CSXY);
344			c_ccvt(c2, C_CSXY);
345			w1 *= c2->cy;
346			w2 *= c1->cy;
347			scale = w1 + w2;
348			if (scale == 0.) {
349			*cres = c_dfcolor;
350			return;
351			}
352			scale = 1. / scale;
353			cres->cx = (w1c1->cx + w2c2->cx) * scale;
354			cres->cy = (w1c1->cy + w2c2->cy) * scale;
355			cres->flags = C_CDXY\|C_CSXY;
356			}
357			}
358
359	greg	3.2	/* multiply two colors -- cres may be c1 or c2 */
360			double
361			c_cmult(C_COLOR cres, C_COLOR c1, double y1, C_COLOR *c2, double y2)
362			{
363			double yres;
364			int i;
365
366			if ((c1->flags\|c2->flags) & C_CDSPEC) {
367			long cmix[C_CNSS], cmax = 0; /* spectral multiply */
368
369			c_ccvt(c1, C_CSSPEC\|C_CSEFF);
370			c_ccvt(c2, C_CSSPEC\|C_CSEFF);
371			for (i = 0; i < C_CNSS; i++) {
372			cmix[i] = c1->ssamp[i] * c2->ssamp[i];
373			if (cmix[i] > cmax)
374			cmax = cmix[i];
375	greg	3.3	else if (cmix[i] < -cmax)
376			cmax = -cmix[i];
377	greg	3.2	}
378			cmax /= C_CMAXV;
379			if (!cmax) {
380			*cres = c_dfcolor;
381			return(0.);
382			}
383			cres->ssum = 0;
384			for (i = 0; i < C_CNSS; i++)
385			cres->ssum += cres->ssamp[i] = cmix[i] / cmax;
386			cres->flags = C_CDSPEC\|C_CSSPEC;
387
388	greg	3.6	c_ccvt(cres, C_CSEFF); /* below is correct */
389	greg	3.5	yres = y1 * y2 * c_y31.ssum * C_CLPWM /
390	greg	3.2	(c1->effc1->ssum c2->effc2->ssum)
391			cres->eff( cres->ssum(double)cmax +
392			C_CNSS/2.0*(cmax-1) );
393			} else {
394			float rgb1[3], rgb2[3]; /* CIE xy multiply */
395
396			c_toSharpRGB(c1, y1, rgb1);
397			c_toSharpRGB(c2, y2, rgb2);
398
399			rgb2[0] = rgb1[0]; rgb2[1] = rgb1[1]; rgb2[2] *= rgb1[2];
400
401			yres = c_fromSharpRGB(rgb2, cres);
402			}
403			return(yres);
404			}
405
406	greg	3.1	#define C1 3.741832e-16 /* W-m^2 */
407			#define C2 1.4388e-2 /* m-K */
408
409			#define bbsp(l,t) (C1/((l)(l)(l)(l)(l)(exp(C2/((t)(l)))-1.)))
410	greg	3.2	#define bblm(t) (C2*0.2/(t))
411	greg	3.1
412			/* set black body spectrum */
413			int
414			c_bbtemp(C_COLOR *clr, double tk)
415			{
416			double sf, wl;
417			int i;
418
419			if (tk < 1000)
420			return(0);
421			wl = bblm(tk); /* scalefactor based on peak */
422			if (wl < C_CMINWL*1e-9)
423			wl = C_CMINWL*1e-9;
424			else if (wl > C_CMAXWL*1e-9)
425			wl = C_CMAXWL*1e-9;
426			sf = C_CMAXV/bbsp(wl,tk);
427			clr->ssum = 0;
428			for (i = 0; i < C_CNSS; i++) {
429			wl = (C_CMINWL + iC_CWLI)1e-9;
430	greg	3.9	clr->ssum += clr->ssamp[i] = sf*bbsp(wl,tk) + frand();
431	greg	3.1	}
432			clr->flags = C_CDSPEC\|C_CSSPEC;
433			return(1);
434			}
435
436			#undef C1
437			#undef C2
438			#undef bbsp
439			#undef bblm
440	greg	3.7
441			#define UV_NORMF 410.
442
443			/* encode (x,y) chromaticity */
444			C_CHROMA
445			c_encodeChroma(C_COLOR *clr)
446			{
447			double df;
448			int ub, vb;
449
450			c_ccvt(clr, C_CSXY);
451			df = UV_NORMF/(-2.clr->cx + 12.clr->cy + 3.);
452	greg	3.9	ub = 4.clr->cxdf + frand();
453			if (ub > 0xff) ub = 0xff;
454	greg	3.10	else ub *= (ub > 0);
455	greg	3.9	vb = 9.clr->cydf + frand();
456			if (vb > 0xff) vb = 0xff;
457	greg	3.10	else vb *= (vb > 0);
458	greg	3.7
459			return(vb<<8 \| ub);
460			}
461
462			/* decode (x,y) chromaticity */
463			void
464			c_decodeChroma(C_COLOR *cres, C_CHROMA ccode)
465			{
466	greg	3.9	double up = (ccode & 0xff)*(1./UV_NORMF);
467			double vp = (ccode>>8 & 0xff)*(1./UV_NORMF);
468	greg	3.7	double df = 1./(6.up - 16.vp + 12.);
469
470			cres->cx = 9.up df;
471			cres->cy = 4.vp df;
472			cres->flags = C_CDXY\|C_CSXY;
473			}
474
475			#undef UV_NORMF