Back to index

radiance  4R0+20100331
data.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: data.c,v 2.29 2009/05/12 16:37:53 greg Exp $";
00003 #endif
00004 /*
00005  *  data.c - routines dealing with interpolated data.
00006  */
00007 
00008 #include "copyright.h"
00009 
00010 #include  <time.h>
00011 
00012 #include  "platform.h"
00013 #include  "paths.h"
00014 #include  "standard.h"
00015 #include  "color.h"
00016 #include  "resolu.h"
00017 #include  "view.h"
00018 #include  "data.h"
00019 
00020                             /* picture memory usage before warning */
00021 #ifndef PSIZWARN
00022 #ifdef SMLMEM
00023 #define PSIZWARN     3000000
00024 #else
00025 #define PSIZWARN     10000000
00026 #endif
00027 #endif
00028 
00029 #ifndef TABSIZ
00030 #define TABSIZ              97            /* table size (prime) */
00031 #endif
00032 
00033 #define hash(s)             (shash(s)%TABSIZ)
00034 
00035 
00036 static DATARRAY       *dtab[TABSIZ];             /* data array list */
00037 
00038 static gethfunc headaspect;
00039 
00040 
00041 extern DATARRAY *
00042 getdata(                           /* get data array dname */
00043        char  *dname
00044 )
00045 {
00046        char  *dfname;
00047        FILE  *fp;
00048        int  asize;
00049        register int  i, j;
00050        register DATARRAY  *dp;
00051                                           /* look for array in list */
00052        for (dp = dtab[hash(dname)]; dp != NULL; dp = dp->next)
00053               if (!strcmp(dname, dp->name))
00054                      return(dp);          /* found! */
00055        /*
00056         *     If we haven't loaded the data already, we will look
00057         *  for it in the directories specified by the library path.
00058         *
00059         *     The file has the following format:
00060         *
00061         *            N
00062         *            beg0   end0   n0
00063         *            beg1   end1   n1
00064         *            . . .
00065         *            begN   endN   nN
00066         *            data, later dimensions changing faster
00067         *            . . .
00068         *
00069         *     For irregularly spaced points, the following can be
00070         *  substituted for begi endi ni:
00071         *
00072         *            0 0 ni p0i p1i .. pni
00073         */
00074 
00075        if ((dfname = getpath(dname, getrlibpath(), R_OK)) == NULL) {
00076               sprintf(errmsg, "cannot find data file \"%s\"", dname);
00077               error(USER, errmsg);
00078        }
00079        if ((fp = fopen(dfname, "r")) == NULL) {
00080               sprintf(errmsg, "cannot open data file \"%s\"", dfname);
00081               error(SYSTEM, errmsg);
00082        }
00083                                                  /* get dimensions */
00084        if (fgetval(fp, 'i', (char *)&asize) <= 0)
00085               goto scanerr;
00086        if ((asize <= 0) | (asize > MAXDDIM)) {
00087               sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
00088               error(USER, errmsg);
00089        }
00090        if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
00091               goto memerr;
00092        dp->name = savestr(dname);
00093        dp->type = DATATY;
00094        dp->nd = asize;
00095        asize = 1;
00096        for (i = 0; i < dp->nd; i++) {
00097               if (fgetval(fp, DATATY, (char *)&dp->dim[i].org) <= 0)
00098                      goto scanerr;
00099               if (fgetval(fp, DATATY, (char *)&dp->dim[i].siz) <= 0)
00100                      goto scanerr;
00101               if (fgetval(fp, 'i', (char *)&dp->dim[i].ne) <= 0)
00102                      goto scanerr;
00103               if (dp->dim[i].ne < 2)
00104                      goto scanerr;
00105               asize *= dp->dim[i].ne;
00106               if ((dp->dim[i].siz -= dp->dim[i].org) == 0) {
00107                      dp->dim[i].p = (DATATYPE *)
00108                                    malloc(dp->dim[i].ne*sizeof(DATATYPE));
00109                      if (dp->dim[i].p == NULL)
00110                             goto memerr;
00111                      for (j = 0; j < dp->dim[i].ne; j++)
00112                             if (fgetval(fp, DATATY,
00113                                           (char *)&dp->dim[i].p[j]) <= 0)
00114                                    goto scanerr;
00115                      for (j = 1; j < dp->dim[i].ne-1; j++)
00116                             if ((dp->dim[i].p[j-1] < dp->dim[i].p[j]) !=
00117                                    (dp->dim[i].p[j] < dp->dim[i].p[j+1]))
00118                                    goto scanerr;
00119                      dp->dim[i].org = dp->dim[i].p[0];
00120                      dp->dim[i].siz = dp->dim[i].p[dp->dim[i].ne-1]
00121                                           - dp->dim[i].p[0];
00122               } else
00123                      dp->dim[i].p = NULL;
00124        }
00125        if ((dp->arr.d = (DATATYPE *)malloc(asize*sizeof(DATATYPE))) == NULL)
00126               goto memerr;
00127        
00128        for (i = 0; i < asize; i++)
00129               if (fgetval(fp, DATATY, (char *)&dp->arr.d[i]) <= 0)
00130                      goto scanerr;
00131        fclose(fp);
00132        i = hash(dname);
00133        dp->next = dtab[i];
00134        return(dtab[i] = dp);
00135 
00136 memerr:
00137        error(SYSTEM, "out of memory in getdata");
00138 scanerr:
00139        sprintf(errmsg, "%s in data file \"%s\"",
00140                      feof(fp) ? "unexpected EOF" : "bad format", dfname);
00141        error(USER, errmsg);
00142        return NULL; /* pro forma return */
00143 }
00144 
00145 
00146 static int
00147 headaspect(                 /* check string for aspect ratio */
00148        char  *s,
00149        void  *iap
00150 )
00151 {
00152        char   fmt[32];
00153 
00154        if (isaspect(s))
00155               *(double*)iap *= aspectval(s);
00156        else if (formatval(fmt, s) && !globmatch(PICFMT, fmt))
00157               *(double*)iap = 0.0;
00158        return(0);
00159 }
00160 
00161 
00162 extern DATARRAY *
00163 getpict(                           /* get picture pname */
00164        char  *pname
00165 )
00166 {
00167        double  inpaspect;
00168        char  *pfname;
00169        FILE  *fp;
00170        COLR  *scanin;
00171        int  sl, ns;
00172        RESOLU inpres;
00173        RREAL  loc[2];
00174        int  y;
00175        register int  x, i;
00176        register DATARRAY  *pp;
00177                                           /* look for array in list */
00178        for (pp = dtab[hash(pname)]; pp != NULL; pp = pp->next)
00179               if (!strcmp(pname, pp->name))
00180                      return(pp);          /* found! */
00181 
00182        if ((pfname = getpath(pname, getrlibpath(), R_OK)) == NULL) {
00183               sprintf(errmsg, "cannot find picture file \"%s\"", pname);
00184               error(USER, errmsg);
00185        }
00186        if ((pp = (DATARRAY *)malloc(3*sizeof(DATARRAY))) == NULL)
00187               goto memerr;
00188 
00189        pp[0].name = savestr(pname);
00190 
00191        if ((fp = fopen(pfname, "r")) == NULL) {
00192               sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
00193               error(SYSTEM, errmsg);
00194        }
00195        SET_FILE_BINARY(fp);
00196                                           /* get dimensions */
00197        inpaspect = 1.0;
00198        getheader(fp, headaspect, &inpaspect);
00199        if (inpaspect <= FTINY || !fgetsresolu(&inpres, fp))
00200               goto readerr;
00201        pp[0].nd = 2;
00202        pp[0].dim[0].ne = inpres.yr;
00203        pp[0].dim[1].ne = inpres.xr;
00204        pp[0].dim[0].org =
00205        pp[0].dim[1].org = 0.0;
00206        if (inpres.xr <= inpres.yr*inpaspect) {
00207               pp[0].dim[0].siz = inpaspect *
00208                                    (double)inpres.yr/inpres.xr;
00209               pp[0].dim[1].siz = 1.0;
00210        } else {
00211               pp[0].dim[0].siz = 1.0;
00212               pp[0].dim[1].siz = (double)inpres.xr/inpres.yr /
00213                                    inpaspect;
00214        }
00215        pp[0].dim[0].p = pp[0].dim[1].p = NULL;
00216        sl = scanlen(&inpres);                           /* allocate array */
00217        ns = numscans(&inpres);
00218        i = ns*sl*sizeof(COLR);
00219 #if PSIZWARN
00220        if (i > PSIZWARN) {                       /* memory warning */
00221               sprintf(errmsg, "picture file \"%s\" using %.1f MB of memory",
00222                             pname, i*(1.0/(1024*1024)));
00223               error(WARNING, errmsg);
00224        }
00225 #endif
00226        if ((pp[0].arr.c = (COLR *)malloc(i)) == NULL)
00227               goto memerr;
00228                                                  /* load picture */
00229        if ((scanin = (COLR *)malloc(sl*sizeof(COLR))) == NULL)
00230               goto memerr;
00231        for (y = 0; y < ns; y++) {
00232               if (freadcolrs(scanin, sl, fp) < 0)
00233                      goto readerr;
00234               for (x = 0; x < sl; x++) {
00235                      pix2loc(loc, &inpres, x, y);
00236                      i = (int)(loc[1]*inpres.yr)*inpres.xr +
00237                                    (int)(loc[0]*inpres.xr);
00238                      copycolr(pp[0].arr.c[i], scanin[x]);
00239               }
00240        }
00241        free((void *)scanin);
00242        fclose(fp);
00243        i = hash(pname);
00244        pp[0].next = dtab[i];              /* link into picture list */
00245        pp[1] = pp[0];
00246        pp[2] = pp[0];
00247        pp[0].type = RED;           /* differentiate RGB records */
00248        pp[1].type = GRN;
00249        pp[2].type = BLU;
00250        return(dtab[i] = pp);
00251 
00252 memerr:
00253        error(SYSTEM, "out of memory in getpict");
00254 readerr:
00255        sprintf(errmsg, "bad picture file \"%s\"", pfname);
00256        error(USER, errmsg);
00257        return NULL; /* pro forma return */
00258 }
00259 
00260 
00261 extern void
00262 freedata(                   /* release data array reference */
00263        DATARRAY  *dta
00264 )
00265 {
00266        DATARRAY  head;
00267        int  hval, nents;
00268        register DATARRAY  *dpl, *dp;
00269        register int  i;
00270 
00271        if (dta == NULL) {                 /* free all if NULL */
00272               hval = 0; nents = TABSIZ;
00273        } else {
00274               hval = hash(dta->name); nents = 1;
00275        }
00276        while (nents--) {
00277               head.next = dtab[hval];
00278               dpl = &head;
00279               while ((dp = dpl->next) != NULL)
00280                      if ((dta == NULL) | (dta == dp)) {
00281                             dpl->next = dp->next;
00282                             if (dp->type == DATATY)
00283                                    free((void *)dp->arr.d);
00284                             else
00285                                    free((void *)dp->arr.c);
00286                             for (i = 0; i < dp->nd; i++)
00287                                    if (dp->dim[i].p != NULL)
00288                                           free((void *)dp->dim[i].p);
00289                             freestr(dp->name);
00290                             free((void *)dp);
00291                      } else
00292                             dpl = dp;
00293               dtab[hval++] = head.next;
00294        }
00295 }
00296 
00297 
00298 extern double
00299 datavalue(           /* interpolate data value at a point */
00300        register DATARRAY  *dp,
00301        double *pt
00302 )
00303 {
00304        DATARRAY  sd;
00305        int  asize;
00306        int  lower, upper;
00307        register int  i;
00308        double x, y0, y1;
00309                                    /* set up dimensions for recursion */
00310        if (dp->nd > 1) {
00311               sd.name = dp->name;
00312               sd.type = dp->type;
00313               sd.nd = dp->nd - 1;
00314               asize = 1;
00315               for (i = 0; i < sd.nd; i++) {
00316                      sd.dim[i].org = dp->dim[i+1].org;
00317                      sd.dim[i].siz = dp->dim[i+1].siz;
00318                      sd.dim[i].p = dp->dim[i+1].p;
00319                      asize *= sd.dim[i].ne = dp->dim[i+1].ne;
00320               }
00321        }
00322                                    /* get independent variable */
00323        if (dp->dim[0].p == NULL) {        /* evenly spaced points */
00324               x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
00325               x *= (double)(dp->dim[0].ne - 1);
00326               i = x;
00327               if (i < 0)
00328                      i = 0;
00329               else if (i > dp->dim[0].ne - 2)
00330                      i = dp->dim[0].ne - 2;
00331        } else {                           /* unevenly spaced points */
00332               if (dp->dim[0].siz > 0.0) {
00333                      lower = 0;
00334                      upper = dp->dim[0].ne;
00335               } else {
00336                      lower = dp->dim[0].ne;
00337                      upper = 0;
00338               }
00339               do {
00340                      i = (lower + upper) >> 1;
00341                      if (pt[0] >= dp->dim[0].p[i])
00342                             lower = i;
00343                      else
00344                             upper = i;
00345               } while (i != (lower + upper) >> 1);
00346               if (i > dp->dim[0].ne - 2)
00347                      i = dp->dim[0].ne - 2;
00348               x = i + (pt[0] - dp->dim[0].p[i]) /
00349                             (dp->dim[0].p[i+1] - dp->dim[0].p[i]);
00350        }
00351                                    /* get dependent variable */
00352        if (dp->nd > 1) {
00353               if (dp->type == DATATY) {
00354                      sd.arr.d = dp->arr.d + i*asize;
00355                      y0 = datavalue(&sd, pt+1);
00356                      sd.arr.d = dp->arr.d + (i+1)*asize;
00357                      y1 = datavalue(&sd, pt+1);
00358               } else {
00359                      sd.arr.c = dp->arr.c + i*asize;
00360                      y0 = datavalue(&sd, pt+1);
00361                      sd.arr.c = dp->arr.c + (i+1)*asize;
00362                      y1 = datavalue(&sd, pt+1);
00363               }
00364        } else {
00365               if (dp->type == DATATY) {
00366                      y0 = dp->arr.d[i];
00367                      y1 = dp->arr.d[i+1];
00368               } else {
00369                      y0 = colrval(dp->arr.c[i],dp->type);
00370                      y1 = colrval(dp->arr.c[i+1],dp->type);
00371               }
00372        }
00373        /*
00374         * Extrapolate as far as one division, then
00375         * taper off harmonically to zero.
00376         */
00377        if (x > i+2)
00378               return( (2*y1-y0)/(x-(i-1)) );
00379 
00380        if (x < i-1)
00381               return( (2*y0-y1)/(i-x) );
00382 
00383        return( y0*((i+1)-x) + y1*(x-i) );
00384 }