src/common/tmap16bit.c

#ifndef lint
static const char RCSid[] = "$Id: tmap16bit.c,v 1.8 2008/05/15 22:15:16 greg Exp $";
#endif
/*
 * Routines for tone-mapping 16-bit/primary pixels
 *
 * Externals declared in tonemap.h
 */

#include "copyright.h"

#include <stdio.h>
#include <math.h>
#include "tmprivat.h"

#define LOGTABBITS      11      /* log table is 1<<LOGTABBITS long */
#define GAMTABBITS      9       /* gamma table is 1<<GAMTABBITS long */
#define MONGAMTSZ       1024    /* monitor gamma table length */

static float    logtab[1<<LOGTABBITS];
static float    gamtab[1<<GAMTABBITS];
static float    gammul[16];
static double   cur_gam = 0.;
static BYTE     mongamtab[MONGAMTSZ];
static double   cur_mongam = 0.;

#define imultpow2(i,s)  ((s)>=0 ? (i)<<(s) : (i)>>-(s))


/* Fill our log table */
static void
mkLogTable()
{       
        int     i;

        if (logtab[0] != 0.f)
                return;
        for (i = 1<<LOGTABBITS; i--; )
                logtab[i] = log(imultpow2(i,15-LOGTABBITS)*(1./(1L<<16))
                                        + .5);
}


/* Fill our input gamma table */
static void
mkGamTable(double gv)
{
        int     i;

        if (gv == cur_gam)
                return;
        for (i = 1<<GAMTABBITS; i--; )
                gamtab[i] = pow((i+.5)*(1./(1<<GAMTABBITS)), gv);
        for (i = 16; i--; )
                gammul[i] = pow((double)(1L<<i), -gv);
        cur_gam = gv;
}


/* Fill our monitor gamma table */
static void
mkMonGamTable(TMstruct *tms)
{
        int     i;
        
        if (tms->mongam == cur_mongam)
                return;
        for (i = MONGAMTSZ; i--; )
                mongamtab[i] = 256.*pow((i+.5)*(1./MONGAMTSZ), 1./tms->mongam);
        cur_mongam = tms->mongam;
}


/* Find normalizing shift value for a 2-byte unsigned integer */
static int
normShift16(int i)
{
        int     s = 0;
        
        if (!i)
                return(-1);
        while (!(i & 0x8000)) {
                i <<= 1;
                ++s;
        }
        return(s);
}


/* Find common normalizing shift for 3 2-byte unsigned integers */
static int
normShift48(uint16 clr48[3])
{
        int     imax = (clr48[1] > clr48[0] ? clr48[1] : clr48[0]);
        if (clr48[2] > imax)
                imax = clr48[2];
        return(normShift16(imax));
}


/* convert at 48-bit tristimulus value to a COLOR */
static void
rgb48_color(COLOR col, uint16 clr48[3], double gv)
{
        int     nshft;

        if (gv == 1.) {                         /* linear case */
                col[0] = clr48[0]*(1./(1L<<16));
                col[1] = clr48[1]*(1./(1L<<16));
                col[2] = clr48[2]*(1./(1L<<16));
                return;
        }
                                                /* non-linear case */
        /* XXX Uncomment if this routine is made public
        if (gv != cur_gam)
                mkGamTable(gv);
        */
        nshft = normShift48(clr48);
        if (nshft < 0) {
                col[0] = col[1] = col[2] = .0f;
                return;
        }
        col[0] = gamtab[imultpow2(clr48[0],GAMTABBITS-16+nshft)] * 
                        gammul[nshft];
        col[1] = gamtab[imultpow2(clr48[1],GAMTABBITS-16+nshft)] *
                        gammul[nshft];
        col[2] = gamtab[imultpow2(clr48[2],GAMTABBITS-16+nshft)] *
                        gammul[nshft];
}


/* Convert 16-bit gray scanline to encoded luminance */
int
tmCvGray16(TMstruct *tms, TMbright *ls, uint16 *scan, int len, double gv)
{
        static const char       funcName[] = "tmCvGray16";
        static double   cur_inpsf = 1.;
        static double   log_inpsf = 0.;
        int             nshft;
        double          d;

        if (tms == NULL)
                returnErr(TM_E_TMINVAL);
        if ((ls == NULL) | (scan == NULL) | (len < 0))
                returnErr(TM_E_ILLEGAL);
        if (gv <= 0.)
                gv = DEFGAM;
                                                /* initialize log table */
        if (logtab[0] == 0.f)
                mkLogTable();
        if (cur_inpsf != tms->inpsf)
                log_inpsf = log(cur_inpsf = tms->inpsf);
                                                /* convert 16-bit grays */
        while (len--) {
                nshft = normShift16(*scan);
                if (nshft < 0) {                /* bogus value */
                        *ls++ = TM_NOBRT;
                        scan++;
                        continue;
                }
                d = logtab[ imultpow2(*scan,LOGTABBITS-15+nshft) &
                                        ((1L<<LOGTABBITS)-1) ];
                d -= M_LN2*nshft;
                d = (double)TM_BRTSCALE * (gv*d + log_inpsf);
                *ls++ = (d>0. ? d+.5 : d-.5);
                scan++;
        }
        returnOK;
}

/* Convert a 48-bit RGB scanline to encoded luminance/chrominance */
int
tmCvRGB48(TMstruct *tms, TMbright *ls, BYTE *cs,
                uint16 (*scan)[3], int len, double gv)
{
        static const char       funcName[] = "tmCvRGB48";
        static double   cur_inpsf = 1.;
        static double   log_inpsf = 0.;
        int             i;

        if (tms == NULL)
                returnErr(TM_E_TMINVAL);
        if ((ls == NULL) | (scan == NULL) | (len < 0))
                returnErr(TM_E_ILLEGAL);
        if (gv <= 0.)
                gv = DEFGAM;
                                                /* sync input gamma table */
        if (gv != cur_gam)
                mkGamTable(gv);
        if (tmNeedMatrix(tms)) {                /* need floating point */
                COLOR   *newscan;
                newscan = (COLOR *)tempbuffer(len*sizeof(COLOR));
                if (newscan == NULL)
                        returnErr(TM_E_NOMEM);
                for (i = len; i--; )
                        rgb48_color(newscan[i], scan[i], gv);
                return(tmCvColors(tms, ls, cs, newscan, len));
        }
                                                /* sync monitor gamma table */
        if (cs != TM_NOCHROM && tms->mongam != cur_mongam)
                mkMonGamTable(tms);
                                                /* initialize log table */
        if (logtab[0] == 0.f)
                mkLogTable();
        if (cur_inpsf != tms->inpsf)
                log_inpsf = log(cur_inpsf = tms->inpsf);
        if (tms->flags & TM_F_MESOPIC)
                tmMkMesofact();
                                                /* convert scanline */
        for (i = len; i--; ) {
                int     nshft = normShift48(scan[i]);
                COLOR   cmon;
                double  lum;
                int     bi;
                
                if (nshft < 0) {
                        bi = TM_NOBRT;                  /* bogus value */
                        lum = 1.;
                        setcolor(cmon, 1., 1., 1.);
                } else {
                        int     j = GAMTABBITS-16+nshft;
                        int     nshft2;
                        double  d;
                                                        /* normalized linear */
                        setcolor(cmon,  gamtab[imultpow2(scan[i][0],j)],
                                        gamtab[imultpow2(scan[i][1],j)],
                                        gamtab[imultpow2(scan[i][2],j)] );
                        lum =   tms->clf[RED]*cmon[RED];
                        lum +=  tms->clf[GRN]*cmon[GRN];
                        lum +=  tms->clf[BLU]*cmon[BLU];
                                                        /* convert to log Y */
                        j = lum * (double)(1L<<16);
                        nshft2 = normShift16(j);
                        d = logtab[ imultpow2(j,LOGTABBITS-15+nshft2) &
                                                ((1L<<LOGTABBITS)-1) ];
                        d -= M_LN2*(gv*nshft + nshft2);
                        d = (double)TM_BRTSCALE*(d + log_inpsf);
                        bi = (int)(d + .5 - (d < 0.));
                }
                                                        /* world luminance */
                ls[i] = bi;
                if (cs == TM_NOCHROM)                   /* no color? */
                        continue;
                                                        /* mesopic adj. */
                if (tms->flags & TM_F_MESOPIC && bi < BMESUPPER) {
                        double  slum = scotlum(cmon);
                        if (bi < BMESLOWER)
                                setcolor(cmon, slum, slum, slum);
                        else {
                                double  pf;
                                pf = (1./256.)*tmMesofact[bi-BMESLOWER];
                                if (tms->flags & TM_F_BW)
                                        cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
                                slum *= 1. - pf;
                                cmon[RED] = slum + pf*cmon[RED];
                                cmon[GRN] = slum + pf*cmon[GRN];
                                cmon[BLU] = slum + pf*cmon[BLU];
                        }
                } else if (tms->flags & TM_F_BW) {
                        cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
                }
                bi = (double)MONGAMTSZ*tms->clf[RED]*cmon[RED]/lum;
                cs[3*i  ] = bi>=MONGAMTSZ ? 255 : mongamtab[bi];
                bi = (double)MONGAMTSZ*tms->clf[GRN]*cmon[GRN]/lum;
                cs[3*i+1] = bi>=MONGAMTSZ ? 255 : mongamtab[bi];
                bi = (double)MONGAMTSZ*tms->clf[BLU]*cmon[BLU]/lum;
                cs[3*i+2] = bi>=MONGAMTSZ ? 255 : mongamtab[bi];
        }
        returnOK;
}
Revision:	1.9
Committed:	Mon Mar 7 20:49:19 2011 UTC (14 years, 8 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.8:	+2 -2 lines
Log Message:	Minor improvement in round-off method
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	1.9	static const char RCSid[] = "$Id: tmap16bit.c,v 1.8 2008/05/15 22:15:16 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* Routines for tone-mapping 16-bit/primary pixels
6			*
7			* Externals declared in tonemap.h
8			*/
9
10			#include "copyright.h"
11
12			#include <stdio.h>
13			#include <math.h>
14			#include "tmprivat.h"
15
16			#define LOGTABBITS 11 /* log table is 1<<LOGTABBITS long */
17			#define GAMTABBITS 9 /* gamma table is 1<<GAMTABBITS long */
18	greg	1.3	#define MONGAMTSZ 1024 /* monitor gamma table length */
19	greg	1.1
20			static float logtab[1<<LOGTABBITS];
21			static float gamtab[1<<GAMTABBITS];
22			static float gammul[16];
23			static double cur_gam = 0.;
24	greg	1.3	static BYTE mongamtab[MONGAMTSZ];
25			static double cur_mongam = 0.;
26	greg	1.1
27			#define imultpow2(i,s) ((s)>=0 ? (i)<<(s) : (i)>>-(s))
28
29
30			/* Fill our log table */
31			static void
32			mkLogTable()
33			{
34			int i;
35
36			if (logtab[0] != 0.f)
37			return;
38			for (i = 1<<LOGTABBITS; i--; )
39			logtab[i] = log(imultpow2(i,15-LOGTABBITS)*(1./(1L<<16))
40			+ .5);
41			}
42
43
44	greg	1.3	/* Fill our input gamma table */
45	greg	1.1	static void
46			mkGamTable(double gv)
47			{
48			int i;
49
50			if (gv == cur_gam)
51			return;
52			for (i = 1<<GAMTABBITS; i--; )
53			gamtab[i] = pow((i+.5)*(1./(1<<GAMTABBITS)), gv);
54			for (i = 16; i--; )
55			gammul[i] = pow((double)(1L<<i), -gv);
56			cur_gam = gv;
57			}
58
59
60	greg	1.3	/* Fill our monitor gamma table */
61			static void
62	greg	1.5	mkMonGamTable(TMstruct *tms)
63	greg	1.3	{
64			int i;
65
66	greg	1.5	if (tms->mongam == cur_mongam)
67	greg	1.3	return;
68			for (i = MONGAMTSZ; i--; )
69	greg	1.5	mongamtab[i] = 256.pow((i+.5)(1./MONGAMTSZ), 1./tms->mongam);
70			cur_mongam = tms->mongam;
71	greg	1.3	}
72
73
74	greg	1.1	/* Find normalizing shift value for a 2-byte unsigned integer */
75			static int
76			normShift16(int i)
77			{
78			int s = 0;
79
80	greg	1.2	if (!i)
81			return(-1);
82	greg	1.1	while (!(i & 0x8000)) {
83			i <<= 1;
84			++s;
85			}
86			return(s);
87			}
88
89
90			/* Find common normalizing shift for 3 2-byte unsigned integers */
91			static int
92			normShift48(uint16 clr48[3])
93			{
94			int imax = (clr48[1] > clr48[0] ? clr48[1] : clr48[0]);
95			if (clr48[2] > imax)
96			imax = clr48[2];
97			return(normShift16(imax));
98			}
99
100
101			/* convert at 48-bit tristimulus value to a COLOR */
102			static void
103	greg	1.2	rgb48_color(COLOR col, uint16 clr48[3], double gv)
104	greg	1.1	{
105			int nshft;
106
107	greg	1.2	if (gv == 1.) { /* linear case */
108	greg	1.1	col[0] = clr48[0]*(1./(1L<<16));
109			col[1] = clr48[1]*(1./(1L<<16));
110			col[2] = clr48[2]*(1./(1L<<16));
111			return;
112			}
113			/* non-linear case */
114	greg	1.3	/* XXX Uncomment if this routine is made public
115	greg	1.2	if (gv != cur_gam)
116			mkGamTable(gv);
117			*/
118	greg	1.1	nshft = normShift48(clr48);
119	greg	1.2	if (nshft < 0) {
120			col[0] = col[1] = col[2] = .0f;
121			return;
122			}
123	greg	1.1	col[0] = gamtab[imultpow2(clr48[0],GAMTABBITS-16+nshft)] *
124			gammul[nshft];
125			col[1] = gamtab[imultpow2(clr48[1],GAMTABBITS-16+nshft)] *
126			gammul[nshft];
127			col[2] = gamtab[imultpow2(clr48[2],GAMTABBITS-16+nshft)] *
128			gammul[nshft];
129			}
130
131
132			/* Convert 16-bit gray scanline to encoded luminance */
133			int
134	greg	1.5	tmCvGray16(TMstruct tms, TMbright ls, uint16 *scan, int len, double gv)
135	greg	1.1	{
136	greg	1.6	static const char funcName[] = "tmCvGray16";
137	greg	1.1	static double cur_inpsf = 1.;
138			static double log_inpsf = 0.;
139			int nshft;
140			double d;
141
142	greg	1.5	if (tms == NULL)
143	greg	1.1	returnErr(TM_E_TMINVAL);
144	schorsch	1.4	if ((ls == NULL) \| (scan == NULL) \| (len < 0))
145	greg	1.1	returnErr(TM_E_ILLEGAL);
146	greg	1.2	if (gv <= 0.)
147			gv = DEFGAM;
148	greg	1.1	/* initialize log table */
149			if (logtab[0] == 0.f)
150			mkLogTable();
151	greg	1.5	if (cur_inpsf != tms->inpsf)
152			log_inpsf = log(cur_inpsf = tms->inpsf);
153	greg	1.1	/* convert 16-bit grays */
154			while (len--) {
155			nshft = normShift16(*scan);
156	greg	1.2	if (nshft < 0) { /* bogus value */
157			*ls++ = TM_NOBRT;
158			scan++;
159			continue;
160			}
161	greg	1.1	d = logtab[ imultpow2(*scan,LOGTABBITS-15+nshft) &
162	greg	1.3	((1L<<LOGTABBITS)-1) ];
163			d -= M_LN2*nshft;
164	greg	1.1	d = (double)TM_BRTSCALE * (gv*d + log_inpsf);
165			*ls++ = (d>0. ? d+.5 : d-.5);
166			scan++;
167			}
168			returnOK;
169			}
170
171			/* Convert a 48-bit RGB scanline to encoded luminance/chrominance */
172			int
173	greg	1.5	tmCvRGB48(TMstruct tms, TMbright ls, BYTE *cs,
174			uint16 (*scan)[3], int len, double gv)
175	greg	1.1	{
176	greg	1.6	static const char funcName[] = "tmCvRGB48";
177	greg	1.1	static double cur_inpsf = 1.;
178			static double log_inpsf = 0.;
179			int i;
180
181	greg	1.5	if (tms == NULL)
182	greg	1.1	returnErr(TM_E_TMINVAL);
183	schorsch	1.4	if ((ls == NULL) \| (scan == NULL) \| (len < 0))
184	greg	1.1	returnErr(TM_E_ILLEGAL);
185	greg	1.2	if (gv <= 0.)
186			gv = DEFGAM;
187	greg	1.3	/* sync input gamma table */
188	greg	1.8	if (gv != cur_gam)
189	greg	1.1	mkGamTable(gv);
190	greg	1.5	if (tmNeedMatrix(tms)) { /* need floating point */
191	greg	1.1	COLOR *newscan;
192			newscan = (COLOR )tempbuffer(lensizeof(COLOR));
193			if (newscan == NULL)
194			returnErr(TM_E_NOMEM);
195			for (i = len; i--; )
196	greg	1.2	rgb48_color(newscan[i], scan[i], gv);
197	greg	1.5	return(tmCvColors(tms, ls, cs, newscan, len));
198	greg	1.1	}
199	greg	1.3	/* sync monitor gamma table */
200	greg	1.5	if (cs != TM_NOCHROM && tms->mongam != cur_mongam)
201			mkMonGamTable(tms);
202	greg	1.1	/* initialize log table */
203			if (logtab[0] == 0.f)
204			mkLogTable();
205	greg	1.5	if (cur_inpsf != tms->inpsf)
206			log_inpsf = log(cur_inpsf = tms->inpsf);
207			if (tms->flags & TM_F_MESOPIC)
208	greg	1.3	tmMkMesofact();
209	greg	1.1	/* convert scanline */
210			for (i = len; i--; ) {
211	greg	1.3	int nshft = normShift48(scan[i]);
212			COLOR cmon;
213			double lum;
214			int bi;
215
216			if (nshft < 0) {
217	greg	1.1	bi = TM_NOBRT; /* bogus value */
218	greg	1.3	lum = 1.;
219	greg	1.7	setcolor(cmon, 1., 1., 1.);
220	greg	1.3	} else {
221			int j = GAMTABBITS-16+nshft;
222			int nshft2;
223			double d;
224			/* normalized linear */
225			setcolor(cmon, gamtab[imultpow2(scan[i][0],j)],
226			gamtab[imultpow2(scan[i][1],j)],
227			gamtab[imultpow2(scan[i][2],j)] );
228	greg	1.5	lum = tms->clf[RED]*cmon[RED];
229			lum += tms->clf[GRN]*cmon[GRN];
230			lum += tms->clf[BLU]*cmon[BLU];
231	greg	1.3	/* convert to log Y */
232			j = lum * (double)(1L<<16);
233			nshft2 = normShift16(j);
234			d = logtab[ imultpow2(j,LOGTABBITS-15+nshft2) &
235			((1L<<LOGTABBITS)-1) ];
236			d -= M_LN2(gvnshft + nshft2);
237			d = (double)TM_BRTSCALE*(d + log_inpsf);
238	greg	1.9	bi = (int)(d + .5 - (d < 0.));
239	greg	1.1	}
240	greg	1.3	/* world luminance */
241	greg	1.1	ls[i] = bi;
242			if (cs == TM_NOCHROM) /* no color? */
243			continue;
244			/* mesopic adj. */
245	greg	1.5	if (tms->flags & TM_F_MESOPIC && bi < BMESUPPER) {
246	greg	1.3	double slum = scotlum(cmon);
247	greg	1.1	if (bi < BMESLOWER)
248	greg	1.3	setcolor(cmon, slum, slum, slum);
249	greg	1.1	else {
250	greg	1.3	double pf;
251			pf = (1./256.)*tmMesofact[bi-BMESLOWER];
252	greg	1.5	if (tms->flags & TM_F_BW)
253	greg	1.3	cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
254			slum *= 1. - pf;
255			cmon[RED] = slum + pf*cmon[RED];
256			cmon[GRN] = slum + pf*cmon[GRN];
257			cmon[BLU] = slum + pf*cmon[BLU];
258	greg	1.1	}
259	greg	1.5	} else if (tms->flags & TM_F_BW) {
260	greg	1.3	cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
261	greg	1.1	}
262	greg	1.5	bi = (double)MONGAMTSZtms->clf[RED]cmon[RED]/lum;
263	greg	1.3	cs[3*i ] = bi>=MONGAMTSZ ? 255 : mongamtab[bi];
264	greg	1.5	bi = (double)MONGAMTSZtms->clf[GRN]cmon[GRN]/lum;
265	greg	1.3	cs[3*i+1] = bi>=MONGAMTSZ ? 255 : mongamtab[bi];
266	greg	1.5	bi = (double)MONGAMTSZtms->clf[BLU]cmon[BLU]/lum;
267	greg	1.3	cs[3*i+2] = bi>=MONGAMTSZ ? 255 : mongamtab[bi];
268	greg	1.1	}
269			returnOK;
270			}