--- ray/src/common/ccolor.c 2012/05/17 17:10:23 3.3 +++ 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.3 2012/05/17 17:10:23 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.3 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,7 +112,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) { @@ -132,6 +134,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; @@ -157,11 +161,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; @@ -216,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; @@ -246,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; } @@ -329,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) );