--- ray/src/common/tonemap.c 2006/03/23 22:01:43 3.22 +++ ray/src/common/tonemap.c 2021/01/07 19:13:57 3.41 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: tonemap.c,v 3.22 2006/03/23 22:01:43 greg Exp $"; +static const char RCSid[] = "$Id: tonemap.c,v 3.41 2021/01/07 19:13:57 greg Exp $"; #endif /* * Tone mapping functions. @@ -12,7 +12,9 @@ static const char RCSid[] = "$Id: tonemap.c,v 3.22 200 #include "copyright.h" #include +#include #include +#include #include "tmprivat.h" #include "tmerrmsg.h" @@ -22,7 +24,12 @@ static const char RCSid[] = "$Id: tonemap.c,v 3.22 200 struct tmPackage *tmPkg[TM_MAXPKG]; int tmNumPkgs = 0; /* number of registered packages */ + /* luminance->brightness lookup */ +static TMbright *tmFloat2BrtLUT = NULL; +#define tmCvLumLUfp(pf) tmFloat2BrtLUT[*(int32 *)(pf) >> 15] + + TMstruct * tmInit( /* initialize new tone mapping */ int flags, @@ -48,7 +55,7 @@ double gamval tmnew->clf[GRN] = rgb2xyzmat[1][1]; tmnew->clf[BLU] = rgb2xyzmat[1][2]; } else { - comprgb2xyzWBmat(cmat, tmnew->monpri=monpri); + comprgb2xyzmat(cmat, tmnew->monpri=monpri); tmnew->clf[RED] = cmat[1][0]; tmnew->clf[GRN] = cmat[1][1]; tmnew->clf[BLU] = cmat[1][2]; @@ -60,7 +67,7 @@ double gamval tmnew->mongam = gamval; /* set color divisors */ for (i = 0; i < 3; i++) - tmnew->cdiv[i] = 256.*pow(tmnew->clf[i], 1./tmnew->mongam); + tmnew->cdiv[i] = TM_BRES*pow(tmnew->clf[i], 1./tmnew->mongam); /* set input transform */ tmnew->inppri = tmnew->monpri; @@ -111,7 +118,7 @@ MEM_PTR dat tms->clf[CIEX] = tms->clf[CIEZ] = 0.; tms->clf[CIEY] = 1.; } else { - comprgb2xyzWBmat(tms->cmat, tms->monpri); + comprgb2xyzmat(tms->cmat, tms->monpri); tms->clf[RED] = tms->cmat[1][0]; tms->clf[GRN] = tms->cmat[1][1]; tms->clf[BLU] = tms->cmat[1][2]; @@ -128,18 +135,16 @@ MEM_PTR dat if (tms->inppri != tms->monpri && PRIMEQ(tms->inppri, tms->monpri)) tms->inppri = tms->monpri; /* no xform */ - comprgb2rgbWBmat(tms->cmat, tms->inppri, tms->monpri); + if (!comprgb2rgbWBmat(tms->cmat, tms->inppri, tms->monpri)) + returnErr(TM_E_ILLEGAL); } for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) tms->cmat[i][j] *= tms->inpsf; /* set color divisors */ for (i = 0; i < 3; i++) - if (tms->clf[i] > .001) - tms->cdiv[i] = - 256.*pow(tms->clf[i], 1./tms->mongam); - else - tms->cdiv[i] = 1; + tms->cdiv[i] = TM_BRES*pow(tms->clf[i] < .001 ? .001 : + tms->clf[i], 1./tms->mongam); /* notify packages */ for (i = tmNumPkgs; i--; ) if (tms->pd[i] != NULL && tmPkg[i]->NewSpace != NULL) @@ -149,7 +154,7 @@ MEM_PTR dat void -tmClearHisto( /* clear current histogram */ +tmClearHisto( /* clear current histogram */ TMstruct *tms ) { @@ -160,28 +165,112 @@ TMstruct *tms } +TMbright +tmCvLuminance( /* convert a single luminance */ +double lum +) +{ + double d; + +#ifdef isfinite + if (!isfinite(lum) || lum <= TM_NOLUM) +#else + if (lum <= TM_NOLUM) +#endif + return(TM_NOBRT); + d = TM_BRTSCALE*log(lum); + return((TMbright)(d + .5 - (d < 0.))); +} + + int +tmCvLums( /* convert luminances using lookup */ +TMbright *ls, +float *scan, +int len +) +{ + if (tmFloat2BrtLUT == NULL) { /* initialize lookup table */ + int32 i; + tmFloat2BrtLUT = (TMbright *)malloc(sizeof(TMbright)*0x10000); + if (tmFloat2BrtLUT == NULL) + return(TM_E_NOMEM); + for (i = 0; i < 0x10000; i++) { + int32 l = (i<<1 | 1) << 14; +#ifndef isfinite + if ((l & 0x7f800000) == 0x7f800000) + tmFloat2BrtLUT[i] = TM_NOBRT; + else +#endif + tmFloat2BrtLUT[i] = tmCvLuminance(*(float *)&l); + } + } + if (len <= 0) + return(TM_E_OK); + if ((ls == NULL) | (scan == NULL)) + return(TM_E_ILLEGAL); + while (len--) { + if (*scan <= TM_NOLUM) { + *ls++ = TM_NOBRT; + ++scan; + continue; + } + *ls++ = tmCvLumLUfp(scan++); + } + return(TM_E_OK); +} + + +int +tmCvGrays( /* convert float gray values */ +TMstruct *tms, +TMbright *ls, +float *scan, +int len +) +{ + static const char funcName[] = "tmCvGrays"; + int i; + + if (tms == NULL) + returnErr(TM_E_TMINVAL); + if ((ls == NULL) | (scan == NULL) | (len < 0)) + returnErr(TM_E_ILLEGAL); + if (tmFloat2BrtLUT == NULL) /* initialize */ + tmCvLums(NULL, NULL, 0); + for (i = len; i--; ) { + float lum = tms->inpsf * scan[i]; + if (lum <= TM_NOLUM) + ls[i] = TM_NOBRT; + else + ls[i] = tmCvLumLUfp(&lum); + } + returnOK; +} + + +int tmCvColors( /* convert float colors */ TMstruct *tms, TMbright *ls, -BYTE *cs, +uby8 *cs, COLOR *scan, int len ) { static const char funcName[] = "tmCvColors"; - static COLOR csmall = {.5*MINLUM, .5*MINLUM, .5*MINLUM}; - static BYTE gamtab[1024]; + static uby8 gamtab[1024]; static double curgam = .0; COLOR cmon; - double lum, slum; - double d; + float lum, slum, d; int i; if (tms == NULL) returnErr(TM_E_TMINVAL); if ((ls == NULL) | (scan == NULL) | (len < 0)) returnErr(TM_E_ILLEGAL); + if (tmFloat2BrtLUT == NULL) /* initialize */ + tmCvLums(NULL, NULL, 0); if (cs != TM_NOCHROM && fabs(tms->mongam - curgam) > .02) { curgam = tms->mongam; /* (re)build table */ for (i = 1024; i--; ) @@ -195,41 +284,38 @@ int len cmon[GRN] = tms->inpsf*scan[i][GRN]; cmon[BLU] = tms->inpsf*scan[i][BLU]; } -#ifdef isnan - if (isnan(cmon[RED]) || isinf(cmon[RED])) cmon[RED] = .0f; - if (isnan(cmon[GRN]) || isinf(cmon[GRN])) cmon[GRN] = .0f; - if (isnan(cmon[BLU]) || isinf(cmon[BLU])) cmon[BLU] = .0f; +#ifdef isfinite + if (!isfinite(cmon[RED]) || cmon[RED] < .0f) cmon[RED] = .0f; + if (!isfinite(cmon[GRN]) || cmon[GRN] < .0f) cmon[GRN] = .0f; + if (!isfinite(cmon[BLU]) || cmon[BLU] < .0f) cmon[BLU] = .0f; +#else + if (cmon[RED] < .0f) cmon[RED] = .0f; + if (cmon[GRN] < .0f) cmon[GRN] = .0f; + if (cmon[BLU] < .0f) cmon[BLU] = .0f; #endif /* world luminance */ lum = tms->clf[RED]*cmon[RED] + tms->clf[GRN]*cmon[GRN] + tms->clf[BLU]*cmon[BLU] ; - /* check range */ - if (clipgamut(cmon, lum, CGAMUT_LOWER, csmall, cwhite)) - lum = tms->clf[RED]*cmon[RED] + - tms->clf[GRN]*cmon[GRN] + - tms->clf[BLU]*cmon[BLU] ; - if (lum < MINLUM) { - ls[i] = MINBRT-1; /* bogus value */ - lum = MINLUM; - } else { - d = TM_BRTSCALE*log(lum); /* encode it */ - ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5); - } + if (lum <= TM_NOLUM) { /* convert brightness */ + lum = cmon[RED] = cmon[GRN] = cmon[BLU] = TM_NOLUM; + ls[i] = TM_NOBRT; + } else + ls[i] = tmCvLumLUfp(&lum); if (cs == TM_NOCHROM) /* no color? */ continue; if (tms->flags & TM_F_MESOPIC && lum < LMESUPPER) { slum = scotlum(cmon); /* mesopic adj. */ - if (lum < LMESLOWER) + if (lum < LMESLOWER) { cmon[RED] = cmon[GRN] = cmon[BLU] = slum; - else { + } else { d = (lum - LMESLOWER)/(LMESUPPER - LMESLOWER); if (tms->flags & TM_F_BW) cmon[RED] = cmon[GRN] = cmon[BLU] = d*lum; else scalecolor(cmon, d); - d = (1.-d)*slum; + d = (1.f-d)*slum; cmon[RED] += d; cmon[GRN] += d; cmon[BLU] += d; @@ -238,60 +324,17 @@ int len cmon[RED] = cmon[GRN] = cmon[BLU] = lum; } d = tms->clf[RED]*cmon[RED]/lum; - cs[3*i ] = d>=.999 ? 255 : gamtab[(int)(1024.*d)]; + cs[3*i ] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)]; d = tms->clf[GRN]*cmon[GRN]/lum; - cs[3*i+1] = d>=.999 ? 255 : gamtab[(int)(1024.*d)]; + cs[3*i+1] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)]; d = tms->clf[BLU]*cmon[BLU]/lum; - cs[3*i+2] = d>=.999 ? 255 : gamtab[(int)(1024.*d)]; + cs[3*i+2] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)]; } returnOK; } -TMbright -tmCvLuminance( /* convert a single luminance */ -double lum -) -{ - double d; - - if (lum <= TM_NOLUM) - return(TM_NOBRT); - d = TM_BRTSCALE*log(lum); - if (d > 0.) - return((TMbright)(d+.5)); - return((TMbright)(d-.5)); -} - - int -tmCvGrays( /* convert float gray values */ -TMstruct *tms, -TMbright *ls, -float *scan, -int len -) -{ - static const char funcName[] = "tmCvGrays"; - double d; - int i; - - if (tms == NULL) - returnErr(TM_E_TMINVAL); - if ((ls == NULL) | (scan == NULL) | (len < 0)) - returnErr(TM_E_ILLEGAL); - for (i = len; i--; ) - if (scan[i] <= TM_NOLUM) { - ls[i] = TM_NOBRT; /* bogus value */ - } else { - d = TM_BRTSCALE*log(scan[i]); /* encode it */ - ls[i] = d>0. ? (int)(d+.5) : (int)(d-.5); - } - returnOK; -} - - -int tmAddHisto( /* add values to histogram */ TMstruct *tms, TMbright *ls, @@ -318,8 +361,8 @@ int wt tms->hbrmin = tms->hbrmax = ls[i]; oldlen = 0; } else { - oldorig = (tms->hbrmin-MINBRT)/HISTEP; - oldlen = (tms->hbrmax-MINBRT)/HISTEP + 1 - oldorig; + oldorig = HISTI(tms->hbrmin); + oldlen = HISTI(tms->hbrmax) + 1 - oldorig; } for (i = len; i--; ) { if ((j = ls[i]) < MINBRT) @@ -329,15 +372,15 @@ int wt else if (j > tms->hbrmax) tms->hbrmax = j; } - horig = (tms->hbrmin-MINBRT)/HISTEP; - hlen = (tms->hbrmax-MINBRT)/HISTEP + 1 - horig; + horig = HISTI(tms->hbrmin); + hlen = HISTI(tms->hbrmax) + 1 - horig; if (hlen > oldlen) { /* (re)allocate histogram */ - int *newhist = (int *)calloc(hlen, sizeof(int)); + HIST_TYP *newhist = (HIST_TYP *)calloc(hlen, sizeof(HIST_TYP)); if (newhist == NULL) returnErr(TM_E_NOMEM); if (oldlen) { /* copy and free old */ - for (i = oldlen, j = i+oldorig-horig; i; ) - newhist[--j] = tms->histo[--i]; + memcpy((MEM_PTR)(newhist+(oldorig-horig)), + (MEM_PTR)tms->histo, oldlen*sizeof(HIST_TYP)); free((MEM_PTR)tms->histo); } tms->histo = newhist; @@ -346,7 +389,7 @@ int wt returnOK; for (i = len; i--; ) /* add in new counts */ if (ls[i] >= MINBRT) - tms->histo[ (ls[i]-MINBRT)/HISTEP - horig ] += wt; + tms->histo[ HISTI(ls[i]) - horig ] += wt; returnOK; } @@ -374,27 +417,6 @@ double La } -static int -tmNewMap( /* allocate new tone-mapping array */ -TMstruct *tms -) -{ - if (tms->lumap != NULL && (tms->mbrmax - tms->mbrmin) != - (tms->hbrmax - tms->hbrmin)) { - free((MEM_PTR)tms->lumap); - tms->lumap = NULL; - } - tms->mbrmin = tms->hbrmin; - tms->mbrmax = tms->hbrmax; - if (tms->mbrmin > tms->mbrmax) - return 0; - if (tms->lumap == NULL) - tms->lumap = (unsigned short *)malloc(sizeof(unsigned short)* - (tms->mbrmax-tms->mbrmin+1)); - return(tms->lumap != NULL); -} - - int tmFixedMapping( /* compute fixed, linear tone-mapping */ TMstruct *tms, @@ -413,10 +435,14 @@ double gamval if (gamval < MINGAM) gamval = tms->mongam; d = log(expmult/tms->inpsf); - for (i = tms->mbrmax-tms->mbrmin+1; i--; ) - tms->lumap[i] = 256. * exp( + for (i = tms->mbrmax-tms->mbrmin+1; i--; ) { + double val = TM_BRES * exp( ( d + (tms->mbrmin+i)*(1./TM_BRTSCALE) ) - / gamval ); + / gamval); + tms->lumap[i] = val; + if (sizeof(TMAP_TYP) == 2 && val >= 0xffff) + tms->lumap[i] = 0xffff; + } returnOK; } @@ -430,11 +456,11 @@ double Ldmax ) { static const char funcName[] = "tmComputeMapping"; - int *histo; + HIST_TYP *histo; float *cumf; - int brt0, histlen, threshold, ceiling, trimmings; + int brt0, histlen; + HIST_TYP threshold, ceiling, trimmings, histot; double logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld; - int32 histot; double sum; double d; int i, j; @@ -449,28 +475,26 @@ double Ldmax Ldmin = Ldmax/Lddyn; logLddyn = log(Lddyn); Ldavg = sqrt(Ldmax*Ldmin); - i = (tms->hbrmin-MINBRT)/HISTEP; - brt0 = MINBRT + HISTEP/2 + i*HISTEP; - histlen = (tms->hbrmax-MINBRT)/HISTEP + 1 - i; + i = HISTI(tms->hbrmin); + brt0 = HISTV(i); + histlen = HISTI(tms->hbrmax) + 1 - i; /* histogram total and mean */ histot = 0; sum = 0; j = brt0 + histlen*HISTEP; for (i = histlen; i--; ) { histot += tms->histo[i]; - sum += (j -= HISTEP) * tms->histo[i]; + sum += (double)(j -= HISTEP) * tms->histo[i]; } threshold = histot*0.005 + .5; - if (threshold < 4) + if (!histot) returnErr(TM_E_TMFAIL); Lwavg = tmLuminance( (double)sum / histot ); - /* allocate space for mapping */ - if (!tmNewMap(tms)) - returnErr(TM_E_NOMEM); /* use linear tone mapping? */ - if (tms->flags & TM_F_LINEAR) + if (tms->flags & TM_F_LINEAR || threshold < 4 || + tms->hbrmax - tms->hbrmin < TM_BRTSCALE*logLddyn) goto linearmap; /* clamp histogram */ - histo = (int *)malloc(histlen*sizeof(int)); + histo = (HIST_TYP *)malloc(histlen*sizeof(HIST_TYP)); cumf = (float *)malloc((histlen+2)*sizeof(float)); if ((histo == NULL) | (cumf == NULL)) returnErr(TM_E_NOMEM); @@ -508,13 +532,16 @@ double Ldmax goto linearmap; } } while (trimmings > threshold); + /* allocate space for mapping */ + if (!tmNewMap(tms)) + returnErr(TM_E_NOMEM); /* assign tone-mapping */ for (i = tms->mbrmax-tms->mbrmin+1; i--; ) { j = d = (double)i/(tms->mbrmax-tms->mbrmin)*histlen; d -= (double)j; Ld = Ldmin*exp(logLddyn*((1.-d)*cumf[j]+d*cumf[j+1])); d = (Ld - Ldmin)/(Ldmax - Ldmin); - tms->lumap[i] = 256.*pow(d, 1./gamval); + tms->lumap[i] = TM_BRES*pow(d, 1./gamval); } free((MEM_PTR)histo); /* clean up and return */ free((MEM_PTR)cumf); @@ -531,29 +558,31 @@ linearmap: /* linear tone-mapping */ int tmMapPixels( /* apply tone-mapping to pixel(s) */ TMstruct *tms, -BYTE *ps, +uby8 *ps, TMbright *ls, -BYTE *cs, +uby8 *cs, int len ) { static const char funcName[] = "tmMapPixels"; - int32 li, pv; + TMbright lv; + TMAP_TYP li; + int pv; if (tms == NULL || tms->lumap == NULL) returnErr(TM_E_TMINVAL); if ((ps == NULL) | (ls == NULL) | (len < 0)) returnErr(TM_E_ILLEGAL); while (len--) { - if ((li = *ls++) < tms->mbrmin) { + if ((lv = *ls++) < tms->mbrmin) { li = 0; } else { - if (li > tms->mbrmax) - li = tms->mbrmax; - li = tms->lumap[li - tms->mbrmin]; + if (lv > tms->mbrmax) + lv = tms->mbrmax; + li = tms->lumap[lv - tms->mbrmin]; } if (cs == TM_NOCHROM) - *ps++ = li>255 ? 255 : li; + *ps++ = li>=TM_BRES ? 255 : (int)(256*li/TM_BRES); else { pv = *cs++ * li / tms->cdiv[RED]; *ps++ = pv>255 ? 255 : pv; @@ -567,8 +596,6 @@ int len } - - TMstruct * tmDup( /* duplicate top tone mapping */ TMstruct *tms @@ -585,20 +612,18 @@ TMstruct *tms return(NULL); *tmnew = *tms; /* copy everything */ if (tmnew->histo != NULL) { /* duplicate histogram */ - len = (tmnew->hbrmax-MINBRT)/HISTEP + 1 - - (tmnew->hbrmin-MINBRT)/HISTEP; - tmnew->histo = (int *)malloc(len*sizeof(int)); + len = HISTI(tmnew->hbrmax) + 1 - HISTI(tmnew->hbrmin); + tmnew->histo = (HIST_TYP *)malloc(len*sizeof(HIST_TYP)); if (tmnew->histo != NULL) - for (i = len; i--; ) - tmnew->histo[i] = tms->histo[i]; + memcpy((MEM_PTR)tmnew->histo, (MEM_PTR)tms->histo, + len*sizeof(HIST_TYP)); } if (tmnew->lumap != NULL) { /* duplicate luminance mapping */ len = tmnew->mbrmax-tmnew->mbrmin+1; - tmnew->lumap = (unsigned short *)malloc( - len*sizeof(unsigned short) ); + tmnew->lumap = (TMAP_TYP *)malloc(len*sizeof(TMAP_TYP)); if (tmnew->lumap != NULL) - for (i = len; i--; ) - tmnew->lumap[i] = tms->lumap[i]; + memcpy((MEM_PTR)tmnew->lumap, (MEM_PTR)tms->lumap, + len*sizeof(TMAP_TYP)); } /* clear package data */ for (i = tmNumPkgs; i--; ) @@ -630,7 +655,7 @@ TMstruct *tms; /******************** Shared but Private library routines *********************/ -BYTE tmMesofact[BMESUPPER-BMESLOWER]; +uby8 tmMesofact[BMESUPPER-BMESLOWER]; void tmMkMesofact() /* build mesopic lookup factor table */ @@ -644,6 +669,27 @@ tmMkMesofact() /* build mesopic lookup factor table tmMesofact[i-BMESLOWER] = 256. * (tmLuminance(i) - LMESLOWER) / (LMESUPPER - LMESLOWER); +} + + +int +tmNewMap( /* allocate new tone-mapping array */ +TMstruct *tms +) +{ + if (tms->lumap != NULL && (tms->mbrmax - tms->mbrmin) != + (tms->hbrmax - tms->hbrmin)) { + free((MEM_PTR)tms->lumap); + tms->lumap = NULL; + } + tms->mbrmin = tms->hbrmin; + tms->mbrmax = tms->hbrmax; + if (tms->mbrmin > tms->mbrmax) + return 0; + if (tms->lumap == NULL) + tms->lumap = (TMAP_TYP *)malloc(sizeof(TMAP_TYP)* + (tms->mbrmax-tms->mbrmin+1)); + return(tms->lumap != NULL); }