source: trunk/gui/src/RpDxToVtk.c @ 3756

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

Fix initial ranges for float->double
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.