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

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: tmapcolrs.c,v 3.34 2019/12/28 18:05:14 greg Exp $";
3 #endif
4 /*
5 * Routines for tone mapping on Radiance RGBE and XYZE pictures.
6 *
7 * Externals declared in tonemap.h
8 */
9
10 #include "copyright.h"
11
12 #include <stdlib.h>
13 #include <math.h>
14
15 #ifdef PCOND
16 #include "paths.h"
17 #endif
18 #include "tmprivat.h"
19 #include "rtio.h"
20 #include "resolu.h"
21
22 #define GAMTSZ 4096
23
24 typedef struct {
25 uby8 gamb[GAMTSZ]; /* gamma lookup table */
26 int clfb[3]; /* encoded tm->clf */
27 int32 cmatb[3][3]; /* encoded color transform */
28 TMbright inpsfb; /* encoded tm->inpsf */
29 } COLRDATA;
30
31 static MEM_PTR colrInit(TMstruct *);
32 static void colrNewSpace(TMstruct *);
33 static gethfunc headline;
34
35 static struct tmPackage colrPkg = { /* our package functions */
36 colrInit, colrNewSpace, free
37 };
38 static int colrReg = -1; /* our package registration number */
39
40 static TMbright logi[256];
41
42
43 int
44 tmCvColrs( /* convert RGBE/XYZE colors */
45 TMstruct *tms,
46 TMbright *ls,
47 uby8 *cs,
48 COLR *scan,
49 int len
50 )
51 {
52 static const char funcName[] = "tmCvColrs";
53 int cmon[4];
54 COLRDATA *cd;
55 int i, j, bi;
56 int32 li, vl;
57
58 if (tms == NULL)
59 returnErr(TM_E_TMINVAL);
60 if ((ls == NULL) | (scan == NULL) | (len < 0))
61 returnErr(TM_E_ILLEGAL);
62 if (colrReg < 0) { /* build tables if necessary */
63 colrReg = tmRegPkg(&colrPkg);
64 if (colrReg < 0)
65 returnErr(TM_E_CODERR1);
66 for (i = 256; i--; )
67 logi[i] = TM_BRTSCALE*log((i+.5)/256.) - .5;
68 tmMkMesofact();
69 }
70 if ((cd = (COLRDATA *)tmPkgData(tms,colrReg)) == NULL)
71 returnErr(TM_E_NOMEM);
72 for (i = len; i--; ) {
73 if (tmNeedMatrix(tms)) { /* apply color xform */
74 bi = 0;
75 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 if (vl < 0)
80 cmon[j] = vl/(int32)0x10000;
81 else if ((cmon[j] = vl>>16) > bi)
82 bi = cmon[j];
83 }
84 cmon[EXP] = scan[i][EXP];
85 while (bi >= 256) { /* handle overflow */
86 cmon[EXP]++;
87 for (j = 3; j--; ) cmon[j] >>= 1;
88 bi >>= 1;
89 }
90 } else
91 copycolr(cmon, scan[i]);
92 /* world luminance */
93 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 else li >>= 12;
98 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 }
105 ls[i] = bi;
106 if (cs == TM_NOCHROM) /* no color? */
107 continue;
108 /* mesopic adj. */
109 if (tms->flags & TM_F_MESOPIC && bi < BMESUPPER) {
110 int pf, sli = normscot(cmon);
111 if (bi < BMESLOWER) {
112 cmon[RED] = cmon[GRN] = cmon[BLU] = sli;
113 } else {
114 if (tms->flags & TM_F_BW)
115 cmon[RED] = cmon[GRN] = cmon[BLU] = li;
116 pf = tmMesofact[bi-BMESLOWER];
117 sli *= 256 - pf;
118 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 }
124 } else if (tms->flags & TM_F_BW) {
125 cmon[RED] = cmon[GRN] = cmon[BLU] = li;
126 } else {
127 for (j = 3; j--; )
128 cmon[j] *= (cmon[j] > 0);
129 }
130 bi = ( (uint32)GAMTSZ*cd->clfb[RED]*cmon[RED]/li ) >> 12;
131 cs[3*i ] = bi>=GAMTSZ ? 255 : cd->gamb[bi];
132 bi = ( (uint32)GAMTSZ*cd->clfb[GRN]*cmon[GRN]/li ) >> 12;
133 cs[3*i+1] = bi>=GAMTSZ ? 255 : cd->gamb[bi];
134 bi = ( (uint32)GAMTSZ*cd->clfb[BLU]*cmon[BLU]/li ) >> 12;
135 cs[3*i+2] = bi>=GAMTSZ ? 255 : cd->gamb[bi];
136 }
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 headline( /* grok a header line */
156 char *s,
157 void *vrh
158 )
159 {
160 char fmt[MAXFMTLEN];
161 struct radhead *rh = vrh;
162
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 return(0);
171 }
172 if (isexpos(s)) {
173 rh->expos *= exposval(s);
174 return(0);
175 }
176 if (isprims(s)) {
177 primsval(rh->mypri, s);
178 rh->primp = rh->mypri;
179 return(0);
180 }
181 return(0);
182 }
183
184
185 int
186 tmLoadPicture( /* convert Radiance picture */
187 TMstruct *tms,
188 TMbright **lpp,
189 uby8 **cpp,
190 int *xp,
191 int *yp,
192 char *fname,
193 FILE *fp
194 )
195 {
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 if (tms == NULL)
204 returnErr(TM_E_TMINVAL);
205 if ((lpp == NULL) | (xp == NULL) | (yp == NULL) |
206 ((fname == NULL) & (fp == TM_GETFILE)))
207 returnErr(TM_E_ILLEGAL);
208 *xp = *yp = 0; /* error precaution */
209 if ((inpf = fp) == TM_GETFILE && (inpf = fopen(fname, "rb")) == NULL)
210 returnErr(TM_E_BADFILE);
211 *lpp = NULL;
212 if (cpp != TM_NOCHROMP) *cpp = NULL;
213 info = rhdefault; /* get our header */
214 getheader(inpf, headline, &info);
215 if ((info.format == FMTBAD) | (info.expos <= 0.) ||
216 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 if ((err = tmSetSpace(tms, info.primp, 1./info.expos, NULL)) != TM_E_OK)
227 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 *cpp = (uby8 *)malloc(3*sizeof(uby8) * *xp * *yp);
234 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 err = tmCvColrs(tms, *lpp + (i * *xp),
246 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 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 returnErr(err);
262 }
263 returnOK;
264 }
265
266
267 #ifdef PCOND
268 static int /* run pcond to map picture */
269 dopcond(psp, xp, yp, flags, monpri, gamval, Lddyn, Ldmax, fname)
270 uby8 **psp;
271 int *xp, *yp;
272 int flags;
273 RGBPRIMP monpri;
274 double gamval, Lddyn, Ldmax;
275 char *fname;
276 {
277 char *funcName = fname;
278 TMstruct *tms = NULL;
279 char cmdbuf[1024];
280 FILE *infp;
281 COLR *scan;
282 uby8 *rp;
283 int y;
284 int x;
285 /* set up gamma correction */
286 if (setcolrcor(pow, 1./gamval) < 0)
287 returnErr(TM_E_NOMEM);
288 /* create command */
289 strcpy(cmdbuf, PCOND);
290 if (flags & TM_F_HCONTR)
291 strcat(cmdbuf, " -s");
292 if (flags & TM_F_MESOPIC)
293 strcat(cmdbuf, " -c");
294 if (flags & TM_F_LINEAR)
295 strcat(cmdbuf, " -l");
296 if (flags & TM_F_ACUITY)
297 strcat(cmdbuf, " -a");
298 if (flags & TM_F_VEIL)
299 strcat(cmdbuf, " -v");
300 if (flags & TM_F_CWEIGHT)
301 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 /* 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 scan = (COLR *)malloc(sizeof(COLR) * *xp);
320 if (flags & TM_F_BW)
321 rp = (uby8 *)malloc(sizeof(uby8) * *xp * *yp);
322 else
323 rp = (uby8 *)malloc(3*sizeof(uby8) * *xp * *yp);
324 if (((*psp = rp) == NULL) | (scan == NULL)) {
325 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 free((MEM_PTR)scan);
333 free((MEM_PTR)*psp);
334 *psp = NULL;
335 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 free((MEM_PTR)scan);
349 pclose(infp);
350 returnOK;
351 }
352 #endif
353
354
355 int /* map a Radiance picture */
356 tmMapPicture(psp, xp, yp, flags, monpri, gamval, Lddyn, Ldmax, fname, fp)
357 uby8 **psp;
358 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 TMstruct *tms = NULL;
367 uby8 *cp;
368 TMbright *lp;
369 int err;
370 /* check arguments */
371 if ((psp == NULL) | (xp == NULL) | (yp == NULL) | (monpri == NULL) |
372 ((fname == NULL) & (fp == TM_GETFILE)))
373 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 #ifdef PCOND
380 /* 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 #endif
385 /* initialize tone mapping */
386 if ((tms = tmInit(flags, monpri, gamval)) == NULL)
387 returnErr(TM_E_NOMEM);
388 /* load & convert picture */
389 err = tmLoadPicture(tms, &lp, (flags&TM_F_BW) ? TM_NOCHROMP : &cp,
390 xp, yp, fname, fp);
391 if (err != TM_E_OK) {
392 tmDone(tms);
393 return(err);
394 }
395 /* allocate space for result */
396 if (flags & TM_F_BW) {
397 *psp = (uby8 *)malloc(sizeof(uby8) * *xp * *yp);
398 if (*psp == NULL) {
399 free((MEM_PTR)lp);
400 tmDone(tms);
401 returnErr(TM_E_NOMEM);
402 }
403 cp = TM_NOCHROM;
404 } else
405 *psp = cp;
406 /* compute color mapping */
407 err = tmAddHisto(tms, lp, *xp * *yp, 1);
408 if (err != TM_E_OK)
409 goto done;
410 err = tmComputeMapping(tms, gamval, Lddyn, Ldmax);
411 if (err != TM_E_OK)
412 goto done;
413 /* map colors */
414 err = tmMapPixels(tms, *psp, lp, cp, *xp * *yp);
415
416 done: /* clean up */
417 free((MEM_PTR)lp);
418 tmDone(tms);
419 if (err != TM_E_OK) { /* free memory on error */
420 free((MEM_PTR)*psp);
421 *psp = NULL;
422 returnErr(err);
423 }
424 returnOK;
425 }
426
427
428 static void
429 colrNewSpace(tms) /* color space changed for tone mapping */
430 TMstruct *tms;
431 {
432 COLRDATA *cd;
433 double d;
434 int i, j;
435
436 cd = (COLRDATA *)tms->pd[colrReg];
437 for (i = 3; i--; )
438 cd->clfb[i] = 0x1000*tms->clf[i] + .5;
439 cd->inpsfb = tmCvLuminance(tms->inpsf);
440 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 }
446
447
448 static MEM_PTR
449 colrInit(tms) /* initialize private data for tone mapping */
450 TMstruct *tms;
451 {
452 COLRDATA *cd;
453 int i;
454 /* 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 }