ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.8
Committed: Mon Aug 17 17:58:47 1998 UTC (25 years, 8 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.7: +8 -6 lines
Log Message:
took "color.h" include out of tonemap.h and put in tmprivat.h
changed routines so zero-length arrays don't return error

File Contents

# Content
1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2
3 #ifndef lint
4 static char SCCSid[] = "$SunId$ SGI";
5 #endif
6
7 /*
8 * Tone mapping functions.
9 * See tonemap.h for detailed function descriptions.
10 */
11
12 #include <stdio.h>
13 #include <math.h>
14 #include "tmprivat.h"
15 #include "tmerrmsg.h"
16
17 #define exp10(x) exp(M_LN10*(x))
18
19 struct tmStruct *tmTop = NULL; /* current tone mapping stack */
20
21 /* our list of conversion packages */
22 struct tmPackage *tmPkg[TM_MAXPKG];
23 int tmNumPkgs = 0; /* number of registered packages */
24
25 int tmLastError; /* last error incurred by library */
26 char *tmLastFunction; /* error-generating function name */
27
28
29 int
30 tmErrorReturn(func, err) /* error return (with message) */
31 char *func;
32 int err;
33 {
34 tmLastFunction = func;
35 tmLastError = err;
36 if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
37 return(err);
38 fputs(func, stderr);
39 fputs(": ", stderr);
40 fputs(tmErrorMessage[err], stderr);
41 fputs("!\n", stderr);
42 return(err);
43 }
44
45
46 struct tmStruct *
47 tmInit(flags, monpri, gamval) /* initialize new tone mapping */
48 int flags;
49 RGBPRIMP monpri;
50 double gamval;
51 {
52 static char funcName[] = "tmInit";
53 COLORMAT cmat;
54 register struct tmStruct *tmnew;
55 register int i;
56 /* allocate structure */
57 tmnew = (struct tmStruct *)malloc(sizeof(struct tmStruct));
58 if (tmnew == NULL)
59 return(NULL);
60
61 tmnew->flags = flags & ~TM_F_UNIMPL;
62 /* set monitor transform */
63 if (monpri == NULL || monpri == stdprims || tmnew->flags & TM_F_BW) {
64 tmnew->monpri = stdprims;
65 tmnew->clf[RED] = rgb2xyzmat[1][0];
66 tmnew->clf[GRN] = rgb2xyzmat[1][1];
67 tmnew->clf[BLU] = rgb2xyzmat[1][2];
68 } else {
69 comprgb2xyzmat(cmat, tmnew->monpri=monpri);
70 tmnew->clf[RED] = cmat[1][0];
71 tmnew->clf[GRN] = cmat[1][1];
72 tmnew->clf[BLU] = cmat[1][2];
73 }
74 /* set gamma value */
75 if (gamval < MINGAM)
76 tmnew->mongam = DEFGAM;
77 else
78 tmnew->mongam = gamval;
79 /* set color divisors */
80 for (i = 0; i < 3; i++)
81 tmnew->cdiv[i] = 256.*pow(tmnew->clf[i], 1./tmnew->mongam);
82
83 /* set input transform */
84 tmnew->inppri = tmnew->monpri;
85 tmnew->cmat[0][0] = tmnew->cmat[1][1] = tmnew->cmat[2][2] =
86 tmnew->inpsf = WHTEFFICACY;
87 tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
88 tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
89 tmnew->hbrmin = tmnew->hbrmax = 0;
90 tmnew->histo = NULL;
91 tmnew->mbrmin = tmnew->mbrmax = 0;
92 tmnew->lumap = NULL;
93 /* zero private data */
94 for (i = TM_MAXPKG; i--; )
95 tmnew->pd[i] = NULL;
96 /* make tmnew current */
97 tmnew->tmprev = tmTop;
98 return(tmTop = tmnew);
99 }
100
101
102 int
103 tmSetSpace(pri, sf) /* set input color space for conversions */
104 RGBPRIMP pri;
105 double sf;
106 {
107 static char funcName[] = "tmSetSpace";
108 register int i, j;
109 /* error check */
110 if (tmTop == NULL)
111 returnErr(TM_E_TMINVAL);
112 if (sf <= 1e-12)
113 returnErr(TM_E_ILLEGAL);
114 /* check if no change */
115 if (pri == tmTop->inppri && FEQ(sf, tmTop->inpsf))
116 returnOK;
117 tmTop->inppri = pri; /* let's set it */
118 tmTop->inpsf = sf;
119
120 if (tmTop->flags & TM_F_BW) { /* color doesn't matter */
121 tmTop->monpri = tmTop->inppri; /* eliminate xform */
122 if (tmTop->inppri == TM_XYZPRIM) {
123 tmTop->clf[CIEX] = tmTop->clf[CIEZ] = 0.;
124 tmTop->clf[CIEY] = 1.;
125 } else {
126 comprgb2xyzmat(tmTop->cmat, tmTop->monpri);
127 tmTop->clf[RED] = tmTop->cmat[1][0];
128 tmTop->clf[GRN] = tmTop->cmat[1][1];
129 tmTop->clf[BLU] = tmTop->cmat[1][2];
130 }
131 tmTop->cmat[0][0] = tmTop->cmat[1][1] = tmTop->cmat[2][2] =
132 tmTop->inpsf;
133 tmTop->cmat[0][1] = tmTop->cmat[0][2] = tmTop->cmat[1][0] =
134 tmTop->cmat[1][2] = tmTop->cmat[2][0] = tmTop->cmat[2][1] = 0.;
135
136 } else if (tmTop->inppri == TM_XYZPRIM) /* input is XYZ */
137 compxyz2rgbmat(tmTop->cmat, tmTop->monpri);
138
139 else { /* input is RGB */
140 if (tmTop->inppri != tmTop->monpri &&
141 PRIMEQ(tmTop->inppri, tmTop->monpri))
142 tmTop->inppri = tmTop->monpri; /* no xform */
143 comprgb2rgbmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
144 }
145 for (i = 0; i < 3; i++)
146 for (j = 0; j < 3; j++)
147 tmTop->cmat[i][j] *= tmTop->inpsf;
148 /* set color divisors */
149 for (i = 0; i < 3; i++)
150 if (tmTop->clf[i] > .001)
151 tmTop->cdiv[i] =
152 256.*pow(tmTop->clf[i], 1./tmTop->mongam);
153 else
154 tmTop->cdiv[i] = 1;
155 /* notify packages */
156 for (i = tmNumPkgs; i--; )
157 if (tmTop->pd[i] != NULL && tmPkg[i]->NewSpace != NULL)
158 (*tmPkg[i]->NewSpace)(tmTop);
159 returnOK;
160 }
161
162
163 void
164 tmClearHisto() /* clear current histogram */
165 {
166 if (tmTop == NULL || tmTop->histo == NULL)
167 return;
168 free((MEM_PTR)tmTop->histo);
169 tmTop->histo = NULL;
170 }
171
172
173 int
174 tmCvColors(ls, cs, scan, len) /* convert float colors */
175 TMbright *ls;
176 BYTE *cs;
177 COLOR *scan;
178 int len;
179 {
180 static char funcName[] = "tmCvColors";
181 static COLOR csmall = {1e-6, 1e-6, 1e-6};
182 COLOR cmon;
183 double lum, slum;
184 register double d;
185 register int i;
186
187 if (tmTop == NULL)
188 returnErr(TM_E_TMINVAL);
189 if (ls == NULL | scan == NULL | len < 0)
190 returnErr(TM_E_ILLEGAL);
191 for (i = len; i--; ) {
192 if (tmNeedMatrix(tmTop)) /* get monitor RGB */
193 colortrans(cmon, tmTop->cmat, scan[i]);
194 else {
195 cmon[RED] = tmTop->inpsf*scan[i][RED];
196 cmon[GRN] = tmTop->inpsf*scan[i][GRN];
197 cmon[BLU] = tmTop->inpsf*scan[i][BLU];
198 }
199 /* world luminance */
200 lum = tmTop->clf[RED]*cmon[RED] +
201 tmTop->clf[GRN]*cmon[GRN] +
202 tmTop->clf[BLU]*cmon[BLU] ;
203 /* check range */
204 if (clipgamut(cmon, lum, CGAMUT_LOWER, csmall, cwhite))
205 lum = tmTop->clf[RED]*cmon[RED] +
206 tmTop->clf[GRN]*cmon[GRN] +
207 tmTop->clf[BLU]*cmon[BLU] ;
208 if (lum < MINLUM) {
209 ls[i] = MINBRT-1; /* bogus value */
210 lum = MINLUM;
211 } else {
212 d = TM_BRTSCALE*log(lum); /* encode it */
213 ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
214 }
215 if (cs == TM_NOCHROM) /* no color? */
216 continue;
217 if (tmTop->flags & TM_F_MESOPIC && lum < LMESUPPER) {
218 slum = scotlum(cmon); /* mesopic adj. */
219 if (lum < LMESLOWER)
220 cmon[RED] = cmon[GRN] = cmon[BLU] = slum;
221 else {
222 d = (lum - LMESLOWER)/(LMESUPPER - LMESLOWER);
223 if (tmTop->flags & TM_F_BW)
224 cmon[RED] = cmon[GRN] =
225 cmon[BLU] = d*lum;
226 else
227 scalecolor(cmon, d);
228 d = (1.-d)*slum;
229 cmon[RED] += d;
230 cmon[GRN] += d;
231 cmon[BLU] += d;
232 }
233 } else if (tmTop->flags & TM_F_BW) {
234 cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
235 }
236 d = tmTop->clf[RED]*cmon[RED]/lum;
237 cs[3*i ] = d>=.999 ? 255 :
238 (int)(256.*pow(d, 1./tmTop->mongam));
239 d = tmTop->clf[GRN]*cmon[GRN]/lum;
240 cs[3*i+1] = d>=.999 ? 255 :
241 (int)(256.*pow(d, 1./tmTop->mongam));
242 d = tmTop->clf[BLU]*cmon[BLU]/lum;
243 cs[3*i+2] = d>=.999 ? 255 :
244 (int)(256.*pow(d, 1./tmTop->mongam));
245 }
246 returnOK;
247 }
248
249
250 int
251 tmAddHisto(ls, len, wt) /* add values to histogram */
252 register TMbright *ls;
253 int len;
254 int wt;
255 {
256 static char funcName[] = "tmAddHisto";
257 int oldorig, oldlen, horig, hlen;
258 register int i, j;
259
260 if (tmTop == NULL)
261 returnErr(TM_E_TMINVAL);
262 if (len < 0)
263 returnErr(TM_E_ILLEGAL);
264 if (len == 0)
265 returnOK;
266 /* first, grow limits */
267 if (tmTop->histo == NULL) {
268 for (i = len; i-- && ls[i] < MINBRT; )
269 ;
270 if (i < 0)
271 returnOK;
272 tmTop->hbrmin = tmTop->hbrmax = ls[i];
273 oldlen = 0;
274 } else {
275 oldorig = (tmTop->hbrmin-MINBRT)/HISTEP;
276 oldlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - oldorig;
277 }
278 for (i = len; i--; ) {
279 if ((j = ls[i]) < MINBRT)
280 continue;
281 if (j < tmTop->hbrmin)
282 tmTop->hbrmin = j;
283 else if (j > tmTop->hbrmax)
284 tmTop->hbrmax = j;
285 }
286 horig = (tmTop->hbrmin-MINBRT)/HISTEP;
287 hlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - horig;
288 if (hlen > oldlen) { /* (re)allocate histogram */
289 register int *newhist = (int *)calloc(hlen, sizeof(int));
290 if (newhist == NULL)
291 returnErr(TM_E_NOMEM);
292 if (oldlen) { /* copy and free old */
293 for (i = oldlen, j = i+oldorig-horig; i; )
294 newhist[--j] = tmTop->histo[--i];
295 free((MEM_PTR)tmTop->histo);
296 }
297 tmTop->histo = newhist;
298 }
299 if (wt == 0)
300 returnOK;
301 for (i = len; i--; ) /* add in new counts */
302 if (ls[i] >= MINBRT)
303 tmTop->histo[ (ls[i]-MINBRT)/HISTEP - horig ] += wt;
304 returnOK;
305 }
306
307
308 static double
309 htcontrs(La) /* human threshold contrast sensitivity, dL(La) */
310 double La;
311 {
312 double l10La, l10dL;
313 /* formula taken from Ferwerda et al. [SG96] */
314 if (La < 1.148e-4)
315 return(1.38e-3);
316 l10La = log10(La);
317 if (l10La < -1.44) /* rod response regime */
318 l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
319 else if (l10La < -.0184)
320 l10dL = l10La - .395;
321 else if (l10La < 1.9) /* cone response regime */
322 l10dL = pow(.249*l10La + .65, 2.7) - .72;
323 else
324 l10dL = l10La - 1.255;
325
326 return(exp10(l10dL));
327 }
328
329
330 int
331 tmComputeMapping(gamval, Lddyn, Ldmax)
332 double gamval;
333 double Lddyn;
334 double Ldmax;
335 {
336 static char funcName[] = "tmComputeMapping";
337 int *histo;
338 float *cumf;
339 int brt0, histlen, histot, threshold, ceiling, trimmings;
340 double logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
341 int4 sum;
342 register double d;
343 register int i, j;
344
345 if (tmTop == NULL || tmTop->histo == NULL)
346 returnErr(TM_E_TMINVAL);
347 /* check arguments */
348 if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
349 if (Ldmax < MINLDMAX) Ldmax = DEFLDMAX;
350 if (gamval < MINGAM) gamval = tmTop->mongam;
351 /* compute handy values */
352 Ldmin = Ldmax/Lddyn;
353 logLddyn = log(Lddyn);
354 Ldavg = sqrt(Ldmax*Ldmin);
355 i = (tmTop->hbrmin-MINBRT)/HISTEP;
356 brt0 = MINBRT + HISTEP/2 + i*HISTEP;
357 histlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - i;
358 /* histogram total and mean */
359 histot = 0; sum = 0;
360 j = brt0 + histlen*HISTEP;
361 for (i = histlen; i--; ) {
362 histot += tmTop->histo[i];
363 sum += (j -= HISTEP) * tmTop->histo[i];
364 }
365 threshold = histot*.025 + .5;
366 if (threshold < 4)
367 returnErr(TM_E_TMFAIL);
368 Lwavg = tmLuminance( (double)sum / histot );
369 if (!(tmTop->flags & TM_F_LINEAR)) { /* clamp histogram */
370 histo = (int *)malloc(histlen*sizeof(int));
371 cumf = (float *)malloc((histlen+1)*sizeof(float));
372 if (histo == NULL | cumf == NULL)
373 returnErr(TM_E_NOMEM);
374 for (i = histlen; i--; ) /* make malleable copy */
375 histo[i] = tmTop->histo[i];
376 do { /* iterate to solution */
377 sum = 0; /* cumulative probability */
378 for (i = 0; i < histlen; i++) {
379 cumf[i] = (double)sum/histot;
380 sum += histo[i];
381 }
382 cumf[i] = 1.;
383 Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
384 ((double)histlen*TM_BRTSCALE) / logLddyn;
385 ceiling = Tr + 1.;
386 trimmings = 0; /* clip to envelope */
387 for (i = histlen; i--; ) {
388 if (tmTop->flags & TM_F_HCONTR) {
389 Lw = tmLuminance(brt0 + i*HISTEP);
390 Ld = Ldmin * exp( logLddyn *
391 .5*(cumf[i]+cumf[i+1]) );
392 ceiling = Tr * (htcontrs(Ld) * Lw) /
393 (htcontrs(Lw) * Ld) + 1.;
394 }
395 if (histo[i] > ceiling) {
396 trimmings += histo[i] - ceiling;
397 histo[i] = ceiling;
398 }
399 }
400 } while ((histot -= trimmings) > threshold &&
401 trimmings > threshold);
402 }
403 /* allocate luminance map */
404 if (tmTop->lumap != NULL)
405 free((MEM_PTR)tmTop->lumap);
406 tmTop->mbrmin = tmTop->hbrmin;
407 tmTop->mbrmax = tmTop->hbrmax;
408 tmTop->lumap = (unsigned short *)malloc(
409 (tmTop->mbrmax-tmTop->mbrmin+1)*sizeof(unsigned short) );
410 if (tmTop->lumap == NULL)
411 returnErr(TM_E_NOMEM);
412 if (tmTop->flags & TM_F_LINEAR || histot <= threshold) {
413 /* linear tone mapping */
414 if (tmTop->flags & TM_F_HCONTR)
415 d = htcontrs(Ldavg) / htcontrs(Lwavg);
416 else
417 d = Ldavg / Lwavg;
418 d = log(d/Ldmax);
419 for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
420 tmTop->lumap[i] = 256. * exp(
421 ( d + (tmTop->mbrmin+i)/(double)TM_BRTSCALE )
422 / gamval );
423 } else {
424 /* histogram adjustment */
425 for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; ) {
426 j = d = (double)i/(tmTop->mbrmax-tmTop->mbrmin)*histlen;
427 d -= (double)j;
428 Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
429 d = (Ld - Ldmin)/(Ldmax - Ldmin);
430 tmTop->lumap[i] = 256.*pow(d, 1./gamval);
431 }
432 }
433 if (!(tmTop->flags & TM_F_LINEAR)) {
434 free((MEM_PTR)histo);
435 free((MEM_PTR)cumf);
436 }
437 returnOK;
438 }
439
440
441 int
442 tmMapPixels(ps, ls, cs, len)
443 register BYTE *ps;
444 TMbright *ls;
445 register BYTE *cs;
446 int len;
447 {
448 static char funcName[] = "tmMapPixels";
449 register int4 li, pv;
450
451 if (tmTop == NULL || tmTop->lumap == NULL)
452 returnErr(TM_E_TMINVAL);
453 if (ps == NULL | ls == NULL | len < 0)
454 returnErr(TM_E_ILLEGAL);
455 while (len--) {
456 if ((li = *ls++) < tmTop->mbrmin)
457 li = tmTop->mbrmin;
458 else if (li > tmTop->mbrmax)
459 li = tmTop->mbrmax;
460 li = tmTop->lumap[li - tmTop->mbrmin];
461 if (cs == TM_NOCHROM)
462 *ps++ = li>255 ? 255 : li;
463 else {
464 pv = *cs++ * li / tmTop->cdiv[RED];
465 *ps++ = pv>255 ? 255 : pv;
466 pv = *cs++ * li / tmTop->cdiv[GRN];
467 *ps++ = pv>255 ? 255 : pv;
468 pv = *cs++ * li / tmTop->cdiv[BLU];
469 *ps++ = pv>255 ? 255 : pv;
470 }
471 }
472 returnOK;
473 }
474
475
476 struct tmStruct *
477 tmPop() /* pop top tone mapping off stack */
478 {
479 register struct tmStruct *tms;
480
481 if ((tms = tmTop) != NULL)
482 tmTop = tms->tmprev;
483 return(tms);
484 }
485
486
487 int
488 tmPull(tms) /* pull a tone mapping from stack */
489 register struct tmStruct *tms;
490 {
491 register struct tmStruct *tms2;
492 /* special cases first */
493 if (tms == NULL | tmTop == NULL)
494 return(0);
495 if (tms == tmTop) {
496 tmTop = tms->tmprev;
497 tms->tmprev = NULL;
498 return(1);
499 }
500 for (tms2 = tmTop; tms2->tmprev != NULL; tms2 = tms2->tmprev)
501 if (tms == tms2->tmprev) { /* remove it */
502 tms2->tmprev = tms->tmprev;
503 tms->tmprev = NULL;
504 return(1);
505 }
506 return(0); /* not found on stack */
507 }
508
509
510 struct tmStruct *
511 tmDup() /* duplicate top tone mapping */
512 {
513 int len;
514 register int i;
515 register struct tmStruct *tmnew;
516
517 if (tmTop == NULL) /* anything to duplicate? */
518 return(NULL);
519 tmnew = (struct tmStruct *)malloc(sizeof(struct tmStruct));
520 if (tmnew == NULL)
521 return(NULL);
522 *tmnew = *tmTop; /* copy everything */
523 if (tmnew->histo != NULL) { /* duplicate histogram */
524 len = (tmnew->hbrmax-MINBRT)/HISTEP + 1 -
525 (tmnew->hbrmin-MINBRT)/HISTEP;
526 tmnew->histo = (int *)malloc(len*sizeof(int));
527 if (tmnew->histo != NULL)
528 for (i = len; i--; )
529 tmnew->histo[i] = tmTop->histo[i];
530 }
531 if (tmnew->lumap != NULL) { /* duplicate luminance mapping */
532 len = tmnew->mbrmax-tmnew->mbrmin+1;
533 tmnew->lumap = (unsigned short *)malloc(
534 len*sizeof(unsigned short) );
535 if (tmnew->lumap != NULL)
536 for (i = len; i--; )
537 tmnew->lumap[i] = tmTop->lumap[i];
538 }
539 /* clear package data */
540 for (i = tmNumPkgs; i--; )
541 tmnew->pd[i] = NULL;
542 tmnew->tmprev = tmTop; /* make copy current */
543 return(tmTop = tmnew);
544 }
545
546
547 int
548 tmPush(tms) /* push tone mapping on top of stack */
549 register struct tmStruct *tms;
550 {
551 static char funcName[] = "tmPush";
552 /* check validity */
553 if (tms == NULL)
554 returnErr(TM_E_ILLEGAL);
555 if (tms == tmTop) /* check necessity */
556 returnOK;
557 /* pull if already in stack */
558 (void)tmPull(tms);
559 /* push it on top */
560 tms->tmprev = tmTop;
561 tmTop = tms;
562 returnOK;
563 }
564
565
566 void
567 tmDone(tms) /* done with tone mapping -- destroy it */
568 register struct tmStruct *tms;
569 {
570 register int i;
571 /* NULL arg. is equiv. to tmTop */
572 if (tms == NULL && (tms = tmTop) == NULL)
573 return;
574 /* take out of stack if present */
575 (void)tmPull(tms);
576 /* free tables */
577 if (tms->histo != NULL)
578 free((MEM_PTR)tms->histo);
579 if (tms->lumap != NULL)
580 free((MEM_PTR)tms->lumap);
581 /* free private data */
582 for (i = tmNumPkgs; i--; )
583 if (tms->pd[i] != NULL)
584 (*tmPkg[i]->Free)(tms->pd[i]);
585 free((MEM_PTR)tms); /* free basic structure */
586 }