ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/tmapcolrs.c
Revision: 3.35
Committed: Tue Oct 13 00:08:46 2020 UTC (3 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.34: +12 -4 lines
Log Message:
fix: Corrected overflow of some color primaries during tone-mapping

File Contents

# User Rev Content
1 greg 3.1 #ifndef lint
2 greg 3.35 static const char RCSid[] = "$Id: tmapcolrs.c,v 3.34 2019/12/28 18:05:14 greg Exp $";
3 greg 3.1 #endif
4     /*
5     * Routines for tone mapping on Radiance RGBE and XYZE pictures.
6 greg 3.10 *
7     * Externals declared in tonemap.h
8     */
9    
10 greg 3.11 #include "copyright.h"
11 greg 3.1
12 greg 3.26 #include <stdlib.h>
13 greg 3.1 #include <math.h>
14 schorsch 3.16
15 greg 3.21 #ifdef PCOND
16 greg 3.31 #include "paths.h"
17 greg 3.10 #endif
18 greg 3.27 #include "tmprivat.h"
19 greg 3.33 #include "rtio.h"
20 greg 3.27 #include "resolu.h"
21 greg 3.1
22 greg 3.28 #define GAMTSZ 4096
23 greg 3.3
24     typedef struct {
25 greg 3.29 uby8 gamb[GAMTSZ]; /* gamma lookup table */
26 greg 3.21 int clfb[3]; /* encoded tm->clf */
27     int32 cmatb[3][3]; /* encoded color transform */
28 greg 3.3 TMbright inpsfb; /* encoded tm->inpsf */
29     } COLRDATA;
30    
31 greg 3.17 static MEM_PTR colrInit(TMstruct *);
32     static void colrNewSpace(TMstruct *);
33 schorsch 3.15 static gethfunc headline;
34    
35 greg 3.3 static struct tmPackage colrPkg = { /* our package functions */
36     colrInit, colrNewSpace, free
37     };
38     static int colrReg = -1; /* our package registration number */
39    
40 greg 3.28 static TMbright logi[256];
41 greg 3.1
42    
43     int
44 greg 3.17 tmCvColrs( /* convert RGBE/XYZE colors */
45     TMstruct *tms,
46     TMbright *ls,
47 greg 3.29 uby8 *cs,
48 greg 3.17 COLR *scan,
49     int len
50     )
51 greg 3.1 {
52 greg 3.18 static const char funcName[] = "tmCvColrs";
53 greg 3.21 int cmon[4];
54 greg 3.28 COLRDATA *cd;
55     int i, j, bi;
56     int32 li, vl;
57 greg 3.1
58 greg 3.17 if (tms == NULL)
59 greg 3.1 returnErr(TM_E_TMINVAL);
60 schorsch 3.14 if ((ls == NULL) | (scan == NULL) | (len < 0))
61 greg 3.1 returnErr(TM_E_ILLEGAL);
62 gwlarson 3.6 if (colrReg < 0) { /* build tables if necessary */
63 greg 3.3 colrReg = tmRegPkg(&colrPkg);
64     if (colrReg < 0)
65     returnErr(TM_E_CODERR1);
66 greg 3.1 for (i = 256; i--; )
67     logi[i] = TM_BRTSCALE*log((i+.5)/256.) - .5;
68 greg 3.10 tmMkMesofact();
69 greg 3.1 }
70 greg 3.17 if ((cd = (COLRDATA *)tmPkgData(tms,colrReg)) == NULL)
71 greg 3.3 returnErr(TM_E_NOMEM);
72 greg 3.1 for (i = len; i--; ) {
73 greg 3.21 if (tmNeedMatrix(tms)) { /* apply color xform */
74 greg 3.35 bi = 0;
75 greg 3.21 for (j = 3; j--; ) {
76     vl = cd->cmatb[j][RED]*(int32)scan[i][RED] +
77     cd->cmatb[j][GRN]*(int32)scan[i][GRN] +
78     cd->cmatb[j][BLU]*(int32)scan[i][BLU] ;
79 greg 3.35 if (vl < 0)
80     cmon[j] = vl/(int32)0x10000;
81     else if ((cmon[j] = vl>>16) > bi)
82     bi = cmon[j];
83 greg 3.21 }
84     cmon[EXP] = scan[i][EXP];
85 greg 3.35 while (bi >= 256) { /* handle overflow */
86     cmon[EXP]++;
87     for (j = 3; j--; ) cmon[j] >>= 1;
88     bi >>= 1;
89     }
90 greg 3.21 } else
91     copycolr(cmon, scan[i]);
92 greg 3.1 /* world luminance */
93 greg 3.28 li = cd->clfb[RED]*(int32)cmon[RED] +
94     cd->clfb[GRN]*(int32)cmon[GRN] +
95     cd->clfb[BLU]*(int32)cmon[BLU] ;
96     if (li >= 1L<<(12+8)) li = 255;
97 greg 3.24 else li >>= 12;
98 greg 3.22 bi = BRT2SCALE(cmon[EXP]-COLXS) + cd->inpsfb;
99     if (li > 0)
100     bi += logi[li];
101     else {
102     bi += logi[0];
103     li = 1; /* avoid /0 */
104 greg 3.1 }
105     ls[i] = bi;
106     if (cs == TM_NOCHROM) /* no color? */
107     continue;
108     /* mesopic adj. */
109 greg 3.17 if (tms->flags & TM_F_MESOPIC && bi < BMESUPPER) {
110 greg 3.21 int pf, sli = normscot(cmon);
111     if (bi < BMESLOWER) {
112 greg 3.1 cmon[RED] = cmon[GRN] = cmon[BLU] = sli;
113 greg 3.21 } else {
114 greg 3.17 if (tms->flags & TM_F_BW)
115 gwlarson 3.8 cmon[RED] = cmon[GRN] = cmon[BLU] = li;
116 greg 3.10 pf = tmMesofact[bi-BMESLOWER];
117 greg 3.1 sli *= 256 - pf;
118 greg 3.21 for (j = 3; j--; ) {
119     cmon[j] = sli + pf*cmon[j];
120     if (cmon[j] <= 0) cmon[j] = 0;
121     else cmon[j] >>= 8;
122     }
123 greg 3.1 }
124 greg 3.17 } else if (tms->flags & TM_F_BW) {
125 gwlarson 3.8 cmon[RED] = cmon[GRN] = cmon[BLU] = li;
126 greg 3.21 } else {
127     for (j = 3; j--; )
128 greg 3.35 cmon[j] *= (cmon[j] > 0);
129 greg 3.1 }
130 greg 3.28 bi = ( (uint32)GAMTSZ*cd->clfb[RED]*cmon[RED]/li ) >> 12;
131 greg 3.3 cs[3*i ] = bi>=GAMTSZ ? 255 : cd->gamb[bi];
132 greg 3.28 bi = ( (uint32)GAMTSZ*cd->clfb[GRN]*cmon[GRN]/li ) >> 12;
133 greg 3.3 cs[3*i+1] = bi>=GAMTSZ ? 255 : cd->gamb[bi];
134 greg 3.28 bi = ( (uint32)GAMTSZ*cd->clfb[BLU]*cmon[BLU]/li ) >> 12;
135 greg 3.3 cs[3*i+2] = bi>=GAMTSZ ? 255 : cd->gamb[bi];
136 greg 3.1 }
137     returnOK;
138     }
139    
140    
141     #define FMTRGB 1 /* Input is RGBE */
142     #define FMTCIE 2 /* Input is CIE XYZE */
143     #define FMTUNK 3 /* Input format is unspecified */
144     #define FMTBAD 4 /* Input is not a recognized format */
145    
146     static struct radhead {
147     int format; /* FMTRGB, FMTCIE, FMTUNK, FMTBAD */
148     double expos; /* input exposure value */
149     RGBPRIMP primp; /* input primaries */
150     RGBPRIMS mypri; /* custom primaries */
151     } rhdefault = {FMTUNK, 1., stdprims, STDPRIMS};
152    
153    
154     static int
155 schorsch 3.15 headline( /* grok a header line */
156 greg 3.28 char *s,
157 schorsch 3.15 void *vrh
158     )
159 greg 3.1 {
160 greg 3.32 char fmt[MAXFMTLEN];
161 greg 3.28 struct radhead *rh = vrh;
162 greg 3.1
163     if (formatval(fmt, s)) {
164     if (!strcmp(fmt, COLRFMT))
165     rh->format = FMTRGB;
166     else if (!strcmp(fmt, CIEFMT))
167     rh->format = FMTCIE;
168     else
169     rh->format = FMTBAD;
170 gwlarson 3.7 return(0);
171 greg 3.1 }
172     if (isexpos(s)) {
173     rh->expos *= exposval(s);
174 gwlarson 3.7 return(0);
175 greg 3.1 }
176     if (isprims(s)) {
177     primsval(rh->mypri, s);
178     rh->primp = rh->mypri;
179 gwlarson 3.7 return(0);
180 greg 3.1 }
181 gwlarson 3.7 return(0);
182 greg 3.1 }
183    
184    
185     int
186 greg 3.17 tmLoadPicture( /* convert Radiance picture */
187     TMstruct *tms,
188     TMbright **lpp,
189 greg 3.29 uby8 **cpp,
190 greg 3.17 int *xp,
191     int *yp,
192     char *fname,
193     FILE *fp
194     )
195 greg 3.1 {
196     char *funcName = fname==NULL ? "tmLoadPicture" : fname;
197     FILE *inpf;
198     struct radhead info;
199     int err;
200     COLR *scanin = NULL;
201     int i;
202     /* check arguments */
203 greg 3.17 if (tms == NULL)
204 greg 3.1 returnErr(TM_E_TMINVAL);
205 schorsch 3.14 if ((lpp == NULL) | (xp == NULL) | (yp == NULL) |
206     ((fname == NULL) & (fp == TM_GETFILE)))
207 greg 3.1 returnErr(TM_E_ILLEGAL);
208     *xp = *yp = 0; /* error precaution */
209 greg 3.20 if ((inpf = fp) == TM_GETFILE && (inpf = fopen(fname, "rb")) == NULL)
210 greg 3.1 returnErr(TM_E_BADFILE);
211 gwlarson 3.6 *lpp = NULL;
212     if (cpp != TM_NOCHROMP) *cpp = NULL;
213 greg 3.1 info = rhdefault; /* get our header */
214 schorsch 3.15 getheader(inpf, headline, &info);
215 schorsch 3.14 if ((info.format == FMTBAD) | (info.expos <= 0.) ||
216 greg 3.1 fgetresolu(xp, yp, inpf) < 0) {
217     err = TM_E_BADFILE; goto done;
218     }
219     if (info.format == FMTUNK) /* assume RGBE format */
220     info.format = FMTRGB;
221     if (info.format == FMTRGB)
222     info.expos /= WHTEFFICACY;
223     else if (info.format == FMTCIE)
224     info.primp = TM_XYZPRIM;
225     /* prepare library */
226 greg 3.19 if ((err = tmSetSpace(tms, info.primp, 1./info.expos, NULL)) != TM_E_OK)
227 greg 3.1 goto done;
228     err = TM_E_NOMEM; /* allocate arrays */
229     *lpp = (TMbright *)malloc(sizeof(TMbright) * *xp * *yp);
230     if (*lpp == NULL)
231     goto done;
232     if (cpp != TM_NOCHROMP) {
233 greg 3.29 *cpp = (uby8 *)malloc(3*sizeof(uby8) * *xp * *yp);
234 greg 3.1 if (*cpp == NULL)
235     goto done;
236     }
237     scanin = (COLR *)malloc(sizeof(COLR) * *xp);
238     if (scanin == NULL)
239     goto done;
240     err = TM_E_BADFILE; /* read & convert scanlines */
241     for (i = 0; i < *yp; i++) {
242     if (freadcolrs(scanin, *xp, inpf) < 0) {
243     err = TM_E_BADFILE; break;
244     }
245 greg 3.17 err = tmCvColrs(tms, *lpp + (i * *xp),
246 greg 3.1 cpp==TM_NOCHROMP ? TM_NOCHROM : *cpp + (i * 3 * *xp),
247     scanin, *xp);
248     if (err != TM_E_OK)
249     break;
250     }
251     done: /* clean up */
252     if (fp == NULL)
253     fclose(inpf);
254     if (scanin != NULL)
255 gwlarson 3.6 free((MEM_PTR)scanin);
256     if (err != TM_E_OK) {
257     if (*lpp != NULL)
258     free((MEM_PTR)*lpp);
259     if (cpp != TM_NOCHROMP && *cpp != NULL)
260     free((MEM_PTR)*cpp);
261 greg 3.1 returnErr(err);
262 gwlarson 3.6 }
263 greg 3.1 returnOK;
264     }
265    
266    
267 gwlarson 3.9 #ifdef PCOND
268 schorsch 3.13 static int /* run pcond to map picture */
269 greg 3.1 dopcond(psp, xp, yp, flags, monpri, gamval, Lddyn, Ldmax, fname)
270 greg 3.29 uby8 **psp;
271 greg 3.1 int *xp, *yp;
272     int flags;
273     RGBPRIMP monpri;
274     double gamval, Lddyn, Ldmax;
275     char *fname;
276     {
277     char *funcName = fname;
278 greg 3.17 TMstruct *tms = NULL;
279 greg 3.10 char cmdbuf[1024];
280 greg 3.1 FILE *infp;
281 greg 3.28 COLR *scan;
282 greg 3.29 uby8 *rp;
283 greg 3.1 int y;
284 greg 3.28 int x;
285 greg 3.1 /* set up gamma correction */
286     if (setcolrcor(pow, 1./gamval) < 0)
287     returnErr(TM_E_NOMEM);
288     /* create command */
289 gwlarson 3.9 strcpy(cmdbuf, PCOND);
290 greg 3.1 if (flags & TM_F_HCONTR)
291 gwlarson 3.9 strcat(cmdbuf, " -s");
292 greg 3.1 if (flags & TM_F_MESOPIC)
293 gwlarson 3.9 strcat(cmdbuf, " -c");
294 greg 3.1 if (flags & TM_F_LINEAR)
295 gwlarson 3.9 strcat(cmdbuf, " -l");
296 greg 3.1 if (flags & TM_F_ACUITY)
297 gwlarson 3.9 strcat(cmdbuf, " -a");
298 greg 3.1 if (flags & TM_F_VEIL)
299 gwlarson 3.9 strcat(cmdbuf, " -v");
300 greg 3.1 if (flags & TM_F_CWEIGHT)
301 gwlarson 3.9 strcat(cmdbuf, " -w");
302     if (monpri != stdprims)
303     sprintf(cmdbuf+strlen(cmdbuf), " -p %f %f %f %f %f %f %f %f",
304     monpri[RED][CIEX], monpri[RED][CIEY],
305     monpri[GRN][CIEX], monpri[GRN][CIEY],
306     monpri[BLU][CIEX], monpri[BLU][CIEY],
307     monpri[WHT][CIEX], monpri[WHT][CIEY]);
308     sprintf(cmdbuf+strlen(cmdbuf), " -d %f -u %f %s", Lddyn, Ldmax, fname);
309 greg 3.1 /* start pcond */
310     if ((infp = popen(cmdbuf, "r")) == NULL)
311     returnErr(TM_E_BADFILE);
312     /* check picture format and size */
313     if (checkheader(infp, COLRFMT, NULL) < 0 ||
314     fgetresolu(xp, yp, infp) < 0) {
315     pclose(infp);
316     returnErr(TM_E_BADFILE);
317     }
318     /* allocate arrays */
319 greg 3.2 scan = (COLR *)malloc(sizeof(COLR) * *xp);
320 greg 3.1 if (flags & TM_F_BW)
321 greg 3.29 rp = (uby8 *)malloc(sizeof(uby8) * *xp * *yp);
322 greg 3.1 else
323 greg 3.29 rp = (uby8 *)malloc(3*sizeof(uby8) * *xp * *yp);
324 schorsch 3.14 if (((*psp = rp) == NULL) | (scan == NULL)) {
325 greg 3.1 pclose(infp);
326     returnErr(TM_E_NOMEM);
327     }
328     /* read and gamma map file */
329     for (y = 0; y < *yp; y++) {
330     if (freadcolrs(scan, *xp, infp) < 0) {
331     pclose(infp);
332 gwlarson 3.6 free((MEM_PTR)scan);
333     free((MEM_PTR)*psp);
334 greg 3.2 *psp = NULL;
335 greg 3.1 returnErr(TM_E_BADFILE);
336     }
337     colrs_gambs(scan, *xp);
338     if (flags & TM_F_BW)
339     for (x = 0; x < *xp; x++)
340     *rp++ = normbright(scan[x]);
341     else
342     for (x = 0; x < *xp; x++) {
343     *rp++ = scan[x][RED];
344     *rp++ = scan[x][GRN];
345     *rp++ = scan[x][BLU];
346     }
347     }
348 gwlarson 3.6 free((MEM_PTR)scan);
349 greg 3.1 pclose(infp);
350     returnOK;
351     }
352 gwlarson 3.9 #endif
353 greg 3.1
354    
355     int /* map a Radiance picture */
356     tmMapPicture(psp, xp, yp, flags, monpri, gamval, Lddyn, Ldmax, fname, fp)
357 greg 3.29 uby8 **psp;
358 greg 3.1 int *xp, *yp;
359     int flags;
360     RGBPRIMP monpri;
361     double gamval, Lddyn, Ldmax;
362     char *fname;
363     FILE *fp;
364     {
365     char *funcName = fname==NULL ? "tmMapPicture" : fname;
366 greg 3.30 TMstruct *tms = NULL;
367 greg 3.29 uby8 *cp;
368 greg 3.1 TMbright *lp;
369     int err;
370     /* check arguments */
371 schorsch 3.14 if ((psp == NULL) | (xp == NULL) | (yp == NULL) | (monpri == NULL) |
372     ((fname == NULL) & (fp == TM_GETFILE)))
373 greg 3.1 returnErr(TM_E_ILLEGAL);
374     /* set defaults */
375     if (gamval < MINGAM) gamval = DEFGAM;
376     if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
377     if (Ldmax < MINLDMAX) Ldmax = DEFLDMAX;
378     if (flags & TM_F_BW) monpri = stdprims;
379 gwlarson 3.9 #ifdef PCOND
380 greg 3.1 /* check for pcond run */
381     if (fp == TM_GETFILE && flags & TM_F_UNIMPL)
382     return( dopcond(psp, xp, yp, flags,
383     monpri, gamval, Lddyn, Ldmax, fname) );
384 gwlarson 3.9 #endif
385 greg 3.1 /* initialize tone mapping */
386 greg 3.17 if ((tms = tmInit(flags, monpri, gamval)) == NULL)
387 greg 3.1 returnErr(TM_E_NOMEM);
388     /* load & convert picture */
389 greg 3.17 err = tmLoadPicture(tms, &lp, (flags&TM_F_BW) ? TM_NOCHROMP : &cp,
390 greg 3.1 xp, yp, fname, fp);
391     if (err != TM_E_OK) {
392 greg 3.17 tmDone(tms);
393 greg 3.1 return(err);
394     }
395     /* allocate space for result */
396     if (flags & TM_F_BW) {
397 greg 3.29 *psp = (uby8 *)malloc(sizeof(uby8) * *xp * *yp);
398 greg 3.2 if (*psp == NULL) {
399 gwlarson 3.6 free((MEM_PTR)lp);
400 greg 3.17 tmDone(tms);
401 greg 3.1 returnErr(TM_E_NOMEM);
402 greg 3.2 }
403 greg 3.1 cp = TM_NOCHROM;
404     } else
405     *psp = cp;
406     /* compute color mapping */
407 greg 3.17 err = tmAddHisto(tms, lp, *xp * *yp, 1);
408 greg 3.1 if (err != TM_E_OK)
409     goto done;
410 greg 3.17 err = tmComputeMapping(tms, gamval, Lddyn, Ldmax);
411 greg 3.1 if (err != TM_E_OK)
412     goto done;
413     /* map colors */
414 greg 3.17 err = tmMapPixels(tms, *psp, lp, cp, *xp * *yp);
415 greg 3.1
416     done: /* clean up */
417 gwlarson 3.6 free((MEM_PTR)lp);
418 greg 3.17 tmDone(tms);
419 greg 3.1 if (err != TM_E_OK) { /* free memory on error */
420 gwlarson 3.6 free((MEM_PTR)*psp);
421 greg 3.1 *psp = NULL;
422     returnErr(err);
423     }
424     returnOK;
425 greg 3.3 }
426    
427    
428     static void
429     colrNewSpace(tms) /* color space changed for tone mapping */
430 greg 3.28 TMstruct *tms;
431 greg 3.3 {
432 greg 3.28 COLRDATA *cd;
433 greg 3.3 double d;
434 greg 3.21 int i, j;
435 greg 3.3
436     cd = (COLRDATA *)tms->pd[colrReg];
437 greg 3.21 for (i = 3; i--; )
438 greg 3.24 cd->clfb[i] = 0x1000*tms->clf[i] + .5;
439 greg 3.23 cd->inpsfb = tmCvLuminance(tms->inpsf);
440 greg 3.21 for (i = 3; i--; )
441     for (j = 3; j--; ) {
442     d = tms->cmat[i][j] / tms->inpsf;
443     cd->cmatb[i][j] = 0x10000*d + (d<0. ? -.5 : .5);
444     }
445 greg 3.3 }
446    
447    
448     static MEM_PTR
449     colrInit(tms) /* initialize private data for tone mapping */
450 greg 3.28 TMstruct *tms;
451 greg 3.3 {
452 greg 3.28 COLRDATA *cd;
453     int i;
454 greg 3.3 /* allocate our data */
455     cd = (COLRDATA *)malloc(sizeof(COLRDATA));
456     if (cd == NULL)
457     return(NULL);
458     tms->pd[colrReg] = (MEM_PTR)cd;
459     /* compute gamma table */
460     for (i = GAMTSZ; i--; )
461     cd->gamb[i] = 256.*pow((i+.5)/GAMTSZ, 1./tms->mongam);
462     /* compute color and scale factors */
463     colrNewSpace(tms);
464     return((MEM_PTR)cd);
465 greg 3.1 }