ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/common/tonemap.c

Comparing ray/src/common/tonemap.c (file contents):
Revision 3.5 by gregl, Wed Nov 12 13:25:48 1997 UTC vs.
Revision 3.13 by schorsch, Sun Jul 27 22:12:01 2003 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1997 Regents of the University of California */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ LBL";
2	>	static const char       RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		* Tone mapping functions.
6		* See tonemap.h for detailed function descriptions.
7	+	* Added von Kries white-balance calculations 10/01 (GW).
8	+	*
9	+	* Externals declared in tonemap.h
10		*/
11
12	+	#include "copyright.h"
13	+
14		#include        <stdio.h>
15		#include        <math.h>
16		#include        "tmprivat.h"
#	Line 26 | Line 28 | int    tmLastError;                    /* last error incurred by library *
28		char    *tmLastFunction;                /* error-generating function name */
29
30
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	–
31		struct tmStruct *
32		tmInit(flags, monpri, gamval)           /* initialize new tone mapping */
33		int     flags;
34		RGBPRIMP        monpri;
35		double  gamval;
36		{
52	–	static char     funcName[] = "tmInit";
37		COLORMAT        cmat;
38		register struct tmStruct        *tmnew;
39		register int    i;
#	Line 66 | Line 50 | double gamval;
50		tmnew->clf[GRN] = rgb2xyzmat[1][1];
51		tmnew->clf[BLU] = rgb2xyzmat[1][2];
52		} else {
53	<	comprgb2xyzmat(cmat, tmnew->monpri=monpri);
53	>	comprgb2xyzWBmat(cmat, tmnew->monpri=monpri);
54		tmnew->clf[RED] = cmat[1][0];
55		tmnew->clf[GRN] = cmat[1][1];
56		tmnew->clf[BLU] = cmat[1][2];
#	Line 86 | Line 70 | double gamval;
70		tmnew->inpsf = WHTEFFICACY;
71		tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
72		tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
73	<	tmnew->brmin = tmnew->brmax = 0;
73	>	tmnew->hbrmin = 10; tmnew->hbrmax = -10;
74		tmnew->histo = NULL;
75	+	tmnew->mbrmin = 10; tmnew->mbrmax = -10;
76		tmnew->lumap = NULL;
77		/* zero private data */
78		for (i = TM_MAXPKG; i--; )
#	Line 122 | Line 107 | double sf;
107		tmTop->clf[CIEX] = tmTop->clf[CIEZ] = 0.;
108		tmTop->clf[CIEY] = 1.;
109		} else {
110	<	comprgb2xyzmat(tmTop->cmat, tmTop->monpri);
110	>	comprgb2xyzWBmat(tmTop->cmat, tmTop->monpri);
111		tmTop->clf[RED] = tmTop->cmat[1][0];
112		tmTop->clf[GRN] = tmTop->cmat[1][1];
113		tmTop->clf[BLU] = tmTop->cmat[1][2];
#	Line 133 | Line 118 | double sf;
118		tmTop->cmat[1][2] = tmTop->cmat[2][0] = tmTop->cmat[2][1] = 0.;
119
120		} else if (tmTop->inppri == TM_XYZPRIM) /* input is XYZ */
121	<	compxyz2rgbmat(tmTop->cmat, tmTop->monpri);
121	>	compxyz2rgbWBmat(tmTop->cmat, tmTop->monpri);
122
123		else {                                  /* input is RGB */
124		if (tmTop->inppri != tmTop->monpri &&
125		PRIMEQ(tmTop->inppri, tmTop->monpri))
126		tmTop->inppri = tmTop->monpri;  /* no xform */
127	<	comprgb2rgbmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
127	>	comprgb2rgbWBmat(tmTop->cmat, tmTop->inppri, tmTop->monpri);
128		}
129		for (i = 0; i < 3; i++)
130		for (j = 0; j < 3; j++)
#	Line 177 | Line 162 | COLOR  *scan;
162		int     len;
163		{
164		static char     funcName[] = "tmCvColors";
165	<	static COLOR    csmall = {1e-6, 1e-6, 1e-6};
165	>	static COLOR    csmall = {.5*MINLUM, .5*MINLUM, .5*MINLUM};
166		COLOR   cmon;
167		double  lum, slum;
168		register double d;
#	Line 185 | Line 170 | int    len;
170
171		if (tmTop == NULL)
172		returnErr(TM_E_TMINVAL);
173	<	if (ls == NULL | scan == NULL | len <= 0)
173	>	if ((ls == NULL) | (scan == NULL) | (len < 0))
174		returnErr(TM_E_ILLEGAL);
175		for (i = len; i--; ) {
176		if (tmNeedMatrix(tmTop))                /* get monitor RGB */
#	Line 247 | Line 232 | int    len;
232
233
234		int
235	+	tmCvGrays(ls, scan, len)                /* convert float gray values */
236	+	TMbright        *ls;
237	+	float   *scan;
238	+	int     len;
239	+	{
240	+	static char     funcName[] = "tmCvGrays";
241	+	register double d;
242	+	register int    i;
243	+
244	+	if (tmTop == NULL)
245	+	returnErr(TM_E_TMINVAL);
246	+	if ((ls == NULL) | (scan == NULL) | (len < 0))
247	+	returnErr(TM_E_ILLEGAL);
248	+	for (i = len; i--; )
249	+	if (scan[i] <= TM_NOLUM) {
250	+	ls[i] = TM_NOBRT;               /* bogus value */
251	+	} else {
252	+	d = TM_BRTSCALE*log(scan[i]);   /* encode it */
253	+	ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5);
254	+	}
255	+	returnOK;
256	+	}
257	+
258	+
259	+	int
260		tmAddHisto(ls, len, wt)                 /* add values to histogram */
261		register TMbright       *ls;
262		int     len;
263		int     wt;
264		{
265		static char     funcName[] = "tmAddHisto";
266	<	int     sum, oldorig, oldlen, horig, hlen;
266	>	int     oldorig=0, oldlen, horig, hlen;
267		register int    i, j;
268
259	–	if (len <= 0)
260	–	returnErr(TM_E_ILLEGAL);
269		if (tmTop == NULL)
270		returnErr(TM_E_TMINVAL);
271	+	if (len < 0)
272	+	returnErr(TM_E_ILLEGAL);
273	+	if (len == 0)
274	+	returnOK;
275		/* first, grow limits */
276		if (tmTop->histo == NULL) {
277		for (i = len; i-- && ls[i] < MINBRT; )
278		;
279		if (i < 0)
280		returnOK;
281	<	tmTop->brmin = tmTop->brmax = ls[i];
281	>	tmTop->hbrmin = tmTop->hbrmax = ls[i];
282		oldlen = 0;
283		} else {
284	<	oldorig = (tmTop->brmin-MINBRT)/HISTEP;
285	<	oldlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - oldorig;
284	>	oldorig = (tmTop->hbrmin-MINBRT)/HISTEP;
285	>	oldlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - oldorig;
286		}
287		for (i = len; i--; ) {
288		if ((j = ls[i]) < MINBRT)
289		continue;
290	<	if (j < tmTop->brmin)
291	<	tmTop->brmin = j;
292	<	else if (j > tmTop->brmax)
293	<	tmTop->brmax = j;
290	>	if (j < tmTop->hbrmin)
291	>	tmTop->hbrmin = j;
292	>	else if (j > tmTop->hbrmax)
293	>	tmTop->hbrmax = j;
294		}
295	<	horig = (tmTop->brmin-MINBRT)/HISTEP;
296	<	hlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - horig;
295	>	horig = (tmTop->hbrmin-MINBRT)/HISTEP;
296	>	hlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - horig;
297		if (hlen > oldlen) {                    /* (re)allocate histogram */
298		register int    *newhist = (int *)calloc(hlen, sizeof(int));
299		if (newhist == NULL)
#	Line 292 | Line 304 | int    wt;
304		free((MEM_PTR)tmTop->histo);
305		}
306		tmTop->histo = newhist;
295	–	if (tmTop->lumap != NULL) {     /* invalid tone map */
296	–	free((MEM_PTR)tmTop->lumap);
297	–	tmTop->lumap = NULL;
298	–	}
307		}
308		if (wt == 0)
309		returnOK;
#	Line 328 | Line 336 | double La;
336		}
337
338
339	+	static int
340	+	tmNewMap()
341	+	{
342	+	if (tmTop->lumap != NULL && (tmTop->mbrmax - tmTop->mbrmin) !=
343	+	(tmTop->hbrmax - tmTop->hbrmin)) {
344	+	free((MEM_PTR)tmTop->lumap);
345	+	tmTop->lumap = NULL;
346	+	}
347	+	tmTop->mbrmin = tmTop->hbrmin;
348	+	tmTop->mbrmax = tmTop->hbrmax;
349	+	if (tmTop->mbrmin > tmTop->mbrmax)
350	+	return 0;
351	+	if (tmTop->lumap == NULL)
352	+	tmTop->lumap = (unsigned short *)malloc(sizeof(unsigned short)*
353	+	(tmTop->mbrmax-tmTop->mbrmin+1));
354	+	return(tmTop->lumap != NULL);
355	+	}
356	+
357	+
358		int
359	+	tmFixedMapping(expmult, gamval)
360	+	double  expmult;
361	+	double  gamval;
362	+	{
363	+	static char     funcName[] = "tmFixedMapping";
364	+	double          d;
365	+	register int    i;
366	+
367	+	if (!tmNewMap())
368	+	returnErr(TM_E_NOMEM);
369	+	if (expmult <= .0)
370	+	expmult = 1.;
371	+	if (gamval < MINGAM)
372	+	gamval = tmTop->mongam;
373	+	d = log(expmult/tmTop->inpsf);
374	+	for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; )
375	+	tmTop->lumap[i] = 256. * exp(
376	+	( d + (tmTop->mbrmin+i)*(1./TM_BRTSCALE) )
377	+	/ gamval );
378	+	returnOK;
379	+	}
380	+
381	+
382	+	int
383		tmComputeMapping(gamval, Lddyn, Ldmax)
384		double  gamval;
385		double  Lddyn;
#	Line 337 | Line 388 | double Ldmax;
388		static char     funcName[] = "tmComputeMapping";
389		int     *histo;
390		float   *cumf;
391	<	int     brt0, histlen, histot, threshold, ceiling, trimmings;
391	>	int     brt0, histlen, threshold, ceiling, trimmings;
392		double  logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
393	<	int4    sum;
393	>	int32   histot;
394	>	double  sum;
395		register double d;
396		register int    i, j;
397
#	Line 353 | Line 405 | double Ldmax;
405		Ldmin = Ldmax/Lddyn;
406		logLddyn = log(Lddyn);
407		Ldavg = sqrt(Ldmax*Ldmin);
408	<	i = (tmTop->brmin-MINBRT)/HISTEP;
408	>	i = (tmTop->hbrmin-MINBRT)/HISTEP;
409		brt0 = MINBRT + HISTEP/2 + i*HISTEP;
410	<	histlen = (tmTop->brmax-MINBRT)/HISTEP + 1 - i;
410	>	histlen = (tmTop->hbrmax-MINBRT)/HISTEP + 1 - i;
411		/* histogram total and mean */
412		histot = 0; sum = 0;
413		j = brt0 + histlen*HISTEP;
#	Line 367 | Line 419 | double Ldmax;
419		if (threshold < 4)
420		returnErr(TM_E_TMFAIL);
421		Lwavg = tmLuminance( (double)sum / histot );
422	<	if (!(tmTop->flags & TM_F_LINEAR)) {    /* clamp histogram */
423	<	histo = (int *)malloc(histlen*sizeof(int));
424	<	cumf = (float *)malloc((histlen+1)*sizeof(float));
425	<	if (histo == NULL | cumf == NULL)
426	<	returnErr(TM_E_NOMEM);
427	<	for (i = histlen; i--; )        /* make malleable copy */
428	<	histo[i] = tmTop->histo[i];
429	<	do {                            /* iterate to solution */
430	<	sum = 0;                /* cumulative probability */
431	<	for (i = 0; i < histlen; i++) {
432	<	cumf[i] = (double)sum/histot;
433	<	sum += histo[i];
422	>	/* allocate space for mapping */
423	>	if (!tmNewMap())
424	>	returnErr(TM_E_NOMEM);
425	>	/* use linear tone mapping? */
426	>	if (tmTop->flags & TM_F_LINEAR)
427	>	goto linearmap;
428	>	/* clamp histogram */
429	>	histo = (int *)malloc(histlen*sizeof(int));
430	>	cumf = (float *)malloc((histlen+2)*sizeof(float));
431	>	if ((histo == NULL) | (cumf == NULL))
432	>	returnErr(TM_E_NOMEM);
433	>	cumf[histlen+1] = 1.;           /* guard for assignment code */
434	>	for (i = histlen; i--; )        /* make malleable copy */
435	>	histo[i] = tmTop->histo[i];
436	>	do {                            /* iterate to solution */
437	>	sum = 0;                /* cumulative probability */
438	>	for (i = 0; i < histlen; i++) {
439	>	cumf[i] = (double)sum/histot;
440	>	sum += histo[i];
441	>	}
442	>	cumf[histlen] = 1.;
443	>	Tr = histot * (double)(tmTop->hbrmax - tmTop->hbrmin) /
444	>	((double)histlen*TM_BRTSCALE) / logLddyn;
445	>	ceiling = Tr + 1.;
446	>	trimmings = 0;          /* clip to envelope */
447	>	for (i = histlen; i--; ) {
448	>	if (tmTop->flags & TM_F_HCONTR) {
449	>	Lw = tmLuminance(brt0 + i*HISTEP);
450	>	Ld = Ldmin * exp( logLddyn *
451	>	.5*(cumf[i]+cumf[i+1]) );
452	>	ceiling = Tr * (htcontrs(Ld) * Lw) /
453	>	(htcontrs(Lw) * Ld) + 1.;
454		}
455	<	cumf[i] = 1.;
456	<	Tr = histot * (double)(tmTop->brmax - tmTop->brmin) /
457	<	((double)histlen*TM_BRTSCALE) / logLddyn;
386	<	ceiling = Tr + 1.;
387	<	trimmings = 0;                  /* clip to envelope */
388	<	for (i = histlen; i--; ) {
389	<	if (tmTop->flags & TM_F_HCONTR) {
390	<	Lw = tmLuminance(brt0 + i*HISTEP);
391	<	Ld = Ldmin * exp( logLddyn *
392	<	.5*(cumf[i]+cumf[i+1]) );
393	<	ceiling = Tr * (htcontrs(Ld) * Lw) /
394	<	(htcontrs(Lw) * Ld) + 1.;
395	<	}
396	<	if (histo[i] > ceiling) {
397	<	trimmings += histo[i] - ceiling;
398	<	histo[i] = ceiling;
399	<	}
455	>	if (histo[i] > ceiling) {
456	>	trimmings += histo[i] - ceiling;
457	>	histo[i] = ceiling;
458		}
401	–	} while ((histot -= trimmings) > threshold &&
402	–	trimmings > threshold);
403	–	}
404	–	/* allocate luminance map */
405	–	if (tmTop->lumap == NULL) {
406	–	tmTop->lumap = (unsigned short *)malloc(
407	–	(tmTop->brmax-tmTop->brmin+1)*sizeof(unsigned short) );
408	–	if (tmTop->lumap == NULL)
409	–	returnErr(TM_E_NOMEM);
410	–	}
411	–	if (tmTop->flags & TM_F_LINEAR || histot <= threshold) {
412	–	/* linear tone mapping */
413	–	if (tmTop->flags & TM_F_HCONTR)
414	–	d = htcontrs(Ldavg) / htcontrs(Lwavg);
415	–	else
416	–	d = Ldavg / Lwavg;
417	–	d = log(d/Ldmax);
418	–	for (i = tmTop->brmax-tmTop->brmin+1; i--; )
419	–	tmTop->lumap[i] = 256. * exp(
420	–	( d + (tmTop->brmin+i)/(double)TM_BRTSCALE )
421	–	/ gamval );
422	–	} else {
423	–	/* histogram adjustment */
424	–	for (i = tmTop->brmax-tmTop->brmin+1; i--; ) {
425	–	j = d = (double)i/(tmTop->brmax-tmTop->brmin)*histlen;
426	–	d -= (double)j;
427	–	Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
428	–	d = (Ld - Ldmin)/(Ldmax - Ldmin);
429	–	tmTop->lumap[i] = 256.*pow(d, 1./gamval);
459		}
460	+	/* check if we're out of data */
461	+	if ((histot -= trimmings) <= threshold) {
462	+	free((MEM_PTR)histo);
463	+	free((MEM_PTR)cumf);
464	+	goto linearmap;
465	+	}
466	+	} while (trimmings > threshold);
467	+	/* assign tone-mapping */
468	+	for (i = tmTop->mbrmax-tmTop->mbrmin+1; i--; ) {
469	+	j = d = (double)i/(tmTop->mbrmax-tmTop->mbrmin)*histlen;
470	+	d -= (double)j;
471	+	Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1]));
472	+	d = (Ld - Ldmin)/(Ldmax - Ldmin);
473	+	tmTop->lumap[i] = 256.*pow(d, 1./gamval);
474		}
475	<	if (!(tmTop->flags & TM_F_LINEAR)) {
476	<	free((MEM_PTR)histo);
434	<	free((MEM_PTR)cumf);
435	<	}
475	>	free((MEM_PTR)histo);           /* clean up and return */
476	>	free((MEM_PTR)cumf);
477		returnOK;
478	+	linearmap:                              /* linear tone-mapping */
479	+	if (tmTop->flags & TM_F_HCONTR)
480	+	d = htcontrs(Ldavg) / htcontrs(Lwavg);
481	+	else
482	+	d = Ldavg / Lwavg;
483	+	return(tmFixedMapping(tmTop->inpsf*d/Ldmax, gamval));
484		}
485
486
#	Line 445 | Line 492 | register BYTE  *cs;
492		int     len;
493		{
494		static char     funcName[] = "tmMapPixels";
495	<	register int4   li, pv;
495	>	register int32  li, pv;
496
497		if (tmTop == NULL || tmTop->lumap == NULL)
498		returnErr(TM_E_TMINVAL);
499	<	if (ps == NULL | ls == NULL | len <= 0)
499	>	if ((ps == NULL) | (ls == NULL) | (len < 0))
500		returnErr(TM_E_ILLEGAL);
501		while (len--) {
502	<	if ((li = *ls++) < tmTop->brmin)
503	<	li = tmTop->brmin;
504	<	else if (li > tmTop->brmax)
505	<	li = tmTop->brmax;
506	<	li = tmTop->lumap[li - tmTop->brmin];
502	>	if ((li = *ls++) < tmTop->mbrmin) {
503	>	li = 0;
504	>	} else {
505	>	if (li > tmTop->mbrmax)
506	>	li = tmTop->mbrmax;
507	>	li = tmTop->lumap[li - tmTop->mbrmin];
508	>	}
509		if (cs == TM_NOCHROM)
510		*ps++ = li>255 ? 255 : li;
511		else {
#	Line 489 | Line 538 | register struct tmStruct       *tms;
538		{
539		register struct tmStruct        *tms2;
540		/* special cases first */
541	<	if (tms == NULL | tmTop == NULL)
541	>	if ((tms == NULL) | (tmTop == NULL))
542		return(0);
543		if (tms == tmTop) {
544		tmTop = tms->tmprev;
#	Line 520 | Line 569 | tmDup()                                /* duplicate top tone mapping */
569		return(NULL);
570		*tmnew = *tmTop;                /* copy everything */
571		if (tmnew->histo != NULL) {     /* duplicate histogram */
572	<	len = (tmnew->brmax-MINBRT)/HISTEP + 1 -
573	<	(tmnew->brmin-MINBRT)/HISTEP;
572	>	len = (tmnew->hbrmax-MINBRT)/HISTEP + 1 -
573	>	(tmnew->hbrmin-MINBRT)/HISTEP;
574		tmnew->histo = (int *)malloc(len*sizeof(int));
575		if (tmnew->histo != NULL)
576		for (i = len; i--; )
577		tmnew->histo[i] = tmTop->histo[i];
578		}
579		if (tmnew->lumap != NULL) {     /* duplicate luminance mapping */
580	<	len = tmnew->brmax-tmnew->brmin+1;
580	>	len = tmnew->mbrmax-tmnew->mbrmin+1;
581		tmnew->lumap = (unsigned short *)malloc(
582		len*sizeof(unsigned short) );
583		if (tmnew->lumap != NULL)
#	Line 582 | Line 631 | register struct tmStruct       *tms;
631		if (tms->pd[i] != NULL)
632		(*tmPkg[i]->Free)(tms->pd[i]);
633		free((MEM_PTR)tms);             /* free basic structure */
634	+	}
635	+
636	+	/******************** Shared but Private library routines *********************/
637	+
638	+	BYTE    tmMesofact[BMESUPPER-BMESLOWER];
639	+
640	+	void
641	+	tmMkMesofact()                          /* build mesopic lookup factor table */
642	+	{
643	+	register int    i;
644	+
645	+	if (tmMesofact[BMESUPPER-BMESLOWER-1])
646	+	return;
647	+
648	+	for (i = BMESLOWER; i < BMESUPPER; i++)
649	+	tmMesofact[i-BMESLOWER] = 256. *
650	+	(tmLuminance(i) - LMESLOWER) /
651	+	(LMESUPPER - LMESLOWER);
652	+	}
653	+
654	+
655	+	int
656	+	tmErrorReturn(func, err)                /* error return (with message) */
657	+	char    *func;
658	+	int     err;
659	+	{
660	+	tmLastFunction = func;
661	+	tmLastError = err;
662	+	if (tmTop != NULL && tmTop->flags & TM_F_NOSTDERR)
663	+	return(err);
664	+	fputs(func, stderr);
665	+	fputs(": ", stderr);
666	+	fputs(tmErrorMessage[err], stderr);
667	+	fputs("!\n", stderr);
668	+	return(err);
669		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)