src/common/color.c

#ifndef lint
static const char       RCSid[] = "$Id: color.c,v 2.38 2024/09/20 17:39:12 greg Exp $";
#endif
/*
 *  color.c - routines for color calculations.
 *
 *  Externals declared in color.h
 */

#include "copyright.h"

#include  <stdio.h>
#include  <stdlib.h>
#include  <math.h>
#include  "color.h"

#ifdef getc_unlocked            /* avoid horrendous overhead of flockfile */
#undef getc
#undef putc
#undef ferror
#define getc    getc_unlocked
#define putc    putc_unlocked
#define ferror  ferror_unlocked
#endif

#define  MINELEN        17      /* minimum scanline length for encoding */
#define  MAXELEN        0x7fff  /* maximum scanline length for encoding */
#define  MINRUN         4       /* minimum run length */


int  CNDX[4] = {0,1,2,3};       /* RGBE indices for SCOLOR, SCOLR */
float  WLPART[4] = {780,588,480,380};   /* RGB wavelength limits+partitions (nm) */


int
setspectrsamp(                  /* assign spectral sampling, 1 if good, -1 if bad */
        int cn[4],              /* input cn[3]=nsamps */
        float wlpt[4]           /* input wlpt[0],wlpt[3]=extrema */
)
{
        static const float      PKWL[3] = {607, 553, 469};
        int                     i, j;

        if (cn[3] < 3)
                return(-1);             /* reject this */

        if (wlpt[0] < wlpt[3]) {
                float   tf = wlpt[0];
                wlpt[0] = wlpt[3]; wlpt[3] = tf;
        }
        if (wlpt[0] - wlpt[3] < 50.f)
                return(-1);             /* also reject */

        if (cn[3] > MAXCSAMP)
                cn[3] = MAXCSAMP;

        if ((wlpt[3] >= PKWL[2]) | (wlpt[0] <= PKWL[0])) {
                wlpt[1] = wlpt[0] + 0.333333f*(wlpt[3]-wlpt[0]);
                wlpt[2] = wlpt[0] + 0.666667f*(wlpt[3]-wlpt[0]);
                cn[0] = 0; cn[1] = cn[3]/3; cn[2] = cn[3]*2/3;
                return(0);              /* unhappy but non-fatal return value */
        }
        wlpt[1] = 588.f;                /* tuned for standard green channel */
        wlpt[2] = 480.f;
        if (cn[3] == 3) {               /* nothing to tune? */
                cn[0] = 0; cn[1] = 1; cn[2] = 2;
        } else {                        /* else find nearest color indices */
                double  curwl[3];
                memset(curwl, 0, sizeof(curwl));
                for (i = cn[3]; i--; ) {
                        const float     cwl = (i+.5f)/cn[3]*(wlpt[3]-wlpt[0]) + wlpt[0];
                        for (j = 3; j--; )
                                if (fabs(cwl - PKWL[j]) < fabs(curwl[j] - PKWL[j])) {
                                        curwl[j] = cwl;
                                        cn[j] = i;
                                }
                }
        }
        return(1);                      /* happy return value */
}


void
setscolor(                      /* assign spectral color from RGB */
        SCOLOR scol,
        double r,
        double g,
        double b
)
{
        double          step, cwl;
        int             i;

        if (NCSAMP == 3) {
                setcolor(scol, r, g, b);
                return;
        }
        step = (WLPART[3] - WLPART[0])/(double)NCSAMP;
        cwl = WLPART[0] + .5*step;
        for (i = 0; i < NCSAMP; i++) {
                if (cwl >= WLPART[1])
                        scol[i] = r;
                else if (cwl >= WLPART[2])
                        scol[i] = g;
                else
                        scol[i] = b;
                cwl += step;
        }
}


void
setscolr(                       /* assign common-exponent spectral from RGB */
        SCOLR sclr,
        double r,
        double g,
        double b
)
{
        SCOLOR  scol;

        setscolor(scol, r, g, b);
        scolor2scolr(sclr, scol, NCSAMP);
}


void
scolor2color(                   /* assign RGB color from spectrum */
        COLOR col,
        const SCOLOR scol,              /* uses average over bands */
        int ncs,
        const float wlpt[4]
)
{
        const double    step = (wlpt[3] - wlpt[0])/(double)ncs;
        double          cwl = wlpt[0] + .5*step;
        int             i, j=0, n=0;

        setcolor(col, 0, 0, 0);
        for (i = 0; i < ncs; i++) {
                if (cwl < wlpt[j+1]) {
                        if (n > 1) col[j] /= (COLORV)n;
                        j++;
                        n = 0;
                }
                col[j] += scol[i];
                n++;
                cwl += step;
        }
        if (n > 1) col[j] /= (COLORV)n;
}


void
scolr2colr(                     /* assign RGBE from common-exponent spectrum */
        COLR clr,
        const SCOLR sclr,
        int ncs,
        const float wlpt[4]
)
#if 0                   /* fancier method seems to be slower(!) */
{
        const double    step = (wlpt[3] - wlpt[0])/(double)ncs;
        double          cwl;
        int             csum[3], cnt[3], eshft;
        int             i, j;

        csum[0] = csum[1] = csum[2] = 0;
        cnt[0] = cnt[1] = cnt[2] = 0;
        cwl = wlpt[j=0] + .5*step;
        for (i = 0; i < ncs; i++) {
                csum[j] += sclr[i];
                ++cnt[j];
                j += ((cwl += step) < wlpt[j+1]);
        }
        eshft = 7;              /* compute exponent shift */
        for (j = 3; (eshft > 0) & (j-- > 0); ) {
                i = 0;
                while (csum[j] < 128*cnt[j] >> i)
                        if (++i >= eshft)
                                break;
                if (eshft > i)
                        eshft = i;
        }
        if (sclr[ncs] <= eshft) {
                clr[RED] = clr[GRN] = clr[BLU] = 0;
                clr[EXP] = 0;
                return;
        }
        for (j = 3; j--; )
                clr[j] = (csum[j]<<eshft)/cnt[j];

        clr[EXP] = sclr[ncs] - eshft;
}
#else
{
        SCOLOR  scol;
        COLOR   col;

        scolr2scolor(scol, sclr, ncs);
        scolor2color(col, scol, ncs, wlpt);
        setcolr(clr, col[RED], col[GRN], col[BLU]);
}
#endif


void
scolor2colr(                    /* assign RGBE from spectral color */
        COLR clr,
        const SCOLOR scol,              /* uses average over bands */
        int ncs,
        const float wlpt[4]
)
{
        COLOR   col;

        scolor2color(col, scol, ncs, wlpt);
        setcolr(clr, col[RED], col[GRN], col[BLU]);
}


void
scolr2color(                    /* assign RGB from common exponent */
        COLOR col,
        const SCOLR sclr,
        int ncs,
        const float wlpt[4]
)
{
        SCOLOR  scol;

        scolr2scolor(scol, sclr, ncs);
        scolor2color(col, scol, ncs, wlpt);
}


void
scolor2scolr(                   /* float spectrum to common exponent */
        SCOLR sclr,
        const SCOLOR scol,
        int ncs
)
{
        int     i = ncs;
        COLORV  p = scol[--i];

        while (i)
                if (scol[--i] > p)
                        p = scol[i];
        if (p <= 1e-32) {
                memset(sclr, 0, ncs+1);
                return;
        }
        p = frexp(p, &i) * 256.0 / p;
        sclr[ncs] = i + COLXS;
        for (i = ncs; i--; )
                sclr[i] = (scol[i] > 0) * (int)(scol[i]*p);
}


void
scolr2scolor(                   /* common exponent to float spectrum */
        SCOLOR scol,
        const SCOLR sclr,
        int ncs
)
{
        double  f;
        int     i;

        if (sclr[ncs] == 0) {
                memset(scol, 0, sizeof(COLORV)*ncs);
                return;
        }
        f = ldexp(1.0, (int)sclr[ncs]-(COLXS+8));

        for (i = ncs; i--; )
                scol[i] = (sclr[i] + 0.5)*f;
}


double
scolor_mean(                    /* compute average for spectral color */
        const SCOLOR  scol
)
{
        int     i = NCSAMP;
        double  sum = 0;

        while (i--)
                sum += scol[i];

        return sum/(double)NCSAMP;
}


double
sintens(                        /* find maximum value from spectrum */
        const SCOLOR  scol
)
{
        int     i = NCSAMP;
        COLORV  peak = scol[--i];

        while (i)
                if (scol[--i] > peak)
                        peak = scol[i];

        return peak;
}


void
convertscolor(                  /* spectrum conversion, zero-fill ends */
        SCOLOR dst,             /* destination spectrum */
        int dnc,                /* destination # of spectral samples/intervals */
        double dwl0,            /* starting destination wavelength (longer) */
        double dwl1,            /* ending destination wavelength (shorter) */
        const COLORV src[],     /* source spectrum array */
        int snc,
        double swl0,            /* long/short wavelengths may be reversed */
        double swl1
)
{
        const int       sdir = 1 - 2*(swl0 < swl1);
        const double    sstp = (swl1 - swl0)/(double)snc;
        const double    dstp = (dwl1 - dwl0)/(double)dnc;
        const double    rdstp = 1./dstp;
        int             si, ssi, di;
        double          wl;

        if ((dnc < 3) | (dwl0 <= dwl1) | (dst == src))
                return;         /* invalid destination */

        if (dnc == snc && (dwl0-swl0)*(dwl0-swl0) + (dwl1-swl1)*(dwl1-swl1) <= .5) {
                memcpy(dst, src, sizeof(COLORV)*dnc);
                return;         /* same spectral sampling */
        }
        memset(dst, 0, sizeof(COLORV)*dnc);
                                /* set starting positions */
        if ((sdir>0 ? swl0 : swl1) <= dwl0) {
                if (sdir > 0) {
                        wl = swl0;
                        ssi = 0;
                } else {
                        wl = swl1;
                        ssi = snc-1;
                }
                si = 0;
                di = (wl - dwl0)*rdstp;
        } else {
                wl = dwl0;
                if (sdir > 0) {
                        ssi = si = (wl - swl0)/sstp;
                } else {
                        si = (wl - swl1)/sstp;
                        ssi = snc-1 - si;
                }
                di = 0;
        }
        swl0 += (sdir < 0)*sstp;
                                /* step through intervals */
        while ((si < snc) & (di < dnc)) {
                double  intvl;
                if (swl0 + (ssi+sdir)*sstp < dwl0 + (di+1)*dstp) {
                        intvl = dwl0 + (di+1)*dstp - wl;
                        dst[di++] += src[ssi]*intvl*rdstp;
                } else {
                        intvl = swl0 + (ssi+sdir)*sstp - wl;
                        dst[di] += src[ssi]*intvl*rdstp;
                        ssi += sdir;
                        si++;
                }
                wl += intvl;
        }
}


void *
tempbuffer(                     /* get a temporary buffer */
        size_t  len
)
{
        static void     *tempbuf = NULL;
        static size_t   tempbuflen = 0;

        if (!len) {             /* call to free */
                if (tempbuflen) {
                        free(tempbuf);
                        tempbuf = NULL;
                        tempbuflen = 0;
                }
                return(NULL);
        }
        if (len <= tempbuflen)  /* big enough already? */
                return(tempbuf);
                                /* else free & reallocate */
        if (tempbuflen)
                free(tempbuf);
        tempbuf = malloc(len);
        tempbuflen = len*(tempbuf != NULL);
        return(tempbuf);
}


int
fwritecolrs(                    /* write out a colr scanline */
        COLR  *scanline,
        int  len,
        FILE  *fp
)
{
        int  i, j, beg, cnt = 1;
        int  c2;
        
        if ((len < MINELEN) | (len > MAXELEN))  /* OOBs, write out flat */
                return(fwrite((char *)scanline,sizeof(COLR),len,fp) - len);
                                        /* put magic header */
        putc(2, fp);
        putc(2, fp);
        putc(len>>8, fp);
        putc(len&0xff, fp);
                                        /* put components seperately */
        for (i = 0; i < 4; i++) {
            for (j = 0; j < len; j += cnt) {    /* find next run */
                for (beg = j; beg < len; beg += cnt) {
                    for (cnt = 1; (cnt < 127) & (beg+cnt < len) &&
                            scanline[beg+cnt][i] == scanline[beg][i]; cnt++)
                        ;
                    if (cnt >= MINRUN)
                        break;                  /* long enough */
                }
                if ((beg-j > 1) & (beg-j < MINRUN)) {
                    c2 = j+1;
                    while (scanline[c2++][i] == scanline[j][i])
                        if (c2 == beg) {        /* short run */
                            putc(128+beg-j, fp);
                            putc(scanline[j][i], fp);
                            j = beg;
                            break;
                        }
                }
                while (j < beg) {               /* write out non-run */
                    if ((c2 = beg-j) > 128) c2 = 128;
                    putc(c2, fp);
                    while (c2--)
                        putc(scanline[j++][i], fp);
                }
                if (cnt >= MINRUN) {            /* write out run */
                    putc(128+cnt, fp);
                    putc(scanline[beg][i], fp);
                } else
                    cnt = 0;
            }
        }
        return(ferror(fp) ? -1 : 0);
}

/*
 * An old-format scanline is either a stream of valid RGBE or XYZE real
 * pixels or at least one real pixel followed by some number of
 * invalid real pixels of the form (1,1,1,n), where n is a count.
 * These can themselves be repeated to create a multibyte repeat
 * count, with the least significant byte first (little-endian order.)
 * Repeat counts are limited by the size of an int; if a repetition
 * leads to an overrun, the rest of the the repetition will be
 * silently ignored.
 */
static int
oldreadcolrs(                   /* read in an old-style colr scanline */
        COLR  *scanline,
        int  len,
        FILE  *fp
)
{
        int  rshift = 0;
        int  i;
        
        while (len > 0) {
                scanline[0][RED] = getc(fp);
                scanline[0][GRN] = getc(fp);
                scanline[0][BLU] = getc(fp);
                scanline[0][EXP] = i = getc(fp);
                if (i == EOF)
                        return(-1);
                if (scanline[0][GRN] == 1 &&
                                (scanline[0][RED] == 1) &
                                (scanline[0][BLU] == 1)) {
                        i = scanline[0][EXP] << rshift;
                        while (i--) {
                                copycolr(scanline[0], scanline[-1]);
                                if (--len <= 0)
                                        return(0);
                                scanline++;
                        }
                        rshift += 8;
                } else {
                        scanline++;
                        len--;
                        rshift = 0;
                }
        }
        return(0);
}

/* 
 * There are two scanline formats: old and new.  The old format
 * compresses runs of RGBE or XYZE four-byte real pixels; the new
 * format breaks the pixels into R, G, B, and E lines (or XYZE lines)
 * which are individually run-length encoded.
 *
 * An old-format scanline always begins with a valid real pixel; at
 * least one of the RGB (or XYZ) values will have its high-order bit
 * set.  A new-format scanline begins with four bytes which are not a
 * valid real pixel: (2, 2, lenhigh, lenlow) where lenhigh is always
 * less than 128 and hence never has a high-order bit set.
 *
 * A new-format scanline is broken into its RGBE or XYZE components.
 * Each is output and run-length encoded separately so that a scanline
 * is broken into four records.  In turn, each record is organized
 * into chunks of up to 128 characters, which begin with a count byte.
 * If the count byte is greater than 128, the following data byte is
 * repeated (count-128) times.  If not, the count byte is followed by
 * that many data bytes.
 */
int
freadcolrs(                     /* read in an encoded colr scanline */
        COLR  *scanline,
        int  len,
        FILE  *fp
)
{
        int  i, j;
        int  code, val;
                                        /* determine scanline type */
        if (len <= 0)
                return(0);
        if ((i = getc(fp)) == EOF)
                return(-1);
        scanline[0][RED] = i;
        scanline[0][GRN] = getc(fp);
        scanline[0][BLU] = getc(fp);
        if ((i = getc(fp)) == EOF)
                return(-1);
        if ((scanline[0][RED] != 2) | (scanline[0][GRN] != 2) |
                        (scanline[0][BLU] & 0x80)) {
                scanline[0][EXP] = i;
                return(oldreadcolrs(scanline+1, len-1, fp));
        }
        if ((scanline[0][BLU]<<8 | i) != len)
                return(-1);             /* length mismatch! */
                                        /* read each component */
        for (i = 0; i < 4; i++)
            for (j = 0; j < len; ) {
                if ((code = getc(fp)) == EOF)
                    return(-1);
                if (code > 128) {       /* run */
                    code &= 127;
                    if ((val = getc(fp)) == EOF)
                        return -1;
                    if (j + code > len)
                        return -1;      /* overrun */
                    while (code--)
                        scanline[j++][i] = val;
                } else {                /* non-run */
                    if (j + code > len)
                        return -1;      /* overrun */
                    while (code--) {
                        if ((val = getc(fp)) == EOF)
                            return -1;
                        scanline[j++][i] = val;
                    }
                }
            }
        return(0);
}


/* read an nc-component common-exponent color scanline */
int
freadscolrs(COLRV *scanline, int nc, int len, FILE *fp)
{
        if (nc < 3)
                return(-1);
        if (nc == 3)
                return(freadcolrs((COLR *)scanline, len, fp));

        if (fread(scanline, nc+1, len, fp) != len)
                return(-1);
        return(0);
}


/* read nc-component common-exponent color scan and convert to COLR's */
int
fread2colrs(COLR *scanline, int len, FILE *fp, int nc, const float wlpt[4])
{
        COLRV  *sclrscan;
        int  n;

        if (nc < 3)
                return(-1);
        if (nc == 3)
                return(freadcolrs(scanline, len, fp));

        sclrscan = (COLRV *)tempbuffer(sizeof(COLRV)*(nc+1)*len);
        if (sclrscan == NULL || freadscolrs(sclrscan, nc, len, fp) < 0)
                return(-1);
        for (n = len; n--; ) {
                scolr2colr(*scanline++, sclrscan, nc, wlpt);
                sclrscan += nc+1;
        }
        return(0);
}


/* write an common-exponent spectral color scanline */
int
fwritescolrs(const COLRV *sscanline, int nc, int len, FILE *fp)
{
        if (nc < 3)
                return(-1);
        if (nc == 3)
                return(fwritecolrs((COLR *)sscanline, len, fp));

        if (fwrite(sscanline, nc+1, len, fp) != len)
                return(-1);
        return(0);
}


int
fwritescan(             /* write out an RGB or XYZ scanline */
        COLOR  *scanline,
        int  len,
        FILE  *fp
)
{
        COLR  *clrscan;
        int  n;
        COLR  *sp;
                                        /* get scanline buffer */
        if ((sp = (COLR *)tempbuffer(len*sizeof(COLR))) == NULL)
                return(-1);
        clrscan = sp;
                                        /* convert scanline */
        n = len;
        while (n-- > 0) {
                setcolr(sp[0], scanline[0][RED],
                                  scanline[0][GRN],
                                  scanline[0][BLU]);
                scanline++;
                sp++;
        }
        return(fwritecolrs(clrscan, len, fp));
}


int
freadscan(              /* read in an RGB or XYZ scanline */
        COLOR  *scanline,
        int  len,
        FILE  *fp
)
{
        COLR  *clrscan;

        if ((clrscan = (COLR *)tempbuffer(len*sizeof(COLR))) == NULL)
                return(-1);
        if (freadcolrs(clrscan, len, fp) < 0)
                return(-1);
                                        /* convert scanline */
        colr_color(scanline[0], clrscan[0]);
        while (--len > 0) {
                scanline++; clrscan++;
                if (clrscan[0][GRN] == clrscan[-1][GRN] &&
                            (clrscan[0][RED] == clrscan[-1][RED]) &
                            (clrscan[0][BLU] == clrscan[-1][BLU]) &
                            (clrscan[0][EXP] == clrscan[-1][EXP]))
                        copycolor(scanline[0], scanline[-1]);
                else
                        colr_color(scanline[0], clrscan[0]);
        }
        return(0);
}


/* read an nc-component color scanline */
int
freadsscan(COLORV *sscanline, int nc, int len, FILE *fp)
{
        COLRV   *tscn = (COLRV *)tempbuffer((nc+1)*len);
        int     i;

        if (tscn == NULL || freadscolrs(tscn, nc, len, fp) < 0)
                return(-1);
        for (i = len; i-- > 0; ) {
                scolr2scolor(sscanline, tscn, nc);
                sscanline += nc;
                tscn += nc+1;
        }
        return(0);
}


/* read an nc-component color scanline and return as RGB */
int
fread2scan(COLOR *scanline, int len, FILE *fp, int nc, const float wlpt[4])
{
        COLRV   *tscn;
        int     i;

        if (nc < 3)
                return(-1);
        if (nc == 3)
                return(freadscan(scanline, len, fp));

        tscn = (COLRV *)tempbuffer((nc+1)*len);
        if (tscn == NULL || freadscolrs(tscn, nc, len, fp) < 0)
                return(-1);
        for (i = len; i-- > 0; ) {
                scolr2color(*scanline++, tscn, nc, wlpt);
                tscn += nc+1;
        }
        return(0);
}


/* write an nc-component spectral color scanline */
int
fwritesscan(const COLORV *sscanline, int nc, int len, FILE *fp)
{
        COLRV   *tscn = (COLRV *)tempbuffer((nc+1)*len);
        int     i;

        if (tscn == NULL)
                return(-1);
        for (i = 0; i < len; i++) {
                scolor2scolr(tscn+i*(nc+1), sscanline, nc);
                sscanline += nc;
        }
        return(fwritescolrs(tscn, nc, len, fp));
}


void
setcolr(                        /* assign a short color value */
        COLR  clr,
        double  r,
        double  g,
        double  b
)
{
        double  d;
        int  e;
        
        d = r > g ? r : g;
        if (b > d) d = b;

        if (d <= 1e-32) {
                clr[RED] = clr[GRN] = clr[BLU] = 0;
                clr[EXP] = 0;
                return;
        }

        d = frexp(d, &e) * 256.0 / d;

        clr[RED] = (r > 0) * (int)(r*d);
        clr[GRN] = (g > 0) * (int)(g*d);
        clr[BLU] = (b > 0) * (int)(b*d);
        clr[EXP] = e + COLXS;
}


void
colr_color(                     /* convert short to float color */
        COLOR  col,
        const COLR  clr
)
{
        double  f;
        
        if (clr[EXP] == 0) {
                col[RED] = col[GRN] = col[BLU] = 0.0;
                return;
        }
        f = ldexp(1.0, (int)clr[EXP]-(COLXS+8));
        col[RED] = (clr[RED] + 0.5)*f;
        col[GRN] = (clr[GRN] + 0.5)*f;
        col[BLU] = (clr[BLU] + 0.5)*f;
}


int
bigdiff(                                /* c1 delta c2 > md? */
        const COLOR  c1,
        const COLOR  c2,
        double  md
)
{
        int  i;

        for (i = 0; i < 3; i++)
                if ((colval(c1,i)-colval(c2,i) > md*colval(c2,i)) |
                                (colval(c2,i)-colval(c1,i) > md*colval(c1,i)))
                        return(1);
        return(0);
}


int
sbigsdiff(                              /* sc1 delta sc2 > md? */
        const SCOLOR  c1,
        const SCOLOR  c2,
        double  md
)
{
        int  i = NCSAMP;

        while (i--)
                if ((c1[i]-c2[i] > md*c2[i]) | (c2[i]-c1[i] > md*c1[i]))
                        return(1);
        return(0);
}
Revision:	2.39
Committed:	Thu Feb 6 02:15:45 2025 UTC (2 months, 3 weeks ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.38:	+8 -3 lines
Log Message:	perf: Short-cut for common case of NCSAMP==3
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.39	static const char RCSid[] = "$Id: color.c,v 2.38 2024/09/20 17:39:12 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* color.c - routines for color calculations.
6			*
7	greg	2.9	* Externals declared in color.h
8			*/
9
10	greg	2.10	#include "copyright.h"
11	greg	1.1
12			#include <stdio.h>
13	greg	2.9	#include <stdlib.h>
14	greg	2.7	#include <math.h>
15	greg	1.1	#include "color.h"
16	greg	2.14
17			#ifdef getc_unlocked /* avoid horrendous overhead of flockfile */
18	greg	2.15	#undef getc
19			#undef putc
20	greg	2.19	#undef ferror
21	greg	2.14	#define getc getc_unlocked
22			#define putc putc_unlocked
23	greg	2.18	#define ferror ferror_unlocked
24	greg	2.14	#endif
25	greg	1.1
26	greg	2.34	#define MINELEN 17 /* minimum scanline length for encoding */
27	greg	2.6	#define MAXELEN 0x7fff /* maximum scanline length for encoding */
28	greg	1.14	#define MINRUN 4 /* minimum run length */
29
30	greg	2.4
31	greg	2.27	int CNDX[4] = {0,1,2,3}; /* RGBE indices for SCOLOR, SCOLR */
32			float WLPART[4] = {780,588,480,380}; /* RGB wavelength limits+partitions (nm) */
33
34
35			int
36			setspectrsamp( /* assign spectral sampling, 1 if good, -1 if bad */
37			int cn[4], /* input cn[3]=nsamps */
38			float wlpt[4] /* input wlpt[0],wlpt[3]=extrema */
39			)
40			{
41			static const float PKWL[3] = {607, 553, 469};
42			int i, j;
43
44			if (cn[3] < 3)
45			return(-1); /* reject this */
46
47			if (wlpt[0] < wlpt[3]) {
48			float tf = wlpt[0];
49			wlpt[0] = wlpt[3]; wlpt[3] = tf;
50			}
51			if (wlpt[0] - wlpt[3] < 50.f)
52			return(-1); /* also reject */
53
54			if (cn[3] > MAXCSAMP)
55			cn[3] = MAXCSAMP;
56
57			if ((wlpt[3] >= PKWL[2]) \| (wlpt[0] <= PKWL[0])) {
58			wlpt[1] = wlpt[0] + 0.333333f*(wlpt[3]-wlpt[0]);
59			wlpt[2] = wlpt[0] + 0.666667f*(wlpt[3]-wlpt[0]);
60			cn[0] = 0; cn[1] = cn[3]/3; cn[2] = cn[3]*2/3;
61			return(0); /* unhappy but non-fatal return value */
62			}
63			wlpt[1] = 588.f; /* tuned for standard green channel */
64			wlpt[2] = 480.f;
65			if (cn[3] == 3) { /* nothing to tune? */
66			cn[0] = 0; cn[1] = 1; cn[2] = 2;
67			} else { /* else find nearest color indices */
68			double curwl[3];
69			memset(curwl, 0, sizeof(curwl));
70			for (i = cn[3]; i--; ) {
71			const float cwl = (i+.5f)/cn[3]*(wlpt[3]-wlpt[0]) + wlpt[0];
72			for (j = 3; j--; )
73			if (fabs(cwl - PKWL[j]) < fabs(curwl[j] - PKWL[j])) {
74			curwl[j] = cwl;
75			cn[j] = i;
76			}
77			}
78			}
79			return(1); /* happy return value */
80			}
81
82
83			void
84			setscolor( /* assign spectral color from RGB */
85			SCOLOR scol,
86			double r,
87			double g,
88			double b
89			)
90			{
91	greg	2.39	double step, cwl;
92	greg	2.27	int i;
93
94	greg	2.39	if (NCSAMP == 3) {
95			setcolor(scol, r, g, b);
96			return;
97			}
98			step = (WLPART[3] - WLPART[0])/(double)NCSAMP;
99			cwl = WLPART[0] + .5*step;
100	greg	2.27	for (i = 0; i < NCSAMP; i++) {
101			if (cwl >= WLPART[1])
102			scol[i] = r;
103			else if (cwl >= WLPART[2])
104			scol[i] = g;
105			else
106			scol[i] = b;
107			cwl += step;
108			}
109			}
110
111
112			void
113	greg	2.38	setscolr( /* assign common-exponent spectral from RGB */
114			SCOLR sclr,
115			double r,
116			double g,
117			double b
118			)
119			{
120			SCOLOR scol;
121
122			setscolor(scol, r, g, b);
123			scolor2scolr(sclr, scol, NCSAMP);
124			}
125
126
127			void
128	greg	2.27	scolor2color( /* assign RGB color from spectrum */
129			COLOR col,
130	greg	2.35	const SCOLOR scol, /* uses average over bands */
131	greg	2.27	int ncs,
132	greg	2.30	const float wlpt[4]
133	greg	2.27	)
134			{
135			const double step = (wlpt[3] - wlpt[0])/(double)ncs;
136			double cwl = wlpt[0] + .5*step;
137			int i, j=0, n=0;
138
139			setcolor(col, 0, 0, 0);
140			for (i = 0; i < ncs; i++) {
141			if (cwl < wlpt[j+1]) {
142			if (n > 1) col[j] /= (COLORV)n;
143			j++;
144			n = 0;
145			}
146			col[j] += scol[i];
147			n++;
148			cwl += step;
149			}
150			if (n > 1) col[j] /= (COLORV)n;
151			}
152
153
154			void
155	greg	2.36	scolr2colr( /* assign RGBE from common-exponent spectrum */
156			COLR clr,
157			const SCOLR sclr,
158			int ncs,
159			const float wlpt[4]
160			)
161	greg	2.37	#if 0 /* fancier method seems to be slower(!) */
162	greg	2.36	{
163			const double step = (wlpt[3] - wlpt[0])/(double)ncs;
164			double cwl;
165			int csum[3], cnt[3], eshft;
166			int i, j;
167
168			csum[0] = csum[1] = csum[2] = 0;
169			cnt[0] = cnt[1] = cnt[2] = 0;
170			cwl = wlpt[j=0] + .5*step;
171			for (i = 0; i < ncs; i++) {
172			csum[j] += sclr[i];
173			++cnt[j];
174			j += ((cwl += step) < wlpt[j+1]);
175			}
176			eshft = 7; /* compute exponent shift */
177			for (j = 3; (eshft > 0) & (j-- > 0); ) {
178			i = 0;
179			while (csum[j] < 128*cnt[j] >> i)
180			if (++i >= eshft)
181			break;
182			if (eshft > i)
183			eshft = i;
184			}
185			if (sclr[ncs] <= eshft) {
186			clr[RED] = clr[GRN] = clr[BLU] = 0;
187			clr[EXP] = 0;
188			return;
189			}
190			for (j = 3; j--; )
191			clr[j] = (csum[j]<<eshft)/cnt[j];
192
193			clr[EXP] = sclr[ncs] - eshft;
194			}
195	greg	2.37	#else
196			{
197			SCOLOR scol;
198			COLOR col;
199
200			scolr2scolor(scol, sclr, ncs);
201			scolor2color(col, scol, ncs, wlpt);
202			setcolr(clr, col[RED], col[GRN], col[BLU]);
203			}
204			#endif
205	greg	2.36
206
207			void
208	greg	2.27	scolor2colr( /* assign RGBE from spectral color */
209			COLR clr,
210	greg	2.35	const SCOLOR scol, /* uses average over bands */
211	greg	2.27	int ncs,
212	greg	2.30	const float wlpt[4]
213	greg	2.27	)
214			{
215			COLOR col;
216
217			scolor2color(col, scol, ncs, wlpt);
218			setcolr(clr, col[RED], col[GRN], col[BLU]);
219			}
220
221
222			void
223	greg	2.33	scolr2color( /* assign RGB from common exponent */
224			COLOR col,
225	greg	2.35	const SCOLR sclr,
226	greg	2.33	int ncs,
227			const float wlpt[4]
228			)
229			{
230			SCOLOR scol;
231
232			scolr2scolor(scol, sclr, ncs);
233			scolor2color(col, scol, ncs, wlpt);
234			}
235
236
237			void
238	greg	2.27	scolor2scolr( /* float spectrum to common exponent */
239			SCOLR sclr,
240	greg	2.35	const SCOLOR scol,
241	greg	2.27	int ncs
242			)
243			{
244			int i = ncs;
245			COLORV p = scol[--i];
246
247			while (i)
248			if (scol[--i] > p)
249			p = scol[i];
250			if (p <= 1e-32) {
251			memset(sclr, 0, ncs+1);
252			return;
253			}
254			p = frexp(p, &i) * 256.0 / p;
255			sclr[ncs] = i + COLXS;
256			for (i = ncs; i--; )
257			sclr[i] = (scol[i] > 0) * (int)(scol[i]*p);
258			}
259
260
261			void
262			scolr2scolor( /* common exponent to float spectrum */
263			SCOLOR scol,
264	greg	2.35	const SCOLR sclr,
265	greg	2.27	int ncs
266			)
267			{
268			double f;
269			int i;
270
271			if (sclr[ncs] == 0) {
272			memset(scol, 0, sizeof(COLORV)*ncs);
273			return;
274			}
275			f = ldexp(1.0, (int)sclr[ncs]-(COLXS+8));
276
277			for (i = ncs; i--; )
278			scol[i] = (sclr[i] + 0.5)*f;
279			}
280
281
282			double
283			scolor_mean( /* compute average for spectral color */
284	greg	2.35	const SCOLOR scol
285	greg	2.27	)
286			{
287			int i = NCSAMP;
288			double sum = 0;
289
290			while (i--)
291			sum += scol[i];
292
293			return sum/(double)NCSAMP;
294			}
295
296
297			double
298			sintens( /* find maximum value from spectrum */
299	greg	2.35	const SCOLOR scol
300	greg	2.27	)
301			{
302			int i = NCSAMP;
303			COLORV peak = scol[--i];
304
305			while (i)
306			if (scol[--i] > peak)
307			peak = scol[i];
308
309			return peak;
310			}
311
312
313			void
314			convertscolor( /* spectrum conversion, zero-fill ends */
315			SCOLOR dst, /* destination spectrum */
316			int dnc, /* destination # of spectral samples/intervals */
317			double dwl0, /* starting destination wavelength (longer) */
318			double dwl1, /* ending destination wavelength (shorter) */
319			const COLORV src[], /* source spectrum array */
320			int snc,
321			double swl0, /* long/short wavelengths may be reversed */
322			double swl1
323			)
324			{
325			const int sdir = 1 - 2*(swl0 < swl1);
326			const double sstp = (swl1 - swl0)/(double)snc;
327			const double dstp = (dwl1 - dwl0)/(double)dnc;
328			const double rdstp = 1./dstp;
329			int si, ssi, di;
330			double wl;
331
332			if ((dnc < 3) \| (dwl0 <= dwl1) \| (dst == src))
333			return; /* invalid destination */
334
335			if (dnc == snc && (dwl0-swl0)(dwl0-swl0) + (dwl1-swl1)(dwl1-swl1) <= .5) {
336			memcpy(dst, src, sizeof(COLORV)*dnc);
337			return; /* same spectral sampling */
338			}
339			memset(dst, 0, sizeof(COLORV)*dnc);
340			/* set starting positions */
341			if ((sdir>0 ? swl0 : swl1) <= dwl0) {
342			if (sdir > 0) {
343			wl = swl0;
344			ssi = 0;
345			} else {
346			wl = swl1;
347			ssi = snc-1;
348			}
349			si = 0;
350			di = (wl - dwl0)*rdstp;
351			} else {
352			wl = dwl0;
353	greg	2.29	if (sdir > 0) {
354			ssi = si = (wl - swl0)/sstp;
355			} else {
356			si = (wl - swl1)/sstp;
357			ssi = snc-1 - si;
358			}
359	greg	2.27	di = 0;
360			}
361			swl0 += (sdir < 0)*sstp;
362			/* step through intervals */
363			while ((si < snc) & (di < dnc)) {
364			double intvl;
365			if (swl0 + (ssi+sdir)sstp < dwl0 + (di+1)dstp) {
366			intvl = dwl0 + (di+1)*dstp - wl;
367			dst[di++] += src[ssi]intvlrdstp;
368			} else {
369			intvl = swl0 + (ssi+sdir)*sstp - wl;
370			dst[di] += src[ssi]intvlrdstp;
371			ssi += sdir;
372			si++;
373			}
374			wl += intvl;
375			}
376			}
377
378
379	greg	2.22	void *
380	greg	2.17	tempbuffer( /* get a temporary buffer */
381	greg	2.27	size_t len
382	greg	2.17	)
383	greg	1.14	{
384	greg	2.27	static void *tempbuf = NULL;
385			static size_t tempbuflen = 0;
386	greg	1.14
387	greg	2.22	if (!len) { /* call to free */
388			if (tempbuflen) {
389	greg	2.20	free(tempbuf);
390	greg	2.22	tempbuf = NULL;
391			tempbuflen = 0;
392			}
393			return(NULL);
394	greg	1.14	}
395	greg	2.22	if (len <= tempbuflen) /* big enough already? */
396			return(tempbuf);
397			/* else free & reallocate */
398			if (tempbuflen)
399			free(tempbuf);
400			tempbuf = malloc(len);
401			tempbuflen = len*(tempbuf != NULL);
402	greg	1.14	return(tempbuf);
403			}
404
405
406	greg	2.9	int
407	greg	2.17	fwritecolrs( /* write out a colr scanline */
408			COLR *scanline,
409			int len,
410			FILE *fp
411			)
412	greg	1.1	{
413	greg	2.17	int i, j, beg, cnt = 1;
414	greg	1.14	int c2;
415	greg	1.1
416	schorsch	2.13	if ((len < MINELEN) \| (len > MAXELEN)) /* OOBs, write out flat */
417	greg	1.14	return(fwrite((char *)scanline,sizeof(COLR),len,fp) - len);
418	greg	2.2	/* put magic header */
419			putc(2, fp);
420	greg	1.14	putc(2, fp);
421			putc(len>>8, fp);
422	greg	2.21	putc(len&0xff, fp);
423	greg	1.14	/* put components seperately */
424			for (i = 0; i < 4; i++) {
425			for (j = 0; j < len; j += cnt) { /* find next run */
426			for (beg = j; beg < len; beg += cnt) {
427	greg	2.20	for (cnt = 1; (cnt < 127) & (beg+cnt < len) &&
428	greg	1.14	scanline[beg+cnt][i] == scanline[beg][i]; cnt++)
429			;
430			if (cnt >= MINRUN)
431			break; /* long enough */
432	greg	1.1	}
433	greg	2.20	if ((beg-j > 1) & (beg-j < MINRUN)) {
434	greg	1.15	c2 = j+1;
435			while (scanline[c2++][i] == scanline[j][i])
436			if (c2 == beg) { /* short run */
437			putc(128+beg-j, fp);
438			putc(scanline[j][i], fp);
439			j = beg;
440			break;
441			}
442			}
443			while (j < beg) { /* write out non-run */
444	greg	1.14	if ((c2 = beg-j) > 128) c2 = 128;
445			putc(c2, fp);
446			while (c2--)
447			putc(scanline[j++][i], fp);
448			}
449			if (cnt >= MINRUN) { /* write out run */
450			putc(128+cnt, fp);
451			putc(scanline[beg][i], fp);
452			} else
453			cnt = 0;
454			}
455	greg	1.1	}
456			return(ferror(fp) ? -1 : 0);
457			}
458
459	greg	2.23	/*
460			* An old-format scanline is either a stream of valid RGBE or XYZE real
461			* pixels or at least one real pixel followed by some number of
462			* invalid real pixels of the form (1,1,1,n), where n is a count.
463			* These can themselves be repeated to create a multibyte repeat
464			* count, with the least significant byte first (little-endian order.)
465			* Repeat counts are limited by the size of an int; if a repetition
466			* leads to an overrun, the rest of the the repetition will be
467			* silently ignored.
468			*/
469	greg	2.9	static int
470	greg	2.19	oldreadcolrs( /* read in an old-style colr scanline */
471	greg	2.17	COLR *scanline,
472			int len,
473			FILE *fp
474			)
475	greg	2.9	{
476	greg	2.25	int rshift = 0;
477	greg	2.17	int i;
478	greg	2.9
479			while (len > 0) {
480			scanline[0][RED] = getc(fp);
481			scanline[0][GRN] = getc(fp);
482			scanline[0][BLU] = getc(fp);
483	greg	2.18	scanline[0][EXP] = i = getc(fp);
484			if (i == EOF)
485	greg	2.9	return(-1);
486	greg	2.20	if (scanline[0][GRN] == 1 &&
487			(scanline[0][RED] == 1) &
488	greg	2.25	(scanline[0][BLU] == 1)) {
489	greg	2.20	i = scanline[0][EXP] << rshift;
490			while (i--) {
491	greg	2.9	copycolr(scanline[0], scanline[-1]);
492	greg	2.20	if (--len <= 0)
493			return(0);
494	greg	2.9	scanline++;
495			}
496			rshift += 8;
497			} else {
498			scanline++;
499			len--;
500			rshift = 0;
501			}
502			}
503			return(0);
504			}
505
506	greg	2.23	/*
507			* There are two scanline formats: old and new. The old format
508			* compresses runs of RGBE or XYZE four-byte real pixels; the new
509			* format breaks the pixels into R, G, B, and E lines (or XYZE lines)
510			* which are individually run-length encoded.
511			*
512			* An old-format scanline always begins with a valid real pixel; at
513			* least one of the RGB (or XYZ) values will have its high-order bit
514			* set. A new-format scanline begins with four bytes which are not a
515			* valid real pixel: (2, 2, lenhigh, lenlow) where lenhigh is always
516			* less than 128 and hence never has a high-order bit set.
517			*
518			* A new-format scanline is broken into its RGBE or XYZE components.
519			* Each is output and run-length encoded separately so that a scanline
520			* is broken into four records. In turn, each record is organized
521			* into chunks of up to 128 characters, which begin with a count byte.
522			* If the count byte is greater than 128, the following data byte is
523			* repeated (count-128) times. If not, the count byte is followed by
524			* that many data bytes.
525			*/
526	greg	2.9	int
527	greg	2.17	freadcolrs( /* read in an encoded colr scanline */
528			COLR *scanline,
529			int len,
530			FILE *fp
531			)
532	greg	1.1	{
533	greg	2.17	int i, j;
534	greg	2.6	int code, val;
535	greg	1.14	/* determine scanline type */
536	greg	2.26	if (len <= 0)
537			return(0);
538	greg	1.14	if ((i = getc(fp)) == EOF)
539			return(-1);
540	greg	2.26	scanline[0][RED] = i;
541	greg	1.14	scanline[0][GRN] = getc(fp);
542			scanline[0][BLU] = getc(fp);
543			if ((i = getc(fp)) == EOF)
544			return(-1);
545	greg	2.26	if ((scanline[0][RED] != 2) \| (scanline[0][GRN] != 2) \|
546			(scanline[0][BLU] & 0x80)) {
547	greg	1.14	scanline[0][EXP] = i;
548			return(oldreadcolrs(scanline+1, len-1, fp));
549			}
550			if ((scanline[0][BLU]<<8 \| i) != len)
551			return(-1); /* length mismatch! */
552			/* read each component */
553			for (i = 0; i < 4; i++)
554			for (j = 0; j < len; ) {
555			if ((code = getc(fp)) == EOF)
556			return(-1);
557			if (code > 128) { /* run */
558	greg	2.6	code &= 127;
559	greg	2.9	if ((val = getc(fp)) == EOF)
560			return -1;
561	greg	2.16	if (j + code > len)
562			return -1; /* overrun */
563	greg	2.6	while (code--)
564			scanline[j++][i] = val;
565	greg	2.16	} else { /* non-run */
566			if (j + code > len)
567			return -1; /* overrun */
568	greg	2.9	while (code--) {
569			if ((val = getc(fp)) == EOF)
570			return -1;
571			scanline[j++][i] = val;
572			}
573	greg	2.16	}
574	greg	1.14	}
575	greg	1.1	return(0);
576			}
577
578
579	greg	2.30	/* read an nc-component common-exponent color scanline */
580			int
581	greg	2.35	freadscolrs(COLRV scanline, int nc, int len, FILE fp)
582	greg	2.30	{
583	greg	2.32	if (nc < 3)
584			return(-1);
585			if (nc == 3)
586			return(freadcolrs((COLR *)scanline, len, fp));
587
588	greg	2.30	if (fread(scanline, nc+1, len, fp) != len)
589			return(-1);
590			return(0);
591			}
592
593
594	greg	2.36	/* read nc-component common-exponent color scan and convert to COLR's */
595			int
596			fread2colrs(COLR scanline, int len, FILE fp, int nc, const float wlpt[4])
597			{
598			COLRV *sclrscan;
599			int n;
600
601			if (nc < 3)
602			return(-1);
603			if (nc == 3)
604			return(freadcolrs(scanline, len, fp));
605
606			sclrscan = (COLRV )tempbuffer(sizeof(COLRV)(nc+1)*len);
607			if (sclrscan == NULL \|\| freadscolrs(sclrscan, nc, len, fp) < 0)
608			return(-1);
609			for (n = len; n--; ) {
610			scolr2colr(*scanline++, sclrscan, nc, wlpt);
611			sclrscan += nc+1;
612			}
613			return(0);
614			}
615
616
617	greg	2.31	/* write an common-exponent spectral color scanline */
618	greg	2.30	int
619	greg	2.35	fwritescolrs(const COLRV sscanline, int nc, int len, FILE fp)
620	greg	2.30	{
621	greg	2.32	if (nc < 3)
622			return(-1);
623			if (nc == 3)
624			return(fwritecolrs((COLR *)sscanline, len, fp));
625
626	greg	2.31	if (fwrite(sscanline, nc+1, len, fp) != len)
627	greg	2.30	return(-1);
628			return(0);
629			}
630
631
632	greg	2.9	int
633	greg	2.32	fwritescan( /* write out an RGB or XYZ scanline */
634	greg	2.17	COLOR *scanline,
635			int len,
636			FILE *fp
637			)
638	greg	1.1	{
639	greg	1.14	COLR *clrscan;
640			int n;
641	greg	2.17	COLR *sp;
642	greg	1.14	/* get scanline buffer */
643			if ((sp = (COLR )tempbuffer(lensizeof(COLR))) == NULL)
644			return(-1);
645			clrscan = sp;
646			/* convert scanline */
647			n = len;
648			while (n-- > 0) {
649			setcolr(sp[0], scanline[0][RED],
650	greg	1.1	scanline[0][GRN],
651			scanline[0][BLU]);
652			scanline++;
653	greg	1.14	sp++;
654	greg	1.1	}
655	greg	1.14	return(fwritecolrs(clrscan, len, fp));
656	greg	1.1	}
657
658
659	greg	2.9	int
660	greg	2.32	freadscan( /* read in an RGB or XYZ scanline */
661	greg	2.17	COLOR *scanline,
662			int len,
663			FILE *fp
664			)
665	greg	1.1	{
666	greg	2.17	COLR *clrscan;
667	greg	1.14
668			if ((clrscan = (COLR )tempbuffer(lensizeof(COLR))) == NULL)
669			return(-1);
670			if (freadcolrs(clrscan, len, fp) < 0)
671			return(-1);
672			/* convert scanline */
673			colr_color(scanline[0], clrscan[0]);
674			while (--len > 0) {
675			scanline++; clrscan++;
676	greg	2.20	if (clrscan[0][GRN] == clrscan[-1][GRN] &&
677			(clrscan[0][RED] == clrscan[-1][RED]) &
678			(clrscan[0][BLU] == clrscan[-1][BLU]) &
679			(clrscan[0][EXP] == clrscan[-1][EXP]))
680	greg	1.14	copycolor(scanline[0], scanline[-1]);
681			else
682			colr_color(scanline[0], clrscan[0]);
683	greg	1.1	}
684			return(0);
685			}
686
687
688	greg	2.30	/* read an nc-component color scanline */
689			int
690			freadsscan(COLORV sscanline, int nc, int len, FILE fp)
691			{
692	greg	2.35	COLRV tscn = (COLRV )tempbuffer((nc+1)*len);
693	greg	2.30	int i;
694
695			if (tscn == NULL \|\| freadscolrs(tscn, nc, len, fp) < 0)
696			return(-1);
697			for (i = len; i-- > 0; ) {
698			scolr2scolor(sscanline, tscn, nc);
699			sscanline += nc;
700			tscn += nc+1;
701			}
702			return(0);
703			}
704
705
706	greg	2.36	/* read an nc-component color scanline and return as RGB */
707			int
708			fread2scan(COLOR scanline, int len, FILE fp, int nc, const float wlpt[4])
709			{
710			COLRV *tscn;
711			int i;
712
713			if (nc < 3)
714			return(-1);
715			if (nc == 3)
716			return(freadscan(scanline, len, fp));
717
718			tscn = (COLRV )tempbuffer((nc+1)len);
719			if (tscn == NULL \|\| freadscolrs(tscn, nc, len, fp) < 0)
720			return(-1);
721			for (i = len; i-- > 0; ) {
722			scolr2color(*scanline++, tscn, nc, wlpt);
723			tscn += nc+1;
724			}
725			return(0);
726			}
727
728
729	greg	2.32	/* write an nc-component spectral color scanline */
730	greg	2.30	int
731	greg	2.35	fwritesscan(const COLORV sscanline, int nc, int len, FILE fp)
732	greg	2.30	{
733	greg	2.35	COLRV tscn = (COLRV )tempbuffer((nc+1)*len);
734	greg	2.30	int i;
735
736			if (tscn == NULL)
737			return(-1);
738			for (i = 0; i < len; i++) {
739	greg	2.31	scolor2scolr(tscn+i*(nc+1), sscanline, nc);
740			sscanline += nc;
741	greg	2.30	}
742	greg	2.31	return(fwritescolrs(tscn, nc, len, fp));
743	greg	2.30	}
744
745
746	greg	2.9	void
747	greg	2.17	setcolr( /* assign a short color value */
748			COLR clr,
749			double r,
750			double g,
751			double b
752			)
753	greg	1.1	{
754			double d;
755			int e;
756
757			d = r > g ? r : g;
758			if (b > d) d = b;
759
760	greg	1.4	if (d <= 1e-32) {
761	greg	1.1	clr[RED] = clr[GRN] = clr[BLU] = 0;
762			clr[EXP] = 0;
763			return;
764			}
765
766	greg	2.21	d = frexp(d, &e) * 256.0 / d;
767	greg	1.1
768	greg	2.27	clr[RED] = (r > 0) * (int)(r*d);
769			clr[GRN] = (g > 0) * (int)(g*d);
770			clr[BLU] = (b > 0) * (int)(b*d);
771	greg	1.1	clr[EXP] = e + COLXS;
772			}
773
774
775	greg	2.9	void
776	greg	2.17	colr_color( /* convert short to float color */
777			COLOR col,
778	greg	2.35	const COLR clr
779	greg	2.17	)
780	greg	1.1	{
781	greg	1.6	double f;
782	greg	1.1
783	greg	2.27	if (clr[EXP] == 0) {
784	greg	1.1	col[RED] = col[GRN] = col[BLU] = 0.0;
785	greg	2.27	return;
786	greg	1.1	}
787	greg	2.27	f = ldexp(1.0, (int)clr[EXP]-(COLXS+8));
788			col[RED] = (clr[RED] + 0.5)*f;
789			col[GRN] = (clr[GRN] + 0.5)*f;
790			col[BLU] = (clr[BLU] + 0.5)*f;
791	greg	1.6	}
792
793
794	greg	2.9	int
795	greg	2.17	bigdiff( /* c1 delta c2 > md? */
796	greg	2.35	const COLOR c1,
797			const COLOR c2,
798	greg	2.17	double md
799			)
800	greg	1.7	{
801	greg	2.17	int i;
802	greg	1.7
803			for (i = 0; i < 3; i++)
804	greg	2.27	if ((colval(c1,i)-colval(c2,i) > md*colval(c2,i)) \|
805			(colval(c2,i)-colval(c1,i) > md*colval(c1,i)))
806			return(1);
807			return(0);
808			}
809
810
811			int
812			sbigsdiff( /* sc1 delta sc2 > md? */
813	greg	2.35	const SCOLOR c1,
814			const SCOLOR c2,
815	greg	2.27	double md
816			)
817			{
818			int i = NCSAMP;
819
820			while (i--)
821			if ((c1[i]-c2[i] > mdc2[i]) \| (c2[i]-c1[i] > mdc1[i]))
822	greg	1.7	return(1);
823			return(0);
824			}