ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/color.c
(Generate patch)

Comparing ray/src/common/color.c (file contents):
Revision 2.20 by greg, Sat Feb 23 19:11:25 2019 UTC vs.
Revision 2.31 by greg, Mon Nov 27 21:00:14 2023 UTC

# Line 28 | Line 28 | static const char      RCSid[] = "$Id$";
28   #define  MINRUN         4       /* minimum run length */
29  
30  
31 < char *
31 > 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 >        const double    step = (WLPART[3] - WLPART[0])/(double)NCSAMP;
92 >        double          cwl = WLPART[0] + .5*step;
93 >        int             i;
94 >
95 >        for (i = 0; i < NCSAMP; i++) {
96 >                if (cwl >= WLPART[1])
97 >                        scol[i] = r;
98 >                else if (cwl >= WLPART[2])
99 >                        scol[i] = g;
100 >                else
101 >                        scol[i] = b;
102 >                cwl += step;
103 >        }
104 > }
105 >
106 >
107 > void
108 > scolor2color(                   /* assign RGB color from spectrum */
109 >        COLOR col,
110 >        SCOLOR scol,            /* uses average over bands */
111 >        int ncs,
112 >        const float wlpt[4]
113 > )
114 > {
115 >        const double    step = (wlpt[3] - wlpt[0])/(double)ncs;
116 >        double          cwl = wlpt[0] + .5*step;
117 >        int             i, j=0, n=0;
118 >
119 >        setcolor(col, 0, 0, 0);
120 >        for (i = 0; i < ncs; i++) {
121 >                if (cwl < wlpt[j+1]) {
122 >                        if (n > 1) col[j] /= (COLORV)n;
123 >                        j++;
124 >                        n = 0;
125 >                }
126 >                col[j] += scol[i];
127 >                n++;
128 >                cwl += step;
129 >        }
130 >        if (n > 1) col[j] /= (COLORV)n;
131 > }
132 >
133 >
134 > void
135 > scolor2colr(                    /* assign RGBE from spectral color */
136 >        COLR clr,
137 >        SCOLOR scol,            /* uses average over bands */
138 >        int ncs,
139 >        const float wlpt[4]
140 > )
141 > {
142 >        COLOR   col;
143 >
144 >        scolor2color(col, scol, ncs, wlpt);
145 >        setcolr(clr, col[RED], col[GRN], col[BLU]);
146 > }
147 >
148 >
149 > void
150 > scolor2scolr(                   /* float spectrum to common exponent */
151 >        SCOLR sclr,
152 >        SCOLOR scol,
153 >        int ncs
154 > )
155 > {
156 >        int     i = ncs;
157 >        COLORV  p = scol[--i];
158 >
159 >        while (i)
160 >                if (scol[--i] > p)
161 >                        p = scol[i];
162 >        if (p <= 1e-32) {
163 >                memset(sclr, 0, ncs+1);
164 >                return;
165 >        }
166 >        p = frexp(p, &i) * 256.0 / p;
167 >        sclr[ncs] = i + COLXS;
168 >        for (i = ncs; i--; )
169 >                sclr[i] = (scol[i] > 0) * (int)(scol[i]*p);
170 > }
171 >
172 >
173 > void
174 > scolr2scolor(                   /* common exponent to float spectrum */
175 >        SCOLOR scol,
176 >        SCOLR sclr,
177 >        int ncs
178 > )
179 > {
180 >        double  f;
181 >        int     i;
182 >
183 >        if (sclr[ncs] == 0) {
184 >                memset(scol, 0, sizeof(COLORV)*ncs);
185 >                return;
186 >        }
187 >        f = ldexp(1.0, (int)sclr[ncs]-(COLXS+8));
188 >
189 >        for (i = ncs; i--; )
190 >                scol[i] = (sclr[i] + 0.5)*f;
191 > }
192 >
193 >
194 > double
195 > scolor_mean(                    /* compute average for spectral color */
196 >        SCOLOR  scol
197 > )
198 > {
199 >        int     i = NCSAMP;
200 >        double  sum = 0;
201 >
202 >        while (i--)
203 >                sum += scol[i];
204 >
205 >        return sum/(double)NCSAMP;
206 > }
207 >
208 >
209 > double
210 > sintens(                        /* find maximum value from spectrum */
211 >        SCOLOR  scol
212 > )
213 > {
214 >        int     i = NCSAMP;
215 >        COLORV  peak = scol[--i];
216 >
217 >        while (i)
218 >                if (scol[--i] > peak)
219 >                        peak = scol[i];
220 >
221 >        return peak;
222 > }
223 >
224 >
225 > void
226 > convertscolor(                  /* spectrum conversion, zero-fill ends */
227 >        SCOLOR dst,             /* destination spectrum */
228 >        int dnc,                /* destination # of spectral samples/intervals */
229 >        double dwl0,            /* starting destination wavelength (longer) */
230 >        double dwl1,            /* ending destination wavelength (shorter) */
231 >        const COLORV src[],     /* source spectrum array */
232 >        int snc,
233 >        double swl0,            /* long/short wavelengths may be reversed */
234 >        double swl1
235 > )
236 > {
237 >        const int       sdir = 1 - 2*(swl0 < swl1);
238 >        const double    sstp = (swl1 - swl0)/(double)snc;
239 >        const double    dstp = (dwl1 - dwl0)/(double)dnc;
240 >        const double    rdstp = 1./dstp;
241 >        int             si, ssi, di;
242 >        double          wl;
243 >
244 >        if ((dnc < 3) | (dwl0 <= dwl1) | (dst == src))
245 >                return;         /* invalid destination */
246 >
247 >        if (dnc == snc && (dwl0-swl0)*(dwl0-swl0) + (dwl1-swl1)*(dwl1-swl1) <= .5) {
248 >                memcpy(dst, src, sizeof(COLORV)*dnc);
249 >                return;         /* same spectral sampling */
250 >        }
251 >        memset(dst, 0, sizeof(COLORV)*dnc);
252 >                                /* set starting positions */
253 >        if ((sdir>0 ? swl0 : swl1) <= dwl0) {
254 >                if (sdir > 0) {
255 >                        wl = swl0;
256 >                        ssi = 0;
257 >                } else {
258 >                        wl = swl1;
259 >                        ssi = snc-1;
260 >                }
261 >                si = 0;
262 >                di = (wl - dwl0)*rdstp;
263 >        } else {
264 >                wl = dwl0;
265 >                if (sdir > 0) {
266 >                        ssi = si = (wl - swl0)/sstp;
267 >                } else {
268 >                        si = (wl - swl1)/sstp;
269 >                        ssi = snc-1 - si;
270 >                }
271 >                di = 0;
272 >        }
273 >        swl0 += (sdir < 0)*sstp;
274 >                                /* step through intervals */
275 >        while ((si < snc) & (di < dnc)) {
276 >                double  intvl;
277 >                if (swl0 + (ssi+sdir)*sstp < dwl0 + (di+1)*dstp) {
278 >                        intvl = dwl0 + (di+1)*dstp - wl;
279 >                        dst[di++] += src[ssi]*intvl*rdstp;
280 >                } else {
281 >                        intvl = swl0 + (ssi+sdir)*sstp - wl;
282 >                        dst[di] += src[ssi]*intvl*rdstp;
283 >                        ssi += sdir;
284 >                        si++;
285 >                }
286 >                wl += intvl;
287 >        }
288 > }
289 >
290 >
291 > void *
292   tempbuffer(                     /* get a temporary buffer */
293 <        unsigned int  len
293 >        size_t  len
294   )
295   {
296 <        static char  *tempbuf = NULL;
297 <        static unsigned  tempbuflen = 0;
296 >        static void     *tempbuf = NULL;
297 >        static size_t   tempbuflen = 0;
298  
299 <        if (!len | (len > tempbuflen)) {
300 <                if (tempbuflen)
299 >        if (!len) {             /* call to free */
300 >                if (tempbuflen) {
301                          free(tempbuf);
302 <                tempbuf = len ? (char *)malloc(len) : (char *)NULL;
303 <                tempbuflen = len*(tempbuf != NULL);
302 >                        tempbuf = NULL;
303 >                        tempbuflen = 0;
304 >                }
305 >                return(NULL);
306          }
307 +        if (len <= tempbuflen)  /* big enough already? */
308 +                return(tempbuf);
309 +                                /* else free & reallocate */
310 +        if (tempbuflen)
311 +                free(tempbuf);
312 +        tempbuf = malloc(len);
313 +        tempbuflen = len*(tempbuf != NULL);
314          return(tempbuf);
315   }
316  
# Line 62 | Line 331 | fwritecolrs(                   /* write out a colr scanline */
331          putc(2, fp);
332          putc(2, fp);
333          putc(len>>8, fp);
334 <        putc(len&255, fp);
334 >        putc(len&0xff, fp);
335                                          /* put components seperately */
336          for (i = 0; i < 4; i++) {
337              for (j = 0; j < len; j += cnt) {    /* find next run */
# Line 99 | Line 368 | fwritecolrs(                   /* write out a colr scanline */
368          return(ferror(fp) ? -1 : 0);
369   }
370  
371 <
371 > /*
372 > * An old-format scanline is either a stream of valid RGBE or XYZE real
373 > * pixels or at least one real pixel followed by some number of
374 > * invalid real pixels of the form (1,1,1,n), where n is a count.
375 > * These can themselves be repeated to create a multibyte repeat
376 > * count, with the least significant byte first (little-endian order.)
377 > * Repeat counts are limited by the size of an int; if a repetition
378 > * leads to an overrun, the rest of the the repetition will be
379 > * silently ignored.
380 > */
381   static int
382   oldreadcolrs(                   /* read in an old-style colr scanline */
383          COLR  *scanline,
# Line 137 | Line 415 | oldreadcolrs(                  /* read in an old-style colr scanline
415          return(0);
416   }
417  
418 <
418 > /*
419 > * There are two scanline formats: old and new.  The old format
420 > * compresses runs of RGBE or XYZE four-byte real pixels; the new
421 > * format breaks the pixels into R, G, B, and E lines (or XYZE lines)
422 > * which are individually run-length encoded.
423 > *
424 > * An old-format scanline always begins with a valid real pixel; at
425 > * least one of the RGB (or XYZ) values will have its high-order bit
426 > * set.  A new-format scanline begins with four bytes which are not a
427 > * valid real pixel: (2, 2, lenhigh, lenlow) where lenhigh is always
428 > * less than 128 and hence never has a high-order bit set.
429 > *
430 > * A new-format scanline is broken into its RGBE or XYZE components.
431 > * Each is output and run-length encoded separately so that a scanline
432 > * is broken into four records.  In turn, each record is organized
433 > * into chunks of up to 128 characters, which begin with a count byte.
434 > * If the count byte is greater than 128, the following data byte is
435 > * repeated (count-128) times.  If not, the count byte is followed by
436 > * that many data bytes.
437 > */
438   int
439   freadcolrs(                     /* read in an encoded colr scanline */
440          COLR  *scanline,
# Line 148 | Line 445 | freadcolrs(                    /* read in an encoded colr scanline */
445          int  i, j;
446          int  code, val;
447                                          /* determine scanline type */
448 <        if ((len < MINELEN) | (len > MAXELEN))
449 <                return(oldreadcolrs(scanline, len, fp));
448 >        if (len <= 0)
449 >                return(0);
450          if ((i = getc(fp)) == EOF)
451                  return(-1);
452 <        if (i != 2) {
156 <                ungetc(i, fp);
157 <                return(oldreadcolrs(scanline, len, fp));
158 <        }
452 >        scanline[0][RED] = i;
453          scanline[0][GRN] = getc(fp);
454          scanline[0][BLU] = getc(fp);
455          if ((i = getc(fp)) == EOF)
456                  return(-1);
457 <        if (scanline[0][GRN] != 2 || scanline[0][BLU] & 128) {
458 <                scanline[0][RED] = 2;
457 >        if ((scanline[0][RED] != 2) | (scanline[0][GRN] != 2) |
458 >                        (scanline[0][BLU] & 0x80)) {
459                  scanline[0][EXP] = i;
460                  return(oldreadcolrs(scanline+1, len-1, fp));
461          }
# Line 194 | Line 488 | freadcolrs(                    /* read in an encoded colr scanline */
488   }
489  
490  
491 + /* read an nc-component common-exponent color scanline */
492   int
493 + freadscolrs(uby8 *scanline, int nc, int len, FILE *fp)
494 + {
495 +        if (fread(scanline, nc+1, len, fp) != len)
496 +                return(-1);
497 +        return(0);
498 + }
499 +
500 +
501 + /* write an common-exponent spectral color scanline */
502 + int
503 + fwritescolrs(uby8 *sscanline, int nc, int len, FILE *fp)
504 + {
505 +        if (fwrite(sscanline, nc+1, len, fp) != len)
506 +                return(-1);
507 +        return(0);
508 + }
509 +
510 +
511 + int
512   fwritescan(                     /* write out a scanline */
513          COLOR  *scanline,
514          int  len,
# Line 250 | Line 564 | freadscan(                     /* read in a scanline */
564   }
565  
566  
567 + /* read an nc-component color scanline */
568 + int
569 + freadsscan(COLORV *sscanline, int nc, int len, FILE *fp)
570 + {
571 +        uby8    *tscn = (uby8 *)tempbuffer((nc+1)*len);
572 +        int     i;
573 +
574 +        if (tscn == NULL || freadscolrs(tscn, nc, len, fp) < 0)
575 +                return(-1);
576 +        for (i = len; i-- > 0; ) {
577 +                scolr2scolor(sscanline, tscn, nc);
578 +                sscanline += nc;
579 +                tscn += nc+1;
580 +        }
581 +        return(0);
582 + }
583 +
584 +
585 + /* write an spectral color scanline (NCSAMP) */
586 + int
587 + fwritesscan(COLORV *sscanline, int nc, int len, FILE *fp)
588 + {
589 +        uby8    *tscn = (uby8 *)tempbuffer((nc+1)*len);
590 +        int     i;
591 +
592 +        if (tscn == NULL)
593 +                return(-1);
594 +        for (i = 0; i < len; i++) {
595 +                scolor2scolr(tscn+i*(nc+1), sscanline, nc);
596 +                sscanline += nc;
597 +        }
598 +        return(fwritescolrs(tscn, nc, len, fp));
599 + }
600 +
601 +
602   void
603   setcolr(                        /* assign a short color value */
604          COLR  clr,
# Line 270 | Line 619 | setcolr(                       /* assign a short color value */
619                  return;
620          }
621  
622 <        d = frexp(d, &e) * 255.9999 / d;
622 >        d = frexp(d, &e) * 256.0 / d;
623  
624 <        if (r > 0.0)
625 <                clr[RED] = r * d;
626 <        else
278 <                clr[RED] = 0;
279 <        if (g > 0.0)
280 <                clr[GRN] = g * d;
281 <        else
282 <                clr[GRN] = 0;
283 <        if (b > 0.0)
284 <                clr[BLU] = b * d;
285 <        else
286 <                clr[BLU] = 0;
287 <
624 >        clr[RED] = (r > 0) * (int)(r*d);
625 >        clr[GRN] = (g > 0) * (int)(g*d);
626 >        clr[BLU] = (b > 0) * (int)(b*d);
627          clr[EXP] = e + COLXS;
628   }
629  
# Line 297 | Line 636 | colr_color(                    /* convert short to float color */
636   {
637          double  f;
638          
639 <        if (clr[EXP] == 0)
639 >        if (clr[EXP] == 0) {
640                  col[RED] = col[GRN] = col[BLU] = 0.0;
641 <        else {
303 <                f = ldexp(1.0, (int)clr[EXP]-(COLXS+8));
304 <                col[RED] = (clr[RED] + 0.5)*f;
305 <                col[GRN] = (clr[GRN] + 0.5)*f;
306 <                col[BLU] = (clr[BLU] + 0.5)*f;
641 >                return;
642          }
643 +        f = ldexp(1.0, (int)clr[EXP]-(COLXS+8));
644 +        col[RED] = (clr[RED] + 0.5)*f;
645 +        col[GRN] = (clr[GRN] + 0.5)*f;
646 +        col[BLU] = (clr[BLU] + 0.5)*f;
647   }
648  
649  
# Line 318 | Line 657 | bigdiff(                               /* c1 delta c2 > md? */
657          int  i;
658  
659          for (i = 0; i < 3; i++)
660 <                if (colval(c1,i)-colval(c2,i) > md*colval(c2,i) ||
661 <                                colval(c2,i)-colval(c1,i) > md*colval(c1,i))
660 >                if ((colval(c1,i)-colval(c2,i) > md*colval(c2,i)) |
661 >                                (colval(c2,i)-colval(c1,i) > md*colval(c1,i)))
662 >                        return(1);
663 >        return(0);
664 > }
665 >
666 >
667 > int
668 > sbigsdiff(                              /* sc1 delta sc2 > md? */
669 >        SCOLOR  c1,
670 >        SCOLOR  c2,
671 >        double  md
672 > )
673 > {
674 >        int  i = NCSAMP;
675 >
676 >        while (i--)
677 >                if ((c1[i]-c2[i] > md*c2[i]) | (c2[i]-c1[i] > md*c1[i]))
678                          return(1);
679          return(0);
680   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines