--- ray/src/common/ccolor.c 2012/05/17 05:47:59 3.2 +++ ray/src/common/ccolor.c 2012/07/17 23:40:24 3.6 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: ccolor.c,v 3.2 2012/05/17 05:47:59 greg Exp $"; +static const char RCSid[] = "$Id: ccolor.c,v 3.6 2012/07/17 23:40:24 greg Exp $"; #endif /* * Spectral color handling routines @@ -9,23 +9,35 @@ static const char RCSid[] = "$Id: ccolor.c,v 3.2 2012/ #include #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; +const 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, +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}, @@ -57,20 +69,9 @@ static const C_COLOR cie_zp = { 1, NULL, C_CDSPEC|C_CS 36057L, .0, .0, }; - /* Sharp primary matrix */ -static const float toSharp[3][3] = { - { 1.2694, -0.0988, -0.1706}, - {-0.8364, 1.8006, 0.0357}, - { 0.0297, -0.0315, 1.0018} -}; - /* inverse Sharp primary matrix */ -static const float fromSharp[3][3] = { - { 0.8156, 0.0472, 0.1372}, - { 0.3791, 0.5769, 0.0440}, - {-0.0123, 0.0167, 0.9955} -}; -static void +/* convert to sharpened RGB color for low-error operations */ +void c_toSharpRGB(C_COLOR *cin, double cieY, float cout[3]) { double xyz[3]; @@ -81,25 +82,26 @@ c_toSharpRGB(C_COLOR *cin, double cieY, float cout[3]) xyz[1] = cieY; xyz[2] = (1. - cin->cx - cin->cy)/cin->cy * cieY; - cout[0] = toSharp[0][0]*xyz[0] + toSharp[0][1]*xyz[1] + - toSharp[0][2]*xyz[2]; - cout[1] = toSharp[1][0]*xyz[0] + toSharp[1][1]*xyz[1] + - toSharp[1][2]*xyz[2]; - cout[2] = toSharp[2][0]*xyz[0] + toSharp[2][1]*xyz[1] + - toSharp[2][2]*xyz[2]; + 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]; } -static double +/* convert back from sharpened RGB color */ +double c_fromSharpRGB(float cin[3], C_COLOR *cout) { double xyz[3], sf; - xyz[0] = fromSharp[0][0]*cin[0] + fromSharp[0][1]*cin[1] + - fromSharp[0][2]*cin[2]; - xyz[1] = fromSharp[1][0]*cin[0] + fromSharp[1][1]*cin[1] + - fromSharp[1][2]*cin[2]; - xyz[2] = fromSharp[2][0]*cin[0] + fromSharp[2][1]*cin[1] + - fromSharp[2][2]*cin[2]; + 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]); @@ -110,6 +112,84 @@ c_fromSharpRGB(float cin[3], C_COLOR *cout) 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) @@ -140,18 +220,19 @@ 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; @@ -170,11 +251,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; } @@ -199,6 +280,8 @@ c_cmix(C_COLOR *cres, double w1, C_COLOR *c1, double w 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; @@ -238,6 +321,8 @@ c_cmult(C_COLOR *cres, C_COLOR *c1, double y1, C_COLOR 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) { @@ -249,8 +334,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) );