src/rt/data.c

/* Copyright (c) 1991 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 *  data.c - routines dealing with interpolated data.
 *
 *     6/4/86
 */

#include  "standard.h"

#include  "color.h"

#include  "data.h"


extern char  *fgetword();
extern double  atof();

extern char  *libpath;                  /* library search path */

static DATARRAY  *dlist = NULL;         /* data array list */

static DATARRAY  *plist = NULL;         /* picture list */


DATARRAY *
getdata(dname)                          /* get data array dname */
char  *dname;
{
        char  word[64];
        char  *dfname;
        FILE  *fp;
        int  asize;
        register int  i, j;
        register DATARRAY  *dp;
                                                /* look for array in list */
        for (dp = dlist; 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 directorys 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, libpath, R_OK)) == NULL) {
                sprintf(errmsg, "cannot find data file \"%s\"", dname);
                error(USER, errmsg);
        }
        if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
                goto memerr;

        dp->name = savestr(dname);

        if ((fp = fopen(dfname, "r")) == NULL) {
                sprintf(errmsg, "cannot open data file \"%s\"", dfname);
                error(SYSTEM, errmsg);
        }
                                                        /* get dimensions */
        if (fgetword(word, sizeof(word), fp) == NULL || !isint(word))
                goto scanerr;
        dp->nd = atoi(word);
        if (dp->nd <= 0 || dp->nd > MAXDDIM) {
                sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
                error(USER, errmsg);
        }
        asize = 1;
        for (i = 0; i < dp->nd; i++) {
                if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
                        goto scanerr;
                dp->dim[i].org = atof(word);
                if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
                        goto scanerr;
                dp->dim[i].siz = atof(word);
                if (fgetword(word, sizeof(word), fp) == NULL || !isint(word))
                        goto scanerr;
                dp->dim[i].ne = atoi(word);
                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 = (double *)malloc(dp->dim[i].ne*sizeof(double));
                        if (dp->dim[i].p == NULL)
                                goto memerr;
                        for (j = 0; j < dp->dim[i].ne; j++) {
                                if (fgetword(word, sizeof(word), fp) == NULL ||
                                                !isflt(word))
                                        goto scanerr;
                                dp->dim[i].p[j] = atof(word);
                        }
                        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 = (DATATYPE *)malloc(asize*sizeof(DATATYPE))) == NULL)
                goto memerr;
        
        for (i = 0; i < asize; i++) {
                if (fgetword(word, sizeof(word), fp) == NULL || !isflt(word))
                        goto scanerr;
                dp->arr[i] = atof(word);
        }
        fclose(fp);
        dp->next = dlist;
        return(dlist = 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);
}


static double  inpaspect;               /* aspect ratio of input picture */

static
headaspect(s)                           /* check string for aspect ratio */
char  *s;
{
        if (isaspect(s))
                inpaspect *= aspectval(s);
}


DATARRAY *
getpict(pname)                          /* get picture pname */
char  *pname;
{
        extern char  *libpath;
        char  *pfname;
        FILE  *fp;
        COLOR  *scanin;
        int  width, height;
        int  x, y;
        register int  i;
        register DATARRAY  *pp;
                                                /* look for array in list */
        for (pp = plist; pp != NULL; pp = pp->next)
                if (!strcmp(pname, pp->name))
                        return(pp);             /* found! */

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

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

        if ((fp = fopen(pfname, "r")) == NULL) {
                sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
                error(SYSTEM, errmsg);
        }
                                                /* get dimensions */
        inpaspect = 1.0;
        getheader(fp, headaspect);
        if (fgetresolu(&width, &height, fp) != (YMAJOR|YDECR))
                goto readerr;
        for (i = 0; i < 3; i++) {
                pp[i].nd = 2;
                pp[i].dim[0].ne = height;
                pp[i].dim[1].ne = width;
                pp[i].dim[0].org =
                pp[i].dim[1].org = 0.0;
                if (width <= height*inpaspect) {
                        pp[i].dim[0].siz = inpaspect*(double)height/width;
                        pp[i].dim[1].siz = 1.0;
                } else {
                        pp[i].dim[0].siz = 1.0;
                        pp[i].dim[1].siz = (double)width/height/inpaspect;
                }
                pp[i].dim[0].p = pp[i].dim[1].p = NULL;
                pp[i].arr = (DATATYPE *)malloc(width*height*sizeof(DATATYPE));
                if (pp[i].arr == NULL)
                        goto memerr;
        }
                                                        /* load picture */
        if ((scanin = (COLOR *)malloc(width*sizeof(COLOR))) == NULL)
                goto memerr;
        for (y = height-1; y >= 0; y--) {
                if (freadscan(scanin, width, fp) < 0)
                        goto readerr;
                for (x = 0; x < width; x++)
                        for (i = 0; i < 3; i++)
                                pp[i].arr[y*width+x] = colval(scanin[x],i);
        }
        free((char *)scanin);
        fclose(fp);
        pp[0].next =
        pp[1].next =
        pp[2].next = plist;
        return(plist = pp);

memerr:
        error(SYSTEM, "out of memory in getpict");
readerr:
        sprintf(errmsg, "bad picture file \"%s\"", pfname);
        error(USER, errmsg);
}


freedata(dname)                 /* free memory associated with dname */
char  *dname;
{
        register DATARRAY  *dp, *dpl;
        register int  i;

        for (dpl = NULL, dp = dlist; dp != NULL; dpl = dp, dp = dp->next)
                if (!strcmp(dname, dp->name)) {
                        if (dpl == NULL)
                                dlist = dp->next;
                        else
                                dpl->next = dp->next;
                        free((char *)dp->arr);
                        for (i = 0; i < dp->nd; i++)
                                if (dp->dim[i].p != NULL)
                                        free((char *)dp->dim[i].p);
                        freestr(dp->name);
                        free((char *)dp);
                        return;
                }
}


freepict(pname)                 /* free memory associated with pname */
char  *pname;
{
        register DATARRAY  *pp, *ppl;

        for (ppl = NULL, pp = plist; pp != NULL; ppl = pp, pp = pp->next)
                if (!strcmp(pname, pp->name)) {
                        if (ppl == NULL)
                                plist = pp->next;
                        else
                                ppl->next = pp->next;
                        free((char *)pp[0].arr);
                        free((char *)pp[1].arr);
                        free((char *)pp[2].arr);
                        freestr(pp[0].name);
                        free((char *)pp);
                        return;
                }
}


double
datavalue(dp, pt)               /* interpolate data value at a point */
register DATARRAY  *dp;
double  *pt;
{
        DATARRAY  sd;
        int  asize;
        int  lower, upper;
        register int  i;
        double  x, y, y0, y1;
                                        /* set up dimensions for recursion */
        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;
        }
                                        /* get independent variable */
        if (dp->dim[0].p == NULL) {             /* evenly spaced points */
                x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
                x = x * (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) {
                y0 = dp->arr[i];
                y1 = dp->arr[i+1];
        } else {
                sd.arr = &dp->arr[i*asize];
                y0 = datavalue(&sd, pt+1);
                sd.arr = &dp->arr[(i+1)*asize];
                y1 = datavalue(&sd, pt+1);
        }
        /*
         * Extrapolate as far as one division, then
         * taper off harmonically to zero.
         */
        if (x > i+2)
                y = (2*y1-y0)/(x-i-1);
        else if (x < i-1)
                y = (2*y0-y1)/(i-x);
        else
                y = y0*((i+1)-x) + y1*(x-i);

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