src/common/spec_rgb.c

#ifndef lint
static const char       RCSid[] = "$Id: spec_rgb.c,v 2.19 2009/11/20 17:51:47 greg Exp $";
#endif
/*
 * Convert colors and spectral ranges.
 * Added von Kries white-balance calculations 10/01 (GW).
 *
 * Externals declared in color.h
 */

#include "copyright.h"

#include <stdio.h>
#include <string.h>
#include "color.h"

#define CEPS    1e-4                    /* color epsilon */

#define CEQ(v1,v2)      ((v1) <= (v2)+CEPS && (v2) <= (v1)+CEPS)

#define XYEQ(c1,c2)     (CEQ((c1)[CIEX],(c2)[CIEX]) && CEQ((c1)[CIEY],(c2)[CIEY]))


RGBPRIMS  stdprims = STDPRIMS;  /* standard primary chromaticities */

COLOR  cblack = BLKCOLOR;               /* global black color */
COLOR  cwhite = WHTCOLOR;               /* global white color */

float  xyneu[2] = {1./3., 1./3.};       /* neutral xy chromaticities */

/*
 *      The following table contains the CIE tristimulus integrals
 *  for X, Y, and Z.  The table is cumulative, so that
 *  each color coordinate integrates to 1.
 */

#define  STARTWL        380             /* starting wavelength (nanometers) */
#define  INCWL          10              /* wavelength increment */
#define  NINC           40              /* # of values */

static BYTE  chroma[3][NINC] = {
        {                                                       /* X */
                0,   0,   0,   2,   6,   13,  22,  30,  36,  41,
                42,  43,  43,  44,  46,  52,  60,  71,  87,  106,
                128, 153, 178, 200, 219, 233, 243, 249, 252, 254,
                255, 255, 255, 255, 255, 255, 255, 255, 255, 255
        }, {                                                    /* Y */
                0,   0,   0,   0,   0,   1,   2,   4,   7,   11,
                17,  24,  34,  48,  64,  84,  105, 127, 148, 169,
                188, 205, 220, 232, 240, 246, 250, 253, 254, 255,
                255, 255, 255, 255, 255, 255, 255, 255, 255, 255
        }, {                                                    /* Z */
                0,   0,   2,   10,  32,  66,  118, 153, 191, 220,
                237, 246, 251, 253, 254, 255, 255, 255, 255, 255,
                255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
                255, 255, 255, 255, 255, 255, 255, 255, 255, 255
        }
};

COLORMAT  xyz2rgbmat = {                /* XYZ to RGB (no white balance) */
        {(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g)/CIE_C_rD,
         (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b)/CIE_C_rD,
         (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g)/CIE_C_rD},
        {(CIE_y_b - CIE_y_r - CIE_y_b*CIE_x_r + CIE_y_r*CIE_x_b)/CIE_C_gD,
         (CIE_x_r - CIE_x_b - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r)/CIE_C_gD,
         (CIE_x_b*CIE_y_r - CIE_x_r*CIE_y_b)/CIE_C_gD},
        {(CIE_y_r - CIE_y_g - CIE_y_r*CIE_x_g + CIE_y_g*CIE_x_r)/CIE_C_bD,
         (CIE_x_g - CIE_x_r - CIE_x_g*CIE_y_r + CIE_x_r*CIE_y_g)/CIE_C_bD,
         (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r)/CIE_C_bD}
};

COLORMAT  rgb2xyzmat = {                /* RGB to XYZ (no white balance) */
        {CIE_x_r*CIE_C_rD/CIE_D,CIE_x_g*CIE_C_gD/CIE_D,CIE_x_b*CIE_C_bD/CIE_D},
        {CIE_y_r*CIE_C_rD/CIE_D,CIE_y_g*CIE_C_gD/CIE_D,CIE_y_b*CIE_C_bD/CIE_D},
        {(1.-CIE_x_r-CIE_y_r)*CIE_C_rD/CIE_D,
         (1.-CIE_x_g-CIE_y_g)*CIE_C_gD/CIE_D,
         (1.-CIE_x_b-CIE_y_b)*CIE_C_bD/CIE_D}
};

COLORMAT  vkmat = {             /* Sharp primary matrix */
        { 1.2694, -0.0988, -0.1706},
        {-0.8364,  1.8006,  0.0357},
        { 0.0297, -0.0315,  1.0018}
};

COLORMAT  ivkmat = {            /* inverse Sharp primary matrix */
        { 0.8156,  0.0472,  0.1372},
        { 0.3791,  0.5769,  0.0440},
        {-0.0123,  0.0167,  0.9955}
};


void
spec_rgb(                       /* compute RGB color from spectral range */
COLOR  col,
int  s,
int  e
)
{
        COLOR  ciecolor;

        spec_cie(ciecolor, s, e);
        cie_rgb(col, ciecolor);
}


void
spec_cie(                       /* compute a color from a spectral range */
COLOR  col,             /* returned color */
int  s,                 /* starting and ending wavelengths */
int  e
)
{
        register int  i, d, r;
        
        s -= STARTWL;
        if (s < 0)
                s = 0;

        e -= STARTWL;
        if (e <= s) {
                col[CIEX] = col[CIEY] = col[CIEZ] = 0.0;
                return;
        }
        if (e >= INCWL*(NINC - 1))
                e = INCWL*(NINC - 1) - 1;

        d = e / INCWL;                  /* interpolate values */
        r = e % INCWL;
        for (i = 0; i < 3; i++)
                col[i] = chroma[i][d]*(INCWL - r) + chroma[i][d + 1]*r;

        d = s / INCWL;
        r = s % INCWL;
        for (i = 0; i < 3; i++)
                col[i] -= chroma[i][d]*(INCWL - r) - chroma[i][d + 1]*r;

        col[CIEX] = (col[CIEX] + 0.5) * (1./(256*INCWL));
        col[CIEY] = (col[CIEY] + 0.5) * (1./(256*INCWL));
        col[CIEZ] = (col[CIEZ] + 0.5) * (1./(256*INCWL));
}


void
cie_rgb(                        /* convert CIE color to standard RGB */
COLOR   rgb,
COLOR  xyz
)
{
        colortrans(rgb, xyz2rgbmat, xyz);
        clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
}


int
clipgamut(                      /* clip to gamut cube */
COLOR  col,
double  brt,
int  gamut,
COLOR  lower,
COLOR  upper
)
{
        int  rflags = 0;
        double  brtmin, brtmax, v, vv;
        COLOR  cgry;
        register int  i;
                                        /* check for no check */
        if (gamut == 0) return(0);
                                        /* check brightness limits */
        brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
        if (gamut & CGAMUT_LOWER && brt < brtmin) {
                copycolor(col, lower);
                return(CGAMUT_LOWER);
        }
        brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
        if (gamut & CGAMUT_UPPER && brt > brtmax) {
                copycolor(col, upper);
                return(CGAMUT_UPPER);
        }
                                        /* compute equivalent grey */
        v = (brt - brtmin)/(brtmax - brtmin);
        for (i = 0; i < 3; i++)
                cgry[i] = v*upper[i] + (1.-v)*lower[i];
        vv = 1.;                        /* check each limit */
        for (i = 0; i < 3; i++)
                if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
                        v = (lower[i] - cgry[i])/(col[i] - cgry[i]);
                        if (v < vv) vv = v;
                        rflags |= CGAMUT_LOWER;
                } else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
                        v = (upper[i] - cgry[i])/(col[i] - cgry[i]);
                        if (v < vv) vv = v;
                        rflags |= CGAMUT_UPPER;
                }
        if (rflags)                     /* desaturate to cube face */
                for (i = 0; i < 3; i++)
                        col[i] = vv*col[i] + (1.-vv)*cgry[i];
        return(rflags);
}


void
colortrans(                     /* convert c1 by mat and put into c2 */
register COLOR  c2,
register COLORMAT  mat,
register COLOR  c1
)
{
        COLOR   cout;

        cout[0] = mat[0][0]*c1[0] + mat[0][1]*c1[1] + mat[0][2]*c1[2];
        cout[1] = mat[1][0]*c1[0] + mat[1][1]*c1[1] + mat[1][2]*c1[2];
        cout[2] = mat[2][0]*c1[0] + mat[2][1]*c1[1] + mat[2][2]*c1[2];

        copycolor(c2, cout);
}


void
multcolormat(                   /* multiply m1 by m2 and put into m3 */
COLORMAT  m3,                   /* m3 can be either m1 or m2 w/o harm */
COLORMAT  m2,
COLORMAT  m1
)
{
        COLORMAT  mt;
        register int  i, j;

        for (i = 0; i < 3; i++)
                for (j = 0; j < 3; j++)
                        mt[i][j] =      m1[i][0]*m2[0][j] +
                                        m1[i][1]*m2[1][j] +
                                        m1[i][2]*m2[2][j] ;
        cpcolormat(m3, mt);
}


int
colorprimsOK(                   /* are color primaries reasonable? */
RGBPRIMS  pr
)
{
        int     i, j;
        
        for (i = 0; i < 4; i++) {
                if ((pr[i][CIEX] <= -.5) | (pr[i][CIEY] <= -.5))
                        return(0);
                if ((pr[i][CIEX] >= 1.5) | (pr[i][CIEY] >= 1.5))
                        return(0);
                if (pr[i][CIEX] + pr[i][CIEY] >= 1.5)
                        return(0);
        }
        for (i = 0; i < 4; i++)
                for (j = i+1; j < 4; j++)
                        if (CEQ(pr[i][CIEX],pr[j][CIEX]) &&
                                        CEQ(pr[i][CIEY],pr[j][CIEY]))
                                return(0);
        return(1);
}


int
compxyz2rgbmat(                 /* compute conversion from CIE to RGB space */
COLORMAT  mat,
register RGBPRIMS  pr
)
{
        double  C_rD, C_gD, C_bD;

        if (pr == stdprims) {   /* can use xyz2rgbmat */
                cpcolormat(mat, xyz2rgbmat);
                return(1);
        }
        if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
                return(0);
        C_rD = (1./pr[WHT][CIEY]) *
                        ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
                          pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
                  pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
        if (CEQ(C_rD,0.))
                return(0);
        C_gD = (1./pr[WHT][CIEY]) *
                        ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
                          pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
                  pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
        if (CEQ(C_gD,0.))
                return(0);
        C_bD = (1./pr[WHT][CIEY]) *
                        ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
                          pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
                  pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
        if (CEQ(C_bD,0.))
                return(0);
        mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
                        pr[BLU][CIEX]*pr[GRN][CIEY] +
                        pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
        mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
                        pr[BLU][CIEX]*pr[GRN][CIEY] +
                        pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
        mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
                        pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
        mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
                        pr[BLU][CIEY]*pr[RED][CIEX] +
                        pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
        mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
                        pr[RED][CIEX]*pr[BLU][CIEY] +
                        pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
        mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
                        pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
        mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
                        pr[RED][CIEY]*pr[GRN][CIEX] +
                        pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
        mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
                        pr[GRN][CIEX]*pr[RED][CIEY] +
                        pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
        mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
                        pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
        return(1);
}


int
comprgb2xyzmat(                 /* compute conversion from RGB to CIE space */
COLORMAT  mat,
register RGBPRIMS  pr
)
{
        double  C_rD, C_gD, C_bD, D;

        if (pr == stdprims) {   /* can use rgb2xyzmat */
                cpcolormat(mat, rgb2xyzmat);
                return(1);
        }
        if (CEQ(pr[WHT][CIEX],0.) | CEQ(pr[WHT][CIEY],0.))
                return(0);
        C_rD = (1./pr[WHT][CIEY]) *
                        ( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
                          pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
                  pr[GRN][CIEX]*pr[BLU][CIEY] - pr[BLU][CIEX]*pr[GRN][CIEY] ) ;
        C_gD = (1./pr[WHT][CIEY]) *
                        ( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
                          pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
                  pr[RED][CIEX]*pr[BLU][CIEY] + pr[BLU][CIEX]*pr[RED][CIEY] ) ;
        C_bD = (1./pr[WHT][CIEY]) *
                        ( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
                          pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
                  pr[RED][CIEX]*pr[GRN][CIEY] - pr[GRN][CIEX]*pr[RED][CIEY] ) ;
        D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
                        pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
                        pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
        if (CEQ(D,0.))
                return(0);
        mat[0][0] = pr[RED][CIEX]*C_rD/D;
        mat[0][1] = pr[GRN][CIEX]*C_gD/D;
        mat[0][2] = pr[BLU][CIEX]*C_bD/D;
        mat[1][0] = pr[RED][CIEY]*C_rD/D;
        mat[1][1] = pr[GRN][CIEY]*C_gD/D;
        mat[1][2] = pr[BLU][CIEY]*C_bD/D;
        mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
        mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
        mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
        return(1);
}


int
comprgb2rgbmat(                 /* compute conversion from RGB1 to RGB2 */
COLORMAT  mat,
RGBPRIMS  pr1,
RGBPRIMS  pr2
)
{
        COLORMAT  pr1toxyz, xyztopr2;

        if (pr1 == pr2) {
                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
                mat[0][1] = mat[0][2] = mat[1][0] =
                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
                return(1);
        }
        if (!comprgb2xyzmat(pr1toxyz, pr1))
                return(0);
        if (!compxyz2rgbmat(xyztopr2, pr2))
                return(0);
                                /* combine transforms */
        multcolormat(mat, pr1toxyz, xyztopr2);
        return(1);
}


int
compxyzWBmat(                   /* CIE von Kries transform from wht1 to wht2 */
COLORMAT  mat,
float  wht1[2],
float  wht2[2]
)
{
        COLOR   cw1, cw2;
        if (XYEQ(wht1,wht2)) {
                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
                mat[0][1] = mat[0][2] = mat[1][0] =
                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
                return(1);
        }
        if (CEQ(wht1[CIEX],0.) | CEQ(wht1[CIEY],0.))
                return(0);
        cw1[RED] = wht1[CIEX]/wht1[CIEY];
        cw1[GRN] = 1.;
        cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
        colortrans(cw1, vkmat, cw1);
        if (CEQ(wht2[CIEX],0.) | CEQ(wht2[CIEY],0.))
                return(0);
        cw2[RED] = wht2[CIEX]/wht2[CIEY];
        cw2[GRN] = 1.;
        cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
        colortrans(cw2, vkmat, cw2);
        if (CEQ(cw1[RED],0.) | CEQ(cw1[GRN],0.) | CEQ(cw1[BLU],0.))
                return(0);
        mat[0][0] = cw2[RED]/cw1[RED];
        mat[1][1] = cw2[GRN]/cw1[GRN];
        mat[2][2] = cw2[BLU]/cw1[BLU];
        mat[0][1] = mat[0][2] = mat[1][0] =
        mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
        multcolormat(mat, vkmat, mat);
        multcolormat(mat, mat, ivkmat);
        return(1);
}


int
compxyz2rgbWBmat(               /* von Kries conversion from CIE to RGB space */
COLORMAT  mat,
RGBPRIMS  pr
)
{
        COLORMAT        wbmat;

        if (!compxyz2rgbmat(mat, pr))
                return(0);
        if (XYEQ(pr[WHT],xyneu))
                return(1);
        if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
                return(0);
        multcolormat(mat, wbmat, mat);
        return(1);
}

int
comprgb2xyzWBmat(               /* von Kries conversion from RGB to CIE space */
COLORMAT  mat,
RGBPRIMS  pr
)
{
        COLORMAT        wbmat;
        
        if (!comprgb2xyzmat(mat, pr))
                return(0);
        if (XYEQ(pr[WHT],xyneu))
                return(1);
        if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
                return(0);
        multcolormat(mat, mat, wbmat);
        return(1);
}

int
comprgb2rgbWBmat(               /* von Kries conversion from RGB1 to RGB2 */
COLORMAT  mat,
RGBPRIMS  pr1,
RGBPRIMS  pr2
)
{
        COLORMAT  pr1toxyz, xyztopr2, wbmat;

        if (pr1 == pr2) {
                mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
                mat[0][1] = mat[0][2] = mat[1][0] =
                mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
                return(1);
        }
        if (!comprgb2xyzmat(pr1toxyz, pr1))
                return(0);
        if (!compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]))
                return(0);
        if (!compxyz2rgbmat(xyztopr2, pr2))
                return(0);
                                /* combine transforms */
        multcolormat(mat, pr1toxyz, wbmat);
        multcolormat(mat, mat, xyztopr2);
        return(1);
}
Revision:	2.20
Committed:	Mon Jan 4 19:40:14 2010 UTC (15 years, 10 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R0
Changes since 2.19:	+7 -2 lines
Log Message:	Improved chromaticity check
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.20	static const char RCSid[] = "$Id: spec_rgb.c,v 2.19 2009/11/20 17:51:47 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* Convert colors and spectral ranges.
6	greg	2.11	* Added von Kries white-balance calculations 10/01 (GW).
7			*
8			* Externals declared in color.h
9			*/
10
11	greg	2.12	#include "copyright.h"
12	greg	1.1
13	greg	2.13	#include <stdio.h>
14			#include <string.h>
15	greg	1.1	#include "color.h"
16	greg	2.11
17			#define CEPS 1e-4 /* color epsilon */
18
19			#define CEQ(v1,v2) ((v1) <= (v2)+CEPS && (v2) <= (v1)+CEPS)
20	greg	1.1
21	greg	2.11	#define XYEQ(c1,c2) (CEQ((c1)[CIEX],(c2)[CIEX]) && CEQ((c1)[CIEY],(c2)[CIEY]))
22	greg	2.5
23	gwlarson	2.10
24	greg	2.11	RGBPRIMS stdprims = STDPRIMS; /* standard primary chromaticities */
25	greg	2.5
26	greg	2.8	COLOR cblack = BLKCOLOR; /* global black color */
27			COLOR cwhite = WHTCOLOR; /* global white color */
28
29	greg	2.11	float xyneu[2] = {1./3., 1./3.}; /* neutral xy chromaticities */
30
31	greg	1.1	/*
32			* The following table contains the CIE tristimulus integrals
33			* for X, Y, and Z. The table is cumulative, so that
34			* each color coordinate integrates to 1.
35			*/
36
37			#define STARTWL 380 /* starting wavelength (nanometers) */
38			#define INCWL 10 /* wavelength increment */
39			#define NINC 40 /* # of values */
40
41			static BYTE chroma[3][NINC] = {
42			{ /* X */
43			0, 0, 0, 2, 6, 13, 22, 30, 36, 41,
44			42, 43, 43, 44, 46, 52, 60, 71, 87, 106,
45			128, 153, 178, 200, 219, 233, 243, 249, 252, 254,
46			255, 255, 255, 255, 255, 255, 255, 255, 255, 255
47			}, { /* Y */
48			0, 0, 0, 0, 0, 1, 2, 4, 7, 11,
49			17, 24, 34, 48, 64, 84, 105, 127, 148, 169,
50			188, 205, 220, 232, 240, 246, 250, 253, 254, 255,
51			255, 255, 255, 255, 255, 255, 255, 255, 255, 255
52			}, { /* Z */
53			0, 0, 2, 10, 32, 66, 118, 153, 191, 220,
54			237, 246, 251, 253, 254, 255, 255, 255, 255, 255,
55			255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
56			255, 255, 255, 255, 255, 255, 255, 255, 255, 255
57			}
58			};
59
60	greg	2.11	COLORMAT xyz2rgbmat = { /* XYZ to RGB (no white balance) */
61	greg	2.3	{(CIE_y_g - CIE_y_b - CIE_x_bCIE_y_g + CIE_y_bCIE_x_g)/CIE_C_rD,
62			(CIE_x_b - CIE_x_g - CIE_x_bCIE_y_g + CIE_x_gCIE_y_b)/CIE_C_rD,
63			(CIE_x_gCIE_y_b - CIE_x_bCIE_y_g)/CIE_C_rD},
64			{(CIE_y_b - CIE_y_r - CIE_y_bCIE_x_r + CIE_y_rCIE_x_b)/CIE_C_gD,
65			(CIE_x_r - CIE_x_b - CIE_x_rCIE_y_b + CIE_x_bCIE_y_r)/CIE_C_gD,
66			(CIE_x_bCIE_y_r - CIE_x_rCIE_y_b)/CIE_C_gD},
67			{(CIE_y_r - CIE_y_g - CIE_y_rCIE_x_g + CIE_y_gCIE_x_r)/CIE_C_bD,
68			(CIE_x_g - CIE_x_r - CIE_x_gCIE_y_r + CIE_x_rCIE_y_g)/CIE_C_bD,
69			(CIE_x_rCIE_y_g - CIE_x_gCIE_y_r)/CIE_C_bD}
70	greg	1.2	};
71	greg	1.1
72	greg	2.11	COLORMAT rgb2xyzmat = { /* RGB to XYZ (no white balance) */
73	greg	2.4	{CIE_x_rCIE_C_rD/CIE_D,CIE_x_gCIE_C_gD/CIE_D,CIE_x_b*CIE_C_bD/CIE_D},
74			{CIE_y_rCIE_C_rD/CIE_D,CIE_y_gCIE_C_gD/CIE_D,CIE_y_b*CIE_C_bD/CIE_D},
75			{(1.-CIE_x_r-CIE_y_r)*CIE_C_rD/CIE_D,
76			(1.-CIE_x_g-CIE_y_g)*CIE_C_gD/CIE_D,
77			(1.-CIE_x_b-CIE_y_b)*CIE_C_bD/CIE_D}
78			};
79	greg	1.2
80	greg	2.11	COLORMAT vkmat = { /* Sharp primary matrix */
81			{ 1.2694, -0.0988, -0.1706},
82			{-0.8364, 1.8006, 0.0357},
83			{ 0.0297, -0.0315, 1.0018}
84			};
85
86			COLORMAT ivkmat = { /* inverse Sharp primary matrix */
87			{ 0.8156, 0.0472, 0.1372},
88			{ 0.3791, 0.5769, 0.0440},
89			{-0.0123, 0.0167, 0.9955}
90			};
91	greg	1.2
92	greg	2.4
93	greg	2.11	void
94	greg	2.15	spec_rgb( /* compute RGB color from spectral range */
95			COLOR col,
96			int s,
97			int e
98			)
99	greg	1.1	{
100			COLOR ciecolor;
101
102			spec_cie(ciecolor, s, e);
103			cie_rgb(col, ciecolor);
104			}
105
106
107	greg	2.11	void
108	greg	2.15	spec_cie( /* compute a color from a spectral range */
109			COLOR col, /* returned color */
110			int s, /* starting and ending wavelengths */
111			int e
112			)
113	greg	1.1	{
114			register int i, d, r;
115
116			s -= STARTWL;
117			if (s < 0)
118			s = 0;
119
120			e -= STARTWL;
121	greg	1.2	if (e <= s) {
122	greg	2.5	col[CIEX] = col[CIEY] = col[CIEZ] = 0.0;
123	greg	1.2	return;
124			}
125	greg	1.1	if (e >= INCWL*(NINC - 1))
126			e = INCWL*(NINC - 1) - 1;
127
128			d = e / INCWL; /* interpolate values */
129			r = e % INCWL;
130			for (i = 0; i < 3; i++)
131			col[i] = chroma[i][d](INCWL - r) + chroma[i][d + 1]r;
132
133			d = s / INCWL;
134			r = s % INCWL;
135			for (i = 0; i < 3; i++)
136			col[i] -= chroma[i][d](INCWL - r) - chroma[i][d + 1]r;
137
138	greg	2.5	col[CIEX] = (col[CIEX] + 0.5) * (1./(256*INCWL));
139			col[CIEY] = (col[CIEY] + 0.5) * (1./(256*INCWL));
140			col[CIEZ] = (col[CIEZ] + 0.5) * (1./(256*INCWL));
141	greg	1.1	}
142
143
144	greg	2.11	void
145	greg	2.15	cie_rgb( /* convert CIE color to standard RGB */
146			COLOR rgb,
147			COLOR xyz
148			)
149	greg	2.8	{
150			colortrans(rgb, xyz2rgbmat, xyz);
151			clipgamut(rgb, xyz[CIEY], CGAMUT_LOWER, cblack, cwhite);
152			}
153
154
155			int
156	greg	2.15	clipgamut( /* clip to gamut cube */
157			COLOR col,
158			double brt,
159			int gamut,
160			COLOR lower,
161			COLOR upper
162			)
163	greg	2.8	{
164			int rflags = 0;
165			double brtmin, brtmax, v, vv;
166			COLOR cgry;
167			register int i;
168			/* check for no check */
169			if (gamut == 0) return(0);
170			/* check brightness limits */
171			brtmin = 1./3.*(lower[0]+lower[1]+lower[2]);
172			if (gamut & CGAMUT_LOWER && brt < brtmin) {
173			copycolor(col, lower);
174			return(CGAMUT_LOWER);
175			}
176			brtmax = 1./3.*(upper[0]+upper[1]+upper[2]);
177			if (gamut & CGAMUT_UPPER && brt > brtmax) {
178			copycolor(col, upper);
179			return(CGAMUT_UPPER);
180			}
181			/* compute equivalent grey */
182			v = (brt - brtmin)/(brtmax - brtmin);
183			for (i = 0; i < 3; i++)
184			cgry[i] = vupper[i] + (1.-v)lower[i];
185			vv = 1.; /* check each limit */
186			for (i = 0; i < 3; i++)
187			if (gamut & CGAMUT_LOWER && col[i] < lower[i]) {
188	greg	2.14	v = (lower[i] - cgry[i])/(col[i] - cgry[i]);
189	greg	2.8	if (v < vv) vv = v;
190			rflags \|= CGAMUT_LOWER;
191			} else if (gamut & CGAMUT_UPPER && col[i] > upper[i]) {
192	greg	2.14	v = (upper[i] - cgry[i])/(col[i] - cgry[i]);
193	greg	2.8	if (v < vv) vv = v;
194			rflags \|= CGAMUT_UPPER;
195			}
196			if (rflags) /* desaturate to cube face */
197			for (i = 0; i < 3; i++)
198			col[i] = vvcol[i] + (1.-vv)cgry[i];
199			return(rflags);
200			}
201
202
203	greg	2.11	void
204	greg	2.15	colortrans( /* convert c1 by mat and put into c2 */
205			register COLOR c2,
206			register COLORMAT mat,
207			register COLOR c1
208			)
209	greg	1.1	{
210	greg	2.9	COLOR cout;
211
212			cout[0] = mat[0][0]c1[0] + mat[0][1]c1[1] + mat[0][2]*c1[2];
213			cout[1] = mat[1][0]c1[0] + mat[1][1]c1[1] + mat[1][2]*c1[2];
214			cout[2] = mat[2][0]c1[0] + mat[2][1]c1[1] + mat[2][2]*c1[2];
215
216			copycolor(c2, cout);
217	greg	2.4	}
218
219
220	greg	2.11	void
221	greg	2.15	multcolormat( /* multiply m1 by m2 and put into m3 */
222			COLORMAT m3, /* m3 can be either m1 or m2 w/o harm */
223			COLORMAT m2,
224			COLORMAT m1
225			)
226	greg	2.4	{
227	greg	2.5	COLORMAT mt;
228			register int i, j;
229	greg	2.4
230			for (i = 0; i < 3; i++)
231	greg	2.5	for (j = 0; j < 3; j++)
232			mt[i][j] = m1[i][0]*m2[0][j] +
233			m1[i][1]*m2[1][j] +
234			m1[i][2]*m2[2][j] ;
235			cpcolormat(m3, mt);
236			}
237
238
239	greg	2.15	int
240			colorprimsOK( /* are color primaries reasonable? */
241			RGBPRIMS pr
242			)
243			{
244	greg	2.20	int i, j;
245	greg	2.15
246			for (i = 0; i < 4; i++) {
247	greg	2.19	if ((pr[i][CIEX] <= -.5) \| (pr[i][CIEY] <= -.5))
248	greg	2.15	return(0);
249	greg	2.18	if ((pr[i][CIEX] >= 1.5) \| (pr[i][CIEY] >= 1.5))
250	greg	2.17	return(0);
251	greg	2.18	if (pr[i][CIEX] + pr[i][CIEY] >= 1.5)
252	greg	2.15	return(0);
253			}
254	greg	2.20	for (i = 0; i < 4; i++)
255			for (j = i+1; j < 4; j++)
256			if (CEQ(pr[i][CIEX],pr[j][CIEX]) &&
257			CEQ(pr[i][CIEY],pr[j][CIEY]))
258			return(0);
259	greg	2.15	return(1);
260			}
261
262
263
264			int
265			compxyz2rgbmat( /* compute conversion from CIE to RGB space */
266			COLORMAT mat,
267			register RGBPRIMS pr
268			)
269	greg	2.5	{
270			double C_rD, C_gD, C_bD;
271
272	greg	2.6	if (pr == stdprims) { /* can use xyz2rgbmat */
273			cpcolormat(mat, xyz2rgbmat);
274	greg	2.16	return(1);
275	greg	2.6	}
276	greg	2.15	if (CEQ(pr[WHT][CIEX],0.) \| CEQ(pr[WHT][CIEY],0.))
277			return(0);
278	greg	2.5	C_rD = (1./pr[WHT][CIEY]) *
279			( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
280			pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
281			pr[GRN][CIEX]pr[BLU][CIEY] - pr[BLU][CIEX]pr[GRN][CIEY] ) ;
282	greg	2.15	if (CEQ(C_rD,0.))
283			return(0);
284	greg	2.5	C_gD = (1./pr[WHT][CIEY]) *
285			( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
286			pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
287			pr[RED][CIEX]pr[BLU][CIEY] + pr[BLU][CIEX]pr[RED][CIEY] ) ;
288	greg	2.15	if (CEQ(C_gD,0.))
289			return(0);
290	greg	2.5	C_bD = (1./pr[WHT][CIEY]) *
291			( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
292			pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
293			pr[RED][CIEX]pr[GRN][CIEY] - pr[GRN][CIEX]pr[RED][CIEY] ) ;
294	greg	2.15	if (CEQ(C_bD,0.))
295			return(0);
296	greg	2.5	mat[0][0] = (pr[GRN][CIEY] - pr[BLU][CIEY] -
297			pr[BLU][CIEX]*pr[GRN][CIEY] +
298			pr[BLU][CIEY]*pr[GRN][CIEX])/C_rD ;
299			mat[0][1] = (pr[BLU][CIEX] - pr[GRN][CIEX] -
300			pr[BLU][CIEX]*pr[GRN][CIEY] +
301			pr[GRN][CIEX]*pr[BLU][CIEY])/C_rD ;
302			mat[0][2] = (pr[GRN][CIEX]*pr[BLU][CIEY] -
303			pr[BLU][CIEX]*pr[GRN][CIEY])/C_rD ;
304			mat[1][0] = (pr[BLU][CIEY] - pr[RED][CIEY] -
305			pr[BLU][CIEY]*pr[RED][CIEX] +
306			pr[RED][CIEY]*pr[BLU][CIEX])/C_gD ;
307			mat[1][1] = (pr[RED][CIEX] - pr[BLU][CIEX] -
308			pr[RED][CIEX]*pr[BLU][CIEY] +
309			pr[BLU][CIEX]*pr[RED][CIEY])/C_gD ;
310			mat[1][2] = (pr[BLU][CIEX]*pr[RED][CIEY] -
311			pr[RED][CIEX]*pr[BLU][CIEY])/C_gD ;
312			mat[2][0] = (pr[RED][CIEY] - pr[GRN][CIEY] -
313			pr[RED][CIEY]*pr[GRN][CIEX] +
314			pr[GRN][CIEY]*pr[RED][CIEX])/C_bD ;
315			mat[2][1] = (pr[GRN][CIEX] - pr[RED][CIEX] -
316			pr[GRN][CIEX]*pr[RED][CIEY] +
317			pr[RED][CIEX]*pr[GRN][CIEY])/C_bD ;
318			mat[2][2] = (pr[RED][CIEX]*pr[GRN][CIEY] -
319			pr[GRN][CIEX]*pr[RED][CIEY])/C_bD ;
320	greg	2.15	return(1);
321	greg	2.5	}
322
323
324	greg	2.15	int
325			comprgb2xyzmat( /* compute conversion from RGB to CIE space */
326			COLORMAT mat,
327			register RGBPRIMS pr
328			)
329	greg	2.5	{
330			double C_rD, C_gD, C_bD, D;
331
332	greg	2.6	if (pr == stdprims) { /* can use rgb2xyzmat */
333			cpcolormat(mat, rgb2xyzmat);
334	greg	2.16	return(1);
335	greg	2.6	}
336	greg	2.15	if (CEQ(pr[WHT][CIEX],0.) \| CEQ(pr[WHT][CIEY],0.))
337			return(0);
338	greg	2.5	C_rD = (1./pr[WHT][CIEY]) *
339			( pr[WHT][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) -
340			pr[WHT][CIEY]*(pr[GRN][CIEX] - pr[BLU][CIEX]) +
341			pr[GRN][CIEX]pr[BLU][CIEY] - pr[BLU][CIEX]pr[GRN][CIEY] ) ;
342			C_gD = (1./pr[WHT][CIEY]) *
343			( pr[WHT][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) -
344			pr[WHT][CIEY]*(pr[BLU][CIEX] - pr[RED][CIEX]) -
345			pr[RED][CIEX]pr[BLU][CIEY] + pr[BLU][CIEX]pr[RED][CIEY] ) ;
346			C_bD = (1./pr[WHT][CIEY]) *
347			( pr[WHT][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) -
348			pr[WHT][CIEY]*(pr[RED][CIEX] - pr[GRN][CIEX]) +
349			pr[RED][CIEX]pr[GRN][CIEY] - pr[GRN][CIEX]pr[RED][CIEY] ) ;
350			D = pr[RED][CIEX]*(pr[GRN][CIEY] - pr[BLU][CIEY]) +
351			pr[GRN][CIEX]*(pr[BLU][CIEY] - pr[RED][CIEY]) +
352			pr[BLU][CIEX]*(pr[RED][CIEY] - pr[GRN][CIEY]) ;
353	greg	2.15	if (CEQ(D,0.))
354			return(0);
355	greg	2.5	mat[0][0] = pr[RED][CIEX]*C_rD/D;
356			mat[0][1] = pr[GRN][CIEX]*C_gD/D;
357			mat[0][2] = pr[BLU][CIEX]*C_bD/D;
358			mat[1][0] = pr[RED][CIEY]*C_rD/D;
359			mat[1][1] = pr[GRN][CIEY]*C_gD/D;
360			mat[1][2] = pr[BLU][CIEY]*C_bD/D;
361			mat[2][0] = (1.-pr[RED][CIEX]-pr[RED][CIEY])*C_rD/D;
362			mat[2][1] = (1.-pr[GRN][CIEX]-pr[GRN][CIEY])*C_gD/D;
363			mat[2][2] = (1.-pr[BLU][CIEX]-pr[BLU][CIEY])*C_bD/D;
364	greg	2.15	return(1);
365	greg	2.5	}
366
367
368	greg	2.15	int
369			comprgb2rgbmat( /* compute conversion from RGB1 to RGB2 */
370			COLORMAT mat,
371			RGBPRIMS pr1,
372			RGBPRIMS pr2
373			)
374	greg	2.5	{
375			COLORMAT pr1toxyz, xyztopr2;
376
377	greg	2.6	if (pr1 == pr2) {
378			mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
379			mat[0][1] = mat[0][2] = mat[1][0] =
380			mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
381	greg	2.16	return(1);
382	greg	2.6	}
383	greg	2.15	if (!comprgb2xyzmat(pr1toxyz, pr1))
384			return(0);
385			if (!compxyz2rgbmat(xyztopr2, pr2))
386			return(0);
387	greg	2.5	/* combine transforms */
388			multcolormat(mat, pr1toxyz, xyztopr2);
389	greg	2.15	return(1);
390	greg	2.11	}
391
392
393	greg	2.15	int
394			compxyzWBmat( /* CIE von Kries transform from wht1 to wht2 */
395			COLORMAT mat,
396			float wht1[2],
397			float wht2[2]
398			)
399	greg	2.11	{
400			COLOR cw1, cw2;
401			if (XYEQ(wht1,wht2)) {
402			mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
403			mat[0][1] = mat[0][2] = mat[1][0] =
404			mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
405	greg	2.15	return(1);
406	greg	2.11	}
407	greg	2.15	if (CEQ(wht1[CIEX],0.) \| CEQ(wht1[CIEY],0.))
408			return(0);
409	greg	2.11	cw1[RED] = wht1[CIEX]/wht1[CIEY];
410			cw1[GRN] = 1.;
411			cw1[BLU] = (1. - wht1[CIEX] - wht1[CIEY])/wht1[CIEY];
412			colortrans(cw1, vkmat, cw1);
413	greg	2.15	if (CEQ(wht2[CIEX],0.) \| CEQ(wht2[CIEY],0.))
414			return(0);
415	greg	2.11	cw2[RED] = wht2[CIEX]/wht2[CIEY];
416			cw2[GRN] = 1.;
417			cw2[BLU] = (1. - wht2[CIEX] - wht2[CIEY])/wht2[CIEY];
418			colortrans(cw2, vkmat, cw2);
419	greg	2.15	if (CEQ(cw1[RED],0.) \| CEQ(cw1[GRN],0.) \| CEQ(cw1[BLU],0.))
420			return(0);
421	greg	2.11	mat[0][0] = cw2[RED]/cw1[RED];
422			mat[1][1] = cw2[GRN]/cw1[GRN];
423			mat[2][2] = cw2[BLU]/cw1[BLU];
424			mat[0][1] = mat[0][2] = mat[1][0] =
425			mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
426			multcolormat(mat, vkmat, mat);
427			multcolormat(mat, mat, ivkmat);
428	greg	2.15	return(1);
429	greg	2.11	}
430
431
432	greg	2.15	int
433			compxyz2rgbWBmat( /* von Kries conversion from CIE to RGB space */
434			COLORMAT mat,
435			RGBPRIMS pr
436			)
437	greg	2.11	{
438			COLORMAT wbmat;
439
440	greg	2.15	if (!compxyz2rgbmat(mat, pr))
441			return(0);
442	greg	2.11	if (XYEQ(pr[WHT],xyneu))
443	greg	2.15	return(1);
444			if (!compxyzWBmat(wbmat, xyneu, pr[WHT]))
445			return(0);
446	greg	2.11	multcolormat(mat, wbmat, mat);
447	greg	2.15	return(1);
448	greg	2.11	}
449
450	greg	2.15	int
451			comprgb2xyzWBmat( /* von Kries conversion from RGB to CIE space */
452			COLORMAT mat,
453			RGBPRIMS pr
454			)
455	greg	2.11	{
456			COLORMAT wbmat;
457
458	greg	2.15	if (!comprgb2xyzmat(mat, pr))
459			return(0);
460	greg	2.11	if (XYEQ(pr[WHT],xyneu))
461	greg	2.15	return(1);
462			if (!compxyzWBmat(wbmat, pr[WHT], xyneu))
463			return(0);
464	greg	2.11	multcolormat(mat, mat, wbmat);
465	greg	2.15	return(1);
466	greg	2.11	}
467
468	greg	2.15	int
469			comprgb2rgbWBmat( /* von Kries conversion from RGB1 to RGB2 */
470			COLORMAT mat,
471			RGBPRIMS pr1,
472			RGBPRIMS pr2
473			)
474	greg	2.11	{
475			COLORMAT pr1toxyz, xyztopr2, wbmat;
476
477			if (pr1 == pr2) {
478			mat[0][0] = mat[1][1] = mat[2][2] = 1.0;
479			mat[0][1] = mat[0][2] = mat[1][0] =
480			mat[1][2] = mat[2][0] = mat[2][1] = 0.0;
481	greg	2.15	return(1);
482	greg	2.11	}
483	greg	2.15	if (!comprgb2xyzmat(pr1toxyz, pr1))
484			return(0);
485			if (!compxyzWBmat(wbmat, pr1[WHT], pr2[WHT]))
486			return(0);
487			if (!compxyz2rgbmat(xyztopr2, pr2))
488			return(0);
489	greg	2.11	/* combine transforms */
490			multcolormat(mat, pr1toxyz, wbmat);
491			multcolormat(mat, mat, xyztopr2);
492	greg	2.15	return(1);
493	greg	1.1	}