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

Last change on this file since 4502 was 4502, checked in by ldelgass, 10 years ago

Cleanups. Remove unused variable, move check for <DX> header to C routine,
don't overwrite _dim written by VTK reader.
File size: 26.5 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 <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; /* Suppress compiler warning. */
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; /* Suppress compiler warning. */
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#ifdef notdef
504                fprintf(stderr, "found class array type shape 3 nPoints=%d\n",
505                        nPoints);
506#endif
507            if (isUniform) {
508                hasVectors = 1;
509                if (GetUniformFieldVectors(interp, nPoints, count, &p, pend, fieldObjPtr)
510                    != TCL_OK) {
511                    return TCL_ERROR;
512                }
513            } else {
514                // This is a 2D point cloud in xy with a uniform zgrid
515                isUniform = 0;
516                nXYPoints = nPoints;
517                if (nXYPoints < 0) {
518                    sprintf(mesg, "bad # points %d", nXYPoints);
519                    Tcl_AppendResult(interp, mesg, (char *)NULL);
520                    return TCL_ERROR;
521                }
522                points = malloc(sizeof(double) * nXYPoints * 2);
523                if (GetPoints(interp, points, nXYPoints, &p, pend) != TCL_OK) {
524                    return TCL_ERROR;
525                }
526            }
527        } else if (sscanf(line, "object %*d class array type %*s rank 1 shape 3"
528                          " items %d data follows", &nPoints) == 1) {
529#ifdef notdef
530                fprintf(stderr, "found class array type shape 3 nPoints=%d\n",
531                        nPoints);
532#endif
533            if (isUniform) {
534                hasVectors = 1;
535                if (GetUniformFieldVectors(interp, nPoints, count, &p, pend, fieldObjPtr)
536                    != TCL_OK) {
537                    return TCL_ERROR;
538                }
539            } else {
540                // This is a 2D point cloud in xy with a uniform zgrid
541                isUniform = 0;
542                nXYPoints = nPoints;
543                if (nXYPoints < 0) {
544                    sprintf(mesg, "bad # points %d", nXYPoints);
545                    Tcl_AppendResult(interp, mesg, (char *)NULL);
546                    return TCL_ERROR;
547                }
548                points = malloc(sizeof(double) * nXYPoints * 2);
549                if (GetPoints(interp, points, nXYPoints, &p, pend) != TCL_OK) {
550                    return TCL_ERROR;
551                }
552            }
553        } else if (sscanf(line, "object %*d class array type %*s rank 0"
554                          " %*s %d data follows", &nPoints) == 1) {
555#ifdef notdef
556            fprintf(stderr, "found class array type rank 0 nPoints=%d\n",
557                nPoints);
558#endif
559            if ((isUniform || isStructuredGrid) &&
560                nPoints != count[0]*count[1]*count[2]) {
561                sprintf(mesg, "inconsistent data: expected %d points"
562                        " but found %d points", count[0]*count[1]*count[2],
563                        nPoints);
564                Tcl_AppendResult(interp, mesg, (char *)NULL);
565                return TCL_ERROR;
566            } else if (!(isUniform || isStructuredGrid) &&
567                       nPoints != nXYPoints * count[2]) {
568                sprintf(mesg, "inconsistent data: expected %d points"
569                        " but found %d points", nXYPoints * count[2],
570                        nPoints);
571                Tcl_AppendResult(interp, mesg, (char *)NULL);
572                return TCL_ERROR;
573            }
574            if (isUniform) {
575                if (GetUniformFieldValues(interp, nPoints, count, &p, pend, fieldObjPtr)
576                    != TCL_OK) {
577                    return TCL_ERROR;
578                }
579            } else if (isStructuredGrid) {
580                if (GetStructuredGridFieldValues(interp, nPoints, &p, pend, fieldObjPtr)
581                    != TCL_OK) {
582                    return TCL_ERROR;
583                }
584            } else {
585                if (GetCloudFieldValues(interp, nXYPoints, count[2], &p, pend, fieldObjPtr)
586                    != TCL_OK) {
587                    return TCL_ERROR;
588                }
589            }
590#ifdef notdef
591        } else {
592            fprintf(stderr, "unknown line (%.80s)\n", line);
593#endif
594        }
595    }
596
597    if (isUniform) {
598        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
599        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
600        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
601        Tcl_AppendToObj(objPtr, "DATASET STRUCTURED_POINTS\n", -1);
602        sprintf(mesg, "DIMENSIONS %d %d %d\n", count[0], count[1], count[2]);
603        Tcl_AppendToObj(objPtr, mesg, -1);
604        sprintf(mesg, "ORIGIN %g %g %g\n", origin[0], origin[1], origin[2]);
605        Tcl_AppendToObj(objPtr, mesg, -1);
606        sprintf(mesg, "SPACING %g %g %g\n", dx, dy, dz);
607        Tcl_AppendToObj(objPtr, mesg, -1);
608        sprintf(mesg, "POINT_DATA %d\n", nPoints);
609        Tcl_AppendToObj(objPtr, mesg, -1);
610        if (hasVectors) {
611            sprintf(mesg, "VECTORS %s double\n", name);
612            Tcl_AppendToObj(objPtr, mesg, -1);
613        } else {
614            sprintf(mesg, "SCALARS %s double 1\n", name);
615            Tcl_AppendToObj(objPtr, mesg, -1);
616            sprintf(mesg, "LOOKUP_TABLE default\n");
617            Tcl_AppendToObj(objPtr, mesg, -1);
618        }
619        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
620    } else if (isStructuredGrid) {
621#ifdef notdef
622        fprintf(stderr, "dv0 %g %g %g\n", dv0[0], dv0[1], dv0[2]);
623        fprintf(stderr, "dv1 %g %g %g\n", dv1[0], dv1[1], dv1[2]);
624        fprintf(stderr, "dv2 %g %g %g\n", dv2[0], dv2[1], dv2[2]);
625#endif
626        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
627        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
628        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
629        Tcl_AppendToObj(objPtr, "DATASET STRUCTURED_GRID\n", -1);
630        sprintf(mesg, "DIMENSIONS %d %d %d\n", count[0], count[1], count[2]);
631        Tcl_AppendToObj(objPtr, mesg, -1);
632        for (iz = 0; iz < count[2]; iz++) {
633            for (iy = 0; iy < count[1]; iy++) {
634                for (ix = 0; ix < count[0]; ix++) {
635                    double x, y, z;
636                    x = origin[0] + dv2[0] * iz + dv1[0] * iy + dv0[0] * ix;
637                    y = origin[1] + dv2[1] * iz + dv1[1] * iy + dv0[1] * ix;
638                    z = origin[2] + dv2[2] * iz + dv1[2] * iy + dv0[2] * ix;
639                    sprintf(mesg, "%g %g %g\n", x, y, z);
640                    Tcl_AppendToObj(pointsObjPtr, mesg, -1);
641                }
642            }
643        }
644        sprintf(mesg, "POINTS %d double\n", nPoints);
645        Tcl_AppendToObj(objPtr, mesg, -1);
646        Tcl_AppendObjToObj(objPtr, pointsObjPtr);
647        sprintf(mesg, "POINT_DATA %d\n", nPoints);
648        Tcl_AppendToObj(objPtr, mesg, -1);
649        sprintf(mesg, "SCALARS %s double 1\n", name);
650        Tcl_AppendToObj(objPtr, mesg, -1);
651        sprintf(mesg, "LOOKUP_TABLE default\n");
652        Tcl_AppendToObj(objPtr, mesg, -1);
653        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
654    } else {
655        /* Fill points.  Have to wait to do this since origin, delta can come after
656         * the point list in the file.
657         */
658        for (iz = 0; iz < count[2]; iz++) {
659            for (i = 0; i < nXYPoints; i++) {
660                sprintf(mesg, "%g %g %g\n", points[i*2], points[i*2+1], (origin[2] + dz * iz));
661                Tcl_AppendToObj(pointsObjPtr, mesg, -1);
662            }
663        }
664
665        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
666        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
667        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
668        Tcl_AppendToObj(objPtr, "DATASET UNSTRUCTURED_GRID\n", -1);
669        sprintf(mesg, "POINTS %d double\n", nPoints);
670        Tcl_AppendToObj(objPtr, mesg, -1);
671        Tcl_AppendObjToObj(objPtr, pointsObjPtr);
672#ifdef DO_WEDGES
673        {
674            double xmin, xmax, ymin, ymax;
675            xmin = xmax = points[0];
676            ymin = ymax = points[1];
677            for (i = 1; i < nXYPoints; i++) {
678                double x, y;
679                x = points[i*2];
680                y = points[i*2+1];
681                if (x < xmin) xmin = x;
682                if (x > xmax) xmax = x;
683                if (y < ymin) ymin = y;
684                if (y > ymax) ymax = y;
685            }
686
687            char fpts[128];
688            sprintf(fpts, "/tmp/tmppts%d", getpid());
689            char fcells[128];
690            sprintf(fcells, "/tmp/tmpcells%d", getpid());
691
692            FILE *ftmp = fopen(fpts, "w");
693            // save corners of bounding box first, to work around meshing
694            // problems in voronoi utility
695            int numBoundaryPoints = 4;
696
697            // Bump out the bounds by an epsilon to avoid problem when
698            // corner points are already nodes
699            double XEPS = (xmax - xmin) / 10.0f;
700            double YEPS = (ymax - ymin) / 10.0f;
701
702            fprintf(ftmp, "%g %g\n", xmin - XEPS, ymin - YEPS);
703            fprintf(ftmp, "%g %g\n", xmax + XEPS, ymin - YEPS);
704            fprintf(ftmp, "%g %g\n", xmax + XEPS, ymax + YEPS);
705            fprintf(ftmp, "%g %g\n", xmin - XEPS, ymax + YEPS);
706
707            for (i = 0; i < nXYPoints; i++) {
708                fprintf(ftmp, "%g %g\n", points[i*2], points[i*2+1]);
709            }
710            fclose(ftmp);
711#ifdef CHECK_WINDINGS
712            double normal[3];
713            normal[0] = normal[1] = 0.0;
714            if (dx >= 0) {
715                normal[2] = 1.0;
716            } else {
717                normal[2] = -1.0;
718            }
719#endif
720            char cmdstr[512];
721            sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
722            if (system(cmdstr) == 0) {
723                int c0, c1, c2;
724                ftmp = fopen(fcells, "r");
725                while (!feof(ftmp)) {
726                    if (fscanf(ftmp, "%d %d %d", &c0, &c1, &c2) == 3) {
727                        c0 -= numBoundaryPoints;
728                        c1 -= numBoundaryPoints;
729                        c2 -= numBoundaryPoints;
730                        // skip boundary points we added
731                        if (c0 >= 0 && c1 >= 0 && c2 >= 0) {
732#ifdef CHECK_WINDINGS
733                            /* Winding of base triangle should point to
734                               outside of cell using right hand rule */
735                            double v1[3], v2[3], c[3];
736                            v1[0] = points[c1*2] - points[c0*2];
737                            v1[1] = points[c1*2+1] - points[c0*2+1];
738                            v1[2] = 0;
739                            v2[0] = points[c2*2] - points[c0*2];
740                            v2[1] = points[c2*2+1] - points[c0*2+1];
741                            v2[2] = 0;
742                            Cross(v1, v2, c);
743                            if (Dot(normal, c) > 0) {
744                                int tmp = c0;
745                                c0 = c2;
746                                c2 = tmp;
747                            }
748#endif
749                            for (iz = 0; iz < count[2] - 1; iz++) {
750                                int base = nXYPoints * iz;
751                                int top = base + nXYPoints;
752                                nCells++;
753                                sprintf(mesg, "%d %d %d %d %d %d %d\n", 6,
754                                        base + c0, base + c1, base + c2,
755                                        top + c0, top + c1, top + c2);
756                                Tcl_AppendToObj(cellsObjPtr, mesg, -1);
757                            }
758                        }
759                    } else {
760                        break;
761                    }
762                }
763                fclose(ftmp);
764             } else {
765                fprintf(stderr, "voronoi mesher failed");
766            }
767            unlink(fpts);
768            unlink(fcells);
769        }
770        sprintf(mesg, "CELLS %d %d\n", nCells, 7*nCells);
771        Tcl_AppendToObj(objPtr, mesg, -1);
772        Tcl_AppendObjToObj(objPtr, cellsObjPtr);
773        sprintf(mesg, "CELL_TYPES %d\n", nCells);
774        Tcl_AppendToObj(objPtr, mesg, -1);
775        sprintf(mesg, "%d\n", 13);
776        for (i = 0; i < nCells; i++) {
777            Tcl_AppendToObj(objPtr, mesg, -1);
778        }
779#endif
780        if (points != NULL) {
781            free(points);
782        }
783        sprintf(mesg, "POINT_DATA %d\n", nPoints);
784        Tcl_AppendToObj(objPtr, mesg, -1);
785        sprintf(mesg, "SCALARS %s double 1\n", name);
786        Tcl_AppendToObj(objPtr, mesg, -1);
787        sprintf(mesg, "LOOKUP_TABLE default\n");
788        Tcl_AppendToObj(objPtr, mesg, -1);
789        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
790    }
791
792    Tcl_DecrRefCount(pointsObjPtr);
793    Tcl_DecrRefCount(cellsObjPtr);
794    Tcl_DecrRefCount(fieldObjPtr);
795    Tcl_SetObjResult(interp, objPtr);
796    return TCL_OK;
797}
798
799/*
800 * ------------------------------------------------------------------------
801 *  RpDxToVtk_Init --
802 *
803 *  Invoked when the Rappture GUI library is being initialized
804 *  to install the "DxToVtk" command into the interpreter.
805 *
806 *  Returns TCL_OK if successful, or TCL_ERROR (along with an error
807 *  message in the interp) if anything goes wrong.
808 * ------------------------------------------------------------------------
809 */
810int
811RpDxToVtk_Init(Tcl_Interp *interp)
812{
813    /* install the widget command */
814    Tcl_CreateObjCommand(interp, "Rappture::DxToVtk", DxToVtkCmd,
815        NULL, NULL);
816    return TCL_OK;
817}
Note: See TracBrowser for help on using the repository browser.