ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.47
Committed: Wed Apr 7 21:13:52 2021 UTC (3 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.46: +6 -6 lines
Log Message:
perf: changed default dynamic range from 32 to 100 to reflect modern displays

File Contents

# User Rev Content
1 greg 3.1 #ifndef lint
2 greg 3.47 static const char RCSid[] = "$Id: tonemap.c,v 3.46 2021/03/02 20:09:14 greg Exp $";
3 greg 3.1 #endif
4     /*
5     * Tone mapping functions.
6     * See tonemap.h for detailed function descriptions.
7 greg 3.9 * Added von Kries white-balance calculations 10/01 (GW).
8     *
9     * Externals declared in tonemap.h
10     */
11    
12 greg 3.10 #include "copyright.h"
13 greg 3.1
14     #include <stdio.h>
15 greg 3.34 #include <stdlib.h>
16 greg 3.1 #include <math.h>
17 greg 3.38 #include <string.h>
18 greg 3.1 #include "tmprivat.h"
19     #include "tmerrmsg.h"
20    
21     #define exp10(x) exp(M_LN10*(x))
22    
23 greg 3.4 /* our list of conversion packages */
24     struct tmPackage *tmPkg[TM_MAXPKG];
25     int tmNumPkgs = 0; /* number of registered packages */
26 greg 3.1
27 greg 3.23 /* luminance->brightness lookup */
28 greg 3.24 static TMbright *tmFloat2BrtLUT = NULL;
29 greg 3.23
30     #define tmCvLumLUfp(pf) tmFloat2BrtLUT[*(int32 *)(pf) >> 15]
31    
32 gregl 3.5
33 greg 3.16 TMstruct *
34     tmInit( /* initialize new tone mapping */
35     int flags,
36     RGBPRIMP monpri,
37     double gamval
38     )
39 greg 3.1 {
40     COLORMAT cmat;
41 greg 3.16 TMstruct *tmnew;
42     int i;
43 greg 3.1 /* allocate structure */
44 greg 3.16 tmnew = (TMstruct *)malloc(sizeof(TMstruct));
45 greg 3.1 if (tmnew == NULL)
46     return(NULL);
47    
48     tmnew->flags = flags & ~TM_F_UNIMPL;
49 greg 3.15 if (tmnew->flags & TM_F_BW)
50     tmnew->flags &= ~TM_F_MESOPIC;
51 greg 3.1 /* set monitor transform */
52     if (monpri == NULL || monpri == stdprims || tmnew->flags & TM_F_BW) {
53     tmnew->monpri = stdprims;
54     tmnew->clf[RED] = rgb2xyzmat[1][0];
55     tmnew->clf[GRN] = rgb2xyzmat[1][1];
56     tmnew->clf[BLU] = rgb2xyzmat[1][2];
57     } else {
58 greg 3.26 comprgb2xyzmat(cmat, tmnew->monpri=monpri);
59 greg 3.1 tmnew->clf[RED] = cmat[1][0];
60     tmnew->clf[GRN] = cmat[1][1];
61     tmnew->clf[BLU] = cmat[1][2];
62     }
63     /* set gamma value */
64     if (gamval < MINGAM)
65     tmnew->mongam = DEFGAM;
66     else
67     tmnew->mongam = gamval;
68 greg 3.4 /* set color divisors */
69     for (i = 0; i < 3; i++)
70 greg 3.40 tmnew->cdiv[i] = TM_BRES*pow(tmnew->clf[i], 1./tmnew->mongam);
71 greg 3.4
72 greg 3.1 /* set input transform */
73     tmnew->inppri = tmnew->monpri;
74     tmnew->cmat[0][0] = tmnew->cmat[1][1] = tmnew->cmat[2][2] =
75     tmnew->inpsf = WHTEFFICACY;
76     tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
77     tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
78 greg 3.18 tmnew->inpdat = NULL;
79 greg 3.9 tmnew->hbrmin = 10; tmnew->hbrmax = -10;
80 greg 3.1 tmnew->histo = NULL;
81 greg 3.9 tmnew->mbrmin = 10; tmnew->mbrmax = -10;
82 greg 3.1 tmnew->lumap = NULL;
83 greg 3.4 /* zero private data */
84     for (i = TM_MAXPKG; i--; )
85     tmnew->pd[i] = NULL;
86 greg 3.17 tmnew->lastError = TM_E_OK;
87     tmnew->lastFunc = "NoErr";
88 greg 3.16 /* return new TMstruct */
89     return(tmnew);
90 greg 3.1 }
91    
92    
93     int
94 greg 3.16 tmSetSpace( /* set input color space for conversions */
95     TMstruct *tms,
96     RGBPRIMP pri,
97 greg 3.18 double sf,
98     MEM_PTR dat
99 greg 3.16 )
100 greg 3.1 {
101 greg 3.17 static const char funcName[] = "tmSetSpace";
102 greg 3.16 int i, j;
103 greg 3.1 /* error check */
104 greg 3.16 if (tms == NULL)
105 greg 3.1 returnErr(TM_E_TMINVAL);
106     if (sf <= 1e-12)
107     returnErr(TM_E_ILLEGAL);
108     /* check if no change */
109 greg 3.18 if (pri == tms->inppri && FEQ(sf, tms->inpsf) && dat == tms->inpdat)
110 greg 3.1 returnOK;
111 greg 3.16 tms->inppri = pri; /* let's set it */
112     tms->inpsf = sf;
113 greg 3.18 tms->inpdat = dat;
114 greg 3.1
115 greg 3.16 if (tms->flags & TM_F_BW) { /* color doesn't matter */
116     tms->monpri = tms->inppri; /* eliminate xform */
117     if (tms->inppri == TM_XYZPRIM) {
118     tms->clf[CIEX] = tms->clf[CIEZ] = 0.;
119     tms->clf[CIEY] = 1.;
120 greg 3.1 } else {
121 greg 3.26 comprgb2xyzmat(tms->cmat, tms->monpri);
122 greg 3.16 tms->clf[RED] = tms->cmat[1][0];
123     tms->clf[GRN] = tms->cmat[1][1];
124     tms->clf[BLU] = tms->cmat[1][2];
125     }
126 greg 3.46 tms->cmat[0][0] = tms->cmat[1][1] = tms->cmat[2][2] = 1.;
127 greg 3.16 tms->cmat[0][1] = tms->cmat[0][2] = tms->cmat[1][0] =
128     tms->cmat[1][2] = tms->cmat[2][0] = tms->cmat[2][1] = 0.;
129 greg 3.1
130 greg 3.46 } else if (tms->inppri == TM_XYZPRIM) { /* input is XYZ */
131 greg 3.16 compxyz2rgbWBmat(tms->cmat, tms->monpri);
132 greg 3.1
133 greg 3.46 } else { /* input is RGB */
134 greg 3.16 if (tms->inppri != tms->monpri &&
135     PRIMEQ(tms->inppri, tms->monpri))
136     tms->inppri = tms->monpri; /* no xform */
137 greg 3.36 if (!comprgb2rgbWBmat(tms->cmat, tms->inppri, tms->monpri))
138     returnErr(TM_E_ILLEGAL);
139 greg 3.1 }
140     for (i = 0; i < 3; i++)
141     for (j = 0; j < 3; j++)
142 greg 3.16 tms->cmat[i][j] *= tms->inpsf;
143 greg 3.4 /* set color divisors */
144     for (i = 0; i < 3; i++)
145 greg 3.40 tms->cdiv[i] = TM_BRES*pow(tms->clf[i] < .001 ? .001 :
146 greg 3.36 tms->clf[i], 1./tms->mongam);
147 greg 3.4 /* notify packages */
148     for (i = tmNumPkgs; i--; )
149 greg 3.16 if (tms->pd[i] != NULL && tmPkg[i]->NewSpace != NULL)
150     (*tmPkg[i]->NewSpace)(tms);
151 greg 3.1 returnOK;
152     }
153    
154    
155     void
156 greg 3.23 tmClearHisto( /* clear current histogram */
157 greg 3.16 TMstruct *tms
158     )
159 greg 3.1 {
160 greg 3.16 if (tms == NULL || tms->histo == NULL)
161 greg 3.1 return;
162 greg 3.46 free(tms->histo);
163 greg 3.16 tms->histo = NULL;
164 greg 3.1 }
165    
166    
167 greg 3.23 TMbright
168     tmCvLuminance( /* convert a single luminance */
169     double lum
170     )
171     {
172     double d;
173    
174     #ifdef isfinite
175     if (!isfinite(lum) || lum <= TM_NOLUM)
176     #else
177     if (lum <= TM_NOLUM)
178     #endif
179     return(TM_NOBRT);
180     d = TM_BRTSCALE*log(lum);
181 greg 3.35 return((TMbright)(d + .5 - (d < 0.)));
182 greg 3.23 }
183    
184    
185     int
186     tmCvLums( /* convert luminances using lookup */
187     TMbright *ls,
188     float *scan,
189     int len
190     )
191     {
192     if (tmFloat2BrtLUT == NULL) { /* initialize lookup table */
193     int32 i;
194     tmFloat2BrtLUT = (TMbright *)malloc(sizeof(TMbright)*0x10000);
195     if (tmFloat2BrtLUT == NULL)
196     return(TM_E_NOMEM);
197     for (i = 0; i < 0x10000; i++) {
198     int32 l = (i<<1 | 1) << 14;
199     #ifndef isfinite
200     if ((l & 0x7f800000) == 0x7f800000)
201     tmFloat2BrtLUT[i] = TM_NOBRT;
202     else
203     #endif
204     tmFloat2BrtLUT[i] = tmCvLuminance(*(float *)&l);
205     }
206     }
207     if (len <= 0)
208     return(TM_E_OK);
209     if ((ls == NULL) | (scan == NULL))
210     return(TM_E_ILLEGAL);
211     while (len--) {
212     if (*scan <= TM_NOLUM) {
213     *ls++ = TM_NOBRT;
214     ++scan;
215     continue;
216     }
217     *ls++ = tmCvLumLUfp(scan++);
218     }
219     return(TM_E_OK);
220     }
221    
222    
223     int
224     tmCvGrays( /* convert float gray values */
225     TMstruct *tms,
226     TMbright *ls,
227     float *scan,
228     int len
229     )
230     {
231     static const char funcName[] = "tmCvGrays";
232     int i;
233    
234     if (tms == NULL)
235     returnErr(TM_E_TMINVAL);
236     if ((ls == NULL) | (scan == NULL) | (len < 0))
237     returnErr(TM_E_ILLEGAL);
238 greg 3.24 if (tmFloat2BrtLUT == NULL) /* initialize */
239     tmCvLums(NULL, NULL, 0);
240 greg 3.23 for (i = len; i--; ) {
241     float lum = tms->inpsf * scan[i];
242     if (lum <= TM_NOLUM)
243     ls[i] = TM_NOBRT;
244     else
245     ls[i] = tmCvLumLUfp(&lum);
246     }
247     returnOK;
248     }
249    
250    
251 greg 3.1 int
252 greg 3.16 tmCvColors( /* convert float colors */
253     TMstruct *tms,
254     TMbright *ls,
255 greg 3.37 uby8 *cs,
256 greg 3.16 COLOR *scan,
257     int len
258     )
259 greg 3.1 {
260 greg 3.17 static const char funcName[] = "tmCvColors";
261 greg 3.37 static uby8 gamtab[1024];
262 greg 3.20 static double curgam = .0;
263 greg 3.1 COLOR cmon;
264 greg 3.23 float lum, slum, d;
265 greg 3.16 int i;
266 greg 3.1
267 greg 3.16 if (tms == NULL)
268 greg 3.1 returnErr(TM_E_TMINVAL);
269 schorsch 3.13 if ((ls == NULL) | (scan == NULL) | (len < 0))
270 greg 3.1 returnErr(TM_E_ILLEGAL);
271 greg 3.24 if (tmFloat2BrtLUT == NULL) /* initialize */
272     tmCvLums(NULL, NULL, 0);
273 greg 3.21 if (cs != TM_NOCHROM && fabs(tms->mongam - curgam) > .02) {
274 greg 3.20 curgam = tms->mongam; /* (re)build table */
275     for (i = 1024; i--; )
276     gamtab[i] = (int)(256.*pow((i+.5)/1024., 1./curgam));
277     }
278 greg 3.1 for (i = len; i--; ) {
279 greg 3.16 if (tmNeedMatrix(tms)) { /* get monitor RGB */
280     colortrans(cmon, tms->cmat, scan[i]);
281 greg 3.15 } else {
282 greg 3.16 cmon[RED] = tms->inpsf*scan[i][RED];
283     cmon[GRN] = tms->inpsf*scan[i][GRN];
284     cmon[BLU] = tms->inpsf*scan[i][BLU];
285 greg 3.1 }
286 greg 3.23 #ifdef isfinite
287     if (!isfinite(cmon[RED]) || cmon[RED] < .0f) cmon[RED] = .0f;
288     if (!isfinite(cmon[GRN]) || cmon[GRN] < .0f) cmon[GRN] = .0f;
289     if (!isfinite(cmon[BLU]) || cmon[BLU] < .0f) cmon[BLU] = .0f;
290     #else
291     if (cmon[RED] < .0f) cmon[RED] = .0f;
292     if (cmon[GRN] < .0f) cmon[GRN] = .0f;
293     if (cmon[BLU] < .0f) cmon[BLU] = .0f;
294 greg 3.22 #endif
295 greg 3.1 /* world luminance */
296 greg 3.16 lum = tms->clf[RED]*cmon[RED] +
297     tms->clf[GRN]*cmon[GRN] +
298     tms->clf[BLU]*cmon[BLU] ;
299 greg 3.25 if (lum <= TM_NOLUM) { /* convert brightness */
300     lum = cmon[RED] = cmon[GRN] = cmon[BLU] = TM_NOLUM;
301 greg 3.23 ls[i] = TM_NOBRT;
302 greg 3.25 } else
303 greg 3.23 ls[i] = tmCvLumLUfp(&lum);
304 greg 3.1 if (cs == TM_NOCHROM) /* no color? */
305     continue;
306 greg 3.16 if (tms->flags & TM_F_MESOPIC && lum < LMESUPPER) {
307 greg 3.1 slum = scotlum(cmon); /* mesopic adj. */
308 greg 3.25 if (lum < LMESLOWER) {
309 greg 3.1 cmon[RED] = cmon[GRN] = cmon[BLU] = slum;
310 greg 3.25 } else {
311 greg 3.1 d = (lum - LMESLOWER)/(LMESUPPER - LMESLOWER);
312 greg 3.16 if (tms->flags & TM_F_BW)
313 greg 3.1 cmon[RED] = cmon[GRN] =
314     cmon[BLU] = d*lum;
315     else
316     scalecolor(cmon, d);
317 greg 3.23 d = (1.f-d)*slum;
318 greg 3.1 cmon[RED] += d;
319     cmon[GRN] += d;
320     cmon[BLU] += d;
321     }
322 greg 3.16 } else if (tms->flags & TM_F_BW) {
323 greg 3.1 cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
324     }
325 greg 3.16 d = tms->clf[RED]*cmon[RED]/lum;
326 greg 3.23 cs[3*i ] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
327 greg 3.16 d = tms->clf[GRN]*cmon[GRN]/lum;
328 greg 3.23 cs[3*i+1] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
329 greg 3.16 d = tms->clf[BLU]*cmon[BLU]/lum;
330 greg 3.23 cs[3*i+2] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
331 greg 3.1 }
332     returnOK;
333     }
334    
335    
336 greg 3.9 int
337 greg 3.16 tmAddHisto( /* add values to histogram */
338     TMstruct *tms,
339     TMbright *ls,
340     int len,
341     int wt
342     )
343 greg 3.1 {
344 greg 3.17 static const char funcName[] = "tmAddHisto";
345 greg 3.9 int oldorig=0, oldlen, horig, hlen;
346 greg 3.16 int i, j;
347 greg 3.1
348 greg 3.16 if (tms == NULL)
349 greg 3.1 returnErr(TM_E_TMINVAL);
350 gwlarson 3.8 if (len < 0)
351     returnErr(TM_E_ILLEGAL);
352     if (len == 0)
353     returnOK;
354 greg 3.1 /* first, grow limits */
355 greg 3.16 if (tms->histo == NULL) {
356 greg 3.1 for (i = len; i-- && ls[i] < MINBRT; )
357     ;
358     if (i < 0)
359     returnOK;
360 greg 3.16 tms->hbrmin = tms->hbrmax = ls[i];
361 greg 3.1 oldlen = 0;
362     } else {
363 greg 3.31 oldorig = HISTI(tms->hbrmin);
364     oldlen = HISTI(tms->hbrmax) + 1 - oldorig;
365 greg 3.1 }
366     for (i = len; i--; ) {
367     if ((j = ls[i]) < MINBRT)
368     continue;
369 greg 3.16 if (j < tms->hbrmin)
370     tms->hbrmin = j;
371     else if (j > tms->hbrmax)
372     tms->hbrmax = j;
373 greg 3.1 }
374 greg 3.31 horig = HISTI(tms->hbrmin);
375     hlen = HISTI(tms->hbrmax) + 1 - horig;
376 greg 3.1 if (hlen > oldlen) { /* (re)allocate histogram */
377 greg 3.38 HIST_TYP *newhist = (HIST_TYP *)calloc(hlen, sizeof(HIST_TYP));
378 greg 3.1 if (newhist == NULL)
379     returnErr(TM_E_NOMEM);
380     if (oldlen) { /* copy and free old */
381 greg 3.46 memcpy(newhist+(oldorig-horig),
382     tms->histo, oldlen*sizeof(HIST_TYP));
383     free(tms->histo);
384 greg 3.1 }
385 greg 3.16 tms->histo = newhist;
386 greg 3.1 }
387     if (wt == 0)
388     returnOK;
389     for (i = len; i--; ) /* add in new counts */
390     if (ls[i] >= MINBRT)
391 greg 3.31 tms->histo[ HISTI(ls[i]) - horig ] += wt;
392 greg 3.1 returnOK;
393     }
394    
395    
396     static double
397 greg 3.16 htcontrs( /* human threshold contrast sensitivity, dL(La) */
398     double La
399     )
400 greg 3.1 {
401     double l10La, l10dL;
402     /* formula taken from Ferwerda et al. [SG96] */
403     if (La < 1.148e-4)
404     return(1.38e-3);
405     l10La = log10(La);
406     if (l10La < -1.44) /* rod response regime */
407     l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
408     else if (l10La < -.0184)
409     l10dL = l10La - .395;
410     else if (l10La < 1.9) /* cone response regime */
411     l10dL = pow(.249*l10La + .65, 2.7) - .72;
412     else
413     l10dL = l10La - 1.255;
414    
415     return(exp10(l10dL));
416     }
417    
418    
419 greg 3.9 int
420 greg 3.16 tmFixedMapping( /* compute fixed, linear tone-mapping */
421     TMstruct *tms,
422     double expmult,
423 greg 3.46 double gamval,
424     double Lddyn
425 greg 3.16 )
426 greg 3.9 {
427 greg 3.17 static const char funcName[] = "tmFixedMapping";
428 greg 3.47 const int maxV = (1L<<(8*sizeof(TMAP_TYP))) - 1;
429     double minD;
430     int i;
431 greg 3.9
432 greg 3.16 if (!tmNewMap(tms))
433 greg 3.9 returnErr(TM_E_NOMEM);
434 greg 3.46 /* check arguments */
435     if (expmult <= .0) expmult = 1.;
436     if (gamval < MINGAM) gamval = tms->mongam;
437     if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
438     minD = 1./Lddyn;
439 greg 3.28 for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
440 greg 3.46 double d;
441     d = expmult/tms->inpsf * tmLuminance(tms->mbrmin + i);
442     if (d <= minD)
443     break; /* map initialized to zeroes */
444     d = (d - minD)/(1. - minD);
445     d = TM_BRES*pow(d, 1./gamval);
446 greg 3.47 tms->lumap[i] = (d > maxV) ? maxV : (int)d;
447 greg 3.28 }
448 greg 3.9 returnOK;
449     }
450    
451    
452 greg 3.1 int
453 greg 3.16 tmComputeMapping( /* compute histogram tone-mapping */
454     TMstruct *tms,
455     double gamval,
456     double Lddyn,
457     double Ldmax
458     )
459 greg 3.1 {
460 greg 3.17 static const char funcName[] = "tmComputeMapping";
461 greg 3.38 HIST_TYP *histo;
462 greg 3.1 float *cumf;
463 greg 3.38 int brt0, histlen;
464 greg 3.39 HIST_TYP threshold, ceiling, trimmings, histot;
465 greg 3.44 double logLddyn, Ldmin, Lwavg, Tr, Lw, Ld;
466 greg 3.9 double sum;
467 greg 3.16 double d;
468     int i, j;
469 greg 3.1
470 greg 3.16 if (tms == NULL || tms->histo == NULL)
471 greg 3.1 returnErr(TM_E_TMINVAL);
472     /* check arguments */
473     if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
474     if (Ldmax < MINLDMAX) Ldmax = DEFLDMAX;
475 greg 3.16 if (gamval < MINGAM) gamval = tms->mongam;
476 greg 3.1 /* compute handy values */
477     Ldmin = Ldmax/Lddyn;
478     logLddyn = log(Lddyn);
479 greg 3.31 i = HISTI(tms->hbrmin);
480 greg 3.32 brt0 = HISTV(i);
481 greg 3.31 histlen = HISTI(tms->hbrmax) + 1 - i;
482 greg 3.1 /* histogram total and mean */
483     histot = 0; sum = 0;
484     j = brt0 + histlen*HISTEP;
485     for (i = histlen; i--; ) {
486 greg 3.16 histot += tms->histo[i];
487 greg 3.30 sum += (double)(j -= HISTEP) * tms->histo[i];
488 greg 3.1 }
489 greg 3.45 threshold = histot*0.002 + .5;
490 greg 3.29 if (!histot)
491 greg 3.1 returnErr(TM_E_TMFAIL);
492     Lwavg = tmLuminance( (double)sum / histot );
493 greg 3.29 /* use linear tone mapping? */
494 greg 3.30 if (tms->flags & TM_F_LINEAR || threshold < 4 ||
495     tms->hbrmax - tms->hbrmin < TM_BRTSCALE*logLddyn)
496 greg 3.29 goto linearmap;
497 greg 3.9 /* clamp histogram */
498 greg 3.38 histo = (HIST_TYP *)malloc(histlen*sizeof(HIST_TYP));
499 greg 3.9 cumf = (float *)malloc((histlen+2)*sizeof(float));
500 schorsch 3.13 if ((histo == NULL) | (cumf == NULL))
501 greg 3.9 returnErr(TM_E_NOMEM);
502     cumf[histlen+1] = 1.; /* guard for assignment code */
503 greg 3.42 /* make malleable copy */
504 greg 3.46 memcpy(histo, tms->histo, histlen*sizeof(HIST_TYP));
505 greg 3.9 do { /* iterate to solution */
506     sum = 0; /* cumulative probability */
507     for (i = 0; i < histlen; i++) {
508 greg 3.45 cumf[i] = sum/histot;
509     sum += (double)histo[i];
510 greg 3.9 }
511     cumf[histlen] = 1.;
512 greg 3.16 Tr = histot * (double)(tms->hbrmax - tms->hbrmin) /
513 greg 3.45 ((double)TM_BRTSCALE*histlen*logLddyn);
514 greg 3.9 ceiling = Tr + 1.;
515     trimmings = 0; /* clip to envelope */
516     for (i = histlen; i--; ) {
517 greg 3.16 if (tms->flags & TM_F_HCONTR) {
518 greg 3.9 Lw = tmLuminance(brt0 + i*HISTEP);
519     Ld = Ldmin * exp( logLddyn *
520     .5*(cumf[i]+cumf[i+1]) );
521     ceiling = Tr * (htcontrs(Ld) * Lw) /
522     (htcontrs(Lw) * Ld) + 1.;
523 greg 3.2 }
524 greg 3.9 if (histo[i] > ceiling) {
525     trimmings += histo[i] - ceiling;
526     histo[i] = ceiling;
527 greg 3.1 }
528     }
529 greg 3.9 /* check if we're out of data */
530     if ((histot -= trimmings) <= threshold) {
531 greg 3.46 free(histo);
532     free(cumf);
533 greg 3.9 goto linearmap;
534     }
535     } while (trimmings > threshold);
536 greg 3.30 /* allocate space for mapping */
537     if (!tmNewMap(tms))
538     returnErr(TM_E_NOMEM);
539 greg 3.9 /* assign tone-mapping */
540 greg 3.16 for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
541     j = d = (double)i/(tms->mbrmax-tms->mbrmin)*histlen;
542 greg 3.9 d -= (double)j;
543     Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
544     d = (Ld - Ldmin)/(Ldmax - Ldmin);
545 greg 3.40 tms->lumap[i] = TM_BRES*pow(d, 1./gamval);
546 greg 3.1 }
547 greg 3.46 free(histo); /* clean up and return */
548     free(cumf);
549 greg 3.1 returnOK;
550 greg 3.9 linearmap: /* linear tone-mapping */
551 greg 3.16 if (tms->flags & TM_F_HCONTR)
552 greg 3.44 d = htcontrs(sqrt(Ldmax*Ldmin)) / htcontrs(Lwavg);
553 greg 3.9 else
554 greg 3.44 d = Ldmax / tmLuminance(tms->hbrmax);
555 greg 3.46 return(tmFixedMapping(tms, tms->inpsf*d/Ldmax, gamval, Lddyn));
556 greg 3.1 }
557    
558    
559     int
560 greg 3.16 tmMapPixels( /* apply tone-mapping to pixel(s) */
561     TMstruct *tms,
562 greg 3.37 uby8 *ps,
563 greg 3.16 TMbright *ls,
564 greg 3.37 uby8 *cs,
565 greg 3.16 int len
566     )
567 greg 3.1 {
568 greg 3.17 static const char funcName[] = "tmMapPixels";
569 greg 3.41 TMbright lv;
570     TMAP_TYP li;
571     int pv;
572 greg 3.1
573 greg 3.16 if (tms == NULL || tms->lumap == NULL)
574 greg 3.1 returnErr(TM_E_TMINVAL);
575 schorsch 3.13 if ((ps == NULL) | (ls == NULL) | (len < 0))
576 greg 3.1 returnErr(TM_E_ILLEGAL);
577     while (len--) {
578 greg 3.41 if ((lv = *ls++) < tms->mbrmin) {
579 greg 3.9 li = 0;
580     } else {
581 greg 3.41 if (lv > tms->mbrmax)
582     lv = tms->mbrmax;
583     li = tms->lumap[lv - tms->mbrmin];
584 greg 3.9 }
585 greg 3.43 if (cs == TM_NOCHROM) {
586     #if !(TM_BRES & 0xff)
587     *ps++ = li>=TM_BRES ? 255 : li/(TM_BRES>>8);
588     #else
589     *ps++ = li>=TM_BRES ? 255 : (li<<8)/TM_BRES;
590     #endif
591     } else {
592 greg 3.16 pv = *cs++ * li / tms->cdiv[RED];
593 greg 3.1 *ps++ = pv>255 ? 255 : pv;
594 greg 3.16 pv = *cs++ * li / tms->cdiv[GRN];
595 greg 3.1 *ps++ = pv>255 ? 255 : pv;
596 greg 3.16 pv = *cs++ * li / tms->cdiv[BLU];
597 greg 3.1 *ps++ = pv>255 ? 255 : pv;
598     }
599     }
600     returnOK;
601     }
602    
603    
604 greg 3.16 TMstruct *
605 greg 3.42 tmDup( /* duplicate tone mapping */
606 greg 3.16 TMstruct *tms
607     )
608 greg 3.3 {
609     int len;
610 greg 3.16 int i;
611     TMstruct *tmnew;
612 greg 3.3
613 greg 3.16 if (tms == NULL) /* anything to duplicate? */
614 greg 3.3 return(NULL);
615 greg 3.16 tmnew = (TMstruct *)malloc(sizeof(TMstruct));
616 greg 3.3 if (tmnew == NULL)
617     return(NULL);
618 greg 3.16 *tmnew = *tms; /* copy everything */
619 greg 3.3 if (tmnew->histo != NULL) { /* duplicate histogram */
620 greg 3.31 len = HISTI(tmnew->hbrmax) + 1 - HISTI(tmnew->hbrmin);
621 greg 3.38 tmnew->histo = (HIST_TYP *)malloc(len*sizeof(HIST_TYP));
622 greg 3.3 if (tmnew->histo != NULL)
623 greg 3.46 memcpy(tmnew->histo, tms->histo, len*sizeof(HIST_TYP));
624 greg 3.3 }
625     if (tmnew->lumap != NULL) { /* duplicate luminance mapping */
626 gregl 3.6 len = tmnew->mbrmax-tmnew->mbrmin+1;
627 greg 3.38 tmnew->lumap = (TMAP_TYP *)malloc(len*sizeof(TMAP_TYP));
628 greg 3.3 if (tmnew->lumap != NULL)
629 greg 3.46 memcpy(tmnew->lumap, tms->lumap, len*sizeof(TMAP_TYP));
630 greg 3.3 }
631 greg 3.4 /* clear package data */
632     for (i = tmNumPkgs; i--; )
633     tmnew->pd[i] = NULL;
634 greg 3.16 /* return copy */
635     return(tmnew);
636 greg 3.1 }
637    
638    
639     void
640     tmDone(tms) /* done with tone mapping -- destroy it */
641 greg 3.16 TMstruct *tms;
642 greg 3.1 {
643 greg 3.16 int i;
644     /* NULL arg. is equiv. to tms */
645     if (tms == NULL)
646 greg 3.1 return;
647     /* free tables */
648     if (tms->histo != NULL)
649 greg 3.46 free(tms->histo);
650 greg 3.1 if (tms->lumap != NULL)
651 greg 3.46 free(tms->lumap);
652 greg 3.4 /* free private data */
653     for (i = tmNumPkgs; i--; )
654     if (tms->pd[i] != NULL)
655     (*tmPkg[i]->Free)(tms->pd[i]);
656 greg 3.46 free(tms); /* free basic structure */
657 greg 3.9 }
658    
659     /******************** Shared but Private library routines *********************/
660    
661 greg 3.37 uby8 tmMesofact[BMESUPPER-BMESLOWER];
662 greg 3.9
663     void
664     tmMkMesofact() /* build mesopic lookup factor table */
665     {
666 greg 3.16 int i;
667 greg 3.9
668     if (tmMesofact[BMESUPPER-BMESLOWER-1])
669     return;
670    
671     for (i = BMESLOWER; i < BMESUPPER; i++)
672     tmMesofact[i-BMESLOWER] = 256. *
673     (tmLuminance(i) - LMESLOWER) /
674     (LMESUPPER - LMESLOWER);
675     }
676    
677    
678     int
679 greg 3.33 tmNewMap( /* allocate new tone-mapping array */
680     TMstruct *tms
681     )
682     {
683     if (tms->lumap != NULL && (tms->mbrmax - tms->mbrmin) !=
684     (tms->hbrmax - tms->hbrmin)) {
685 greg 3.46 free(tms->lumap);
686 greg 3.33 tms->lumap = NULL;
687     }
688     tms->mbrmin = tms->hbrmin;
689     tms->mbrmax = tms->hbrmax;
690     if (tms->mbrmin > tms->mbrmax)
691 greg 3.47 return(0);
692 greg 3.33 if (tms->lumap == NULL)
693 greg 3.46 tms->lumap = (TMAP_TYP *)calloc(tms->mbrmax-tms->mbrmin+1,
694     sizeof(TMAP_TYP));
695     else
696     memset(tms->lumap, 0, (tms->mbrmax-tms->mbrmin+1)*sizeof(TMAP_TYP));
697    
698 greg 3.33 return(tms->lumap != NULL);
699     }
700    
701    
702     int
703 greg 3.16 tmErrorReturn( /* error return (with message) */
704 greg 3.17 const char *func,
705 greg 3.16 TMstruct *tms,
706     int err
707     )
708 greg 3.9 {
709 greg 3.17 if (tms != NULL) {
710     tms->lastFunc = func;
711     tms->lastError = err;
712     if (tms->flags & TM_F_NOSTDERR)
713     return(err);
714     }
715 greg 3.9 fputs(func, stderr);
716     fputs(": ", stderr);
717     fputs(tmErrorMessage[err], stderr);
718     fputs("!\n", stderr);
719     return(err);
720 greg 3.1 }