source: branches/1.3/gui/src/RpDxToVtk.c @ 3844

Last change on this file since 3844 was 3844, checked in by ldelgass, 11 years ago

Sync with trunk. Branch now differs only from trunk by r3722 (branch is version
1.3, trunk is version 1.4)
File size: 14.6 KB
Line 
1
2/*
3 * ----------------------------------------------------------------------
4 *  RpDxToVtk -
5 *
6 * ======================================================================
7 *  AUTHOR:  Michael McLennan, Purdue University
8 *  Copyright (c) 2004-2012  HUBzero Foundation, LLC
9 *
10 *  See the file "license.terms" for information on usage and
11 *  redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES.
12 * ======================================================================
13 */
14#include <string.h>
15#include <stdlib.h>
16#include <ctype.h>
17#include <math.h>
18#include <limits.h>
19#include <float.h>
20#include "tcl.h"
21
22static INLINE char *
23SkipSpaces(char *string)
24{
25    while (isspace(*string)) {
26        string++;
27    }
28    return string;
29}
30
31static INLINE char *
32GetLine(char **stringPtr, const char *endPtr)
33{
34    char *line, *p;
35
36    line = SkipSpaces(*stringPtr);
37    for (p = line; p < endPtr; p++) {
38        if (*p == '\n') {
39            *stringPtr = p + 1;
40            return line;
41        }
42    }
43    *stringPtr = p;
44    return line;
45}
46
47static int
48GetUniformFieldValues(Tcl_Interp *interp, int nPoints, int *counts, char **stringPtr,
49                      const char *endPtr, Tcl_Obj *objPtr)
50{
51    int i;
52    const char *p;
53    char mesg[2000];
54    double *array, scale, vmin, vmax;
55    int iX, iY, iZ;
56
57    p = *stringPtr;
58    array = malloc(sizeof(double) * nPoints);
59    if (array == NULL) {
60        return TCL_ERROR;
61    }
62    vmin = DBL_MAX, vmax = -DBL_MAX;
63    iX = iY = iZ = 0;
64    for (i = 0; i < nPoints; i++) {
65        double value;
66        char *nextPtr;
67        int loc;
68
69        if (p >= endPtr) {
70            Tcl_AppendResult(interp, "unexpected EOF in reading points",
71                             (char *)NULL);
72            return TCL_ERROR;
73        }
74        value = strtod(p, &nextPtr);
75        if (nextPtr == p) {
76            Tcl_AppendResult(interp, "bad value found in reading points",
77                             (char *)NULL);
78            return TCL_ERROR;
79        }
80        p = nextPtr;
81        loc = iZ*counts[0]*counts[1] + iY*counts[0] + iX;
82        if (++iZ >= counts[2]) {
83            iZ = 0;
84            if (++iY >= counts[1]) {
85                iY = 0;
86                ++iX;
87            }
88        }
89        array[loc] = value;
90        if (value < vmin) {
91            vmin = value;
92        }
93        if (value > vmax) {
94            vmax = value;
95        }
96    }
97    scale = 1.0 / (vmax - vmin);
98    for (i = 0; i < nPoints; i++) {
99#ifdef notdef
100        sprintf(mesg, "%g\n", (array[i] - vmin) * scale);
101#endif
102        sprintf(mesg, "%g\n", array[i]);
103        Tcl_AppendToObj(objPtr, mesg, -1);
104    }
105    free(array);
106    *stringPtr = (char *)p;
107    return TCL_OK;
108}
109
110static int
111GetCloudFieldValues(Tcl_Interp *interp, int nXYPoints, int nZPoints, char **stringPtr,
112                    const char *endPtr, Tcl_Obj *objPtr)
113{
114    int i;
115    const char *p;
116    char mesg[2000];
117    double *array, scale, vmin, vmax;
118    int iXY, iZ;
119    int nPoints;
120
121    nPoints = nXYPoints * nZPoints;
122
123    p = *stringPtr;
124    array = malloc(sizeof(double) * nPoints);
125    if (array == NULL) {
126        return TCL_ERROR;
127    }
128    vmin = DBL_MAX, vmax = -DBL_MAX;
129    iXY = iZ = 0;
130    for (i = 0; i < nPoints; i++) {
131        double value;
132        char *nextPtr;
133        int loc;
134
135        if (p >= endPtr) {
136            Tcl_AppendResult(interp, "unexpected EOF in reading points",
137                             (char *)NULL);
138            return TCL_ERROR;
139        }
140        value = strtod(p, &nextPtr);
141        if (nextPtr == p) {
142            Tcl_AppendResult(interp, "bad value found in reading points",
143                             (char *)NULL);
144            return TCL_ERROR;
145        }
146        p = nextPtr;
147        loc = nXYPoints * iZ + iXY;
148        if (++iZ >= nZPoints) {
149            iZ = 0;
150            ++iXY;
151        }
152        array[loc] = value;
153        if (value < vmin) {
154            vmin = value;
155        }
156        if (value > vmax) {
157            vmax = value;
158        }
159    }
160    scale = 1.0 / (vmax - vmin);
161    for (i = 0; i < nPoints; i++) {
162#ifdef notdef
163        sprintf(mesg, "%g\n", (array[i] - vmin) * scale);
164#endif
165        sprintf(mesg, "%g\n", array[i]);
166        Tcl_AppendToObj(objPtr, mesg, -1);
167    }
168    free(array);
169    *stringPtr = (char *)p;
170    return TCL_OK;
171}
172
173static int
174GetPoints(Tcl_Interp *interp, double *array, int nXYPoints,
175          char **stringPtr, const char *endPtr)
176{
177    int i;
178    const char *p;
179
180    p = *stringPtr;
181    if (array == NULL) {
182        return TCL_ERROR;
183    }
184    for (i = 0; i < nXYPoints; i++) {
185        double x, y, z;
186        char *nextPtr;
187
188        if (p >= endPtr) {
189            Tcl_AppendResult(interp, "unexpected EOF in reading points",
190                             (char *)NULL);
191            return TCL_ERROR;
192        }
193        x = strtod(p, &nextPtr);
194        if (nextPtr == p) {
195            Tcl_AppendResult(interp, "bad value found in reading points",
196                             (char *)NULL);
197            return TCL_ERROR;
198        }
199        p = nextPtr;
200        y = strtod(p, &nextPtr);
201        if (nextPtr == p) {
202            Tcl_AppendResult(interp, "bad value found in reading points",
203                             (char *)NULL);
204            return TCL_ERROR;
205        }
206        p = nextPtr;
207        z = strtod(p, &nextPtr);
208        if (nextPtr == p) {
209            Tcl_AppendResult(interp, "bad value found in reading points",
210                             (char *)NULL);
211            return TCL_ERROR;
212        }
213        p = nextPtr;
214
215        array[i*2  ] = x;
216        array[i*2+1] = y;
217    }
218
219    *stringPtr = (char *)p;
220    return TCL_OK;
221}
222
223/*
224 *  DxToVtk string
225 *
226 * In DX format:
227 *  rank 0 means scalars,
228 *  rank 1 means vectors,
229 *  rank 2 means matrices,
230 *  rank 3 means tensors
231 *
232 *  For rank 1, shape is a single number equal to the number of dimensions.
233 *  e.g. rank 1 shape 3 means a 3-component vector field
234 *
235 */
236
237static int
238DxToVtkCmd(ClientData clientData, Tcl_Interp *interp, int objc,
239           Tcl_Obj *const *objv)
240{
241    double *points;
242    Tcl_Obj *objPtr, *pointsObjPtr, *fieldObjPtr;
243    char *p, *pend;
244    char *string;
245    char mesg[2000];
246    double dx, dy, dz;
247    double origin[3];
248    int count[3];
249    int length, nAxes, nPoints, nXYPoints;
250    char *name;
251    int isUniform;
252    int i, iz;
253
254    name = "myScalars";
255    nAxes = nPoints = nXYPoints = 0;
256    dx = dy = dz = 0.0; /* Suppress compiler warning. */
257    origin[0] = origin[1] = origin[2] = 0.0; /* May not have an origin line. */
258    count[0] = count[1] = count[2] = 0; /* Suppress compiler warning. */
259    isUniform = 1; /* Default to expecting uniform grid */
260   
261    if (objc != 2) {
262        Tcl_AppendResult(interp, "wrong # arguments: should be \"",
263                         Tcl_GetString(objv[0]), " string\"", (char *)NULL);
264        return TCL_ERROR;
265    }
266    string = Tcl_GetStringFromObj(objv[1], &length);
267    if (strncmp("<ODX>", string, 5) == 0) {
268        string += 5;
269        length -= 5;
270    }
271    pointsObjPtr = Tcl_NewStringObj("", -1);
272    fieldObjPtr = Tcl_NewStringObj("", -1);
273    for (p = string, pend = p + length; p < pend; /*empty*/) {
274        char *line;
275        double ddx, ddy, ddz;
276
277        line = GetLine(&p, pend);
278        if (line >= pend) {
279            break;                        /* EOF */
280        }
281        if ((line[0] == '#') || (line[0] == '\n')) {
282            continue;                        /* Skip blank or comment lines. */
283        }
284        if (sscanf(line, "object %*d class gridpositions counts %d %d %d",
285                   count, count + 1, count + 2) == 3) {
286            isUniform = 1;
287            if ((count[0] < 0) || (count[1] < 0) || (count[2] < 0)) {
288                sprintf(mesg, "invalid grid size: x=%d, y=%d, z=%d",
289                        count[0], count[1], count[2]);
290                Tcl_AppendResult(interp, mesg, (char *)NULL);
291                return TCL_ERROR;
292            }
293#ifdef notdef
294            fprintf(stderr, "found gridpositions counts %d %d %d\n",
295                    count[0], count[1], count[2]);
296#endif
297        } else if (sscanf(line, "origin %lg %lg %lg", origin, origin + 1,
298                origin + 2) == 3) {
299            /* Found origin. */
300#ifdef notdef
301            fprintf(stderr, "found origin %g %g %g\n",
302                    origin[0], origin[1], origin[2]);
303#endif
304        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
305            if (nAxes == 3) {
306                Tcl_AppendResult(interp, "too many delta statements",
307                        (char *)NULL);
308                return TCL_ERROR;
309            }
310            if (ddx != 0.0) {
311                if (ddy != 0.0 || ddz != 0.0) {
312                    Tcl_AppendResult(interp, "invalid delta statement",
313                                     (char *)NULL);
314                    return TCL_ERROR;
315                }
316                dx = ddx;
317            } else if (ddy != 0.0) {
318                if (ddx != 0.0 || ddz != 0.0) {
319                    Tcl_AppendResult(interp, "invalid delta statement",
320                                     (char *)NULL);
321                    return TCL_ERROR;
322                }
323                dy = ddy;
324            } else if (ddz != 0.0) {
325                if (ddx != 0.0 || ddy != 0.0) {
326                    Tcl_AppendResult(interp, "invalid delta statement",
327                                     (char *)NULL);
328                    return TCL_ERROR;
329                }
330                dz = ddz;
331            }
332            nAxes++;
333#ifdef notdef
334            fprintf(stderr, "found delta %g %g %g\n", ddx, ddy, ddx);
335#endif
336        } else if (sscanf(line, "object %*d class regulararray count %d",
337                          &count[2]) == 1) {
338            // Z grid
339        } else if (sscanf(line, "object %*d class array type %*s rank 1 shape 3"
340                          " items %d data follows", &nXYPoints) == 1) {
341            // This is a 2D point cloud in xy with a uniform zgrid
342            isUniform = 0;
343#ifdef notdef
344            fprintf(stderr, "found class array type shape 3 nPoints=%d\n",
345                    nPoints);
346#endif
347            if (nXYPoints < 0) {
348                sprintf(mesg, "bad # points %d", nXYPoints);
349                Tcl_AppendResult(interp, mesg, (char *)NULL);
350                return TCL_ERROR;
351            }
352            points = malloc(sizeof(double) * nXYPoints * 2);
353            if (GetPoints(interp, points, nXYPoints, &p, pend) != TCL_OK) {
354                return TCL_ERROR;
355            }
356        } else if (sscanf(line, "object %*d class array type %*s rank 0"
357                          " %*s %d data follows", &nPoints) == 1) {
358#ifdef notdef
359            fprintf(stderr, "found class array type rank 0 nPoints=%d\n",
360                nPoints);
361#endif
362            if (isUniform && nPoints != count[0]*count[1]*count[2]) {
363                sprintf(mesg, "inconsistent data: expected %d points"
364                        " but found %d points", count[0]*count[1]*count[2],
365                        nPoints);
366                Tcl_AppendResult(interp, mesg, (char *)NULL);
367                return TCL_ERROR;
368            } else if (!isUniform && nPoints != nXYPoints * count[2]) {
369                sprintf(mesg, "inconsistent data: expected %d points"
370                        " but found %d points", nXYPoints * count[2],
371                        nPoints);
372                Tcl_AppendResult(interp, mesg, (char *)NULL);
373                return TCL_ERROR;
374            }
375            if (isUniform) {
376                if (GetUniformFieldValues(interp, nPoints, count, &p, pend, fieldObjPtr)
377                    != TCL_OK) {
378                    return TCL_ERROR;
379                }
380            } else {
381                if (GetCloudFieldValues(interp, nXYPoints, count[2], &p, pend, fieldObjPtr)
382                    != TCL_OK) {
383                    return TCL_ERROR;
384                }
385            }
386#ifdef notdef
387        } else {
388            fprintf(stderr, "unknown line (%.80s)\n", line);
389#endif
390        }
391    }
392
393    if (isUniform) {
394        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
395        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
396        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
397        Tcl_AppendToObj(objPtr, "DATASET STRUCTURED_POINTS\n", -1);
398        sprintf(mesg, "DIMENSIONS %d %d %d\n", count[0], count[1], count[2]);
399        Tcl_AppendToObj(objPtr, mesg, -1);
400        sprintf(mesg, "ORIGIN %g %g %g\n", origin[0], origin[1], origin[2]);
401        Tcl_AppendToObj(objPtr, mesg, -1);
402        sprintf(mesg, "SPACING %g %g %g\n", dx, dy, dz);
403        Tcl_AppendToObj(objPtr, mesg, -1);
404        sprintf(mesg, "POINT_DATA %d\n", nPoints);
405        Tcl_AppendToObj(objPtr, mesg, -1);
406        sprintf(mesg, "SCALARS %s double 1\n", name);
407        Tcl_AppendToObj(objPtr, mesg, -1);
408        sprintf(mesg, "LOOKUP_TABLE default\n");
409        Tcl_AppendToObj(objPtr, mesg, -1);
410        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
411    } else {
412        /* Fill points.  Have to wait to do this since origin, delta can come after
413         * the point list in the file.
414         */
415        for (iz = 0; iz < count[2]; iz++) {
416            for (i = 0; i < nXYPoints; i++) {
417                sprintf(mesg, "%g %g %g\n", points[i*2], points[i*2+1], origin[2] + dz * iz);
418                Tcl_AppendToObj(pointsObjPtr, mesg, -1);
419            }
420        }
421        free(points);
422
423        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
424        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
425        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
426        Tcl_AppendToObj(objPtr, "DATASET UNSTRUCTURED_GRID\n", -1);
427        sprintf(mesg, "POINTS %d double\n", nPoints);
428        Tcl_AppendToObj(objPtr, mesg, -1);
429        Tcl_AppendObjToObj(objPtr, pointsObjPtr);
430        sprintf(mesg, "POINT_DATA %d\n", nPoints);
431        Tcl_AppendToObj(objPtr, mesg, -1);
432        sprintf(mesg, "SCALARS %s double 1\n", name);
433        Tcl_AppendToObj(objPtr, mesg, -1);
434        sprintf(mesg, "LOOKUP_TABLE default\n");
435        Tcl_AppendToObj(objPtr, mesg, -1);
436        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
437    }
438
439    Tcl_DecrRefCount(pointsObjPtr);
440    Tcl_DecrRefCount(fieldObjPtr);
441    Tcl_SetObjResult(interp, objPtr);
442    return TCL_OK;
443}
444
445/*
446 * ------------------------------------------------------------------------
447 *  RpDxToVtk_Init --
448 *
449 *  Invoked when the Rappture GUI library is being initialized
450 *  to install the "DxToVtk" command into the interpreter.
451 *
452 *  Returns TCL_OK if successful, or TCL_ERROR (along with an error
453 *  message in the interp) if anything goes wrong.
454 * ------------------------------------------------------------------------
455 */
456int
457RpDxToVtk_Init(Tcl_Interp *interp)
458{
459    /* install the widget command */
460    Tcl_CreateObjCommand(interp, "Rappture::DxToVtk", DxToVtkCmd,
461        NULL, NULL);
462    return TCL_OK;
463}
Note: See TracBrowser for help on using the repository browser.