home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / SRC / RT / DATA.C < prev    next >
C/C++ Source or Header  |  1993-10-07  |  9KB  |  363 lines
  1. /* Copyright (c) 1991 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)data.c 2.2 12/19/91 LBL";
  5. #endif
  6.  
  7. /*
  8.  *  data.c - routines dealing with interpolated data.
  9.  *
  10.  *     6/4/86
  11.  */
  12.  
  13. #include  "standard.h"
  14.  
  15. #include  "color.h"
  16.  
  17. #include  "resolu.h"
  18.  
  19. #include  "data.h"
  20.  
  21.  
  22. extern char  *fgetword();
  23.  
  24. extern char  *libpath;                  /* library search path */
  25.  
  26. static DATARRAY  *dlist = NULL;         /* data array list */
  27.  
  28. static DATARRAY  *plist = NULL;         /* picture list */
  29.  
  30.  
  31. DATARRAY *
  32. getdata(dname)                          /* get data array dname */
  33. char  *dname;
  34. {
  35.     char  word[64];
  36.     char  *dfname;
  37.     FILE  *fp;
  38.     int  asize;
  39.     register int  i, j;
  40.     register DATARRAY  *dp;
  41.                         /* look for array in list */
  42.     for (dp = dlist; dp != NULL; dp = dp->next)
  43.         if (!strcmp(dname, dp->name))
  44.             return(dp);             /* found! */
  45.  
  46.     /*
  47.      *      If we haven't loaded the data already, we will look
  48.      *  for it in the directorys specified by the library path.
  49.      *
  50.      *      The file has the following format:
  51.      *
  52.      *              N
  53.      *              beg0    end0    n0
  54.      *              beg1    end1    n1
  55.      *              . . .
  56.      *              begN    endN    nN
  57.      *              data, later dimensions changing faster
  58.      *              . . .
  59.      *
  60.      *      For irregularly spaced points, the following can be
  61.      *  substituted for begi endi ni:
  62.      *
  63.      *              0 0 ni p0i p1i .. pni
  64.      */
  65.  
  66.     if ((dfname = getpath(dname, libpath, R_OK)) == NULL) {
  67.         sprintf(errmsg, "cannot find data file \"%s\"", dname);
  68.         error(USER, errmsg);
  69.     }
  70.     if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
  71.         goto memerr;
  72.  
  73.     dp->name = savestr(dname);
  74.  
  75.     if ((fp = fopen(dfname, "r")) == NULL) {
  76.         sprintf(errmsg, "cannot open data file \"%s\"", dfname);
  77.         error(SYSTEM, errmsg);
  78.     }
  79.                             /* get dimensions */
  80.     if (fgetword(word, sizeof(word), fp) == NULL || !isint(word))
  81.         goto scanerr;
  82.     dp->nd = atoi(word);
  83.     if (dp->nd <= 0 || dp->nd > MAXDDIM) {
  84.         sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
  85.         error(USER, errmsg);
  86.     }
  87.     asize = 1;
  88.     for (i = 0; i < dp->nd; i++) {
  89.         if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
  90.             goto scanerr;
  91.         dp->dim[i].org = atof(word);
  92.         if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
  93.             goto scanerr;
  94.         dp->dim[i].siz = atof(word);
  95.         if (fgetword(word, sizeof(word), fp) == NULL || !isint(word))
  96.             goto scanerr;
  97.         dp->dim[i].ne = atoi(word);
  98.         if (dp->dim[i].ne < 2)
  99.             goto scanerr;
  100.         asize *= dp->dim[i].ne;
  101.         if ((dp->dim[i].siz -= dp->dim[i].org) == 0) {
  102.             dp->dim[i].p = (double *)malloc(dp->dim[i].ne*sizeof(double));
  103.             if (dp->dim[i].p == NULL)
  104.                 goto memerr;
  105.             for (j = 0; j < dp->dim[i].ne; j++) {
  106.                 if (fgetword(word, sizeof(word), fp) == NULL ||
  107.                         !isflt(word))
  108.                     goto scanerr;
  109.                 dp->dim[i].p[j] = atof(word);
  110.             }
  111.             for (j = 1; j < dp->dim[i].ne-1; j++)
  112.                 if ((dp->dim[i].p[j-1] < dp->dim[i].p[j]) !=
  113.                     (dp->dim[i].p[j] < dp->dim[i].p[j+1]))
  114.                     goto scanerr;
  115.             dp->dim[i].org = dp->dim[i].p[0];
  116.             dp->dim[i].siz = dp->dim[i].p[dp->dim[i].ne-1]
  117.                         - dp->dim[i].p[0];
  118.         } else
  119.             dp->dim[i].p = NULL;
  120.     }
  121.     if ((dp->arr = (DATATYPE *)malloc(asize*sizeof(DATATYPE))) == NULL)
  122.         goto memerr;
  123.     
  124.     for (i = 0; i < asize; i++) {
  125.         if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
  126.             goto scanerr;
  127.         dp->arr[i] = atof(word);
  128.     }
  129.     fclose(fp);
  130.     dp->next = dlist;
  131.     return(dlist = dp);
  132.  
  133. memerr:
  134.     error(SYSTEM, "out of memory in getdata");
  135. scanerr:
  136.     sprintf(errmsg, "%s in data file \"%s\"",
  137.             feof(fp) ? "unexpected EOF" : "bad format", dfname);
  138.     error(USER, errmsg);
  139. }
  140.  
  141.  
  142. static double  inpaspect;               /* aspect ratio of input picture */
  143.  
  144. static
  145. headaspect(s)                           /* check string for aspect ratio */
  146. char  *s;
  147. {
  148.     if (isaspect(s))
  149.         inpaspect *= aspectval(s);
  150. }
  151.  
  152.  
  153. DATARRAY *
  154. getpict(pname)                          /* get picture pname */
  155. char  *pname;
  156. {
  157.     extern char  *libpath;
  158.     char  *pfname;
  159.     FILE  *fp;
  160.     COLOR  *scanin;
  161.     int  sl, ns;
  162.     RESOLU  inpres;
  163.     FLOAT  loc[2];
  164.     int  y;
  165.     register int  x, i;
  166.     register DATARRAY  *pp;
  167.                         /* look for array in list */
  168.     for (pp = plist; pp != NULL; pp = pp->next)
  169.         if (!strcmp(pname, pp->name))
  170.             return(pp);             /* found! */
  171.  
  172.     if ((pfname = getpath(pname, libpath, R_OK)) == NULL) {
  173.         sprintf(errmsg, "cannot find picture file \"%s\"", pname);
  174.         error(USER, errmsg);
  175.     }
  176.     if ((pp = (DATARRAY *)calloc(3, sizeof(DATARRAY))) == NULL)
  177.         goto memerr;
  178.  
  179.     pp[0].name = 
  180.     pp[1].name =
  181.     pp[2].name = savestr(pname);
  182.  
  183.     if ((fp = fopen(pfname, "r")) == NULL) {
  184.         sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
  185.         error(SYSTEM, errmsg);
  186.     }
  187. #ifdef MSDOS
  188.     setmode(fileno(fp), O_BINARY);
  189. #endif
  190.                         /* get dimensions */
  191.     inpaspect = 1.0;
  192.     getheader(fp, headaspect);
  193.     if (!fgetsresolu(&inpres, fp))
  194.         goto readerr;
  195.     for (i = 0; i < 3; i++) {
  196.         pp[i].nd = 2;
  197.         pp[i].dim[0].ne = inpres.yr;
  198.         pp[i].dim[1].ne = inpres.xr;
  199.         pp[i].dim[0].org =
  200.         pp[i].dim[1].org = 0.0;
  201.         if (inpres.xr <= inpres.yr*inpaspect) {
  202.             pp[i].dim[0].siz = inpaspect *
  203.                         (double)inpres.yr/inpres.xr;
  204.             pp[i].dim[1].siz = 1.0;
  205.         } else {
  206.             pp[i].dim[0].siz = 1.0;
  207.             pp[i].dim[1].siz = (double)inpres.xr/inpres.yr /
  208.                         inpaspect;
  209.         }
  210.         pp[i].dim[0].p = pp[i].dim[1].p = NULL;
  211.         pp[i].arr = (DATATYPE *)
  212.                 malloc(inpres.xr*inpres.yr*sizeof(DATATYPE));
  213.         if (pp[i].arr == NULL)
  214.             goto memerr;
  215.     }
  216.                             /* load picture */
  217.     sl = scanlen(&inpres);
  218.     ns = numscans(&inpres);
  219.     if ((scanin = (COLOR *)malloc(sl*sizeof(COLOR))) == NULL)
  220.         goto memerr;
  221.     for (y = 0; y < ns; y++) {
  222.         if (freadscan(scanin, sl, fp) < 0)
  223.             goto readerr;
  224.         for (x = 0; x < sl; x++) {
  225.             pix2loc(loc, &inpres, x, y);
  226.             i = (int)(loc[1]*inpres.yr)*inpres.xr +
  227.                     (int)(loc[0]*inpres.xr);
  228.             pp[0].arr[i] = colval(scanin[x],RED);
  229.             pp[1].arr[i] = colval(scanin[x],GRN);
  230.             pp[2].arr[i] = colval(scanin[x],BLU);
  231.         }
  232.     }
  233.     free((char *)scanin);
  234.     fclose(fp);
  235.     pp[0].next =
  236.     pp[1].next =
  237.     pp[2].next = plist;
  238.     return(plist = pp);
  239.  
  240. memerr:
  241.     error(SYSTEM, "out of memory in getpict");
  242. readerr:
  243.     sprintf(errmsg, "bad picture file \"%s\"", pfname);
  244.     error(USER, errmsg);
  245. }
  246.  
  247.  
  248. freedata(dname)                 /* free memory associated with dname */
  249. char  *dname;
  250. {
  251.     register DATARRAY  *dp, *dpl;
  252.     register int  i;
  253.  
  254.     for (dpl = NULL, dp = dlist; dp != NULL; dpl = dp, dp = dp->next)
  255.         if (!strcmp(dname, dp->name)) {
  256.             if (dpl == NULL)
  257.                 dlist = dp->next;
  258.             else
  259.                 dpl->next = dp->next;
  260.             free((char *)dp->arr);
  261.             for (i = 0; i < dp->nd; i++)
  262.                 if (dp->dim[i].p != NULL)
  263.                     free((char *)dp->dim[i].p);
  264.             freestr(dp->name);
  265.             free((char *)dp);
  266.             return;
  267.         }
  268. }
  269.  
  270.  
  271. freepict(pname)                 /* free memory associated with pname */
  272. char  *pname;
  273. {
  274.     register DATARRAY  *pp, *ppl;
  275.  
  276.     for (ppl = NULL, pp = plist; pp != NULL; ppl = pp, pp = pp->next)
  277.         if (!strcmp(pname, pp->name)) {
  278.             if (ppl == NULL)
  279.                 plist = pp->next;
  280.             else
  281.                 ppl->next = pp->next;
  282.             free((char *)pp[0].arr);
  283.             free((char *)pp[1].arr);
  284.             free((char *)pp[2].arr);
  285.             freestr(pp[0].name);
  286.             free((char *)pp);
  287.             return;
  288.         }
  289. }
  290.  
  291.  
  292. double
  293. datavalue(dp, pt)               /* interpolate data value at a point */
  294. register DATARRAY  *dp;
  295. double  *pt;
  296. {
  297.     DATARRAY  sd;
  298.     int  asize;
  299.     int  lower, upper;
  300.     register int  i;
  301.     double  x, y, y0, y1;
  302.                     /* set up dimensions for recursion */
  303.     sd.nd = dp->nd - 1;
  304.     asize = 1;
  305.     for (i = 0; i < sd.nd; i++) {
  306.         sd.dim[i].org = dp->dim[i+1].org;
  307.         sd.dim[i].siz = dp->dim[i+1].siz;
  308.         sd.dim[i].p = dp->dim[i+1].p;
  309.         asize *= sd.dim[i].ne = dp->dim[i+1].ne;
  310.     }
  311.                     /* get independent variable */
  312.     if (dp->dim[0].p == NULL) {             /* evenly spaced points */
  313.         x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
  314.         x = x * (dp->dim[0].ne - 1);
  315.         i = x;
  316.         if (i < 0)
  317.             i = 0;
  318.         else if (i > dp->dim[0].ne - 2)
  319.             i = dp->dim[0].ne - 2;
  320.     } else {                                /* unevenly spaced points */
  321.         if (dp->dim[0].siz > 0.0) {
  322.             lower = 0;
  323.             upper = dp->dim[0].ne;
  324.         } else {
  325.             lower = dp->dim[0].ne;
  326.             upper = 0;
  327.         }
  328.         do {
  329.             i = (lower + upper) >> 1;
  330.             if (pt[0] >= dp->dim[0].p[i])
  331.                 lower = i;
  332.             else
  333.                 upper = i;
  334.         } while (i != (lower + upper) >> 1);
  335.         if (i > dp->dim[0].ne - 2)
  336.             i = dp->dim[0].ne - 2;
  337.         x = i + (pt[0] - dp->dim[0].p[i]) /
  338.                 (dp->dim[0].p[i+1] - dp->dim[0].p[i]);
  339.     }
  340.                     /* get dependent variable */
  341.     if (dp->nd == 1) {
  342.         y0 = dp->arr[i];
  343.         y1 = dp->arr[i+1];
  344.     } else {
  345.         sd.arr = &dp->arr[i*asize];
  346.         y0 = datavalue(&sd, pt+1);
  347.         sd.arr = &dp->arr[(i+1)*asize];
  348.         y1 = datavalue(&sd, pt+1);
  349.     }
  350.     /*
  351.      * Extrapolate as far as one division, then
  352.      * taper off harmonically to zero.
  353.      */
  354.     if (x > i+2)
  355.         y = (2*y1-y0)/(x-i-1);
  356.     else if (x < i-1)
  357.         y = (2*y0-y1)/(i-x);
  358.     else
  359.         y = y0*((i+1)-x) + y1*(x-i);
  360.  
  361.     return(y);
  362. }
  363.