src/rt/data.c

#ifndef lint
static const char       RCSid[] = "$Id: data.c,v 2.36 2023/12/13 23:26:16 greg Exp $";
#endif
/*
 *  data.c - routines dealing with interpolated data.
 */

#include "copyright.h"

#include  <time.h>

#include  "platform.h"
#include  "paths.h"
#include  "standard.h"
#include  "color.h"
#include  "resolu.h"
#include  "view.h"
#include  "data.h"

                                /* picture memory usage before warning */
#ifndef PSIZWARN
#ifdef SMLMEM
#define PSIZWARN        3000000
#else
#define PSIZWARN        50000000
#endif
#endif

#ifndef TABSIZ
#define TABSIZ          997             /* table size (prime) */
#endif

#define hash(s)         (shash(s)%TABSIZ)


static DATARRAY  *dtab[TABSIZ];         /* data array list */

static gethfunc headaspect;


DATARRAY *
getdata(                                /* get data array dname */
        char  *dname
)
{
        char  *dfname;
        FILE  *fp;
        int  asize=0;
        int  i, j;
        DATARRAY  *dp;
                                                /* look for array in list */
        for (dp = dtab[hash(dname)]; dp != NULL; dp = dp->next)
                if (!strcmp(dname, dp->name))
                        return(dp);             /* found! */
        /*
         *      If we haven't loaded the data already, we will look
         *  for it in the directories specified by the library path.
         *
         *      The file has the following format:
         *
         *              N
         *              beg0    end0    n0
         *              beg1    end1    n1
         *              . . .
         *              begN    endN    nN
         *              data, later dimensions changing faster
         *              . . .
         *
         *      For irregularly spaced points, the following can be
         *  substituted for begi endi ni:
         *
         *              0 0 ni p0i p1i .. pni
         */

        if ((dfname = getpath(dname, getrlibpath(), R_OK)) == NULL) {
                sprintf(errmsg, "cannot find data file \"%s\"", dname);
                error(SYSTEM, errmsg);
        }
        if ((fp = fopen(dfname, "r")) == NULL) {
                sprintf(errmsg, "cannot open data file \"%s\"", dfname);
                error(SYSTEM, errmsg);
        }
                                                        /* get dimensions */
        if (fgetval(fp, 'i', &asize) <= 0)
                goto scanerr;
        if ((asize <= 0) | (asize > MAXDDIM)) {
                sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
                error(USER, errmsg);
        }
        if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
                goto memerr;
        dp->name = savestr(dname);
        dp->type = DATATY;
        dp->nd = asize;
        asize = 1;
        for (i = 0; i < dp->nd; i++) {
                if (fgetval(fp, DATATY, &dp->dim[i].org) <= 0)
                        goto scanerr;
                if (fgetval(fp, DATATY, &dp->dim[i].siz) <= 0)
                        goto scanerr;
                if (fgetval(fp, 'i', &dp->dim[i].ne) <= 0)
                        goto scanerr;
                if (dp->dim[i].ne < 2)
                        goto scanerr;
                asize *= dp->dim[i].ne;
                if ((dp->dim[i].siz -= dp->dim[i].org) == 0) {
                        dp->dim[i].p = (DATATYPE *)
                                        malloc(dp->dim[i].ne*sizeof(DATATYPE));
                        if (dp->dim[i].p == NULL)
                                goto memerr;
                        for (j = 0; j < dp->dim[i].ne; j++)
                                if (fgetval(fp, DATATY, &dp->dim[i].p[j]) <= 0)
                                        goto scanerr;
                        for (j = 1; j < dp->dim[i].ne-1; j++)
                                if ((dp->dim[i].p[j-1] < dp->dim[i].p[j]) !=
                                        (dp->dim[i].p[j] < dp->dim[i].p[j+1]))
                                        goto scanerr;
                        dp->dim[i].org = dp->dim[i].p[0];
                        dp->dim[i].siz = dp->dim[i].p[dp->dim[i].ne-1]
                                                - dp->dim[i].p[0];
                } else
                        dp->dim[i].p = NULL;
        }
        if ((dp->arr.d = (DATATYPE *)malloc(asize*sizeof(DATATYPE))) == NULL)
                goto memerr;
        
        for (i = 0; i < asize; i++)
                if (fgetval(fp, DATATY, &dp->arr.d[i]) <= 0)
                        goto scanerr;
        fclose(fp);
        i = hash(dname);
        dp->next = dtab[i];
        return(dtab[i] = dp);
memerr:
        error(SYSTEM, "out of memory in getdata");
scanerr:
        sprintf(errmsg, "%s in data file \"%s\"",
                        feof(fp) ? "unexpected EOF" : "bad format", dfname);
        error(USER, errmsg);
        return NULL;    /* pro forma return */
}


static int
headaspect(                     /* check string for aspect ratio */
        char  *s,
        void  *iap
)
{
        char    fmt[MAXFMTLEN];

        if (isaspect(s))
                *(double*)iap *= aspectval(s);
        else if (formatval(fmt, s) && strcmp(fmt, COLRFMT))
                *(double*)iap = 0.0;
        return(0);
}

DATARRAY *
getpict(                                /* get picture pname */
        char  *pname
)
{
        double  inpaspect;
        char  *pfname;
        FILE  *fp;
        COLR  *scanin;
        int  sl, ns;
        RESOLU  inpres;
        RREAL  loc[2];
        int  y;
        int  x, i;
        DATARRAY  *pp;
                                                /* look for array in list */
        for (pp = dtab[hash(pname)]; pp != NULL; pp = pp->next)
                if (!strcmp(pname, pp->name))
                        return(pp);             /* found! */

        if ((pfname = getpath(pname, getrlibpath(), R_OK)) == NULL) {
                sprintf(errmsg, "cannot find picture file \"%s\"", pname);
                error(SYSTEM, errmsg);
        }
        if ((pp = (DATARRAY *)malloc(3*sizeof(DATARRAY))) == NULL)
                goto memerr;

        pp[0].name = savestr(pname);

        if ((fp = fopen(pfname, "rb")) == NULL) {
                sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
                error(SYSTEM, errmsg);
        }
                                                /* get dimensions */
        inpaspect = 1.0;
        getheader(fp, headaspect, &inpaspect);
        if (inpaspect <= FTINY || !fgetsresolu(&inpres, fp))
                goto readerr;
        pp[0].nd = 2;
        pp[0].dim[0].ne = inpres.yr;
        pp[0].dim[1].ne = inpres.xr;
        pp[0].dim[0].org =
        pp[0].dim[1].org = 0.0;
        if (inpres.xr <= inpres.yr*inpaspect) {
                pp[0].dim[0].siz = inpaspect *
                                        (double)inpres.yr/inpres.xr;
                pp[0].dim[1].siz = 1.0;
        } else {
                pp[0].dim[0].siz = 1.0;
                pp[0].dim[1].siz = (double)inpres.xr/inpres.yr /
                                        inpaspect;
        }
        pp[0].dim[0].p = pp[0].dim[1].p = NULL;
        sl = scanlen(&inpres);                          /* allocate array */
        ns = numscans(&inpres);
        i = ns*sl*sizeof(COLR);
#if PSIZWARN
        if (i > PSIZWARN) {                             /* memory warning */
                sprintf(errmsg, "picture file \"%s\" using %.1f MB of memory",
                                pname, i*(1.0/(1024*1024)));
                error(WARNING, errmsg);
        }
#endif
        if ((pp[0].arr.c = (COLR *)malloc(i)) == NULL)
                goto memerr;
                                                        /* load picture */
        if ((scanin = (COLR *)malloc(sl*sizeof(COLR))) == NULL)
                goto memerr;
        for (y = 0; y < ns; y++) {
                if (freadcolrs(scanin, sl, fp) < 0)
                        goto readerr;
                for (x = 0; x < sl; x++) {
                        pix2loc(loc, &inpres, x, y);
                        i = (int)(loc[1]*inpres.yr)*inpres.xr +
                                        (int)(loc[0]*inpres.xr);
                        copycolr(pp[0].arr.c[i], scanin[x]);
                }
        }
        free(scanin);
        fclose(fp);
        i = hash(pname);
        pp[0].next = dtab[i];           /* link into picture list */
        pp[1] = pp[0];
        pp[2] = pp[0];
        pp[0].type = RED;               /* differentiate RGB records */
        pp[1].type = GRN;
        pp[2].type = BLU;
        return(dtab[i] = pp);
memerr:
        error(SYSTEM, "out of memory in getpict");
readerr:
        sprintf(errmsg, "bad picture file \"%s\"", pfname);
        error(USER, errmsg);
        return NULL;    /* pro forma return */
}


/* header info type for hyperspectral image */
typedef struct {
        float   wlpart[4];      /* wavelength partitions */
        int     nc;             /* number of components */
        double  inpaspect;      /* pixel aspect ratio */
} SPECINFO;

static int
specheadline(                           /* get info for spectral image */
        char  *s,
        void  *cdp
)
{
        SPECINFO        *sip = (SPECINFO *)cdp;
        char            fmt[MAXFMTLEN];

        if (isaspect(s))
                sip->inpaspect *= aspectval(s);
        else if (isncomp(s))
                sip->nc = ncompval(s);
        else if (iswlsplit(s))
                wlsplitval(sip->wlpart, s);
        else if (formatval(fmt, s) && strcmp(fmt, SPECFMT))
                return(-1);
        return(0);
}

DATARRAY *
getspec(                /* load hyperspectral image as data */
        char *sname
)
{
        SPECINFO        si;
        char            *pfname;
        FILE            *fp;
        int             sl, ns;
        int             y, i;
        DATARRAY        *pp;
                                                /* look for array in list */
        for (pp = dtab[hash(sname)]; pp != NULL; pp = pp->next)
                if (!strcmp(sname, pp->name))
                        return(pp);             /* found! */

        if ((pfname = getpath(sname, getrlibpath(), R_OK)) == NULL) {
                sprintf(errmsg, "cannot find hyperspectral image \"%s\"", sname);
                error(SYSTEM, errmsg);
        }
        if ((fp = fopen(pfname, "rb")) == NULL) {
                sprintf(errmsg, "cannot open hyperspectral image \"%s\"", pfname);
                error(SYSTEM, errmsg);
        }
        si.wlpart[3] = 0;
        si.nc = 0;
        si.inpaspect = 1.0;
        if (getheader(fp, specheadline, &si) < 0 ||
                        (si.nc <= 3) | (si.nc > MAXCSAMP) | (si.wlpart[3] < 1) ||
                        !fscnresolu(&sl, &ns, fp))
                goto readerr;

        if ((pp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
                goto memerr;

        pp->name = savestr(sname);
        pp->type = SPECTY;
        pp->nd = 3;
        pp->dim[0].ne = ns;
        pp->dim[1].ne = sl;
        pp->dim[0].org =
        pp->dim[1].org = 0.0;
        if (sl <= ns*si.inpaspect) {
                pp->dim[0].siz = si.inpaspect * (double)ns/sl;
                pp->dim[1].siz = 1.0;
        } else {
                pp->dim[0].siz = 1.0;
                pp->dim[1].siz = (double)sl/ns / si.inpaspect;
        }
        pp->dim[2].ne = si.nc;
        pp->dim[2].siz = si.wlpart[3] - si.wlpart[0];
        pp->dim[2].org = si.wlpart[0] + 0.5*pp->dim[2].siz/si.nc;
        pp->dim[2].siz *= (si.nc - 1.0)/si.nc;
        pp->dim[0].p = pp->dim[1].p = pp->dim[2].p = NULL;
        i = ns*sl*(si.nc+1);
#if PSIZWARN
        if (i > PSIZWARN) {                     /* memory warning */
                sprintf(errmsg, "hyperspectral image \"%s\" using %.1f MB of memory",
                                sname, i*(1.0/(1024*1024)));
                error(WARNING, errmsg);
        }
#endif
        if ((pp->arr.s = (uby8 *)malloc(i)) == NULL)
                goto memerr;
        for (y = 0; y < ns; y++)                /* read each scanline */
                if (freadscolrs(pp->arr.s + y*sl*(si.nc+1), si.nc, sl, fp) < 0)
                        goto readerr;
        fclose(fp);
        i = hash(sname);                        /* insert in hash table */
        pp->next = dtab[i];
        return(dtab[i] = pp);
memerr:
        error(SYSTEM, "out of memory in getspec");
readerr:
        sprintf(errmsg, "bad hyperspectral image \"%s\"", pfname);
        error(USER, errmsg);
        return NULL;    /* pro forma return */
}


void
freedata(                       /* release data array reference */
        DATARRAY  *dta
)
{
        DATARRAY  head;
        int  hval, nents;
        DATARRAY  *dpl, *dp;
        int  i;

        if (dta == NULL) {                      /* free all if NULL */
                hval = 0; nents = TABSIZ;
        } else {
                hval = hash(dta->name); nents = 1;
                if (!*dta->name) {              /* not a data file? */
                        dta->next = dtab[hval];
                        dtab[hval] = dta;       /* ...fake position */
                }
        }
        while (nents--) {
                head.next = dtab[hval];
                dpl = &head;
                while ((dp = dpl->next) != NULL)
                        if ((dta == NULL) | (dta == dp)) {
                                dpl->next = dp->next;
                                free(dp->arr.p);
                                for (i = 0; i < dp->nd; i++)
                                        if (dp->dim[i].p != NULL)
                                                free(dp->dim[i].p);
                                freestr(dp->name);
                                free(dp);
                        } else
                                dpl = dp;
                dtab[hval++] = head.next;
        }
}


/* internal call to interpolate data value or vector */
static double
data_interp(DATARRAY *dp, double *pt, double coef, DATATYPE *rvec)
{
        DATARRAY        sd;
        int             stride, i;
        double          x, c0, c1, y0, y1;
                                        /* set up dimensions for recursion */
        if (dp->nd > 1) {
                sd.name = dp->name;
                sd.type = dp->type;
                sd.nd = dp->nd - 1;
                memcpy(sd.dim, dp->dim+1, sd.nd*sizeof(struct dadim));
                stride = sd.dim[i = sd.nd-1].ne + (sd.type==SPECTY);
                while (i-- > 0)
                        stride *= sd.dim[i].ne;
        }
                                        /* get independent variable */
        if (dp->dim[0].p == NULL) {             /* evenly spaced points */
                x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
                x *= (double)(dp->dim[0].ne - 1);
                i = x;
                if (i < 0)
                        i = 0;
                else if (i > dp->dim[0].ne - 2)
                        i = dp->dim[0].ne - 2;
        } else {                                /* unevenly spaced points */
                int     lower, upper;
                if (dp->dim[0].siz > 0.0) {
                        lower = 0;
                        upper = dp->dim[0].ne;
                } else {
                        lower = dp->dim[0].ne;
                        upper = 0;
                }
                do {
                        i = (lower + upper) >> 1;
                        if (pt[0] >= dp->dim[0].p[i])
                                lower = i;
                        else
                                upper = i;
                } while (i != (lower + upper) >> 1);

                if (i > dp->dim[0].ne - 2)
                        i = dp->dim[0].ne - 2;

                x = i + (pt[0] - dp->dim[0].p[i]) /
                                (dp->dim[0].p[i+1] - dp->dim[0].p[i]);
        }
        /*
         * Compute interpolation coefficients:
         * extrapolate as far as one division, then
         * taper off harmonically to zero.
         */
        if (x > i+2) {
                c0 = 1./(i-1 - x);
                c1 = -2.*c0;
        } else if (x < i-1) {
                c1 = 1./(i - x);
                c0 = -2.*c1;
        } else {
                c0 = i+1 - x;
                c1 = x - i;
        }
        c0 *= coef;
        c1 *= coef;
                                        /* check if vector interp */
        if ((dp->nd == 2) & (rvec != NULL)) {
                if (dp->type == DATATY) {
                        sd.arr.d = dp->arr.d + i*stride;
                        for (i = sd.dim[0].ne; i--; )
                                rvec[i] += c0*sd.arr.d[i]
                                        + c1*sd.arr.d[i+stride];
                        return(0.);
                }
                if (dp->type == SPECTY) {
                        double  f;
                        sd.arr.s = dp->arr.s + i*stride;
                        f = ldexp(1.0, (int)sd.arr.s[sd.dim[0].ne]
                                        - (COLXS+8));
                        for (i = sd.dim[0].ne; i--; )
                                rvec[i] += c0*f*(sd.arr.s[i] + 0.5);
                        sd.arr.s += stride;
                        f = ldexp(1.0, (int)sd.arr.s[sd.dim[0].ne]
                                        - (COLXS+8));
                        for (i = sd.dim[0].ne; i--; )
                                rvec[i] += c1*f*(sd.arr.s[i] + 0.5);
                        return(0.);
                }
                sd.arr.c = dp->arr.c + i*stride;
                for (i = sd.dim[0].ne; i--; )
                        rvec[i] += c0*colrval(sd.arr.c[i],sd.type)
                                + c1*colrval(sd.arr.c[i+stride],sd.type);
                return(0.);
        }
                                        /* get dependent variable */
        if (dp->nd > 1) {
                if (dp->type == DATATY) {
                        sd.arr.d = dp->arr.d + i*stride;
                        y0 = data_interp(&sd, pt+1, c0, rvec);
                        sd.arr.d += stride;
                } else if (dp->type == SPECTY) {
                        sd.arr.s = dp->arr.s + i*stride;
                        y0 = data_interp(&sd, pt+1, c0, rvec);
                        sd.arr.s += stride;
                } else {
                        sd.arr.c = dp->arr.c + i*stride;
                        y0 = data_interp(&sd, pt+1, c0, rvec);
                        sd.arr.c += stride;
                }
                y1 = data_interp(&sd, pt+1, c1, rvec);
        } else {                        /* end of recursion */
                if (dp->type == DATATY) {
                        y0 = dp->arr.d[i];
                        y1 = dp->arr.d[i+1];
                } else if (dp->type == SPECTY) {
                        if (dp->arr.s[dp->dim[0].ne]) {
                                double  f = ldexp(1.0, -(COLXS+8) +
                                                (int)dp->arr.s[dp->dim[0].ne]);
                                y0 = (dp->arr.s[i] + 0.5)*f;
                                y1 = (dp->arr.s[i+1] + 0.5)*f;
                        } else
                                y0 = y1 = 0.0;
                } else {
                        y0 = colrval(dp->arr.c[i],dp->type);
                        y1 = colrval(dp->arr.c[i+1],dp->type);
                }
                y0 *= c0;
                y1 *= c1;
        }
        return(y0 + y1);        /* coefficients already applied */
}


double
datavalue(              /* interpolate data value at a point */
        DATARRAY  *dp,
        double  *pt
)
{
        return(data_interp(dp, pt, 1., NULL));
}


/* Interpolate final vector corresponding to last dimension in data array */
DATARRAY *
datavector(DATARRAY *dp, double *pt)
{
        DATARRAY        *newdp;

        if (dp->nd < 2)
                error(INTERNAL, "datavector() called with 1-D array");
                                        /* create vector array */
        newdp = (DATARRAY *)malloc(sizeof(DATARRAY) -
                                (MAXDDIM-1)*sizeof(struct dadim) +
                                sizeof(DATATYPE)*dp->dim[dp->nd-1].ne);
        if (newdp == NULL)
                error(SYSTEM, "out of memory in datavector");
        newdp->next = NULL;
        newdp->name = dp->name;
        newdp->type = DATATY;
        newdp->nd = 1;                  /* vector data goes here */
        newdp->dim[0] = dp->dim[dp->nd-1];
        newdp->arr.d = (DATATYPE *)(newdp->dim + 1);
        memset(newdp->arr.d, 0, sizeof(DATATYPE)*newdp->dim[0].ne);

        (void)data_interp(dp, pt, 1., newdp->arr.d);

        return(newdp);                  /* will be free'd using free() */
}


#if 0
double
datavalue(              /* interpolate data value at a point */
        DATARRAY  *dp,
        double  *pt
)
{
        DATARRAY  sd;
        int  asize;
        int  lower, upper;
        int  i;
        double  x, y0, y1;
                                        /* set up dimensions for recursion */
        if (dp->nd > 1) {
                sd.name = dp->name;
                sd.type = dp->type;
                sd.nd = dp->nd - 1;
                asize = 1;
                for (i = 0; i < sd.nd; i++) {
                        sd.dim[i].org = dp->dim[i+1].org;
                        sd.dim[i].siz = dp->dim[i+1].siz;
                        sd.dim[i].p = dp->dim[i+1].p;
                        asize *= (sd.dim[i].ne = dp->dim[i+1].ne) +
                                ((sd.type==SPECTY) & (i==sd.nd-1));
                }
        }
                                        /* get independent variable */
        if (dp->dim[0].p == NULL) {             /* evenly spaced points */
                x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
                x *= (double)(dp->dim[0].ne - 1);
                i = x;
                if (i < 0)
                        i = 0;
                else if (i > dp->dim[0].ne - 2)
                        i = dp->dim[0].ne - 2;
        } else {                                /* unevenly spaced points */
                if (dp->dim[0].siz > 0.0) {
                        lower = 0;
                        upper = dp->dim[0].ne;
                } else {
                        lower = dp->dim[0].ne;
                        upper = 0;
                }
                do {
                        i = (lower + upper) >> 1;
                        if (pt[0] >= dp->dim[0].p[i])
                                lower = i;
                        else
                                upper = i;
                } while (i != (lower + upper) >> 1);

                if (i > dp->dim[0].ne - 2)
                        i = dp->dim[0].ne - 2;

                x = i + (pt[0] - dp->dim[0].p[i]) /
                                (dp->dim[0].p[i+1] - dp->dim[0].p[i]);
        }
                                        /* get dependent variable */
        if (dp->nd > 1) {
                if (dp->type == DATATY) {
                        sd.arr.d = dp->arr.d + i*asize;
                        y0 = datavalue(&sd, pt+1);
                        sd.arr.d += asize;
                        y1 = datavalue(&sd, pt+1);
                } else if (dp->type == SPECTY) {
                        sd.arr.s = dp->arr.s + i*asize;
                        y0 = datavalue(&sd, pt+1);
                        sd.arr.s += asize;
                        y1 = datavalue(&sd, pt+1);
                } else {
                        sd.arr.c = dp->arr.c + i*asize;
                        y0 = datavalue(&sd, pt+1);
                        sd.arr.c += asize;
                        y1 = datavalue(&sd, pt+1);
                }
        } else {
                if (dp->type == DATATY) {
                        y0 = dp->arr.d[i];
                        y1 = dp->arr.d[i+1];
                } else if (dp->type == SPECTY) {
                        if (dp->arr.s[dp->dim[0].ne]) {
                                double  f = ldexp(1.0, -(COLXS+8) +
                                                (int)dp->arr.s[dp->dim[0].ne]);
                                y0 = (dp->arr.s[i] + 0.5)*f;
                                y1 = (dp->arr.s[i+1] + 0.5)*f;
                        } else
                                y0 = y1 = 0.0;
                } else {
                        y0 = colrval(dp->arr.c[i],dp->type);
                        y1 = colrval(dp->arr.c[i+1],dp->type);
                }
        }
        /*
         * Extrapolate as far as one division, then
         * taper off harmonically to zero.
         */
        if (x > i+2)
                return( (2*y1-y0)/(x-(i-1)) );

        if (x < i-1)
                return( (2*y0-y1)/(i-x) );

        return( y0*((i+1)-x) + y1*(x-i) );
}
#endif
Revision:	2.37
Committed:	Tue Mar 12 16:54:51 2024 UTC (19 months, 1 week ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.36:	+179 -1 lines
Log Message:	perf: added datavector() call for quicker spectral interpolation
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.37	static const char RCSid[] = "$Id: data.c,v 2.36 2023/12/13 23:26:16 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* data.c - routines dealing with interpolated data.
6			*/
7
8	greg	2.16	#include "copyright.h"
9	greg	2.22
10			#include <time.h>
11	greg	2.15
12	greg	2.25	#include "platform.h"
13	schorsch	2.28	#include "paths.h"
14	greg	1.1	#include "standard.h"
15			#include "color.h"
16	greg	1.11	#include "resolu.h"
17	schorsch	2.27	#include "view.h"
18	greg	1.1	#include "data.h"
19
20	greg	2.7	/* picture memory usage before warning */
21			#ifndef PSIZWARN
22	greg	2.18	#ifdef SMLMEM
23	greg	2.29	#define PSIZWARN 3000000
24	greg	2.18	#else
25	greg	2.34	#define PSIZWARN 50000000
26	greg	2.7	#endif
27			#endif
28	greg	1.1
29	greg	2.7	#ifndef TABSIZ
30	greg	2.35	#define TABSIZ 997 /* table size (prime) */
31	greg	2.7	#endif
32	greg	2.6
33			#define hash(s) (shash(s)%TABSIZ)
34
35
36			static DATARRAY dtab[TABSIZ]; / data array list */
37	greg	1.1
38	schorsch	2.26	static gethfunc headaspect;
39
40	greg	1.1
41	greg	2.31	DATARRAY *
42	schorsch	2.27	getdata( /* get data array dname */
43			char *dname
44			)
45	greg	1.1	{
46			char *dfname;
47			FILE *fp;
48	greg	2.30	int asize=0;
49	greg	2.31	int i, j;
50			DATARRAY *dp;
51	greg	1.1	/* look for array in list */
52	greg	2.6	for (dp = dtab[hash(dname)]; dp != NULL; dp = dp->next)
53	greg	1.1	if (!strcmp(dname, dp->name))
54			return(dp); /* found! */
55			/*
56			* If we haven't loaded the data already, we will look
57	greg	2.5	* for it in the directories specified by the library path.
58	greg	1.1	*
59			* The file has the following format:
60			*
61	greg	1.4	* N
62	greg	1.1	* beg0 end0 n0
63			* beg1 end1 n1
64			* . . .
65	greg	1.4	* begN endN nN
66	greg	1.1	* data, later dimensions changing faster
67			* . . .
68			*
69	greg	1.4	* For irregularly spaced points, the following can be
70			* substituted for begi endi ni:
71			*
72	greg	1.6	* 0 0 ni p0i p1i .. pni
73	greg	1.1	*/
74
75	greg	2.17	if ((dfname = getpath(dname, getrlibpath(), R_OK)) == NULL) {
76	greg	1.1	sprintf(errmsg, "cannot find data file \"%s\"", dname);
77	greg	2.32	error(SYSTEM, errmsg);
78	greg	1.1	}
79			if ((fp = fopen(dfname, "r")) == NULL) {
80			sprintf(errmsg, "cannot open data file \"%s\"", dfname);
81			error(SYSTEM, errmsg);
82			}
83			/* get dimensions */
84	greg	2.36	if (fgetval(fp, 'i', &asize) <= 0)
85	greg	1.1	goto scanerr;
86	schorsch	2.24	if ((asize <= 0) \| (asize > MAXDDIM)) {
87	greg	1.1	sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
88			error(USER, errmsg);
89			}
90	greg	2.12	if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
91			goto memerr;
92			dp->name = savestr(dname);
93			dp->type = DATATY;
94			dp->nd = asize;
95	greg	1.1	asize = 1;
96			for (i = 0; i < dp->nd; i++) {
97	greg	2.36	if (fgetval(fp, DATATY, &dp->dim[i].org) <= 0)
98	greg	1.6	goto scanerr;
99	greg	2.36	if (fgetval(fp, DATATY, &dp->dim[i].siz) <= 0)
100	greg	1.10	goto scanerr;
101	greg	2.36	if (fgetval(fp, 'i', &dp->dim[i].ne) <= 0)
102	greg	1.10	goto scanerr;
103	greg	1.6	if (dp->dim[i].ne < 2)
104			goto scanerr;
105			asize *= dp->dim[i].ne;
106			if ((dp->dim[i].siz -= dp->dim[i].org) == 0) {
107	greg	2.12	dp->dim[i].p = (DATATYPE *)
108			malloc(dp->dim[i].ne*sizeof(DATATYPE));
109	greg	1.4	if (dp->dim[i].p == NULL)
110			goto memerr;
111	greg	2.10	for (j = 0; j < dp->dim[i].ne; j++)
112	greg	2.36	if (fgetval(fp, DATATY, &dp->dim[i].p[j]) <= 0)
113	greg	1.4	goto scanerr;
114			for (j = 1; j < dp->dim[i].ne-1; j++)
115			if ((dp->dim[i].p[j-1] < dp->dim[i].p[j]) !=
116			(dp->dim[i].p[j] < dp->dim[i].p[j+1]))
117			goto scanerr;
118			dp->dim[i].org = dp->dim[i].p[0];
119			dp->dim[i].siz = dp->dim[i].p[dp->dim[i].ne-1]
120			- dp->dim[i].p[0];
121			} else
122	greg	1.6	dp->dim[i].p = NULL;
123	greg	1.1	}
124	greg	2.12	if ((dp->arr.d = (DATATYPE )malloc(asizesizeof(DATATYPE))) == NULL)
125	greg	1.1	goto memerr;
126
127	greg	2.10	for (i = 0; i < asize; i++)
128	greg	2.36	if (fgetval(fp, DATATY, &dp->arr.d[i]) <= 0)
129	greg	1.1	goto scanerr;
130			fclose(fp);
131	greg	2.6	i = hash(dname);
132			dp->next = dtab[i];
133			return(dtab[i] = dp);
134	greg	1.1	memerr:
135			error(SYSTEM, "out of memory in getdata");
136			scanerr:
137			sprintf(errmsg, "%s in data file \"%s\"",
138			feof(fp) ? "unexpected EOF" : "bad format", dfname);
139			error(USER, errmsg);
140	greg	2.36	return NULL; /* pro forma return */
141	greg	1.1	}
142
143
144	greg	2.15	static int
145	schorsch	2.26	headaspect( /* check string for aspect ratio */
146			char *s,
147			void *iap
148			)
149	greg	1.5	{
150	greg	2.33	char fmt[MAXFMTLEN];
151	gwlarson	2.14
152	greg	1.5	if (isaspect(s))
153	schorsch	2.26	(double)iap *= aspectval(s);
154	greg	2.36	else if (formatval(fmt, s) && strcmp(fmt, COLRFMT))
155	schorsch	2.26	(double)iap = 0.0;
156	gwlarson	2.14	return(0);
157	greg	1.5	}
158
159	greg	2.31	DATARRAY *
160	schorsch	2.27	getpict( /* get picture pname */
161			char *pname
162			)
163	greg	1.1	{
164	greg	2.4	double inpaspect;
165	greg	1.1	char *pfname;
166			FILE *fp;
167	greg	2.12	COLR *scanin;
168	greg	1.11	int sl, ns;
169	greg	2.3	RESOLU inpres;
170	schorsch	2.21	RREAL loc[2];
171	greg	1.11	int y;
172	greg	2.31	int x, i;
173			DATARRAY *pp;
174	greg	1.1	/* look for array in list */
175	greg	2.12	for (pp = dtab[hash(pname)]; pp != NULL; pp = pp->next)
176	greg	1.1	if (!strcmp(pname, pp->name))
177			return(pp); /* found! */
178
179	greg	2.17	if ((pfname = getpath(pname, getrlibpath(), R_OK)) == NULL) {
180	greg	1.1	sprintf(errmsg, "cannot find picture file \"%s\"", pname);
181	greg	2.32	error(SYSTEM, errmsg);
182	greg	1.1	}
183	greg	2.12	if ((pp = (DATARRAY )malloc(3sizeof(DATARRAY))) == NULL)
184	greg	1.1	goto memerr;
185
186	greg	2.12	pp[0].name = savestr(pname);
187	greg	1.1
188	greg	2.36	if ((fp = fopen(pfname, "rb")) == NULL) {
189	greg	1.1	sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
190			error(SYSTEM, errmsg);
191			}
192			/* get dimensions */
193	greg	1.5	inpaspect = 1.0;
194	schorsch	2.26	getheader(fp, headaspect, &inpaspect);
195	gwlarson	2.14	if (inpaspect <= FTINY \|\| !fgetsresolu(&inpres, fp))
196	greg	1.1	goto readerr;
197	greg	2.12	pp[0].nd = 2;
198			pp[0].dim[0].ne = inpres.yr;
199			pp[0].dim[1].ne = inpres.xr;
200			pp[0].dim[0].org =
201			pp[0].dim[1].org = 0.0;
202			if (inpres.xr <= inpres.yr*inpaspect) {
203			pp[0].dim[0].siz = inpaspect *
204			(double)inpres.yr/inpres.xr;
205			pp[0].dim[1].siz = 1.0;
206			} else {
207			pp[0].dim[0].siz = 1.0;
208			pp[0].dim[1].siz = (double)inpres.xr/inpres.yr /
209			inpaspect;
210			}
211			pp[0].dim[0].p = pp[0].dim[1].p = NULL;
212			sl = scanlen(&inpres); /* allocate array */
213			ns = numscans(&inpres);
214			i = nsslsizeof(COLR);
215	greg	2.7	#if PSIZWARN
216	greg	2.12	if (i > PSIZWARN) { /* memory warning */
217	greg	2.29	sprintf(errmsg, "picture file \"%s\" using %.1f MB of memory",
218			pname, i(1.0/(10241024)));
219	greg	2.7	error(WARNING, errmsg);
220			}
221			#endif
222	greg	2.12	if ((pp[0].arr.c = (COLR *)malloc(i)) == NULL)
223			goto memerr;
224	greg	1.1	/* load picture */
225	greg	2.12	if ((scanin = (COLR )malloc(slsizeof(COLR))) == NULL)
226	greg	1.1	goto memerr;
227	greg	1.11	for (y = 0; y < ns; y++) {
228	greg	2.12	if (freadcolrs(scanin, sl, fp) < 0)
229	greg	1.1	goto readerr;
230	greg	1.11	for (x = 0; x < sl; x++) {
231			pix2loc(loc, &inpres, x, y);
232			i = (int)(loc[1]inpres.yr)inpres.xr +
233			(int)(loc[0]*inpres.xr);
234	greg	2.12	copycolr(pp[0].arr.c[i], scanin[x]);
235	greg	1.11	}
236	greg	1.1	}
237	greg	2.36	free(scanin);
238	greg	1.1	fclose(fp);
239	greg	2.6	i = hash(pname);
240	greg	2.12	pp[0].next = dtab[i]; /* link into picture list */
241	schorsch	2.23	pp[1] = pp[0];
242			pp[2] = pp[0];
243	greg	2.12	pp[0].type = RED; /* differentiate RGB records */
244			pp[1].type = GRN;
245			pp[2].type = BLU;
246			return(dtab[i] = pp);
247	greg	1.1	memerr:
248			error(SYSTEM, "out of memory in getpict");
249			readerr:
250			sprintf(errmsg, "bad picture file \"%s\"", pfname);
251			error(USER, errmsg);
252	greg	2.36	return NULL; /* pro forma return */
253			}
254
255
256			/* header info type for hyperspectral image */
257			typedef struct {
258			float wlpart[4]; /* wavelength partitions */
259			int nc; /* number of components */
260			double inpaspect; /* pixel aspect ratio */
261			} SPECINFO;
262
263			static int
264			specheadline( /* get info for spectral image */
265			char *s,
266			void *cdp
267			)
268			{
269			SPECINFO sip = (SPECINFO )cdp;
270			char fmt[MAXFMTLEN];
271
272			if (isaspect(s))
273			sip->inpaspect *= aspectval(s);
274			else if (isncomp(s))
275			sip->nc = ncompval(s);
276			else if (iswlsplit(s))
277			wlsplitval(sip->wlpart, s);
278			else if (formatval(fmt, s) && strcmp(fmt, SPECFMT))
279			return(-1);
280			return(0);
281			}
282
283			DATARRAY *
284			getspec( /* load hyperspectral image as data */
285			char *sname
286			)
287			{
288			SPECINFO si;
289			char *pfname;
290			FILE *fp;
291			int sl, ns;
292			int y, i;
293			DATARRAY *pp;
294			/* look for array in list */
295			for (pp = dtab[hash(sname)]; pp != NULL; pp = pp->next)
296			if (!strcmp(sname, pp->name))
297			return(pp); /* found! */
298
299			if ((pfname = getpath(sname, getrlibpath(), R_OK)) == NULL) {
300			sprintf(errmsg, "cannot find hyperspectral image \"%s\"", sname);
301			error(SYSTEM, errmsg);
302			}
303			if ((fp = fopen(pfname, "rb")) == NULL) {
304			sprintf(errmsg, "cannot open hyperspectral image \"%s\"", pfname);
305			error(SYSTEM, errmsg);
306			}
307			si.wlpart[3] = 0;
308			si.nc = 0;
309			si.inpaspect = 1.0;
310			if (getheader(fp, specheadline, &si) < 0 \|\|
311			(si.nc <= 3) \| (si.nc > MAXCSAMP) \| (si.wlpart[3] < 1) \|\|
312			!fscnresolu(&sl, &ns, fp))
313			goto readerr;
314
315			if ((pp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
316			goto memerr;
317
318			pp->name = savestr(sname);
319			pp->type = SPECTY;
320			pp->nd = 3;
321			pp->dim[0].ne = ns;
322			pp->dim[1].ne = sl;
323			pp->dim[0].org =
324			pp->dim[1].org = 0.0;
325			if (sl <= ns*si.inpaspect) {
326			pp->dim[0].siz = si.inpaspect * (double)ns/sl;
327			pp->dim[1].siz = 1.0;
328			} else {
329			pp->dim[0].siz = 1.0;
330			pp->dim[1].siz = (double)sl/ns / si.inpaspect;
331			}
332			pp->dim[2].ne = si.nc;
333			pp->dim[2].siz = si.wlpart[3] - si.wlpart[0];
334			pp->dim[2].org = si.wlpart[0] + 0.5*pp->dim[2].siz/si.nc;
335			pp->dim[2].siz *= (si.nc - 1.0)/si.nc;
336			pp->dim[0].p = pp->dim[1].p = pp->dim[2].p = NULL;
337			i = nssl(si.nc+1);
338			#if PSIZWARN
339			if (i > PSIZWARN) { /* memory warning */
340			sprintf(errmsg, "hyperspectral image \"%s\" using %.1f MB of memory",
341			sname, i(1.0/(10241024)));
342			error(WARNING, errmsg);
343			}
344			#endif
345			if ((pp->arr.s = (uby8 *)malloc(i)) == NULL)
346			goto memerr;
347			for (y = 0; y < ns; y++) /* read each scanline */
348			if (freadscolrs(pp->arr.s + ysl(si.nc+1), si.nc, sl, fp) < 0)
349			goto readerr;
350			fclose(fp);
351			i = hash(sname); /* insert in hash table */
352			pp->next = dtab[i];
353			return(dtab[i] = pp);
354			memerr:
355			error(SYSTEM, "out of memory in getspec");
356			readerr:
357			sprintf(errmsg, "bad hyperspectral image \"%s\"", pfname);
358			error(USER, errmsg);
359			return NULL; /* pro forma return */
360	greg	1.1	}
361
362
363	greg	2.31	void
364	schorsch	2.27	freedata( /* release data array reference */
365			DATARRAY *dta
366			)
367	greg	1.1	{
368	greg	2.5	DATARRAY head;
369	greg	2.6	int hval, nents;
370	greg	2.31	DATARRAY dpl, dp;
371			int i;
372	greg	1.1
373	greg	2.15	if (dta == NULL) { /* free all if NULL */
374	greg	2.6	hval = 0; nents = TABSIZ;
375			} else {
376	greg	2.15	hval = hash(dta->name); nents = 1;
377	greg	2.37	if (!dta->name) { / not a data file? */
378			dta->next = dtab[hval];
379			dtab[hval] = dta; /* ...fake position */
380			}
381	greg	2.6	}
382			while (nents--) {
383			head.next = dtab[hval];
384			dpl = &head;
385			while ((dp = dpl->next) != NULL)
386	schorsch	2.24	if ((dta == NULL) \| (dta == dp)) {
387	greg	2.6	dpl->next = dp->next;
388	greg	2.36	free(dp->arr.p);
389	greg	2.6	for (i = 0; i < dp->nd; i++)
390			if (dp->dim[i].p != NULL)
391	greg	2.36	free(dp->dim[i].p);
392	greg	2.6	freestr(dp->name);
393	greg	2.36	free(dp);
394	greg	2.6	} else
395			dpl = dp;
396			dtab[hval++] = head.next;
397			}
398	greg	1.1	}
399
400
401	greg	2.37	/* internal call to interpolate data value or vector */
402			static double
403			data_interp(DATARRAY dp, double pt, double coef, DATATYPE *rvec)
404			{
405			DATARRAY sd;
406			int stride, i;
407			double x, c0, c1, y0, y1;
408			/* set up dimensions for recursion */
409			if (dp->nd > 1) {
410			sd.name = dp->name;
411			sd.type = dp->type;
412			sd.nd = dp->nd - 1;
413			memcpy(sd.dim, dp->dim+1, sd.nd*sizeof(struct dadim));
414			stride = sd.dim[i = sd.nd-1].ne + (sd.type==SPECTY);
415			while (i-- > 0)
416			stride *= sd.dim[i].ne;
417			}
418			/* get independent variable */
419			if (dp->dim[0].p == NULL) { /* evenly spaced points */
420			x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
421			x *= (double)(dp->dim[0].ne - 1);
422			i = x;
423			if (i < 0)
424			i = 0;
425			else if (i > dp->dim[0].ne - 2)
426			i = dp->dim[0].ne - 2;
427			} else { /* unevenly spaced points */
428			int lower, upper;
429			if (dp->dim[0].siz > 0.0) {
430			lower = 0;
431			upper = dp->dim[0].ne;
432			} else {
433			lower = dp->dim[0].ne;
434			upper = 0;
435			}
436			do {
437			i = (lower + upper) >> 1;
438			if (pt[0] >= dp->dim[0].p[i])
439			lower = i;
440			else
441			upper = i;
442			} while (i != (lower + upper) >> 1);
443
444			if (i > dp->dim[0].ne - 2)
445			i = dp->dim[0].ne - 2;
446
447			x = i + (pt[0] - dp->dim[0].p[i]) /
448			(dp->dim[0].p[i+1] - dp->dim[0].p[i]);
449			}
450			/*
451			* Compute interpolation coefficients:
452			* extrapolate as far as one division, then
453			* taper off harmonically to zero.
454			*/
455			if (x > i+2) {
456			c0 = 1./(i-1 - x);
457			c1 = -2.*c0;
458			} else if (x < i-1) {
459			c1 = 1./(i - x);
460			c0 = -2.*c1;
461			} else {
462			c0 = i+1 - x;
463			c1 = x - i;
464			}
465			c0 *= coef;
466			c1 *= coef;
467			/* check if vector interp */
468			if ((dp->nd == 2) & (rvec != NULL)) {
469			if (dp->type == DATATY) {
470			sd.arr.d = dp->arr.d + i*stride;
471			for (i = sd.dim[0].ne; i--; )
472			rvec[i] += c0*sd.arr.d[i]
473			+ c1*sd.arr.d[i+stride];
474			return(0.);
475			}
476			if (dp->type == SPECTY) {
477			double f;
478			sd.arr.s = dp->arr.s + i*stride;
479			f = ldexp(1.0, (int)sd.arr.s[sd.dim[0].ne]
480			- (COLXS+8));
481			for (i = sd.dim[0].ne; i--; )
482			rvec[i] += c0f(sd.arr.s[i] + 0.5);
483			sd.arr.s += stride;
484			f = ldexp(1.0, (int)sd.arr.s[sd.dim[0].ne]
485			- (COLXS+8));
486			for (i = sd.dim[0].ne; i--; )
487			rvec[i] += c1f(sd.arr.s[i] + 0.5);
488			return(0.);
489			}
490			sd.arr.c = dp->arr.c + i*stride;
491			for (i = sd.dim[0].ne; i--; )
492			rvec[i] += c0*colrval(sd.arr.c[i],sd.type)
493			+ c1*colrval(sd.arr.c[i+stride],sd.type);
494			return(0.);
495			}
496			/* get dependent variable */
497			if (dp->nd > 1) {
498			if (dp->type == DATATY) {
499			sd.arr.d = dp->arr.d + i*stride;
500			y0 = data_interp(&sd, pt+1, c0, rvec);
501			sd.arr.d += stride;
502			} else if (dp->type == SPECTY) {
503			sd.arr.s = dp->arr.s + i*stride;
504			y0 = data_interp(&sd, pt+1, c0, rvec);
505			sd.arr.s += stride;
506			} else {
507			sd.arr.c = dp->arr.c + i*stride;
508			y0 = data_interp(&sd, pt+1, c0, rvec);
509			sd.arr.c += stride;
510			}
511			y1 = data_interp(&sd, pt+1, c1, rvec);
512			} else { /* end of recursion */
513			if (dp->type == DATATY) {
514			y0 = dp->arr.d[i];
515			y1 = dp->arr.d[i+1];
516			} else if (dp->type == SPECTY) {
517			if (dp->arr.s[dp->dim[0].ne]) {
518			double f = ldexp(1.0, -(COLXS+8) +
519			(int)dp->arr.s[dp->dim[0].ne]);
520			y0 = (dp->arr.s[i] + 0.5)*f;
521			y1 = (dp->arr.s[i+1] + 0.5)*f;
522			} else
523			y0 = y1 = 0.0;
524			} else {
525			y0 = colrval(dp->arr.c[i],dp->type);
526			y1 = colrval(dp->arr.c[i+1],dp->type);
527			}
528			y0 *= c0;
529			y1 *= c1;
530			}
531			return(y0 + y1); /* coefficients already applied */
532			}
533
534
535			double
536			datavalue( /* interpolate data value at a point */
537			DATARRAY *dp,
538			double *pt
539			)
540			{
541			return(data_interp(dp, pt, 1., NULL));
542			}
543
544
545			/* Interpolate final vector corresponding to last dimension in data array */
546			DATARRAY *
547			datavector(DATARRAY dp, double pt)
548			{
549			DATARRAY *newdp;
550
551			if (dp->nd < 2)
552			error(INTERNAL, "datavector() called with 1-D array");
553			/* create vector array */
554			newdp = (DATARRAY *)malloc(sizeof(DATARRAY) -
555			(MAXDDIM-1)*sizeof(struct dadim) +
556			sizeof(DATATYPE)*dp->dim[dp->nd-1].ne);
557			if (newdp == NULL)
558			error(SYSTEM, "out of memory in datavector");
559			newdp->next = NULL;
560			newdp->name = dp->name;
561			newdp->type = DATATY;
562			newdp->nd = 1; /* vector data goes here */
563			newdp->dim[0] = dp->dim[dp->nd-1];
564			newdp->arr.d = (DATATYPE *)(newdp->dim + 1);
565			memset(newdp->arr.d, 0, sizeof(DATATYPE)*newdp->dim[0].ne);
566
567			(void)data_interp(dp, pt, 1., newdp->arr.d);
568
569			return(newdp); /* will be free'd using free() */
570			}
571
572
573			#if 0
574	greg	2.31	double
575	schorsch	2.27	datavalue( /* interpolate data value at a point */
576	greg	2.31	DATARRAY *dp,
577	schorsch	2.27	double *pt
578			)
579	greg	1.1	{
580			DATARRAY sd;
581			int asize;
582	greg	1.6	int lower, upper;
583	greg	2.31	int i;
584	greg	2.9	double x, y0, y1;
585	greg	1.4	/* set up dimensions for recursion */
586	greg	2.12	if (dp->nd > 1) {
587			sd.name = dp->name;
588			sd.type = dp->type;
589			sd.nd = dp->nd - 1;
590			asize = 1;
591			for (i = 0; i < sd.nd; i++) {
592			sd.dim[i].org = dp->dim[i+1].org;
593			sd.dim[i].siz = dp->dim[i+1].siz;
594			sd.dim[i].p = dp->dim[i+1].p;
595	greg	2.36	asize *= (sd.dim[i].ne = dp->dim[i+1].ne) +
596			((sd.type==SPECTY) & (i==sd.nd-1));
597	greg	2.12	}
598	greg	1.1	}
599	greg	1.4	/* get independent variable */
600			if (dp->dim[0].p == NULL) { /* evenly spaced points */
601			x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
602	greg	2.9	x *= (double)(dp->dim[0].ne - 1);
603	greg	1.4	i = x;
604			if (i < 0)
605			i = 0;
606			else if (i > dp->dim[0].ne - 2)
607			i = dp->dim[0].ne - 2;
608			} else { /* unevenly spaced points */
609	greg	1.6	if (dp->dim[0].siz > 0.0) {
610			lower = 0;
611			upper = dp->dim[0].ne;
612			} else {
613			lower = dp->dim[0].ne;
614			upper = 0;
615			}
616			do {
617			i = (lower + upper) >> 1;
618			if (pt[0] >= dp->dim[0].p[i])
619			lower = i;
620	greg	1.8	else
621	greg	1.6	upper = i;
622			} while (i != (lower + upper) >> 1);
623	greg	2.36
624	greg	1.8	if (i > dp->dim[0].ne - 2)
625	greg	1.4	i = dp->dim[0].ne - 2;
626	greg	2.36
627	greg	1.4	x = i + (pt[0] - dp->dim[0].p[i]) /
628			(dp->dim[0].p[i+1] - dp->dim[0].p[i]);
629			}
630			/* get dependent variable */
631	greg	2.12	if (dp->nd > 1) {
632			if (dp->type == DATATY) {
633			sd.arr.d = dp->arr.d + i*asize;
634			y0 = datavalue(&sd, pt+1);
635	greg	2.36	sd.arr.d += asize;
636			y1 = datavalue(&sd, pt+1);
637			} else if (dp->type == SPECTY) {
638			sd.arr.s = dp->arr.s + i*asize;
639			y0 = datavalue(&sd, pt+1);
640			sd.arr.s += asize;
641	greg	2.12	y1 = datavalue(&sd, pt+1);
642			} else {
643			sd.arr.c = dp->arr.c + i*asize;
644			y0 = datavalue(&sd, pt+1);
645	greg	2.36	sd.arr.c += asize;
646	greg	2.12	y1 = datavalue(&sd, pt+1);
647			}
648	greg	1.1	} else {
649	greg	2.12	if (dp->type == DATATY) {
650			y0 = dp->arr.d[i];
651			y1 = dp->arr.d[i+1];
652	greg	2.36	} else if (dp->type == SPECTY) {
653			if (dp->arr.s[dp->dim[0].ne]) {
654			double f = ldexp(1.0, -(COLXS+8) +
655			(int)dp->arr.s[dp->dim[0].ne]);
656			y0 = (dp->arr.s[i] + 0.5)*f;
657			y1 = (dp->arr.s[i+1] + 0.5)*f;
658			} else
659			y0 = y1 = 0.0;
660	greg	2.12	} else {
661			y0 = colrval(dp->arr.c[i],dp->type);
662			y1 = colrval(dp->arr.c[i+1],dp->type);
663			}
664	greg	1.1	}
665			/*
666			* Extrapolate as far as one division, then
667			* taper off harmonically to zero.
668			*/
669			if (x > i+2)
670	greg	2.9	return( (2*y1-y0)/(x-(i-1)) );
671	greg	1.1
672	greg	2.9	if (x < i-1)
673			return( (2*y0-y1)/(i-x) );
674
675			return( y0((i+1)-x) + y1(x-i) );
676	greg	1.1	}
677	greg	2.37	#endif