source: trunk/gui/src/RpDxToVtk.c

Last change on this file was 5673, checked in by ldelgass, 9 years ago

Fix line endings, set eol-style to native on all C/C++ sources.
  • Property svn:eol-style set to native
File size: 26.9 KB
Line 
1 /* -*- mode: c; c-basic-offset: 4; indent-tabs-mode: nil -*- */
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 <stdio.h>
15#include <string.h>
16#include <stdlib.h>
17#include <sys/types.h>
18#include <unistd.h>
19#include <ctype.h>
20#include <math.h>
21#include <limits.h>
22#include <float.h>
23#include "tcl.h"
24
25#define DO_WEDGES
26//#define CHECK_WINDINGS
27
28static INLINE char *
29SkipSpaces(char *string)
30{
31    while (isspace(*string)) {
32        string++;
33    }
34    return string;
35}
36
37static INLINE char *
38GetLine(char **stringPtr, const char *endPtr)
39{
40    char *line, *p;
41
42    line = SkipSpaces(*stringPtr);
43    for (p = line; p < endPtr; p++) {
44        if (*p == '\n') {
45            *stringPtr = p + 1;
46            return line;
47        }
48    }
49    *stringPtr = p;
50    return line;
51}
52
53static int
54GetUniformFieldValues(Tcl_Interp *interp, int nPoints, int *counts, char **stringPtr,
55                      const char *endPtr, Tcl_Obj *objPtr)
56{
57    int i;
58    const char *p;
59    char mesg[2000];
60    double *array;
61    int iX, iY, iZ;
62
63    p = *stringPtr;
64    array = malloc(sizeof(double) * nPoints);
65    if (array == NULL) {
66        return TCL_ERROR;
67    }
68    iX = iY = iZ = 0;
69    for (i = 0; i < nPoints; i++) {
70        double value;
71        char *nextPtr;
72        int loc;
73
74        if (p >= endPtr) {
75            Tcl_AppendResult(interp, "unexpected EOF in reading field values",
76                             (char *)NULL);
77            free(array);
78            return TCL_ERROR;
79        }
80        value = strtod(p, &nextPtr);
81        if (nextPtr == p) {
82            Tcl_AppendResult(interp, "bad value found in reading field values",
83                             (char *)NULL);
84            free(array);
85            return TCL_ERROR;
86        }
87        p = nextPtr;
88        loc = iZ*counts[0]*counts[1] + iY*counts[0] + iX;
89        if (++iZ >= counts[2]) {
90            iZ = 0;
91            if (++iY >= counts[1]) {
92                iY = 0;
93                ++iX;
94            }
95        }
96        array[loc] = value;
97    }
98    for (i = 0; i < nPoints; i++) {
99        sprintf(mesg, "%g\n", array[i]);
100        Tcl_AppendToObj(objPtr, mesg, -1);
101    }
102    free(array);
103    *stringPtr = (char *)p;
104    return TCL_OK;
105}
106
107static int
108GetStructuredGridFieldValues(Tcl_Interp *interp, int nPoints, char **stringPtr,
109                             const char *endPtr, Tcl_Obj *objPtr)
110{
111    int i;
112    const char *p;
113    char mesg[2000];
114    double *array;
115
116    p = *stringPtr;
117    array = malloc(sizeof(double) * nPoints);
118    if (array == NULL) {
119        return TCL_ERROR;
120    }
121    for (i = 0; i < nPoints; i++) {
122        double value;
123        char *nextPtr;
124
125        if (p >= endPtr) {
126            Tcl_AppendResult(interp, "unexpected EOF in reading field values",
127                             (char *)NULL);
128            free(array);
129            return TCL_ERROR;
130        }
131        value = strtod(p, &nextPtr);
132        if (nextPtr == p) {
133            Tcl_AppendResult(interp, "bad value found in reading field values",
134                             (char *)NULL);
135            free(array);
136            return TCL_ERROR;
137        }
138        p = nextPtr;
139        array[i] = value;
140    }
141    for (i = 0; i < nPoints; i++) {
142        sprintf(mesg, "%g\n", array[i]);
143        Tcl_AppendToObj(objPtr, mesg, -1);
144    }
145    free(array);
146    *stringPtr = (char *)p;
147    return TCL_OK;
148}
149
150static int
151GetCloudFieldValues(Tcl_Interp *interp, int nXYPoints, int nZPoints, char **stringPtr,
152                    const char *endPtr, Tcl_Obj *objPtr)
153{
154    int i;
155    const char *p;
156    char mesg[2000];
157    double *array;
158    int iXY, iZ;
159    int nPoints;
160
161    nPoints = nXYPoints * nZPoints;
162
163    p = *stringPtr;
164    array = malloc(sizeof(double) * nPoints);
165    if (array == NULL) {
166        return TCL_ERROR;
167    }
168    iXY = iZ = 0;
169    for (i = 0; i < nPoints; i++) {
170        double value;
171        char *nextPtr;
172        int loc;
173
174        if (p >= endPtr) {
175            Tcl_AppendResult(interp, "unexpected EOF in reading field values",
176                             (char *)NULL);
177            free(array);
178            return TCL_ERROR;
179        }
180        value = strtod(p, &nextPtr);
181        if (nextPtr == p) {
182            Tcl_AppendResult(interp, "bad value found in reading field values",
183                             (char *)NULL);
184            free(array);
185            return TCL_ERROR;
186        }
187        p = nextPtr;
188        loc = nXYPoints * iZ + iXY;
189        if (++iZ >= nZPoints) {
190            iZ = 0;
191            ++iXY;
192        }
193        array[loc] = value;
194    }
195    for (i = 0; i < nPoints; i++) {
196        sprintf(mesg, "%g\n", array[i]);
197        Tcl_AppendToObj(objPtr, mesg, -1);
198    }
199    free(array);
200    *stringPtr = (char *)p;
201    return TCL_OK;
202}
203
204static int
205GetPoints(Tcl_Interp *interp, double *array, int nXYPoints,
206          char **stringPtr, const char *endPtr)
207{
208    int i;
209    const char *p;
210
211    p = *stringPtr;
212    if (array == NULL) {
213        return TCL_ERROR;
214    }
215    for (i = 0; i < nXYPoints; i++) {
216        double x, y;
217        char *nextPtr;
218
219        if (p >= endPtr) {
220            Tcl_AppendResult(interp, "unexpected EOF in reading points",
221                             (char *)NULL);
222            return TCL_ERROR;
223        }
224        x = strtod(p, &nextPtr);
225        if (nextPtr == p) {
226            Tcl_AppendResult(interp, "bad value found in reading points",
227                             (char *)NULL);
228            return TCL_ERROR;
229        }
230        p = nextPtr;
231        y = strtod(p, &nextPtr);
232        if (nextPtr == p) {
233            Tcl_AppendResult(interp, "bad value found in reading points",
234                             (char *)NULL);
235            return TCL_ERROR;
236        }
237        p = nextPtr;
238        /* z is unused */
239        strtod(p, &nextPtr);
240        if (nextPtr == p) {
241            Tcl_AppendResult(interp, "bad value found in reading points",
242                             (char *)NULL);
243            return TCL_ERROR;
244        }
245        p = nextPtr;
246
247        array[i*2  ] = x;
248        array[i*2+1] = y;
249    }
250
251    *stringPtr = (char *)p;
252    return TCL_OK;
253}
254
255static int
256GetUniformFieldVectors(Tcl_Interp *interp, int nPoints, int *counts,
257                       char **stringPtr, const char *endPtr, Tcl_Obj *objPtr)
258{
259    int i;
260    const char *p;
261    char mesg[2000];
262    double *array;
263    int iX, iY, iZ;
264
265    p = *stringPtr;
266    array = malloc(sizeof(double) * nPoints * 3);
267    if (array == NULL) {
268        return TCL_ERROR;
269    }
270    iX = iY = iZ = 0;
271    for (i = 0; i < nPoints; i++) {
272        double x, y, z;
273        char *nextPtr;
274        int loc;
275
276        if (p >= endPtr) {
277            Tcl_AppendResult(interp, "unexpected EOF in reading vectors",
278                             (char *)NULL);
279            free(array);
280            return TCL_ERROR;
281        }
282        x = strtod(p, &nextPtr);
283        if (nextPtr == p) {
284            Tcl_AppendResult(interp, "bad value found in reading vectors",
285                             (char *)NULL);
286            free(array);
287            return TCL_ERROR;
288        }
289        p = nextPtr;
290        y = strtod(p, &nextPtr);
291        if (nextPtr == p) {
292            Tcl_AppendResult(interp, "bad value found in reading vectors",
293                             (char *)NULL);
294            free(array);
295            return TCL_ERROR;
296        }
297        p = nextPtr;
298        z = strtod(p, &nextPtr);
299        if (nextPtr == p) {
300            Tcl_AppendResult(interp, "bad value found in reading vectors",
301                             (char *)NULL);
302            free(array);
303            return TCL_ERROR;
304        }
305        p = nextPtr;
306        loc = iZ*counts[0]*counts[1] + iY*counts[0] + iX;
307        if (++iZ >= counts[2]) {
308            iZ = 0;
309            if (++iY >= counts[1]) {
310                iY = 0;
311                ++iX;
312            }
313        }
314        array[loc*3  ] = x;
315        array[loc*3+1] = y;
316        array[loc*3+2] = z;
317    }
318    for (i = 0; i < nPoints; i++) {
319        sprintf(mesg, "%g %g %g\n", array[i*3], array[i*3+1], array[i*3+2]);
320        Tcl_AppendToObj(objPtr, mesg, -1);
321    }
322    free(array);
323    *stringPtr = (char *)p;
324    return TCL_OK;
325}
326
327#if defined(DO_WEDGES) && defined(CHECK_WINDINGS)
328static void
329Normalize(double v[3])
330{
331    double len = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
332    v[0] /= len;
333    v[1] /= len;
334    v[2] /= len;
335}
336
337static double
338Dot(double v1[3], double v2[3])
339{
340    return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
341}
342
343static void
344Cross(double v1[3], double v2[3], double vout[3])
345{
346    vout[0] = v1[1]*v2[2] - v1[2]*v2[1];
347    vout[1] = v1[2]*v2[0] - v1[0]*v2[2];
348    vout[2] = v1[0]*v2[1] - v1[1]*v2[0];
349}
350#endif
351/*
352 *  DxToVtk string
353 *
354 * In DX format:
355 *  rank 0 means scalars,
356 *  rank 1 means vectors,
357 *  rank 2 means matrices,
358 *  rank 3 means tensors
359 *
360 *  For rank 1, shape is a single number equal to the number of dimensions.
361 *  e.g. rank 1 shape 3 means a 3-component vector field
362 *
363 */
364
365static int
366DxToVtkCmd(ClientData clientData, Tcl_Interp *interp, int objc,
367           Tcl_Obj *const *objv)
368{
369    double *points;
370    Tcl_Obj *objPtr, *pointsObjPtr, *fieldObjPtr, *cellsObjPtr;
371    char *p, *pend;
372    char *string;
373    char mesg[2000];
374    double dv0[3], dv1[3], dv2[3];
375    double dx, dy, dz;
376    double origin[3];
377    int count[3];
378    int length, nAxes, nPoints, nXYPoints, nCells;
379    char *name;
380    int isUniform;
381    int isStructuredGrid;
382    int hasVectors;
383    int i, ix, iy, iz;
384
385    name = "component";
386    points = NULL;
387    nAxes = nPoints = nXYPoints = nCells = 0;
388    dx = dy = dz = 0.0;
389    origin[0] = origin[1] = origin[2] = 0.0; /* May not have an origin line. */
390    dv0[0] = dv0[1] = dv0[2] = 0.0;
391    dv1[0] = dv1[1] = dv1[2] = 0.0;
392    dv2[0] = dv2[1] = dv2[2] = 0.0;
393    count[0] = count[1] = count[2] = 0;
394    isUniform = 0;
395    isStructuredGrid = 0;
396    hasVectors = 0;
397
398    if (objc != 2) {
399        Tcl_AppendResult(interp, "wrong # arguments: should be \"",
400                         Tcl_GetString(objv[0]), " string\"", (char *)NULL);
401        return TCL_ERROR;
402    }
403    string = Tcl_GetStringFromObj(objv[1], &length);
404    if (strncmp("<ODX>", string, 5) == 0) {
405        string += 5;
406        length -= 5;
407    } else if (strncmp("<DX>", string, 4) == 0) {
408        string += 4;
409        length -= 4;
410    }
411    pointsObjPtr = Tcl_NewStringObj("", -1);
412    cellsObjPtr = Tcl_NewStringObj("", -1);
413    fieldObjPtr = Tcl_NewStringObj("", -1);
414    for (p = string, pend = p + length; p < pend; /*empty*/) {
415        char *line;
416        double ddx, ddy, ddz;
417
418        line = GetLine(&p, pend);
419        if (line >= pend) {
420            break;                        /* EOF */
421        }
422        if ((line[0] == '#') || (line[0] == '\n')) {
423            continue;                        /* Skip blank or comment lines. */
424        }
425        if (sscanf(line, "object %*d class gridpositions counts %d %d %d",
426                   count, count + 1, count + 2) == 3) {
427            isUniform = 1;
428            if ((count[0] < 0) || (count[1] < 0) || (count[2] < 0)) {
429                sprintf(mesg, "Invalid grid size: x=%d, y=%d, z=%d",
430                        count[0], count[1], count[2]);
431                Tcl_AppendResult(interp, mesg, (char *)NULL);
432                return TCL_ERROR;
433            }
434#ifdef notdef
435            fprintf(stderr, "found gridpositions counts %d %d %d\n",
436                    count[0], count[1], count[2]);
437#endif
438        } else if (sscanf(line, "origin %lg %lg %lg", origin, origin + 1,
439                origin + 2) == 3) {
440            /* Found origin. */
441#ifdef notdef
442            fprintf(stderr, "found origin %g %g %g\n",
443                    origin[0], origin[1], origin[2]);
444#endif
445        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
446            if (nAxes == 3) {
447                Tcl_AppendResult(interp, "too many delta statements",
448                                 (char *)NULL);
449                return TCL_ERROR;
450            }
451            switch (nAxes) {
452            case 0:
453                dv0[0] = ddx;
454                dv0[1] = ddy;
455                dv0[2] = ddz;
456                break;
457            case 1:
458                dv1[0] = ddx;
459                dv1[1] = ddy;
460                dv1[2] = ddz;
461                break;
462            case 2:
463                dv2[0] = ddx;
464                dv2[1] = ddy;
465                dv2[2] = ddz;
466                break;
467            default:
468                break;
469            }
470
471            if (ddx != 0.0) {
472                if (ddy != 0.0 || ddz != 0.0) {
473                    /* Not axis aligned or skewed grid */
474                    isUniform = 0;
475                    isStructuredGrid = 1;
476                }
477                dx = ddx;
478            } else if (ddy != 0.0) {
479                if (ddx != 0.0 || ddz != 0.0) {
480                    /* Not axis aligned or skewed grid */
481                    isUniform = 0;
482                    isStructuredGrid = 1;
483                }
484                dy = ddy;
485            } else if (ddz != 0.0) {
486                if (ddx != 0.0 || ddy != 0.0) {
487                    /* Not axis aligned or skewed grid */
488                    isUniform = 0;
489                    isStructuredGrid = 1;
490                }
491                dz = ddz;
492            }
493            nAxes++;
494#ifdef notdef
495            fprintf(stderr, "found delta %g %g %g\n", ddx, ddy, ddx);
496#endif
497        } else if (sscanf(line, "object %*d class regulararray count %d",
498                          &count[2]) == 1) {
499            // Z grid
500            isUniform = 0;
501        } else if (sscanf(line, "object %*d class array type %*s shape 3 rank 1"
502                          " items %d data follows", &nPoints) == 1) {
503            // XXX: Deprecated, invalid ordering of keywords
504            Tcl_AppendResult(interp, "Invalid DX: 'rank' keyword should precede 'shape'.",
505                             (char *)NULL);
506            return TCL_ERROR;
507        } else if (sscanf(line, "object %*d class array type %*s rank 1 shape 3"
508                          " items %d data follows", &nPoints) == 1) {
509#ifdef notdef
510                fprintf(stderr, "found class array type shape 3 nPoints=%d\n",
511                        nPoints);
512#endif
513            if (isUniform) {
514                hasVectors = 1;
515                if (GetUniformFieldVectors(interp, nPoints, count, &p, pend, fieldObjPtr)
516                    != TCL_OK) {
517                    return TCL_ERROR;
518                }
519            } else {
520                // This is a 2D point cloud in xy with a uniform zgrid
521                isUniform = 0;
522                nXYPoints = nPoints;
523                if (nXYPoints < 0) {
524                    sprintf(mesg, "bad # points %d", nXYPoints);
525                    Tcl_AppendResult(interp, mesg, (char *)NULL);
526                    return TCL_ERROR;
527                }
528                points = malloc(sizeof(double) * nXYPoints * 2);
529                if (GetPoints(interp, points, nXYPoints, &p, pend) != TCL_OK) {
530                    return TCL_ERROR;
531                }
532            }
533        } else if (sscanf(line, "object %*d class array type %*s rank 0"
534                          " items %d data follows", &nPoints) == 1) {
535#ifdef notdef
536            fprintf(stderr, "found class array type rank 0 nPoints=%d\n",
537                nPoints);
538#endif
539            if ((isUniform || isStructuredGrid) &&
540                nPoints != count[0]*count[1]*count[2]) {
541                sprintf(mesg, "inconsistent data: expected %d points"
542                        " but found %d points", count[0]*count[1]*count[2],
543                        nPoints);
544                Tcl_AppendResult(interp, mesg, (char *)NULL);
545                return TCL_ERROR;
546            } else if (!(isUniform || isStructuredGrid) &&
547                       nPoints != nXYPoints * count[2]) {
548                sprintf(mesg, "inconsistent data: expected %d points"
549                        " but found %d points", nXYPoints * count[2],
550                        nPoints);
551                Tcl_AppendResult(interp, mesg, (char *)NULL);
552                return TCL_ERROR;
553            }
554            if (isUniform) {
555                if (GetUniformFieldValues(interp, nPoints, count, &p, pend, fieldObjPtr)
556                    != TCL_OK) {
557                    return TCL_ERROR;
558                }
559            } else if (isStructuredGrid) {
560                if (GetStructuredGridFieldValues(interp, nPoints, &p, pend, fieldObjPtr)
561                    != TCL_OK) {
562                    return TCL_ERROR;
563                }
564            } else {
565                if (GetCloudFieldValues(interp, nXYPoints, count[2], &p, pend, fieldObjPtr)
566                    != TCL_OK) {
567                    return TCL_ERROR;
568                }
569            }
570        } else if (sscanf(line, "object %*d class array type %*s rank 0"
571                          " times %d data follows", &nPoints) == 1) {
572            /* Handle incorrect keyword: Do any tools make this mistake? */
573            Tcl_AppendResult(interp, "Invalid DX: found 'times' where 'items' keyword was expected.",
574                             (char *)NULL);
575            return TCL_ERROR;
576#ifdef notdef
577        } else {
578            fprintf(stderr, "unknown line (%.80s)\n", line);
579#endif
580        }
581    }
582
583    if (isUniform) {
584        if (nPoints > 1 && nAxes == 0) {
585            Tcl_AppendResult(interp, "Invalid DX file: uniform grid with no deltas found",
586                             (char *)NULL);
587            return TCL_ERROR;
588        }
589        if (nPoints > 1 && ((dx == dy) && (dx == dz) && (dx == 0.0))) {
590            sprintf(mesg, "Invalid deltas in DX file: %g %g %g", dx, dy, dz);
591            Tcl_AppendResult(interp, mesg, (char *)NULL);
592            return TCL_ERROR;
593        }
594        if (dx < 0.0 || dy < 0.0 || dz < 0.0) {
595            sprintf(mesg, "Negative deltas not supported in DX file: %g %g %g",
596                    dx, dy, dz);
597            Tcl_AppendResult(interp, mesg, (char *)NULL);
598            return TCL_ERROR;
599        }
600        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
601        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
602        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
603        Tcl_AppendToObj(objPtr, "DATASET STRUCTURED_POINTS\n", -1);
604        sprintf(mesg, "DIMENSIONS %d %d %d\n", count[0], count[1], count[2]);
605        Tcl_AppendToObj(objPtr, mesg, -1);
606        sprintf(mesg, "ORIGIN %g %g %g\n", origin[0], origin[1], origin[2]);
607        Tcl_AppendToObj(objPtr, mesg, -1);
608        sprintf(mesg, "SPACING %g %g %g\n", dx, dy, dz);
609        Tcl_AppendToObj(objPtr, mesg, -1);
610        sprintf(mesg, "POINT_DATA %d\n", nPoints);
611        Tcl_AppendToObj(objPtr, mesg, -1);
612        if (hasVectors) {
613            sprintf(mesg, "VECTORS %s double\n", name);
614            Tcl_AppendToObj(objPtr, mesg, -1);
615        } else {
616            sprintf(mesg, "SCALARS %s double 1\n", name);
617            Tcl_AppendToObj(objPtr, mesg, -1);
618            sprintf(mesg, "LOOKUP_TABLE default\n");
619            Tcl_AppendToObj(objPtr, mesg, -1);
620        }
621        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
622    } else if (isStructuredGrid) {
623#ifdef notdef
624        fprintf(stderr, "dv0 %g %g %g\n", dv0[0], dv0[1], dv0[2]);
625        fprintf(stderr, "dv1 %g %g %g\n", dv1[0], dv1[1], dv1[2]);
626        fprintf(stderr, "dv2 %g %g %g\n", dv2[0], dv2[1], dv2[2]);
627#endif
628        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
629        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
630        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
631        Tcl_AppendToObj(objPtr, "DATASET STRUCTURED_GRID\n", -1);
632        sprintf(mesg, "DIMENSIONS %d %d %d\n", count[0], count[1], count[2]);
633        Tcl_AppendToObj(objPtr, mesg, -1);
634        for (iz = 0; iz < count[2]; iz++) {
635            for (iy = 0; iy < count[1]; iy++) {
636                for (ix = 0; ix < count[0]; ix++) {
637                    double x, y, z;
638                    x = origin[0] + dv2[0] * iz + dv1[0] * iy + dv0[0] * ix;
639                    y = origin[1] + dv2[1] * iz + dv1[1] * iy + dv0[1] * ix;
640                    z = origin[2] + dv2[2] * iz + dv1[2] * iy + dv0[2] * ix;
641                    sprintf(mesg, "%g %g %g\n", x, y, z);
642                    Tcl_AppendToObj(pointsObjPtr, mesg, -1);
643                }
644            }
645        }
646        sprintf(mesg, "POINTS %d double\n", nPoints);
647        Tcl_AppendToObj(objPtr, mesg, -1);
648        Tcl_AppendObjToObj(objPtr, pointsObjPtr);
649        sprintf(mesg, "POINT_DATA %d\n", nPoints);
650        Tcl_AppendToObj(objPtr, mesg, -1);
651        sprintf(mesg, "SCALARS %s double 1\n", name);
652        Tcl_AppendToObj(objPtr, mesg, -1);
653        sprintf(mesg, "LOOKUP_TABLE default\n");
654        Tcl_AppendToObj(objPtr, mesg, -1);
655        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
656    } else {
657        /* Fill points.  Have to wait to do this since origin, delta can come after
658         * the point list in the file.
659         */
660        for (iz = 0; iz < count[2]; iz++) {
661            for (i = 0; i < nXYPoints; i++) {
662                sprintf(mesg, "%g %g %g\n", points[i*2], points[i*2+1], (origin[2] + dz * iz));
663                Tcl_AppendToObj(pointsObjPtr, mesg, -1);
664            }
665        }
666
667        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
668        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
669        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
670        Tcl_AppendToObj(objPtr, "DATASET UNSTRUCTURED_GRID\n", -1);
671        sprintf(mesg, "POINTS %d double\n", nPoints);
672        Tcl_AppendToObj(objPtr, mesg, -1);
673        Tcl_AppendObjToObj(objPtr, pointsObjPtr);
674#ifdef DO_WEDGES
675        {
676            double xmin, xmax, ymin, ymax;
677            xmin = xmax = points[0];
678            ymin = ymax = points[1];
679            for (i = 1; i < nXYPoints; i++) {
680                double x, y;
681                x = points[i*2];
682                y = points[i*2+1];
683                if (x < xmin) xmin = x;
684                if (x > xmax) xmax = x;
685                if (y < ymin) ymin = y;
686                if (y > ymax) ymax = y;
687            }
688
689            char fpts[128];
690            sprintf(fpts, "/tmp/tmppts%d", getpid());
691            char fcells[128];
692            sprintf(fcells, "/tmp/tmpcells%d", getpid());
693
694            FILE *ftmp = fopen(fpts, "w");
695            // save corners of bounding box first, to work around meshing
696            // problems in voronoi utility
697            int numBoundaryPoints = 4;
698
699            // Bump out the bounds by an epsilon to avoid problem when
700            // corner points are already nodes
701            double XEPS = (xmax - xmin) / 10.0f;
702            double YEPS = (ymax - ymin) / 10.0f;
703
704            fprintf(ftmp, "%g %g\n", xmin - XEPS, ymin - YEPS);
705            fprintf(ftmp, "%g %g\n", xmax + XEPS, ymin - YEPS);
706            fprintf(ftmp, "%g %g\n", xmax + XEPS, ymax + YEPS);
707            fprintf(ftmp, "%g %g\n", xmin - XEPS, ymax + YEPS);
708
709            for (i = 0; i < nXYPoints; i++) {
710                fprintf(ftmp, "%g %g\n", points[i*2], points[i*2+1]);
711            }
712            fclose(ftmp);
713#ifdef CHECK_WINDINGS
714            double normal[3];
715            normal[0] = normal[1] = 0.0;
716            if (dx >= 0) {
717                normal[2] = 1.0;
718            } else {
719                normal[2] = -1.0;
720            }
721#endif
722            char cmdstr[512];
723            sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
724            if (system(cmdstr) == 0) {
725                int c0, c1, c2;
726                ftmp = fopen(fcells, "r");
727                while (!feof(ftmp)) {
728                    if (fscanf(ftmp, "%d %d %d", &c0, &c1, &c2) == 3) {
729                        c0 -= numBoundaryPoints;
730                        c1 -= numBoundaryPoints;
731                        c2 -= numBoundaryPoints;
732                        // skip boundary points we added
733                        if (c0 >= 0 && c1 >= 0 && c2 >= 0) {
734#ifdef CHECK_WINDINGS
735                            /* Winding of base triangle should point to
736                               outside of cell using right hand rule */
737                            double v1[3], v2[3], c[3];
738                            v1[0] = points[c1*2] - points[c0*2];
739                            v1[1] = points[c1*2+1] - points[c0*2+1];
740                            v1[2] = 0;
741                            v2[0] = points[c2*2] - points[c0*2];
742                            v2[1] = points[c2*2+1] - points[c0*2+1];
743                            v2[2] = 0;
744                            Cross(v1, v2, c);
745                            if (Dot(normal, c) > 0) {
746                                int tmp = c0;
747                                c0 = c2;
748                                c2 = tmp;
749                            }
750#endif
751                            for (iz = 0; iz < count[2] - 1; iz++) {
752                                int base = nXYPoints * iz;
753                                int top = base + nXYPoints;
754                                nCells++;
755                                sprintf(mesg, "%d %d %d %d %d %d %d\n", 6,
756                                        base + c0, base + c1, base + c2,
757                                        top + c0, top + c1, top + c2);
758                                Tcl_AppendToObj(cellsObjPtr, mesg, -1);
759                            }
760                        }
761                    } else {
762                        break;
763                    }
764                }
765                fclose(ftmp);
766             } else {
767                fprintf(stderr, "voronoi mesher failed");
768            }
769            unlink(fpts);
770            unlink(fcells);
771        }
772        sprintf(mesg, "CELLS %d %d\n", nCells, 7*nCells);
773        Tcl_AppendToObj(objPtr, mesg, -1);
774        Tcl_AppendObjToObj(objPtr, cellsObjPtr);
775        sprintf(mesg, "CELL_TYPES %d\n", nCells);
776        Tcl_AppendToObj(objPtr, mesg, -1);
777        sprintf(mesg, "%d\n", 13);
778        for (i = 0; i < nCells; i++) {
779            Tcl_AppendToObj(objPtr, mesg, -1);
780        }
781#endif
782        if (points != NULL) {
783            free(points);
784        }
785        sprintf(mesg, "POINT_DATA %d\n", nPoints);
786        Tcl_AppendToObj(objPtr, mesg, -1);
787        sprintf(mesg, "SCALARS %s double 1\n", name);
788        Tcl_AppendToObj(objPtr, mesg, -1);
789        sprintf(mesg, "LOOKUP_TABLE default\n");
790        Tcl_AppendToObj(objPtr, mesg, -1);
791        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
792    }
793
794    Tcl_DecrRefCount(pointsObjPtr);
795    Tcl_DecrRefCount(cellsObjPtr);
796    Tcl_DecrRefCount(fieldObjPtr);
797    Tcl_SetObjResult(interp, objPtr);
798    return TCL_OK;
799}
800
801/*
802 * ------------------------------------------------------------------------
803 *  RpDxToVtk_Init --
804 *
805 *  Invoked when the Rappture GUI library is being initialized
806 *  to install the "DxToVtk" command into the interpreter.
807 *
808 *  Returns TCL_OK if successful, or TCL_ERROR (along with an error
809 *  message in the interp) if anything goes wrong.
810 * ------------------------------------------------------------------------
811 */
812int
813RpDxToVtk_Init(Tcl_Interp *interp)
814{
815    /* install the widget command */
816    Tcl_CreateObjCommand(interp, "Rappture::DxToVtk", DxToVtkCmd,
817        NULL, NULL);
818    return TCL_OK;
819}
Note: See TracBrowser for help on using the repository browser.