src/common/tonemap.c

#ifndef lint
static const char       RCSid[] = "$Id: tonemap.c,v 3.40 2021/01/07 02:13:49 greg Exp $";
#endif
/*
 * Tone mapping functions.
 * See tonemap.h for detailed function descriptions.
 * Added von Kries white-balance calculations 10/01 (GW).
 *
 * Externals declared in tonemap.h
 */

#include "copyright.h"

#include        <stdio.h>
#include        <stdlib.h>
#include        <math.h>
#include        <string.h>
#include        "tmprivat.h"
#include        "tmerrmsg.h"

#define exp10(x)        exp(M_LN10*(x))

                                        /* our list of conversion packages */
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,
RGBPRIMP        monpri,
double  gamval
)
{
        COLORMAT        cmat;
        TMstruct        *tmnew;
        int     i;
                                                /* allocate structure */
        tmnew = (TMstruct *)malloc(sizeof(TMstruct));
        if (tmnew == NULL)
                return(NULL);

        tmnew->flags = flags & ~TM_F_UNIMPL;
        if (tmnew->flags & TM_F_BW)
                tmnew->flags &= ~TM_F_MESOPIC;
                                                /* set monitor transform */
        if (monpri == NULL || monpri == stdprims || tmnew->flags & TM_F_BW) {
                tmnew->monpri = stdprims;
                tmnew->clf[RED] = rgb2xyzmat[1][0];
                tmnew->clf[GRN] = rgb2xyzmat[1][1];
                tmnew->clf[BLU] = rgb2xyzmat[1][2];
        } else {
                comprgb2xyzmat(cmat, tmnew->monpri=monpri);
                tmnew->clf[RED] = cmat[1][0];
                tmnew->clf[GRN] = cmat[1][1];
                tmnew->clf[BLU] = cmat[1][2];
        }
                                                /* set gamma value */
        if (gamval < MINGAM)
                tmnew->mongam = DEFGAM;
        else
                tmnew->mongam = gamval;
                                                /* set color divisors */
        for (i = 0; i < 3; i++)
                tmnew->cdiv[i] = TM_BRES*pow(tmnew->clf[i], 1./tmnew->mongam);

                                                /* set input transform */
        tmnew->inppri = tmnew->monpri;
        tmnew->cmat[0][0] = tmnew->cmat[1][1] = tmnew->cmat[2][2] =
                        tmnew->inpsf = WHTEFFICACY;
        tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
        tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
        tmnew->inpdat = NULL;
        tmnew->hbrmin = 10; tmnew->hbrmax = -10;
        tmnew->histo = NULL;
        tmnew->mbrmin = 10; tmnew->mbrmax = -10;
        tmnew->lumap = NULL;
                                                /* zero private data */
        for (i = TM_MAXPKG; i--; )
                tmnew->pd[i] = NULL;
        tmnew->lastError = TM_E_OK;
        tmnew->lastFunc = "NoErr";
                                                /* return new TMstruct */
        return(tmnew);
}


int
tmSetSpace(                     /* set input color space for conversions */
TMstruct        *tms,
RGBPRIMP        pri,
double  sf,
MEM_PTR dat
)
{
        static const char funcName[] = "tmSetSpace";
        int     i, j;
                                                /* error check */
        if (tms == NULL)
                returnErr(TM_E_TMINVAL);
        if (sf <= 1e-12)
                returnErr(TM_E_ILLEGAL);
                                                /* check if no change */
        if (pri == tms->inppri && FEQ(sf, tms->inpsf) && dat == tms->inpdat)
                returnOK;
        tms->inppri = pri;                      /* let's set it */
        tms->inpsf = sf;
        tms->inpdat = dat;

        if (tms->flags & TM_F_BW) {             /* color doesn't matter */
                tms->monpri = tms->inppri;              /* eliminate xform */
                if (tms->inppri == TM_XYZPRIM) {
                        tms->clf[CIEX] = tms->clf[CIEZ] = 0.;
                        tms->clf[CIEY] = 1.;
                } else {
                        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];
                }
                tms->cmat[0][0] = tms->cmat[1][1] = tms->cmat[2][2] =
                                tms->inpsf;
                tms->cmat[0][1] = tms->cmat[0][2] = tms->cmat[1][0] =
                tms->cmat[1][2] = tms->cmat[2][0] = tms->cmat[2][1] = 0.;

        } else if (tms->inppri == TM_XYZPRIM)   /* input is XYZ */
                compxyz2rgbWBmat(tms->cmat, tms->monpri);

        else {                                  /* input is RGB */
                if (tms->inppri != tms->monpri &&
                                PRIMEQ(tms->inppri, tms->monpri))
                        tms->inppri = tms->monpri;      /* no xform */
                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++)
                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)
                        (*tmPkg[i]->NewSpace)(tms);
        returnOK;
}


void
tmClearHisto(                           /* clear current histogram */
TMstruct        *tms
)
{
        if (tms == NULL || tms->histo == NULL)
                return;
        free((MEM_PTR)tms->histo);
        tms->histo = NULL;
}


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,
uby8    *cs,
COLOR   *scan,
int     len
)
{
        static const char funcName[] = "tmCvColors";
        static uby8     gamtab[1024];
        static double   curgam = .0;
        COLOR   cmon;
        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--; )
                        gamtab[i] = (int)(256.*pow((i+.5)/1024., 1./curgam));
        }
        for (i = len; i--; ) {
                if (tmNeedMatrix(tms)) {                /* get monitor RGB */
                        colortrans(cmon, tms->cmat, scan[i]);
                } else {
                        cmon[RED] = tms->inpsf*scan[i][RED];
                        cmon[GRN] = tms->inpsf*scan[i][GRN];
                        cmon[BLU] = tms->inpsf*scan[i][BLU];
                }
#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] ;
                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) {
                                cmon[RED] = cmon[GRN] = cmon[BLU] = slum;
                        } 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.f-d)*slum;
                                cmon[RED] += d;
                                cmon[GRN] += d;
                                cmon[BLU] += d;
                        }
                } else if (tms->flags & TM_F_BW) {
                        cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
                }
                d = tms->clf[RED]*cmon[RED]/lum;
                cs[3*i  ] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
                d = tms->clf[GRN]*cmon[GRN]/lum;
                cs[3*i+1] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
                d = tms->clf[BLU]*cmon[BLU]/lum;
                cs[3*i+2] = d>=.999f ? 255 : gamtab[(int)(1024.f*d)];
        }
        returnOK;
}


int
tmAddHisto(                             /* add values to histogram */
TMstruct        *tms,
TMbright        *ls,
int     len,
int     wt
)
{
        static const char funcName[] = "tmAddHisto";
        int     oldorig=0, oldlen, horig, hlen;
        int     i, j;

        if (tms == NULL)
                returnErr(TM_E_TMINVAL);
        if (len < 0)
                returnErr(TM_E_ILLEGAL);
        if (len == 0)
                returnOK;
                                                /* first, grow limits */
        if (tms->histo == NULL) {
                for (i = len; i-- && ls[i] < MINBRT; )
                        ;
                if (i < 0)
                        returnOK;
                tms->hbrmin = tms->hbrmax = ls[i];
                oldlen = 0;
        } else {
                oldorig = HISTI(tms->hbrmin);
                oldlen = HISTI(tms->hbrmax) + 1 - oldorig;
        }
        for (i = len; i--; ) {
                if ((j = ls[i]) < MINBRT)
                        continue;
                if (j < tms->hbrmin)
                        tms->hbrmin = j;
                else if (j > tms->hbrmax)
                        tms->hbrmax = j;
        }
        horig = HISTI(tms->hbrmin);
        hlen = HISTI(tms->hbrmax) + 1 - horig;
        if (hlen > oldlen) {                    /* (re)allocate histogram */
                HIST_TYP  *newhist = (HIST_TYP *)calloc(hlen, sizeof(HIST_TYP));
                if (newhist == NULL)
                        returnErr(TM_E_NOMEM);
                if (oldlen) {                   /* copy and free old */
                        memcpy((MEM_PTR)(newhist+(oldorig-horig)),
                                (MEM_PTR)tms->histo, oldlen*sizeof(HIST_TYP));
                        free((MEM_PTR)tms->histo);
                }
                tms->histo = newhist;
        }
        if (wt == 0)
                returnOK;
        for (i = len; i--; )                    /* add in new counts */
                if (ls[i] >= MINBRT)
                        tms->histo[ HISTI(ls[i]) - horig ] += wt;
        returnOK;
}


static double
htcontrs(               /* human threshold contrast sensitivity, dL(La) */
double  La
)
{
        double  l10La, l10dL;
                                /* formula taken from Ferwerda et al. [SG96] */
        if (La < 1.148e-4)
                return(1.38e-3);
        l10La = log10(La);
        if (l10La < -1.44)              /* rod response regime */
                l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
        else if (l10La < -.0184)
                l10dL = l10La - .395;
        else if (l10La < 1.9)           /* cone response regime */
                l10dL = pow(.249*l10La + .65, 2.7) - .72;
        else
                l10dL = l10La - 1.255;

        return(exp10(l10dL));
}


int
tmFixedMapping(                 /* compute fixed, linear tone-mapping */
TMstruct        *tms,
double  expmult,
double  gamval
)
{
        static const char funcName[] = "tmFixedMapping";
        double          d;
        int     i;
        
        if (!tmNewMap(tms))
                returnErr(TM_E_NOMEM);
        if (expmult <= .0)
                expmult = 1.;
        if (gamval < MINGAM)
                gamval = tms->mongam;
        d = log(expmult/tms->inpsf);
        for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
                double  val = TM_BRES * exp(
                        ( d + (tms->mbrmin+i)*(1./TM_BRTSCALE) )
                        / gamval);
                tms->lumap[i] = val;
                if (sizeof(TMAP_TYP) == 2 && val >= 0xffff)
                        tms->lumap[i] = 0xffff;
        }
        returnOK;
}


int
tmComputeMapping(                       /* compute histogram tone-mapping */
TMstruct        *tms,
double  gamval,
double  Lddyn,
double  Ldmax
)
{
        static const char funcName[] = "tmComputeMapping";
        HIST_TYP        *histo;
        float   *cumf;
        int     brt0, histlen;
        HIST_TYP        threshold, ceiling, trimmings, histot;
        double  logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
        double  sum;
        double  d;
        int     i, j;

        if (tms == NULL || tms->histo == NULL)
                returnErr(TM_E_TMINVAL);
                                        /* check arguments */
        if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
        if (Ldmax < MINLDMAX) Ldmax = DEFLDMAX;
        if (gamval < MINGAM) gamval = tms->mongam;
                                        /* compute handy values */
        Ldmin = Ldmax/Lddyn;
        logLddyn = log(Lddyn);
        Ldavg = sqrt(Ldmax*Ldmin);
        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 += (double)(j -= HISTEP) * tms->histo[i];
        }
        threshold = histot*0.005 + .5;
        if (!histot)
                returnErr(TM_E_TMFAIL);
        Lwavg = tmLuminance( (double)sum / histot );
                                        /* use linear tone mapping? */
        if (tms->flags & TM_F_LINEAR || threshold < 4 ||
                        tms->hbrmax - tms->hbrmin < TM_BRTSCALE*logLddyn)
                goto linearmap;
                                        /* clamp histogram */
        histo = (HIST_TYP *)malloc(histlen*sizeof(HIST_TYP));
        cumf = (float *)malloc((histlen+2)*sizeof(float));
        if ((histo == NULL) | (cumf == NULL))
                returnErr(TM_E_NOMEM);
        cumf[histlen+1] = 1.;           /* guard for assignment code */
        for (i = histlen; i--; )        /* make malleable copy */
                histo[i] = tms->histo[i];
        do {                            /* iterate to solution */
                sum = 0;                /* cumulative probability */
                for (i = 0; i < histlen; i++) {
                        cumf[i] = (double)sum/histot;
                        sum += histo[i];
                }
                cumf[histlen] = 1.;
                Tr = histot * (double)(tms->hbrmax - tms->hbrmin) /
                        ((double)histlen*TM_BRTSCALE) / logLddyn;
                ceiling = Tr + 1.;
                trimmings = 0;          /* clip to envelope */
                for (i = histlen; i--; ) {
                        if (tms->flags & TM_F_HCONTR) {
                                Lw = tmLuminance(brt0 + i*HISTEP);
                                Ld = Ldmin * exp( logLddyn *
                                        .5*(cumf[i]+cumf[i+1]) );
                                ceiling = Tr * (htcontrs(Ld) * Lw) /
                                        (htcontrs(Lw) * Ld) + 1.;
                        }
                        if (histo[i] > ceiling) {
                                trimmings += histo[i] - ceiling;
                                histo[i] = ceiling;
                        }
                }
                                        /* check if we're out of data */
                if ((histot -= trimmings) <= threshold) {
                        free((MEM_PTR)histo);
                        free((MEM_PTR)cumf);
                        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] = TM_BRES*pow(d, 1./gamval);
        }
        free((MEM_PTR)histo);           /* clean up and return */
        free((MEM_PTR)cumf);
        returnOK;
linearmap:                              /* linear tone-mapping */
        if (tms->flags & TM_F_HCONTR)
                d = htcontrs(Ldavg) / htcontrs(Lwavg);
        else
                d = Ldavg / Lwavg;
        return(tmFixedMapping(tms, tms->inpsf*d/Ldmax, gamval));
}


int
tmMapPixels(                    /* apply tone-mapping to pixel(s) */
TMstruct        *tms,
uby8    *ps,
TMbright        *ls,
uby8    *cs,
int     len
)
{
        static const char funcName[] = "tmMapPixels";
        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 ((lv = *ls++) < tms->mbrmin) {
                        li = 0;
                } else {
                        if (lv > tms->mbrmax)
                                lv = tms->mbrmax;
                        li = tms->lumap[lv - tms->mbrmin];
                }
                if (cs == TM_NOCHROM)
                        *ps++ = li>=TM_BRES ? 255 : (int)(256*li/TM_BRES);
                else {
                        pv = *cs++ * li / tms->cdiv[RED];
                        *ps++ = pv>255 ? 255 : pv;
                        pv = *cs++ * li / tms->cdiv[GRN];
                        *ps++ = pv>255 ? 255 : pv;
                        pv = *cs++ * li / tms->cdiv[BLU];
                        *ps++ = pv>255 ? 255 : pv;
                }
        }
        returnOK;
}


TMstruct *
tmDup(                          /* duplicate top tone mapping */
TMstruct        *tms
)
{
        int     len;
        int     i;
        TMstruct        *tmnew;

        if (tms == NULL)                /* anything to duplicate? */
                return(NULL);
        tmnew = (TMstruct *)malloc(sizeof(TMstruct));
        if (tmnew == NULL)
                return(NULL);
        *tmnew = *tms;          /* copy everything */
        if (tmnew->histo != NULL) {     /* duplicate histogram */
                len = HISTI(tmnew->hbrmax) + 1 - HISTI(tmnew->hbrmin);
                tmnew->histo = (HIST_TYP *)malloc(len*sizeof(HIST_TYP));
                if (tmnew->histo != NULL)
                        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 = (TMAP_TYP *)malloc(len*sizeof(TMAP_TYP));
                if (tmnew->lumap != NULL)
                        memcpy((MEM_PTR)tmnew->lumap, (MEM_PTR)tms->lumap,
                                        len*sizeof(TMAP_TYP));
        }
                                        /* clear package data */
        for (i = tmNumPkgs; i--; )
                tmnew->pd[i] = NULL;
                                        /* return copy */
        return(tmnew);
}


void
tmDone(tms)                     /* done with tone mapping -- destroy it */
TMstruct        *tms;
{
        int     i;
                                        /* NULL arg. is equiv. to tms */
        if (tms == NULL)
                return;
                                        /* free tables */
        if (tms->histo != NULL)
                free((MEM_PTR)tms->histo);
        if (tms->lumap != NULL)
                free((MEM_PTR)tms->lumap);
                                        /* free private data */
        for (i = tmNumPkgs; i--; )
                if (tms->pd[i] != NULL)
                        (*tmPkg[i]->Free)(tms->pd[i]);
        free((MEM_PTR)tms);             /* free basic structure */
}

/******************** Shared but Private library routines *********************/

uby8    tmMesofact[BMESUPPER-BMESLOWER];

void
tmMkMesofact()                          /* build mesopic lookup factor table */
{
        int     i;

        if (tmMesofact[BMESUPPER-BMESLOWER-1])
                return;

        for (i = BMESLOWER; i < BMESUPPER; i++)
                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);
}


int
tmErrorReturn(                          /* error return (with message) */
const char      *func,
TMstruct        *tms,
int     err
)
{
        if (tms != NULL) {
                tms->lastFunc = func;
                tms->lastError = err;
                if (tms->flags & TM_F_NOSTDERR)
                        return(err);
        }
        fputs(func, stderr);
        fputs(": ", stderr);
        fputs(tmErrorMessage[err], stderr);
        fputs("!\n", stderr);
        return(err);
}
Revision:	3.41
Committed:	Thu Jan 7 19:13:57 2021 UTC (3 years, 4 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.40:	+12 -8 lines
Log Message:	fix: fixed issue with last change and made it so TMAP_TYP can be float
#	User	Rev	Content
1	greg	3.1	#ifndef lint
2	greg	3.41	static const char RCSid[] = "$Id: tonemap.c,v 3.40 2021/01/07 02:13:49 greg Exp $";
3	greg	3.1	#endif
4			/*
5			* Tone mapping functions.
6			* See tonemap.h for detailed function descriptions.
7	greg	3.9	* Added von Kries white-balance calculations 10/01 (GW).
8			*
9			* Externals declared in tonemap.h
10			*/
11
12	greg	3.10	#include "copyright.h"
13	greg	3.1
14			#include <stdio.h>
15	greg	3.34	#include <stdlib.h>
16	greg	3.1	#include <math.h>
17	greg	3.38	#include <string.h>
18	greg	3.1	#include "tmprivat.h"
19			#include "tmerrmsg.h"
20
21			#define exp10(x) exp(M_LN10*(x))
22
23	greg	3.4	/* our list of conversion packages */
24			struct tmPackage *tmPkg[TM_MAXPKG];
25			int tmNumPkgs = 0; /* number of registered packages */
26	greg	3.1
27	greg	3.23	/* luminance->brightness lookup */
28	greg	3.24	static TMbright *tmFloat2BrtLUT = NULL;
29	greg	3.23
30			#define tmCvLumLUfp(pf) tmFloat2BrtLUT[(int32 )(pf) >> 15]
31
32	gregl	3.5
33	greg	3.16	TMstruct *
34			tmInit( /* initialize new tone mapping */
35			int flags,
36			RGBPRIMP monpri,
37			double gamval
38			)
39	greg	3.1	{
40			COLORMAT cmat;
41	greg	3.16	TMstruct *tmnew;
42			int i;
43	greg	3.1	/* allocate structure */
44	greg	3.16	tmnew = (TMstruct *)malloc(sizeof(TMstruct));
45	greg	3.1	if (tmnew == NULL)
46			return(NULL);
47
48			tmnew->flags = flags & ~TM_F_UNIMPL;
49	greg	3.15	if (tmnew->flags & TM_F_BW)
50			tmnew->flags &= ~TM_F_MESOPIC;
51	greg	3.1	/* set monitor transform */
52			if (monpri == NULL \|\| monpri == stdprims \|\| tmnew->flags & TM_F_BW) {
53			tmnew->monpri = stdprims;
54			tmnew->clf[RED] = rgb2xyzmat[1][0];
55			tmnew->clf[GRN] = rgb2xyzmat[1][1];
56			tmnew->clf[BLU] = rgb2xyzmat[1][2];
57			} else {
58	greg	3.26	comprgb2xyzmat(cmat, tmnew->monpri=monpri);
59	greg	3.1	tmnew->clf[RED] = cmat[1][0];
60			tmnew->clf[GRN] = cmat[1][1];
61			tmnew->clf[BLU] = cmat[1][2];
62			}
63			/* set gamma value */
64			if (gamval < MINGAM)
65			tmnew->mongam = DEFGAM;
66			else
67			tmnew->mongam = gamval;
68	greg	3.4	/* set color divisors */
69			for (i = 0; i < 3; i++)
70	greg	3.40	tmnew->cdiv[i] = TM_BRES*pow(tmnew->clf[i], 1./tmnew->mongam);
71	greg	3.4
72	greg	3.1	/* set input transform */
73			tmnew->inppri = tmnew->monpri;
74			tmnew->cmat[0][0] = tmnew->cmat[1][1] = tmnew->cmat[2][2] =
75			tmnew->inpsf = WHTEFFICACY;
76			tmnew->cmat[0][1] = tmnew->cmat[0][2] = tmnew->cmat[1][0] =
77			tmnew->cmat[1][2] = tmnew->cmat[2][0] = tmnew->cmat[2][1] = 0.;
78	greg	3.18	tmnew->inpdat = NULL;
79	greg	3.9	tmnew->hbrmin = 10; tmnew->hbrmax = -10;
80	greg	3.1	tmnew->histo = NULL;
81	greg	3.9	tmnew->mbrmin = 10; tmnew->mbrmax = -10;
82	greg	3.1	tmnew->lumap = NULL;
83	greg	3.4	/* zero private data */
84			for (i = TM_MAXPKG; i--; )
85			tmnew->pd[i] = NULL;
86	greg	3.17	tmnew->lastError = TM_E_OK;
87			tmnew->lastFunc = "NoErr";
88	greg	3.16	/* return new TMstruct */
89			return(tmnew);
90	greg	3.1	}
91
92
93			int
94	greg	3.16	tmSetSpace( /* set input color space for conversions */
95			TMstruct *tms,
96			RGBPRIMP pri,
97	greg	3.18	double sf,
98			MEM_PTR dat
99	greg	3.16	)
100	greg	3.1	{
101	greg	3.17	static const char funcName[] = "tmSetSpace";
102	greg	3.16	int i, j;
103	greg	3.1	/* error check */
104	greg	3.16	if (tms == NULL)
105	greg	3.1	returnErr(TM_E_TMINVAL);
106			if (sf <= 1e-12)
107			returnErr(TM_E_ILLEGAL);
108			/* check if no change */
109	greg	3.18	if (pri == tms->inppri && FEQ(sf, tms->inpsf) && dat == tms->inpdat)
110	greg	3.1	returnOK;
111	greg	3.16	tms->inppri = pri; /* let's set it */
112			tms->inpsf = sf;
113	greg	3.18	tms->inpdat = dat;
114	greg	3.1
115	greg	3.16	if (tms->flags & TM_F_BW) { /* color doesn't matter */
116			tms->monpri = tms->inppri; /* eliminate xform */
117			if (tms->inppri == TM_XYZPRIM) {
118			tms->clf[CIEX] = tms->clf[CIEZ] = 0.;
119			tms->clf[CIEY] = 1.;
120	greg	3.1	} else {
121	greg	3.26	comprgb2xyzmat(tms->cmat, tms->monpri);
122	greg	3.16	tms->clf[RED] = tms->cmat[1][0];
123			tms->clf[GRN] = tms->cmat[1][1];
124			tms->clf[BLU] = tms->cmat[1][2];
125			}
126			tms->cmat[0][0] = tms->cmat[1][1] = tms->cmat[2][2] =
127			tms->inpsf;
128			tms->cmat[0][1] = tms->cmat[0][2] = tms->cmat[1][0] =
129			tms->cmat[1][2] = tms->cmat[2][0] = tms->cmat[2][1] = 0.;
130	greg	3.1
131	greg	3.16	} else if (tms->inppri == TM_XYZPRIM) /* input is XYZ */
132			compxyz2rgbWBmat(tms->cmat, tms->monpri);
133	greg	3.1
134			else { /* input is RGB */
135	greg	3.16	if (tms->inppri != tms->monpri &&
136			PRIMEQ(tms->inppri, tms->monpri))
137			tms->inppri = tms->monpri; /* no xform */
138	greg	3.36	if (!comprgb2rgbWBmat(tms->cmat, tms->inppri, tms->monpri))
139			returnErr(TM_E_ILLEGAL);
140	greg	3.1	}
141			for (i = 0; i < 3; i++)
142			for (j = 0; j < 3; j++)
143	greg	3.16	tms->cmat[i][j] *= tms->inpsf;
144	greg	3.4	/* set color divisors */
145			for (i = 0; i < 3; i++)
146	greg	3.40	tms->cdiv[i] = TM_BRES*pow(tms->clf[i] < .001 ? .001 :
147	greg	3.36	tms->clf[i], 1./tms->mongam);
148	greg	3.4	/* notify packages */
149			for (i = tmNumPkgs; i--; )
150	greg	3.16	if (tms->pd[i] != NULL && tmPkg[i]->NewSpace != NULL)
151			(*tmPkg[i]->NewSpace)(tms);
152	greg	3.1	returnOK;
153			}
154
155
156			void
157	greg	3.23	tmClearHisto( /* clear current histogram */
158	greg	3.16	TMstruct *tms
159			)
160	greg	3.1	{
161	greg	3.16	if (tms == NULL \|\| tms->histo == NULL)
162	greg	3.1	return;
163	greg	3.16	free((MEM_PTR)tms->histo);
164			tms->histo = NULL;
165	greg	3.1	}
166
167
168	greg	3.23	TMbright
169			tmCvLuminance( /* convert a single luminance */
170			double lum
171			)
172			{
173			double d;
174
175			#ifdef isfinite
176			if (!isfinite(lum) \|\| lum <= TM_NOLUM)
177			#else
178			if (lum <= TM_NOLUM)
179			#endif
180			return(TM_NOBRT);
181			d = TM_BRTSCALE*log(lum);
182	greg	3.35	return((TMbright)(d + .5 - (d < 0.)));
183	greg	3.23	}
184
185
186			int
187			tmCvLums( /* convert luminances using lookup */
188			TMbright *ls,
189			float *scan,
190			int len
191			)
192			{
193			if (tmFloat2BrtLUT == NULL) { /* initialize lookup table */
194			int32 i;
195			tmFloat2BrtLUT = (TMbright )malloc(sizeof(TMbright)0x10000);
196			if (tmFloat2BrtLUT == NULL)
197			return(TM_E_NOMEM);
198			for (i = 0; i < 0x10000; i++) {
199			int32 l = (i<<1 \| 1) << 14;
200			#ifndef isfinite
201			if ((l & 0x7f800000) == 0x7f800000)
202			tmFloat2BrtLUT[i] = TM_NOBRT;
203			else
204			#endif
205			tmFloat2BrtLUT[i] = tmCvLuminance((float )&l);
206			}
207			}
208			if (len <= 0)
209			return(TM_E_OK);
210			if ((ls == NULL) \| (scan == NULL))
211			return(TM_E_ILLEGAL);
212			while (len--) {
213			if (*scan <= TM_NOLUM) {
214			*ls++ = TM_NOBRT;
215			++scan;
216			continue;
217			}
218			*ls++ = tmCvLumLUfp(scan++);
219			}
220			return(TM_E_OK);
221			}
222
223
224			int
225			tmCvGrays( /* convert float gray values */
226			TMstruct *tms,
227			TMbright *ls,
228			float *scan,
229			int len
230			)
231			{
232			static const char funcName[] = "tmCvGrays";
233			int i;
234
235			if (tms == NULL)
236			returnErr(TM_E_TMINVAL);
237			if ((ls == NULL) \| (scan == NULL) \| (len < 0))
238			returnErr(TM_E_ILLEGAL);
239	greg	3.24	if (tmFloat2BrtLUT == NULL) /* initialize */
240			tmCvLums(NULL, NULL, 0);
241	greg	3.23	for (i = len; i--; ) {
242			float lum = tms->inpsf * scan[i];
243			if (lum <= TM_NOLUM)
244			ls[i] = TM_NOBRT;
245			else
246			ls[i] = tmCvLumLUfp(&lum);
247			}
248			returnOK;
249			}
250
251
252	greg	3.1	int
253	greg	3.16	tmCvColors( /* convert float colors */
254			TMstruct *tms,
255			TMbright *ls,
256	greg	3.37	uby8 *cs,
257	greg	3.16	COLOR *scan,
258			int len
259			)
260	greg	3.1	{
261	greg	3.17	static const char funcName[] = "tmCvColors";
262	greg	3.37	static uby8 gamtab[1024];
263	greg	3.20	static double curgam = .0;
264	greg	3.1	COLOR cmon;
265	greg	3.23	float lum, slum, d;
266	greg	3.16	int i;
267	greg	3.1
268	greg	3.16	if (tms == NULL)
269	greg	3.1	returnErr(TM_E_TMINVAL);
270	schorsch	3.13	if ((ls == NULL) \| (scan == NULL) \| (len < 0))
271	greg	3.1	returnErr(TM_E_ILLEGAL);
272	greg	3.24	if (tmFloat2BrtLUT == NULL) /* initialize */
273			tmCvLums(NULL, NULL, 0);
274	greg	3.21	if (cs != TM_NOCHROM && fabs(tms->mongam - curgam) > .02) {
275	greg	3.20	curgam = tms->mongam; /* (re)build table */
276			for (i = 1024; i--; )
277			gamtab[i] = (int)(256.*pow((i+.5)/1024., 1./curgam));
278			}
279	greg	3.1	for (i = len; i--; ) {
280	greg	3.16	if (tmNeedMatrix(tms)) { /* get monitor RGB */
281			colortrans(cmon, tms->cmat, scan[i]);
282	greg	3.15	} else {
283	greg	3.16	cmon[RED] = tms->inpsf*scan[i][RED];
284			cmon[GRN] = tms->inpsf*scan[i][GRN];
285			cmon[BLU] = tms->inpsf*scan[i][BLU];
286	greg	3.1	}
287	greg	3.23	#ifdef isfinite
288			if (!isfinite(cmon[RED]) \|\| cmon[RED] < .0f) cmon[RED] = .0f;
289			if (!isfinite(cmon[GRN]) \|\| cmon[GRN] < .0f) cmon[GRN] = .0f;
290			if (!isfinite(cmon[BLU]) \|\| cmon[BLU] < .0f) cmon[BLU] = .0f;
291			#else
292			if (cmon[RED] < .0f) cmon[RED] = .0f;
293			if (cmon[GRN] < .0f) cmon[GRN] = .0f;
294			if (cmon[BLU] < .0f) cmon[BLU] = .0f;
295	greg	3.22	#endif
296	greg	3.1	/* world luminance */
297	greg	3.16	lum = tms->clf[RED]*cmon[RED] +
298			tms->clf[GRN]*cmon[GRN] +
299			tms->clf[BLU]*cmon[BLU] ;
300	greg	3.25	if (lum <= TM_NOLUM) { /* convert brightness */
301			lum = cmon[RED] = cmon[GRN] = cmon[BLU] = TM_NOLUM;
302	greg	3.23	ls[i] = TM_NOBRT;
303	greg	3.25	} else
304	greg	3.23	ls[i] = tmCvLumLUfp(&lum);
305	greg	3.1	if (cs == TM_NOCHROM) /* no color? */
306			continue;
307	greg	3.16	if (tms->flags & TM_F_MESOPIC && lum < LMESUPPER) {
308	greg	3.1	slum = scotlum(cmon); /* mesopic adj. */
309	greg	3.25	if (lum < LMESLOWER) {
310	greg	3.1	cmon[RED] = cmon[GRN] = cmon[BLU] = slum;
311	greg	3.25	} else {
312	greg	3.1	d = (lum - LMESLOWER)/(LMESUPPER - LMESLOWER);
313	greg	3.16	if (tms->flags & TM_F_BW)
314	greg	3.1	cmon[RED] = cmon[GRN] =
315			cmon[BLU] = d*lum;
316			else
317			scalecolor(cmon, d);
318	greg	3.23	d = (1.f-d)*slum;
319	greg	3.1	cmon[RED] += d;
320			cmon[GRN] += d;
321			cmon[BLU] += d;
322			}
323	greg	3.16	} else if (tms->flags & TM_F_BW) {
324	greg	3.1	cmon[RED] = cmon[GRN] = cmon[BLU] = lum;
325			}
326	greg	3.16	d = tms->clf[RED]*cmon[RED]/lum;
327	greg	3.23	cs[3i ] = d>=.999f ? 255 : gamtab[(int)(1024.fd)];
328	greg	3.16	d = tms->clf[GRN]*cmon[GRN]/lum;
329	greg	3.23	cs[3i+1] = d>=.999f ? 255 : gamtab[(int)(1024.fd)];
330	greg	3.16	d = tms->clf[BLU]*cmon[BLU]/lum;
331	greg	3.23	cs[3i+2] = d>=.999f ? 255 : gamtab[(int)(1024.fd)];
332	greg	3.1	}
333			returnOK;
334			}
335
336
337	greg	3.9	int
338	greg	3.16	tmAddHisto( /* add values to histogram */
339			TMstruct *tms,
340			TMbright *ls,
341			int len,
342			int wt
343			)
344	greg	3.1	{
345	greg	3.17	static const char funcName[] = "tmAddHisto";
346	greg	3.9	int oldorig=0, oldlen, horig, hlen;
347	greg	3.16	int i, j;
348	greg	3.1
349	greg	3.16	if (tms == NULL)
350	greg	3.1	returnErr(TM_E_TMINVAL);
351	gwlarson	3.8	if (len < 0)
352			returnErr(TM_E_ILLEGAL);
353			if (len == 0)
354			returnOK;
355	greg	3.1	/* first, grow limits */
356	greg	3.16	if (tms->histo == NULL) {
357	greg	3.1	for (i = len; i-- && ls[i] < MINBRT; )
358			;
359			if (i < 0)
360			returnOK;
361	greg	3.16	tms->hbrmin = tms->hbrmax = ls[i];
362	greg	3.1	oldlen = 0;
363			} else {
364	greg	3.31	oldorig = HISTI(tms->hbrmin);
365			oldlen = HISTI(tms->hbrmax) + 1 - oldorig;
366	greg	3.1	}
367			for (i = len; i--; ) {
368			if ((j = ls[i]) < MINBRT)
369			continue;
370	greg	3.16	if (j < tms->hbrmin)
371			tms->hbrmin = j;
372			else if (j > tms->hbrmax)
373			tms->hbrmax = j;
374	greg	3.1	}
375	greg	3.31	horig = HISTI(tms->hbrmin);
376			hlen = HISTI(tms->hbrmax) + 1 - horig;
377	greg	3.1	if (hlen > oldlen) { /* (re)allocate histogram */
378	greg	3.38	HIST_TYP newhist = (HIST_TYP )calloc(hlen, sizeof(HIST_TYP));
379	greg	3.1	if (newhist == NULL)
380			returnErr(TM_E_NOMEM);
381			if (oldlen) { /* copy and free old */
382	greg	3.39	memcpy((MEM_PTR)(newhist+(oldorig-horig)),
383	greg	3.38	(MEM_PTR)tms->histo, oldlen*sizeof(HIST_TYP));
384	greg	3.16	free((MEM_PTR)tms->histo);
385	greg	3.1	}
386	greg	3.16	tms->histo = newhist;
387	greg	3.1	}
388			if (wt == 0)
389			returnOK;
390			for (i = len; i--; ) /* add in new counts */
391			if (ls[i] >= MINBRT)
392	greg	3.31	tms->histo[ HISTI(ls[i]) - horig ] += wt;
393	greg	3.1	returnOK;
394			}
395
396
397			static double
398	greg	3.16	htcontrs( /* human threshold contrast sensitivity, dL(La) */
399			double La
400			)
401	greg	3.1	{
402			double l10La, l10dL;
403			/* formula taken from Ferwerda et al. [SG96] */
404			if (La < 1.148e-4)
405			return(1.38e-3);
406			l10La = log10(La);
407			if (l10La < -1.44) /* rod response regime */
408			l10dL = pow(.405*l10La + 1.6, 2.18) - 2.86;
409			else if (l10La < -.0184)
410			l10dL = l10La - .395;
411			else if (l10La < 1.9) /* cone response regime */
412			l10dL = pow(.249*l10La + .65, 2.7) - .72;
413			else
414			l10dL = l10La - 1.255;
415
416			return(exp10(l10dL));
417			}
418
419
420	greg	3.9	int
421	greg	3.16	tmFixedMapping( /* compute fixed, linear tone-mapping */
422			TMstruct *tms,
423			double expmult,
424			double gamval
425			)
426	greg	3.9	{
427	greg	3.17	static const char funcName[] = "tmFixedMapping";
428	greg	3.9	double d;
429	greg	3.16	int i;
430	greg	3.9
431	greg	3.16	if (!tmNewMap(tms))
432	greg	3.9	returnErr(TM_E_NOMEM);
433			if (expmult <= .0)
434			expmult = 1.;
435			if (gamval < MINGAM)
436	greg	3.16	gamval = tms->mongam;
437			d = log(expmult/tms->inpsf);
438	greg	3.28	for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
439	greg	3.40	double val = TM_BRES * exp(
440	greg	3.16	( d + (tms->mbrmin+i)*(1./TM_BRTSCALE) )
441	greg	3.28	/ gamval);
442	greg	3.41	tms->lumap[i] = val;
443			if (sizeof(TMAP_TYP) == 2 && val >= 0xffff)
444			tms->lumap[i] = 0xffff;
445	greg	3.28	}
446	greg	3.9	returnOK;
447			}
448
449
450	greg	3.1	int
451	greg	3.16	tmComputeMapping( /* compute histogram tone-mapping */
452			TMstruct *tms,
453			double gamval,
454			double Lddyn,
455			double Ldmax
456			)
457	greg	3.1	{
458	greg	3.17	static const char funcName[] = "tmComputeMapping";
459	greg	3.38	HIST_TYP *histo;
460	greg	3.1	float *cumf;
461	greg	3.38	int brt0, histlen;
462	greg	3.39	HIST_TYP threshold, ceiling, trimmings, histot;
463	greg	3.1	double logLddyn, Ldmin, Ldavg, Lwavg, Tr, Lw, Ld;
464	greg	3.9	double sum;
465	greg	3.16	double d;
466			int i, j;
467	greg	3.1
468	greg	3.16	if (tms == NULL \|\| tms->histo == NULL)
469	greg	3.1	returnErr(TM_E_TMINVAL);
470			/* check arguments */
471			if (Lddyn < MINLDDYN) Lddyn = DEFLDDYN;
472			if (Ldmax < MINLDMAX) Ldmax = DEFLDMAX;
473	greg	3.16	if (gamval < MINGAM) gamval = tms->mongam;
474	greg	3.1	/* compute handy values */
475			Ldmin = Ldmax/Lddyn;
476			logLddyn = log(Lddyn);
477			Ldavg = sqrt(Ldmax*Ldmin);
478	greg	3.31	i = HISTI(tms->hbrmin);
479	greg	3.32	brt0 = HISTV(i);
480	greg	3.31	histlen = HISTI(tms->hbrmax) + 1 - i;
481	greg	3.1	/* histogram total and mean */
482			histot = 0; sum = 0;
483			j = brt0 + histlen*HISTEP;
484			for (i = histlen; i--; ) {
485	greg	3.16	histot += tms->histo[i];
486	greg	3.30	sum += (double)(j -= HISTEP) * tms->histo[i];
487	greg	3.1	}
488	greg	3.14	threshold = histot*0.005 + .5;
489	greg	3.29	if (!histot)
490	greg	3.1	returnErr(TM_E_TMFAIL);
491			Lwavg = tmLuminance( (double)sum / histot );
492	greg	3.29	/* use linear tone mapping? */
493	greg	3.30	if (tms->flags & TM_F_LINEAR \|\| threshold < 4 \|\|
494			tms->hbrmax - tms->hbrmin < TM_BRTSCALE*logLddyn)
495	greg	3.29	goto linearmap;
496	greg	3.9	/* clamp histogram */
497	greg	3.38	histo = (HIST_TYP )malloc(histlensizeof(HIST_TYP));
498	greg	3.9	cumf = (float )malloc((histlen+2)sizeof(float));
499	schorsch	3.13	if ((histo == NULL) \| (cumf == NULL))
500	greg	3.9	returnErr(TM_E_NOMEM);
501			cumf[histlen+1] = 1.; /* guard for assignment code */
502			for (i = histlen; i--; ) /* make malleable copy */
503	greg	3.16	histo[i] = tms->histo[i];
504	greg	3.9	do { /* iterate to solution */
505			sum = 0; /* cumulative probability */
506			for (i = 0; i < histlen; i++) {
507			cumf[i] = (double)sum/histot;
508			sum += histo[i];
509			}
510			cumf[histlen] = 1.;
511	greg	3.16	Tr = histot * (double)(tms->hbrmax - tms->hbrmin) /
512	greg	3.9	((double)histlen*TM_BRTSCALE) / logLddyn;
513			ceiling = Tr + 1.;
514			trimmings = 0; /* clip to envelope */
515			for (i = histlen; i--; ) {
516	greg	3.16	if (tms->flags & TM_F_HCONTR) {
517	greg	3.9	Lw = tmLuminance(brt0 + i*HISTEP);
518			Ld = Ldmin * exp( logLddyn *
519			.5*(cumf[i]+cumf[i+1]) );
520			ceiling = Tr * (htcontrs(Ld) * Lw) /
521			(htcontrs(Lw) * Ld) + 1.;
522	greg	3.2	}
523	greg	3.9	if (histo[i] > ceiling) {
524			trimmings += histo[i] - ceiling;
525			histo[i] = ceiling;
526	greg	3.1	}
527			}
528	greg	3.9	/* check if we're out of data */
529			if ((histot -= trimmings) <= threshold) {
530			free((MEM_PTR)histo);
531			free((MEM_PTR)cumf);
532			goto linearmap;
533			}
534			} while (trimmings > threshold);
535	greg	3.30	/* allocate space for mapping */
536			if (!tmNewMap(tms))
537			returnErr(TM_E_NOMEM);
538	greg	3.9	/* assign tone-mapping */
539	greg	3.16	for (i = tms->mbrmax-tms->mbrmin+1; i--; ) {
540			j = d = (double)i/(tms->mbrmax-tms->mbrmin)*histlen;
541	greg	3.9	d -= (double)j;
542			Ld = Ldminexp(logLddyn((1.-d)cumf[j]+dcumf[j+1]));
543			d = (Ld - Ldmin)/(Ldmax - Ldmin);
544	greg	3.40	tms->lumap[i] = TM_BRES*pow(d, 1./gamval);
545	greg	3.1	}
546	greg	3.9	free((MEM_PTR)histo); /* clean up and return */
547			free((MEM_PTR)cumf);
548	greg	3.1	returnOK;
549	greg	3.9	linearmap: /* linear tone-mapping */
550	greg	3.16	if (tms->flags & TM_F_HCONTR)
551	greg	3.9	d = htcontrs(Ldavg) / htcontrs(Lwavg);
552			else
553			d = Ldavg / Lwavg;
554	greg	3.16	return(tmFixedMapping(tms, tms->inpsf*d/Ldmax, gamval));
555	greg	3.1	}
556
557
558			int
559	greg	3.16	tmMapPixels( /* apply tone-mapping to pixel(s) */
560			TMstruct *tms,
561	greg	3.37	uby8 *ps,
562	greg	3.16	TMbright *ls,
563	greg	3.37	uby8 *cs,
564	greg	3.16	int len
565			)
566	greg	3.1	{
567	greg	3.17	static const char funcName[] = "tmMapPixels";
568	greg	3.41	TMbright lv;
569			TMAP_TYP li;
570			int pv;
571	greg	3.1
572	greg	3.16	if (tms == NULL \|\| tms->lumap == NULL)
573	greg	3.1	returnErr(TM_E_TMINVAL);
574	schorsch	3.13	if ((ps == NULL) \| (ls == NULL) \| (len < 0))
575	greg	3.1	returnErr(TM_E_ILLEGAL);
576			while (len--) {
577	greg	3.41	if ((lv = *ls++) < tms->mbrmin) {
578	greg	3.9	li = 0;
579			} else {
580	greg	3.41	if (lv > tms->mbrmax)
581			lv = tms->mbrmax;
582			li = tms->lumap[lv - tms->mbrmin];
583	greg	3.9	}
584	greg	3.1	if (cs == TM_NOCHROM)
585	greg	3.41	ps++ = li>=TM_BRES ? 255 : (int)(256li/TM_BRES);
586	greg	3.1	else {
587	greg	3.16	pv = cs++ li / tms->cdiv[RED];
588	greg	3.1	*ps++ = pv>255 ? 255 : pv;
589	greg	3.16	pv = cs++ li / tms->cdiv[GRN];
590	greg	3.1	*ps++ = pv>255 ? 255 : pv;
591	greg	3.16	pv = cs++ li / tms->cdiv[BLU];
592	greg	3.1	*ps++ = pv>255 ? 255 : pv;
593			}
594			}
595			returnOK;
596			}
597
598
599	greg	3.16	TMstruct *
600			tmDup( /* duplicate top tone mapping */
601			TMstruct *tms
602			)
603	greg	3.3	{
604			int len;
605	greg	3.16	int i;
606			TMstruct *tmnew;
607	greg	3.3
608	greg	3.16	if (tms == NULL) /* anything to duplicate? */
609	greg	3.3	return(NULL);
610	greg	3.16	tmnew = (TMstruct *)malloc(sizeof(TMstruct));
611	greg	3.3	if (tmnew == NULL)
612			return(NULL);
613	greg	3.16	tmnew = tms; /* copy everything */
614	greg	3.3	if (tmnew->histo != NULL) { /* duplicate histogram */
615	greg	3.31	len = HISTI(tmnew->hbrmax) + 1 - HISTI(tmnew->hbrmin);
616	greg	3.38	tmnew->histo = (HIST_TYP )malloc(lensizeof(HIST_TYP));
617	greg	3.3	if (tmnew->histo != NULL)
618	greg	3.38	memcpy((MEM_PTR)tmnew->histo, (MEM_PTR)tms->histo,
619			len*sizeof(HIST_TYP));
620	greg	3.3	}
621			if (tmnew->lumap != NULL) { /* duplicate luminance mapping */
622	gregl	3.6	len = tmnew->mbrmax-tmnew->mbrmin+1;
623	greg	3.38	tmnew->lumap = (TMAP_TYP )malloc(lensizeof(TMAP_TYP));
624	greg	3.3	if (tmnew->lumap != NULL)
625	greg	3.38	memcpy((MEM_PTR)tmnew->lumap, (MEM_PTR)tms->lumap,
626			len*sizeof(TMAP_TYP));
627	greg	3.3	}
628	greg	3.4	/* clear package data */
629			for (i = tmNumPkgs; i--; )
630			tmnew->pd[i] = NULL;
631	greg	3.16	/* return copy */
632			return(tmnew);
633	greg	3.1	}
634
635
636			void
637			tmDone(tms) /* done with tone mapping -- destroy it */
638	greg	3.16	TMstruct *tms;
639	greg	3.1	{
640	greg	3.16	int i;
641			/* NULL arg. is equiv. to tms */
642			if (tms == NULL)
643	greg	3.1	return;
644			/* free tables */
645			if (tms->histo != NULL)
646	greg	3.4	free((MEM_PTR)tms->histo);
647	greg	3.1	if (tms->lumap != NULL)
648	greg	3.4	free((MEM_PTR)tms->lumap);
649			/* free private data */
650			for (i = tmNumPkgs; i--; )
651			if (tms->pd[i] != NULL)
652			(*tmPkg[i]->Free)(tms->pd[i]);
653			free((MEM_PTR)tms); /* free basic structure */
654	greg	3.9	}
655
656			/****************** Shared but Private library routines *******************/
657
658	greg	3.37	uby8 tmMesofact[BMESUPPER-BMESLOWER];
659	greg	3.9
660			void
661			tmMkMesofact() /* build mesopic lookup factor table */
662			{
663	greg	3.16	int i;
664	greg	3.9
665			if (tmMesofact[BMESUPPER-BMESLOWER-1])
666			return;
667
668			for (i = BMESLOWER; i < BMESUPPER; i++)
669			tmMesofact[i-BMESLOWER] = 256. *
670			(tmLuminance(i) - LMESLOWER) /
671			(LMESUPPER - LMESLOWER);
672			}
673
674
675			int
676	greg	3.33	tmNewMap( /* allocate new tone-mapping array */
677			TMstruct *tms
678			)
679			{
680			if (tms->lumap != NULL && (tms->mbrmax - tms->mbrmin) !=
681			(tms->hbrmax - tms->hbrmin)) {
682			free((MEM_PTR)tms->lumap);
683			tms->lumap = NULL;
684			}
685			tms->mbrmin = tms->hbrmin;
686			tms->mbrmax = tms->hbrmax;
687			if (tms->mbrmin > tms->mbrmax)
688			return 0;
689			if (tms->lumap == NULL)
690	greg	3.38	tms->lumap = (TMAP_TYP )malloc(sizeof(TMAP_TYP)
691	greg	3.33	(tms->mbrmax-tms->mbrmin+1));
692			return(tms->lumap != NULL);
693			}
694
695
696			int
697	greg	3.16	tmErrorReturn( /* error return (with message) */
698	greg	3.17	const char *func,
699	greg	3.16	TMstruct *tms,
700			int err
701			)
702	greg	3.9	{
703	greg	3.17	if (tms != NULL) {
704			tms->lastFunc = func;
705			tms->lastError = err;
706			if (tms->flags & TM_F_NOSTDERR)
707			return(err);
708			}
709	greg	3.9	fputs(func, stderr);
710			fputs(": ", stderr);
711			fputs(tmErrorMessage[err], stderr);
712			fputs("!\n", stderr);
713			return(err);
714	greg	3.1	}