ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tonemap.c
Revision: 3.52
Committed: Wed May 26 17:50:26 2021 UTC (2 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.51: +3 -3 lines
Log Message:
fix: corner case where much of histogram gets truncated

File Contents

# User Rev Content
1 greg 3.1 #ifndef lint
2 greg 3.52 static const char RCSid[] = "$Id: tonemap.c,v 3.51 2021/05/04 21:50:54 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.51 tms->hbrmin = 10; tms->hbrmax = -10;
164 greg 3.16 tms->histo = NULL;
165 greg 3.1 }
166    
167    
168 greg 3.23 TMbright
169     tmCvLuminance( /* convert a single luminance */
170     double lum
171     )
172     {
173     double d;
174    
175     #ifdef isfinite
176     if (!isfinite(lum) || lum <= TM_NOLUM)
177     #else
178     if (lum <= TM_NOLUM)
179     #endif
180     return(TM_NOBRT);
181     d = TM_BRTSCALE*log(lum);
182 greg 3.35 return((TMbright)(d + .5 - (d < 0.)));
183 greg 3.23 }
184    
185    
186     int
187     tmCvLums( /* convert luminances using lookup */
188     TMbright *ls,
189     float *scan,
190     int len
191     )
192     {
193     if (tmFloat2BrtLUT == NULL) { /* initialize lookup table */
194     int32 i;
195     tmFloat2BrtLUT = (TMbright *)malloc(sizeof(TMbright)*0x10000);
196     if (tmFloat2BrtLUT == NULL)
197     return(TM_E_NOMEM);
198     for (i = 0; i < 0x10000; i++) {
199     int32 l = (i<<1 | 1) << 14;
200     #ifndef isfinite
201     if ((l & 0x7f800000) == 0x7f800000)
202     tmFloat2BrtLUT[i] = TM_NOBRT;
203     else
204     #endif
205     tmFloat2BrtLUT[i] = tmCvLuminance(*(float *)&l);
206     }
207     }
208     if (len <= 0)
209     return(TM_E_OK);
210     if ((ls == NULL) | (scan == NULL))
211     return(TM_E_ILLEGAL);
212     while (len--) {
213     if (*scan <= TM_NOLUM) {
214     *ls++ = TM_NOBRT;
215     ++scan;
216     continue;
217     }
218     *ls++ = tmCvLumLUfp(scan++);
219     }
220     return(TM_E_OK);
221     }
222    
223    
224     int
225     tmCvGrays( /* convert float gray values */
226     TMstruct *tms,
227     TMbright *ls,
228     float *scan,
229     int len
230     )
231     {
232     static const char funcName[] = "tmCvGrays";
233     int i;
234    
235     if (tms == NULL)
236     returnErr(TM_E_TMINVAL);
237     if ((ls == NULL) | (scan == NULL) | (len < 0))
238     returnErr(TM_E_ILLEGAL);
239 greg 3.24 if (tmFloat2BrtLUT == NULL) /* initialize */
240     tmCvLums(NULL, NULL, 0);
241 greg 3.23 for (i = len; i--; ) {
242     float lum = tms->inpsf * scan[i];
243     if (lum <= TM_NOLUM)
244     ls[i] = TM_NOBRT;
245     else
246     ls[i] = tmCvLumLUfp(&lum);
247     }
248     returnOK;
249     }
250    
251    
252 greg 3.1 int
253 greg 3.16 tmCvColors( /* convert float colors */
254     TMstruct *tms,
255     TMbright *ls,
256 greg 3.37 uby8 *cs,
257 greg 3.16 COLOR *scan,
258     int len
259     )
260 greg 3.1 {
261 greg 3.17 static const char funcName[] = "tmCvColors";
262 greg 3.37 static uby8 gamtab[1024];
263 greg 3.20 static double curgam = .0;
264 greg 3.1 COLOR cmon;
265 greg 3.23 float lum, slum, d;
266 greg 3.16 int i;
267 greg 3.1
268 greg 3.16 if (tms == NULL)
269 greg 3.1 returnErr(TM_E_TMINVAL);
270 schorsch 3.13 if ((ls == NULL) | (scan == NULL) | (len < 0))
271 greg 3.1 returnErr(TM_E_ILLEGAL);
272 greg 3.24 if (tmFloat2BrtLUT == NULL) /* initialize */
273     tmCvLums(NULL, NULL, 0);
274 greg 3.21 if (cs != TM_NOCHROM && fabs(tms->mongam - curgam) > .02) {
275 greg 3.20 curgam = tms->mongam; /* (re)build table */
276     for (i = 1024; i--; )
277     gamtab[i] = (int)(256.*pow((i+.5)/1024., 1./curgam));
278     }
279 greg 3.1 for (i = len; i--; ) {
280 greg 3.16 if (tmNeedMatrix(tms)) { /* get monitor RGB */
281     colortrans(cmon, tms->cmat, scan[i]);
282 greg 3.15 } else {
283 greg 3.16 cmon[RED] = tms->inpsf*scan[i][RED];
284     cmon[GRN] = tms->inpsf*scan[i][GRN];
285     cmon[BLU] = tms->inpsf*scan[i][BLU];
286 greg 3.1 }
287 greg 3.23 #ifdef isfinite
288     if (!isfinite(cmon[RED]) || cmon[RED] < .0f) cmon[RED] = .0f;
289     if (!isfinite(cmon[GRN]) || cmon[GRN] < .0f) cmon[GRN] = .0f;
290     if (!isfinite(cmon[BLU]) || cmon[BLU] < .0f) cmon[BLU] = .0f;
291     #else
292     if (cmon[RED] < .0f) cmon[RED] = .0f;
293     if (cmon[GRN] < .0f) cmon[GRN] = .0f;
294     if (cmon[BLU] < .0f) cmon[BLU] = .0f;
295 greg 3.22 #endif
296 greg 3.1 /* world luminance */
297 greg 3.16 lum = tms->clf[RED]*cmon[RED] +
298     tms->clf[GRN]*cmon[GRN] +
299     tms->clf[BLU]*cmon[BLU] ;
300 greg 3.25 if (lum <= TM_NOLUM) { /* convert brightness */
301     lum = cmon[RED] = cmon[GRN] = cmon[BLU] = TM_NOLUM;
302 greg 3.23 ls[i] = TM_NOBRT;
303 greg 3.25 } else
304 greg 3.23 ls[i] = tmCvLumLUfp(&lum);
305 greg 3.1 if (cs == TM_NOCHROM) /* no color? */
306     continue;
307 greg 3.16 if (tms->flags & TM_F_MESOPIC && lum < LMESUPPER) {
308 greg 3.1 slum = scotlum(cmon); /* mesopic adj. */
309 greg 3.25 if (lum < LMESLOWER) {
310 greg 3.1 cmon[RED] = cmon[GRN] = cmon[BLU] = slum;
311 greg 3.25 } else {
312 greg 3.1 d = (lum - LMESLOWER)/(LMESUPPER - LMESLOWER);
313 greg 3.16 if (tms->flags & TM_F_BW)
314 greg 3.1 cmon[RED] = cmon[GRN] =
315     cmon[BLU] = d*lum;
316     else
317     scalecolor(cmon, d);
318 greg 3.23 d = (1.f-d)*slum;
319 greg 3.1 cmon[RED] += d;
320     cmon[GRN] += d;
321     cmon[BLU] += d;
322     }
323 greg 3.16 } else if (tms->flags & TM_F_BW) {
324 greg 3.1 cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
325     }
326 greg 3.16 d = tms->clf[RED]*cmon[RED]/lum;
327 greg 3.23 cs[3*i ] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
328 greg 3.16 d = tms->clf[GRN]*cmon[GRN]/lum;
329 greg 3.23 cs[3*i+1] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
330 greg 3.16 d = tms->clf[BLU]*cmon[BLU]/lum;
331 greg 3.23 cs[3*i+2] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
332 greg 3.1 }
333     returnOK;
334     }
335    
336    
337 greg 3.9 int
338 greg 3.16 tmAddHisto( /* add values to histogram */
339     TMstruct *tms,
340     TMbright *ls,
341     int len,
342     int wt
343     )
344 greg 3.1 {
345 greg 3.17 static const char funcName[] = "tmAddHisto";
346 greg 3.9 int oldorig=0, oldlen, horig, hlen;
347 greg 3.51 int i;
348 greg 3.1
349 greg 3.16 if (tms == NULL)
350 greg 3.1 returnErr(TM_E_TMINVAL);
351 gwlarson 3.8 if (len < 0)
352     returnErr(TM_E_ILLEGAL);
353     if (len == 0)
354     returnOK;
355 greg 3.1 /* first, grow limits */
356 greg 3.16 if (tms->histo == NULL) {
357 greg 3.1 for (i = len; i-- && ls[i] < MINBRT; )
358     ;
359     if (i < 0)
360     returnOK;
361 greg 3.16 tms->hbrmin = tms->hbrmax = ls[i];
362 greg 3.1 oldlen = 0;
363     } else {
364 greg 3.31 oldorig = HISTI(tms->hbrmin);
365     oldlen = HISTI(tms->hbrmax) + 1 - oldorig;
366 greg 3.1 }
367     for (i = len; i--; ) {
368 greg 3.51 if (ls[i] < MINBRT)
369 greg 3.1 continue;
370 greg 3.51 if (ls[i] < tms->hbrmin)
371     tms->hbrmin = ls[i];
372     else if (ls[i] > tms->hbrmax)
373     tms->hbrmax = ls[i];
374 greg 3.1 }
375 greg 3.31 horig = HISTI(tms->hbrmin);
376     hlen = HISTI(tms->hbrmax) + 1 - horig;
377 greg 3.1 if (hlen > oldlen) { /* (re)allocate histogram */
378 greg 3.38 HIST_TYP *newhist = (HIST_TYP *)calloc(hlen, sizeof(HIST_TYP));
379 greg 3.1 if (newhist == NULL)
380     returnErr(TM_E_NOMEM);
381     if (oldlen) { /* copy and free old */
382 greg 3.46 memcpy(newhist+(oldorig-horig),
383     tms->histo, oldlen*sizeof(HIST_TYP));
384     free(tms->histo);
385 greg 3.1 }
386 greg 3.16 tms->histo = newhist;
387 greg 3.1 }
388     if (wt == 0)
389     returnOK;
390     for (i = len; i--; ) /* add in new counts */
391     if (ls[i] >= MINBRT)
392 greg 3.31 tms->histo[ HISTI(ls[i]) - horig ] += wt;
393 greg 3.1 returnOK;
394     }
395    
396    
397     static double
398 greg 3.16 htcontrs( /* human threshold contrast sensitivity, dL(La) */
399     double La
400     )
401 greg 3.1 {
402     double l10La, l10dL;
403     /* formula taken from Ferwerda et al. [SG96] */
404     if (La < 1.148e-4)
405     return(1.38e-3);
406     l10La = log10(La);
407     if (l10La < -1.44) /* rod response regime */
408     l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
409     else if (l10La < -.0184)
410     l10dL = l10La - .395;
411     else if (l10La < 1.9) /* cone response regime */
412     l10dL = pow(.249*l10La + .65, 2.7) - .72;
413     else
414     l10dL = l10La - 1.255;
415    
416     return(exp10(l10dL));
417     }
418    
419    
420 greg 3.9 int
421 greg 3.16 tmFixedMapping( /* compute fixed, linear tone-mapping */
422     TMstruct *tms,
423     double expmult,
424 greg 3.46 double gamval,
425     double Lddyn
426 greg 3.16 )
427 greg 3.9 {
428 greg 3.17 static const char funcName[] = "tmFixedMapping";
429 greg 3.47 const int maxV = (1L<<(8*sizeof(TMAP_TYP))) - 1;
430     double minD;
431     int i;
432 greg 3.9
433 greg 3.16 if (!tmNewMap(tms))
434 greg 3.9 returnErr(TM_E_NOMEM);
435 greg 3.46 /* check arguments */
436     if (expmult <= .0) expmult = 1.;
437     if (gamval < MINGAM) gamval = tms->mongam;
438     if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
439     minD = 1./Lddyn;
440 greg 3.28 for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
441 greg 3.46 double d;
442     d = expmult/tms->inpsf * tmLuminance(tms->mbrmin + i);
443 greg 3.48 if (d >= 2.*minD)
444     d -= minD;
445     else /* soft black crushing */
446     d *= d/(4.*minD);
447     d /= 1. - minD;
448 greg 3.46 d = TM_BRES*pow(d, 1./gamval);
449 greg 3.47 tms->lumap[i] = (d > maxV) ? maxV : (int)d;
450 greg 3.28 }
451 greg 3.9 returnOK;
452     }
453    
454    
455 greg 3.1 int
456 greg 3.16 tmComputeMapping( /* compute histogram tone-mapping */
457     TMstruct *tms,
458     double gamval,
459     double Lddyn,
460     double Ldmax
461     )
462 greg 3.1 {
463 greg 3.17 static const char funcName[] = "tmComputeMapping";
464 greg 3.38 HIST_TYP *histo;
465 greg 3.1 float *cumf;
466 greg 3.38 int brt0, histlen;
467 greg 3.39 HIST_TYP threshold, ceiling, trimmings, histot;
468 greg 3.44 double logLddyn, Ldmin, Lwavg, Tr, Lw, Ld;
469 greg 3.9 double sum;
470 greg 3.16 double d;
471     int i, j;
472 greg 3.1
473 greg 3.16 if (tms == NULL || tms->histo == NULL)
474 greg 3.1 returnErr(TM_E_TMINVAL);
475     /* check arguments */
476     if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
477     if (Ldmax < MINLDMAX) Ldmax = DEFLDMAX;
478 greg 3.16 if (gamval < MINGAM) gamval = tms->mongam;
479 greg 3.1 /* compute handy values */
480     Ldmin = Ldmax/Lddyn;
481     logLddyn = log(Lddyn);
482 greg 3.31 i = HISTI(tms->hbrmin);
483 greg 3.32 brt0 = HISTV(i);
484 greg 3.31 histlen = HISTI(tms->hbrmax) + 1 - i;
485 greg 3.1 /* histogram total and mean */
486     histot = 0; sum = 0;
487     j = brt0 + histlen*HISTEP;
488     for (i = histlen; i--; ) {
489 greg 3.16 histot += tms->histo[i];
490 greg 3.30 sum += (double)(j -= HISTEP) * tms->histo[i];
491 greg 3.1 }
492 greg 3.29 if (!histot)
493 greg 3.1 returnErr(TM_E_TMFAIL);
494 greg 3.52 threshold = histot/500 + 1;
495 greg 3.1 Lwavg = tmLuminance( (double)sum / histot );
496 greg 3.29 /* use linear tone mapping? */
497 greg 3.52 if (tms->flags & TM_F_LINEAR ||
498 greg 3.30 tms->hbrmax - tms->hbrmin < TM_BRTSCALE*logLddyn)
499 greg 3.29 goto linearmap;
500 greg 3.9 /* clamp histogram */
501 greg 3.38 histo = (HIST_TYP *)malloc(histlen*sizeof(HIST_TYP));
502 greg 3.9 cumf = (float *)malloc((histlen+2)*sizeof(float));
503 schorsch 3.13 if ((histo == NULL) | (cumf == NULL))
504 greg 3.9 returnErr(TM_E_NOMEM);
505     cumf[histlen+1] = 1.; /* guard for assignment code */
506 greg 3.42 /* make malleable copy */
507 greg 3.46 memcpy(histo, tms->histo, histlen*sizeof(HIST_TYP));
508 greg 3.9 do { /* iterate to solution */
509     sum = 0; /* cumulative probability */
510     for (i = 0; i < histlen; i++) {
511 greg 3.45 cumf[i] = sum/histot;
512     sum += (double)histo[i];
513 greg 3.9 }
514     cumf[histlen] = 1.;
515 greg 3.16 Tr = histot * (double)(tms->hbrmax - tms->hbrmin) /
516 greg 3.45 ((double)TM_BRTSCALE*histlen*logLddyn);
517 greg 3.9 ceiling = Tr + 1.;
518     trimmings = 0; /* clip to envelope */
519     for (i = histlen; i--; ) {
520 greg 3.16 if (tms->flags & TM_F_HCONTR) {
521 greg 3.9 Lw = tmLuminance(brt0 + i*HISTEP);
522     Ld = Ldmin * exp( logLddyn *
523     .5*(cumf[i]+cumf[i+1]) );
524     ceiling = Tr * (htcontrs(Ld) * Lw) /
525     (htcontrs(Lw) * Ld) + 1.;
526 greg 3.2 }
527 greg 3.9 if (histo[i] > ceiling) {
528     trimmings += histo[i] - ceiling;
529     histo[i] = ceiling;
530 greg 3.1 }
531     }
532 greg 3.9 /* check if we're out of data */
533     if ((histot -= trimmings) <= threshold) {
534 greg 3.46 free(histo);
535     free(cumf);
536 greg 3.9 goto linearmap;
537     }
538 greg 3.50 } while (40*trimmings > histot);
539 greg 3.30 /* allocate space for mapping */
540     if (!tmNewMap(tms))
541     returnErr(TM_E_NOMEM);
542 greg 3.9 /* assign tone-mapping */
543 greg 3.16 for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
544     j = d = (double)i/(tms->mbrmax-tms->mbrmin)*histlen;
545 greg 3.9 d -= (double)j;
546     Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
547     d = (Ld - Ldmin)/(Ldmax - Ldmin);
548 greg 3.40 tms->lumap[i] = TM_BRES*pow(d, 1./gamval);
549 greg 3.1 }
550 greg 3.46 free(histo); /* clean up and return */
551     free(cumf);
552 greg 3.1 returnOK;
553 greg 3.9 linearmap: /* linear tone-mapping */
554 greg 3.16 if (tms->flags & TM_F_HCONTR)
555 greg 3.44 d = htcontrs(sqrt(Ldmax*Ldmin)) / htcontrs(Lwavg);
556 greg 3.9 else
557 greg 3.44 d = Ldmax / tmLuminance(tms->hbrmax);
558 greg 3.46 return(tmFixedMapping(tms, tms->inpsf*d/Ldmax, gamval, Lddyn));
559 greg 3.1 }
560    
561    
562     int
563 greg 3.16 tmMapPixels( /* apply tone-mapping to pixel(s) */
564     TMstruct *tms,
565 greg 3.37 uby8 *ps,
566 greg 3.16 TMbright *ls,
567 greg 3.37 uby8 *cs,
568 greg 3.16 int len
569     )
570 greg 3.1 {
571 greg 3.17 static const char funcName[] = "tmMapPixels";
572 greg 3.41 TMbright lv;
573     TMAP_TYP li;
574     int pv;
575 greg 3.1
576 greg 3.16 if (tms == NULL || tms->lumap == NULL)
577 greg 3.1 returnErr(TM_E_TMINVAL);
578 schorsch 3.13 if ((ps == NULL) | (ls == NULL) | (len < 0))
579 greg 3.1 returnErr(TM_E_ILLEGAL);
580     while (len--) {
581 greg 3.41 if ((lv = *ls++) < tms->mbrmin) {
582 greg 3.9 li = 0;
583     } else {
584 greg 3.41 if (lv > tms->mbrmax)
585     lv = tms->mbrmax;
586     li = tms->lumap[lv - tms->mbrmin];
587 greg 3.9 }
588 greg 3.43 if (cs == TM_NOCHROM) {
589     #if !(TM_BRES & 0xff)
590     *ps++ = li>=TM_BRES ? 255 : li/(TM_BRES>>8);
591     #else
592     *ps++ = li>=TM_BRES ? 255 : (li<<8)/TM_BRES;
593     #endif
594     } else {
595 greg 3.16 pv = *cs++ * li / tms->cdiv[RED];
596 greg 3.1 *ps++ = pv>255 ? 255 : pv;
597 greg 3.16 pv = *cs++ * li / tms->cdiv[GRN];
598 greg 3.1 *ps++ = pv>255 ? 255 : pv;
599 greg 3.16 pv = *cs++ * li / tms->cdiv[BLU];
600 greg 3.1 *ps++ = pv>255 ? 255 : pv;
601     }
602     }
603     returnOK;
604     }
605    
606    
607 greg 3.16 TMstruct *
608 greg 3.42 tmDup( /* duplicate tone mapping */
609 greg 3.16 TMstruct *tms
610     )
611 greg 3.3 {
612     int len;
613 greg 3.16 int i;
614     TMstruct *tmnew;
615 greg 3.3
616 greg 3.16 if (tms == NULL) /* anything to duplicate? */
617 greg 3.3 return(NULL);
618 greg 3.16 tmnew = (TMstruct *)malloc(sizeof(TMstruct));
619 greg 3.3 if (tmnew == NULL)
620     return(NULL);
621 greg 3.16 *tmnew = *tms; /* copy everything */
622 greg 3.3 if (tmnew->histo != NULL) { /* duplicate histogram */
623 greg 3.31 len = HISTI(tmnew->hbrmax) + 1 - HISTI(tmnew->hbrmin);
624 greg 3.38 tmnew->histo = (HIST_TYP *)malloc(len*sizeof(HIST_TYP));
625 greg 3.3 if (tmnew->histo != NULL)
626 greg 3.46 memcpy(tmnew->histo, tms->histo, len*sizeof(HIST_TYP));
627 greg 3.3 }
628     if (tmnew->lumap != NULL) { /* duplicate luminance mapping */
629 gregl 3.6 len = tmnew->mbrmax-tmnew->mbrmin+1;
630 greg 3.38 tmnew->lumap = (TMAP_TYP *)malloc(len*sizeof(TMAP_TYP));
631 greg 3.3 if (tmnew->lumap != NULL)
632 greg 3.46 memcpy(tmnew->lumap, tms->lumap, len*sizeof(TMAP_TYP));
633 greg 3.3 }
634 greg 3.4 /* clear package data */
635     for (i = tmNumPkgs; i--; )
636     tmnew->pd[i] = NULL;
637 greg 3.16 /* return copy */
638     return(tmnew);
639 greg 3.1 }
640    
641    
642     void
643     tmDone(tms) /* done with tone mapping -- destroy it */
644 greg 3.16 TMstruct *tms;
645 greg 3.1 {
646 greg 3.16 int i;
647     /* NULL arg. is equiv. to tms */
648     if (tms == NULL)
649 greg 3.1 return;
650     /* free tables */
651     if (tms->histo != NULL)
652 greg 3.46 free(tms->histo);
653 greg 3.1 if (tms->lumap != NULL)
654 greg 3.46 free(tms->lumap);
655 greg 3.4 /* free private data */
656     for (i = tmNumPkgs; i--; )
657     if (tms->pd[i] != NULL)
658     (*tmPkg[i]->Free)(tms->pd[i]);
659 greg 3.46 free(tms); /* free basic structure */
660 greg 3.9 }
661    
662     /******************** Shared but Private library routines *********************/
663    
664 greg 3.37 uby8 tmMesofact[BMESUPPER-BMESLOWER];
665 greg 3.9
666     void
667     tmMkMesofact() /* build mesopic lookup factor table */
668     {
669 greg 3.16 int i;
670 greg 3.9
671     if (tmMesofact[BMESUPPER-BMESLOWER-1])
672     return;
673    
674     for (i = BMESLOWER; i < BMESUPPER; i++)
675     tmMesofact[i-BMESLOWER] = 256. *
676     (tmLuminance(i) - LMESLOWER) /
677     (LMESUPPER - LMESLOWER);
678     }
679    
680    
681     int
682 greg 3.33 tmNewMap( /* allocate new tone-mapping array */
683     TMstruct *tms
684     )
685     {
686     if (tms->lumap != NULL && (tms->mbrmax - tms->mbrmin) !=
687     (tms->hbrmax - tms->hbrmin)) {
688 greg 3.46 free(tms->lumap);
689 greg 3.33 tms->lumap = NULL;
690     }
691     tms->mbrmin = tms->hbrmin;
692     tms->mbrmax = tms->hbrmax;
693     if (tms->mbrmin > tms->mbrmax)
694 greg 3.47 return(0);
695 greg 3.33 if (tms->lumap == NULL)
696 greg 3.46 tms->lumap = (TMAP_TYP *)calloc(tms->mbrmax-tms->mbrmin+1,
697     sizeof(TMAP_TYP));
698     else
699     memset(tms->lumap, 0, (tms->mbrmax-tms->mbrmin+1)*sizeof(TMAP_TYP));
700    
701 greg 3.33 return(tms->lumap != NULL);
702     }
703    
704    
705     int
706 greg 3.16 tmErrorReturn( /* error return (with message) */
707 greg 3.17 const char *func,
708 greg 3.16 TMstruct *tms,
709     int err
710     )
711 greg 3.9 {
712 greg 3.17 if (tms != NULL) {
713     tms->lastFunc = func;
714     tms->lastError = err;
715     if (tms->flags & TM_F_NOSTDERR)
716     return(err);
717     }
718 greg 3.9 fputs(func, stderr);
719     fputs(": ", stderr);
720     fputs(tmErrorMessage[err], stderr);
721     fputs("!\n", stderr);
722     return(err);
723 greg 3.1 }