ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/color.c
Revision: 2.33
Committed: Wed Jan 17 16:06:56 2024 UTC (3 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.32: +16 -1 lines
Log Message:
feat: added missing call for scolr2color()

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: color.c,v 2.32 2023/12/07 23:16:58 greg Exp $";
3 #endif
4 /*
5 * color.c - routines for color calculations.
6 *
7 * Externals declared in color.h
8 */
9
10 #include "copyright.h"
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <math.h>
15 #include "color.h"
16
17 #ifdef getc_unlocked /* avoid horrendous overhead of flockfile */
18 #undef getc
19 #undef putc
20 #undef ferror
21 #define getc getc_unlocked
22 #define putc putc_unlocked
23 #define ferror ferror_unlocked
24 #endif
25
26 #define MINELEN 8 /* minimum scanline length for encoding */
27 #define MAXELEN 0x7fff /* maximum scanline length for encoding */
28 #define MINRUN 4 /* minimum run length */
29
30
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 scolr2color( /* assign RGB from common exponent */
151 COLOR col,
152 SCOLR sclr,
153 int ncs,
154 const float wlpt[4]
155 )
156 {
157 SCOLOR scol;
158
159 scolr2scolor(scol, sclr, ncs);
160 scolor2color(col, scol, ncs, wlpt);
161 }
162
163
164 void
165 scolor2scolr( /* float spectrum to common exponent */
166 SCOLR sclr,
167 SCOLOR scol,
168 int ncs
169 )
170 {
171 int i = ncs;
172 COLORV p = scol[--i];
173
174 while (i)
175 if (scol[--i] > p)
176 p = scol[i];
177 if (p <= 1e-32) {
178 memset(sclr, 0, ncs+1);
179 return;
180 }
181 p = frexp(p, &i) * 256.0 / p;
182 sclr[ncs] = i + COLXS;
183 for (i = ncs; i--; )
184 sclr[i] = (scol[i] > 0) * (int)(scol[i]*p);
185 }
186
187
188 void
189 scolr2scolor( /* common exponent to float spectrum */
190 SCOLOR scol,
191 SCOLR sclr,
192 int ncs
193 )
194 {
195 double f;
196 int i;
197
198 if (sclr[ncs] == 0) {
199 memset(scol, 0, sizeof(COLORV)*ncs);
200 return;
201 }
202 f = ldexp(1.0, (int)sclr[ncs]-(COLXS+8));
203
204 for (i = ncs; i--; )
205 scol[i] = (sclr[i] + 0.5)*f;
206 }
207
208
209 double
210 scolor_mean( /* compute average for spectral color */
211 SCOLOR scol
212 )
213 {
214 int i = NCSAMP;
215 double sum = 0;
216
217 while (i--)
218 sum += scol[i];
219
220 return sum/(double)NCSAMP;
221 }
222
223
224 double
225 sintens( /* find maximum value from spectrum */
226 SCOLOR scol
227 )
228 {
229 int i = NCSAMP;
230 COLORV peak = scol[--i];
231
232 while (i)
233 if (scol[--i] > peak)
234 peak = scol[i];
235
236 return peak;
237 }
238
239
240 void
241 convertscolor( /* spectrum conversion, zero-fill ends */
242 SCOLOR dst, /* destination spectrum */
243 int dnc, /* destination # of spectral samples/intervals */
244 double dwl0, /* starting destination wavelength (longer) */
245 double dwl1, /* ending destination wavelength (shorter) */
246 const COLORV src[], /* source spectrum array */
247 int snc,
248 double swl0, /* long/short wavelengths may be reversed */
249 double swl1
250 )
251 {
252 const int sdir = 1 - 2*(swl0 < swl1);
253 const double sstp = (swl1 - swl0)/(double)snc;
254 const double dstp = (dwl1 - dwl0)/(double)dnc;
255 const double rdstp = 1./dstp;
256 int si, ssi, di;
257 double wl;
258
259 if ((dnc < 3) | (dwl0 <= dwl1) | (dst == src))
260 return; /* invalid destination */
261
262 if (dnc == snc && (dwl0-swl0)*(dwl0-swl0) + (dwl1-swl1)*(dwl1-swl1) <= .5) {
263 memcpy(dst, src, sizeof(COLORV)*dnc);
264 return; /* same spectral sampling */
265 }
266 memset(dst, 0, sizeof(COLORV)*dnc);
267 /* set starting positions */
268 if ((sdir>0 ? swl0 : swl1) <= dwl0) {
269 if (sdir > 0) {
270 wl = swl0;
271 ssi = 0;
272 } else {
273 wl = swl1;
274 ssi = snc-1;
275 }
276 si = 0;
277 di = (wl - dwl0)*rdstp;
278 } else {
279 wl = dwl0;
280 if (sdir > 0) {
281 ssi = si = (wl - swl0)/sstp;
282 } else {
283 si = (wl - swl1)/sstp;
284 ssi = snc-1 - si;
285 }
286 di = 0;
287 }
288 swl0 += (sdir < 0)*sstp;
289 /* step through intervals */
290 while ((si < snc) & (di < dnc)) {
291 double intvl;
292 if (swl0 + (ssi+sdir)*sstp < dwl0 + (di+1)*dstp) {
293 intvl = dwl0 + (di+1)*dstp - wl;
294 dst[di++] += src[ssi]*intvl*rdstp;
295 } else {
296 intvl = swl0 + (ssi+sdir)*sstp - wl;
297 dst[di] += src[ssi]*intvl*rdstp;
298 ssi += sdir;
299 si++;
300 }
301 wl += intvl;
302 }
303 }
304
305
306 void *
307 tempbuffer( /* get a temporary buffer */
308 size_t len
309 )
310 {
311 static void *tempbuf = NULL;
312 static size_t tempbuflen = 0;
313
314 if (!len) { /* call to free */
315 if (tempbuflen) {
316 free(tempbuf);
317 tempbuf = NULL;
318 tempbuflen = 0;
319 }
320 return(NULL);
321 }
322 if (len <= tempbuflen) /* big enough already? */
323 return(tempbuf);
324 /* else free & reallocate */
325 if (tempbuflen)
326 free(tempbuf);
327 tempbuf = malloc(len);
328 tempbuflen = len*(tempbuf != NULL);
329 return(tempbuf);
330 }
331
332
333 int
334 fwritecolrs( /* write out a colr scanline */
335 COLR *scanline,
336 int len,
337 FILE *fp
338 )
339 {
340 int i, j, beg, cnt = 1;
341 int c2;
342
343 if ((len < MINELEN) | (len > MAXELEN)) /* OOBs, write out flat */
344 return(fwrite((char *)scanline,sizeof(COLR),len,fp) - len);
345 /* put magic header */
346 putc(2, fp);
347 putc(2, fp);
348 putc(len>>8, fp);
349 putc(len&0xff, fp);
350 /* put components seperately */
351 for (i = 0; i < 4; i++) {
352 for (j = 0; j < len; j += cnt) { /* find next run */
353 for (beg = j; beg < len; beg += cnt) {
354 for (cnt = 1; (cnt < 127) & (beg+cnt < len) &&
355 scanline[beg+cnt][i] == scanline[beg][i]; cnt++)
356 ;
357 if (cnt >= MINRUN)
358 break; /* long enough */
359 }
360 if ((beg-j > 1) & (beg-j < MINRUN)) {
361 c2 = j+1;
362 while (scanline[c2++][i] == scanline[j][i])
363 if (c2 == beg) { /* short run */
364 putc(128+beg-j, fp);
365 putc(scanline[j][i], fp);
366 j = beg;
367 break;
368 }
369 }
370 while (j < beg) { /* write out non-run */
371 if ((c2 = beg-j) > 128) c2 = 128;
372 putc(c2, fp);
373 while (c2--)
374 putc(scanline[j++][i], fp);
375 }
376 if (cnt >= MINRUN) { /* write out run */
377 putc(128+cnt, fp);
378 putc(scanline[beg][i], fp);
379 } else
380 cnt = 0;
381 }
382 }
383 return(ferror(fp) ? -1 : 0);
384 }
385
386 /*
387 * An old-format scanline is either a stream of valid RGBE or XYZE real
388 * pixels or at least one real pixel followed by some number of
389 * invalid real pixels of the form (1,1,1,n), where n is a count.
390 * These can themselves be repeated to create a multibyte repeat
391 * count, with the least significant byte first (little-endian order.)
392 * Repeat counts are limited by the size of an int; if a repetition
393 * leads to an overrun, the rest of the the repetition will be
394 * silently ignored.
395 */
396 static int
397 oldreadcolrs( /* read in an old-style colr scanline */
398 COLR *scanline,
399 int len,
400 FILE *fp
401 )
402 {
403 int rshift = 0;
404 int i;
405
406 while (len > 0) {
407 scanline[0][RED] = getc(fp);
408 scanline[0][GRN] = getc(fp);
409 scanline[0][BLU] = getc(fp);
410 scanline[0][EXP] = i = getc(fp);
411 if (i == EOF)
412 return(-1);
413 if (scanline[0][GRN] == 1 &&
414 (scanline[0][RED] == 1) &
415 (scanline[0][BLU] == 1)) {
416 i = scanline[0][EXP] << rshift;
417 while (i--) {
418 copycolr(scanline[0], scanline[-1]);
419 if (--len <= 0)
420 return(0);
421 scanline++;
422 }
423 rshift += 8;
424 } else {
425 scanline++;
426 len--;
427 rshift = 0;
428 }
429 }
430 return(0);
431 }
432
433 /*
434 * There are two scanline formats: old and new. The old format
435 * compresses runs of RGBE or XYZE four-byte real pixels; the new
436 * format breaks the pixels into R, G, B, and E lines (or XYZE lines)
437 * which are individually run-length encoded.
438 *
439 * An old-format scanline always begins with a valid real pixel; at
440 * least one of the RGB (or XYZ) values will have its high-order bit
441 * set. A new-format scanline begins with four bytes which are not a
442 * valid real pixel: (2, 2, lenhigh, lenlow) where lenhigh is always
443 * less than 128 and hence never has a high-order bit set.
444 *
445 * A new-format scanline is broken into its RGBE or XYZE components.
446 * Each is output and run-length encoded separately so that a scanline
447 * is broken into four records. In turn, each record is organized
448 * into chunks of up to 128 characters, which begin with a count byte.
449 * If the count byte is greater than 128, the following data byte is
450 * repeated (count-128) times. If not, the count byte is followed by
451 * that many data bytes.
452 */
453 int
454 freadcolrs( /* read in an encoded colr scanline */
455 COLR *scanline,
456 int len,
457 FILE *fp
458 )
459 {
460 int i, j;
461 int code, val;
462 /* determine scanline type */
463 if (len <= 0)
464 return(0);
465 if ((i = getc(fp)) == EOF)
466 return(-1);
467 scanline[0][RED] = i;
468 scanline[0][GRN] = getc(fp);
469 scanline[0][BLU] = getc(fp);
470 if ((i = getc(fp)) == EOF)
471 return(-1);
472 if ((scanline[0][RED] != 2) | (scanline[0][GRN] != 2) |
473 (scanline[0][BLU] & 0x80)) {
474 scanline[0][EXP] = i;
475 return(oldreadcolrs(scanline+1, len-1, fp));
476 }
477 if ((scanline[0][BLU]<<8 | i) != len)
478 return(-1); /* length mismatch! */
479 /* read each component */
480 for (i = 0; i < 4; i++)
481 for (j = 0; j < len; ) {
482 if ((code = getc(fp)) == EOF)
483 return(-1);
484 if (code > 128) { /* run */
485 code &= 127;
486 if ((val = getc(fp)) == EOF)
487 return -1;
488 if (j + code > len)
489 return -1; /* overrun */
490 while (code--)
491 scanline[j++][i] = val;
492 } else { /* non-run */
493 if (j + code > len)
494 return -1; /* overrun */
495 while (code--) {
496 if ((val = getc(fp)) == EOF)
497 return -1;
498 scanline[j++][i] = val;
499 }
500 }
501 }
502 return(0);
503 }
504
505
506 /* read an nc-component common-exponent color scanline */
507 int
508 freadscolrs(uby8 *scanline, int nc, int len, FILE *fp)
509 {
510 if (nc < 3)
511 return(-1);
512 if (nc == 3)
513 return(freadcolrs((COLR *)scanline, len, fp));
514
515 if (fread(scanline, nc+1, len, fp) != len)
516 return(-1);
517 return(0);
518 }
519
520
521 /* write an common-exponent spectral color scanline */
522 int
523 fwritescolrs(uby8 *sscanline, int nc, int len, FILE *fp)
524 {
525 if (nc < 3)
526 return(-1);
527 if (nc == 3)
528 return(fwritecolrs((COLR *)sscanline, len, fp));
529
530 if (fwrite(sscanline, nc+1, len, fp) != len)
531 return(-1);
532 return(0);
533 }
534
535
536 int
537 fwritescan( /* write out an RGB or XYZ scanline */
538 COLOR *scanline,
539 int len,
540 FILE *fp
541 )
542 {
543 COLR *clrscan;
544 int n;
545 COLR *sp;
546 /* get scanline buffer */
547 if ((sp = (COLR *)tempbuffer(len*sizeof(COLR))) == NULL)
548 return(-1);
549 clrscan = sp;
550 /* convert scanline */
551 n = len;
552 while (n-- > 0) {
553 setcolr(sp[0], scanline[0][RED],
554 scanline[0][GRN],
555 scanline[0][BLU]);
556 scanline++;
557 sp++;
558 }
559 return(fwritecolrs(clrscan, len, fp));
560 }
561
562
563 int
564 freadscan( /* read in an RGB or XYZ scanline */
565 COLOR *scanline,
566 int len,
567 FILE *fp
568 )
569 {
570 COLR *clrscan;
571
572 if ((clrscan = (COLR *)tempbuffer(len*sizeof(COLR))) == NULL)
573 return(-1);
574 if (freadcolrs(clrscan, len, fp) < 0)
575 return(-1);
576 /* convert scanline */
577 colr_color(scanline[0], clrscan[0]);
578 while (--len > 0) {
579 scanline++; clrscan++;
580 if (clrscan[0][GRN] == clrscan[-1][GRN] &&
581 (clrscan[0][RED] == clrscan[-1][RED]) &
582 (clrscan[0][BLU] == clrscan[-1][BLU]) &
583 (clrscan[0][EXP] == clrscan[-1][EXP]))
584 copycolor(scanline[0], scanline[-1]);
585 else
586 colr_color(scanline[0], clrscan[0]);
587 }
588 return(0);
589 }
590
591
592 /* read an nc-component color scanline */
593 int
594 freadsscan(COLORV *sscanline, int nc, int len, FILE *fp)
595 {
596 uby8 *tscn = (uby8 *)tempbuffer((nc+1)*len);
597 int i;
598
599 if (tscn == NULL || freadscolrs(tscn, nc, len, fp) < 0)
600 return(-1);
601 for (i = len; i-- > 0; ) {
602 scolr2scolor(sscanline, tscn, nc);
603 sscanline += nc;
604 tscn += nc+1;
605 }
606 return(0);
607 }
608
609
610 /* write an nc-component spectral color scanline */
611 int
612 fwritesscan(COLORV *sscanline, int nc, int len, FILE *fp)
613 {
614 uby8 *tscn = (uby8 *)tempbuffer((nc+1)*len);
615 int i;
616
617 if (tscn == NULL)
618 return(-1);
619 for (i = 0; i < len; i++) {
620 scolor2scolr(tscn+i*(nc+1), sscanline, nc);
621 sscanline += nc;
622 }
623 return(fwritescolrs(tscn, nc, len, fp));
624 }
625
626
627 void
628 setcolr( /* assign a short color value */
629 COLR clr,
630 double r,
631 double g,
632 double b
633 )
634 {
635 double d;
636 int e;
637
638 d = r > g ? r : g;
639 if (b > d) d = b;
640
641 if (d <= 1e-32) {
642 clr[RED] = clr[GRN] = clr[BLU] = 0;
643 clr[EXP] = 0;
644 return;
645 }
646
647 d = frexp(d, &e) * 256.0 / d;
648
649 clr[RED] = (r > 0) * (int)(r*d);
650 clr[GRN] = (g > 0) * (int)(g*d);
651 clr[BLU] = (b > 0) * (int)(b*d);
652 clr[EXP] = e + COLXS;
653 }
654
655
656 void
657 colr_color( /* convert short to float color */
658 COLOR col,
659 COLR clr
660 )
661 {
662 double f;
663
664 if (clr[EXP] == 0) {
665 col[RED] = col[GRN] = col[BLU] = 0.0;
666 return;
667 }
668 f = ldexp(1.0, (int)clr[EXP]-(COLXS+8));
669 col[RED] = (clr[RED] + 0.5)*f;
670 col[GRN] = (clr[GRN] + 0.5)*f;
671 col[BLU] = (clr[BLU] + 0.5)*f;
672 }
673
674
675 int
676 bigdiff( /* c1 delta c2 > md? */
677 COLOR c1,
678 COLOR c2,
679 double md
680 )
681 {
682 int i;
683
684 for (i = 0; i < 3; i++)
685 if ((colval(c1,i)-colval(c2,i) > md*colval(c2,i)) |
686 (colval(c2,i)-colval(c1,i) > md*colval(c1,i)))
687 return(1);
688 return(0);
689 }
690
691
692 int
693 sbigsdiff( /* sc1 delta sc2 > md? */
694 SCOLOR c1,
695 SCOLOR c2,
696 double md
697 )
698 {
699 int i = NCSAMP;
700
701 while (i--)
702 if ((c1[i]-c2[i] > md*c2[i]) | (c2[i]-c1[i] > md*c1[i]))
703 return(1);
704 return(0);
705 }