--- ray/src/common/ccolor.c 2012/05/18 20:43:13 3.4 +++ ray/src/common/ccolor.c 2017/04/05 00:54:50 3.10 @@ -1,14 +1,18 @@ #ifndef lint -static const char RCSid[] = "$Id: ccolor.c,v 3.4 2012/05/18 20:43:13 greg Exp $"; +static const char RCSid[] = "$Id: ccolor.c,v 3.10 2017/04/05 00:54:50 greg Exp $"; #endif /* * Spectral color handling routines */ #include +#include #include #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}, @@ -22,22 +26,24 @@ float XYZfromSharp[3][3] = { {-0.0123, 0.0167, 0.9955} }; -C_COLOR c_dfcolor = C_DEFCOLOR; +const C_COLOR c_dfcolor = C_DEFCOLOR; +const C_CHROMA c_dfchroma = 49750; /* c_encodeChroma(&c_dfcolor) */ + /* CIE 1931 Standard Observer curves */ -static const C_COLOR cie_xf = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF, +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 }; -static const C_COLOR cie_yf = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF, +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 }; -static const C_COLOR cie_zf = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF, +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}, @@ -96,13 +102,17 @@ 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]; + 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; @@ -112,7 +122,7 @@ c_fromSharpRGB(float cin[3], C_COLOR *cout) return(xyz[1]); } -/* assign arbitrary spectrum */ +/* assign arbitrary spectrum and return Y value */ double c_sset(C_COLOR *clr, double wlmin, double wlmax, const float spec[], int nwl) { @@ -134,6 +144,8 @@ c_sset(C_COLOR *clr, double wlmin, double wlmax, const wlmax -= wlstep; --nwl; } + if ((nwl <= 1) | (wlmin >= wlmax)) + return(0.); imax = nwl; /* box filter if necessary */ boxpos = 0; boxstep = 1; @@ -159,11 +171,11 @@ c_sset(C_COLOR *clr, double wlmin, double wlmax, const scale = va[i]; else if (va[i] < -scale) scale = -va[i]; - yval += va[i] * cie_yf.ssamp[i]; + yval += va[i] * c_y31.ssamp[i]; } if (scale <= 1e-7) return(0.); - yval /= (double)cie_yf.ssum; + yval /= (double)c_y31.ssum; scale = C_CMAXV / scale; clr->ssum = 0; /* convert to our spacing */ wl0 = wlmin; @@ -177,9 +189,9 @@ c_sset(C_COLOR *clr, double wlmin, double wlmax, const pos++; } if ((wl+1e-7 >= wl0) & (wl-1e-7 <= wl0)) - clr->ssamp[i] = scale*va[pos] + .5; + clr->ssamp[i] = scale*va[pos] + frand(); else /* interpolate if necessary */ - clr->ssamp[i] = .5 + scale / wlstep * + clr->ssamp[i] = frand() + scale / wlstep * ( va[pos]*(wl0+wlstep - wl) + va[pos+1]*(wl - wl0) ); clr->ssum += clr->ssamp[i]; @@ -218,25 +230,26 @@ c_ccvt(C_COLOR *clr, int fl) 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 += 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)cie_xf.ssum; - y /= (double)cie_yf.ssum; - z /= (double)cie_zf.ssum; + 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; - } else if (fl & C_CSSPEC) { /* cxy -> cspec */ + } + 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; + + z*cie_zp.ssamp[i] + frand(); if (clr->ssamp[i] < 0) /* out of gamut! */ clr->ssamp[i] = 0; else @@ -248,11 +261,11 @@ c_ccvt(C_COLOR *clr, int fl) if (clr->flags & C_CSSPEC) { /* from spectrum */ y = 0.; for (i = 0; i < C_CNSS; i++) - y += cie_yf.ssamp[i] * clr->ssamp[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*cie_xf.eff + clr->cy*cie_yf.eff + - (1. - clr->cx - clr->cy)*cie_zf.eff; + clr->eff = clr->cx*c_x31.eff + clr->cy*c_y31.eff + + (1. - clr->cx - clr->cy)*c_z31.eff; } clr->flags |= C_CSEFF; } @@ -283,7 +296,7 @@ c_cmix(C_COLOR *cres, double w1, C_COLOR *c1, double w scale = C_CMAXV / scale; cres->ssum = 0; for (i = 0; i < C_CNSS; i++) - cres->ssum += cres->ssamp[i] = scale*cmix[i] + .5; + cres->ssum += cres->ssamp[i] = scale*cmix[i] + frand(); cres->flags = C_CDSPEC|C_CSSPEC; } else { /* CIE xy mixing */ c_ccvt(c1, C_CSXY); @@ -331,8 +344,8 @@ c_cmult(C_COLOR *cres, C_COLOR *c1, double y1, C_COLOR 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 / + 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) ); @@ -373,7 +386,7 @@ c_bbtemp(C_COLOR *clr, double 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->ssum += clr->ssamp[i] = sf*bbsp(wl,tk) + frand(); } clr->flags = C_CDSPEC|C_CSSPEC; return(1); @@ -383,3 +396,39 @@ c_bbtemp(C_COLOR *clr, double tk) #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