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

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

Add support for converting DX with vectors to VTK
File size: 26.4 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    }
408    pointsObjPtr = Tcl_NewStringObj("", -1);
409    cellsObjPtr = Tcl_NewStringObj("", -1);
410    fieldObjPtr = Tcl_NewStringObj("", -1);
411    for (p = string, pend = p + length; p < pend; /*empty*/) {
412        char *line;
413        double ddx, ddy, ddz;
414
415        line = GetLine(&p, pend);
416        if (line >= pend) {
417            break;                        /* EOF */
418        }
419        if ((line[0] == '#') || (line[0] == '\n')) {
420            continue;                        /* Skip blank or comment lines. */
421        }
422        if (sscanf(line, "object %*d class gridpositions counts %d %d %d",
423                   count, count + 1, count + 2) == 3) {
424            isUniform = 1;
425            if ((count[0] < 0) || (count[1] < 0) || (count[2] < 0)) {
426                sprintf(mesg, "invalid grid size: x=%d, y=%d, z=%d",
427                        count[0], count[1], count[2]);
428                Tcl_AppendResult(interp, mesg, (char *)NULL);
429                return TCL_ERROR;
430            }
431#ifdef notdef
432            fprintf(stderr, "found gridpositions counts %d %d %d\n",
433                    count[0], count[1], count[2]);
434#endif
435        } else if (sscanf(line, "origin %lg %lg %lg", origin, origin + 1,
436                origin + 2) == 3) {
437            /* Found origin. */
438#ifdef notdef
439            fprintf(stderr, "found origin %g %g %g\n",
440                    origin[0], origin[1], origin[2]);
441#endif
442        } else if (sscanf(line, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
443            if (nAxes == 3) {
444                Tcl_AppendResult(interp, "too many delta statements",
445                        (char *)NULL);
446                return TCL_ERROR;
447            }
448            switch (nAxes) {
449            case 0:
450                dv0[0] = ddx;
451                dv0[1] = ddy;
452                dv0[2] = ddz;
453                break;
454            case 1:
455                dv1[0] = ddx;
456                dv1[1] = ddy;
457                dv1[2] = ddz;
458                break;
459            case 2:
460                dv2[0] = ddx;
461                dv2[1] = ddy;
462                dv2[2] = ddz;
463                break;
464            default:
465                break;
466            }
467
468            if (ddx != 0.0) {
469                if (ddy != 0.0 || ddz != 0.0) {
470                    /* Not axis aligned or skewed grid */
471                    isUniform = 0;
472                    isStructuredGrid = 1;
473                }
474                dx = ddx;
475            } else if (ddy != 0.0) {
476                if (ddx != 0.0 || ddz != 0.0) {
477                    /* Not axis aligned or skewed grid */
478                    isUniform = 0;
479                    isStructuredGrid = 1;
480                }
481                dy = ddy;
482            } else if (ddz != 0.0) {
483                if (ddx != 0.0 || ddy != 0.0) {
484                    /* Not axis aligned or skewed grid */
485                    isUniform = 0;
486                    isStructuredGrid = 1;
487                }
488                dz = ddz;
489            }
490            nAxes++;
491#ifdef notdef
492            fprintf(stderr, "found delta %g %g %g\n", ddx, ddy, ddx);
493#endif
494        } else if (sscanf(line, "object %*d class regulararray count %d",
495                          &count[2]) == 1) {
496            // Z grid
497            isUniform = 0;
498        } else if (sscanf(line, "object %*d class array type %*s shape 3 rank 1"
499                          " items %d data follows", &nPoints) == 1) {
500#ifdef notdef
501                fprintf(stderr, "found class array type shape 3 nPoints=%d\n",
502                        nPoints);
503#endif
504            if (isUniform) {
505                hasVectors = 1;
506                if (GetUniformFieldVectors(interp, nPoints, count, &p, pend, fieldObjPtr)
507                    != TCL_OK) {
508                    return TCL_ERROR;
509                }
510            } else {
511                // This is a 2D point cloud in xy with a uniform zgrid
512                isUniform = 0;
513                nXYPoints = nPoints;
514                if (nXYPoints < 0) {
515                    sprintf(mesg, "bad # points %d", nXYPoints);
516                    Tcl_AppendResult(interp, mesg, (char *)NULL);
517                    return TCL_ERROR;
518                }
519                points = malloc(sizeof(double) * nXYPoints * 2);
520                if (GetPoints(interp, points, nXYPoints, &p, pend) != TCL_OK) {
521                    return TCL_ERROR;
522                }
523            }
524        } else if (sscanf(line, "object %*d class array type %*s rank 1 shape 3"
525                          " items %d data follows", &nPoints) == 1) {
526#ifdef notdef
527                fprintf(stderr, "found class array type shape 3 nPoints=%d\n",
528                        nPoints);
529#endif
530            if (isUniform) {
531                hasVectors = 1;
532                if (GetUniformFieldVectors(interp, nPoints, count, &p, pend, fieldObjPtr)
533                    != TCL_OK) {
534                    return TCL_ERROR;
535                }
536            } else {
537                // This is a 2D point cloud in xy with a uniform zgrid
538                isUniform = 0;
539                nXYPoints = nPoints;
540                if (nXYPoints < 0) {
541                    sprintf(mesg, "bad # points %d", nXYPoints);
542                    Tcl_AppendResult(interp, mesg, (char *)NULL);
543                    return TCL_ERROR;
544                }
545                points = malloc(sizeof(double) * nXYPoints * 2);
546                if (GetPoints(interp, points, nXYPoints, &p, pend) != TCL_OK) {
547                    return TCL_ERROR;
548                }
549            }
550        } else if (sscanf(line, "object %*d class array type %*s rank 0"
551                          " %*s %d data follows", &nPoints) == 1) {
552#ifdef notdef
553            fprintf(stderr, "found class array type rank 0 nPoints=%d\n",
554                nPoints);
555#endif
556            if ((isUniform || isStructuredGrid) &&
557                nPoints != count[0]*count[1]*count[2]) {
558                sprintf(mesg, "inconsistent data: expected %d points"
559                        " but found %d points", count[0]*count[1]*count[2],
560                        nPoints);
561                Tcl_AppendResult(interp, mesg, (char *)NULL);
562                return TCL_ERROR;
563            } else if (!(isUniform || isStructuredGrid) &&
564                       nPoints != nXYPoints * count[2]) {
565                sprintf(mesg, "inconsistent data: expected %d points"
566                        " but found %d points", nXYPoints * count[2],
567                        nPoints);
568                Tcl_AppendResult(interp, mesg, (char *)NULL);
569                return TCL_ERROR;
570            }
571            if (isUniform) {
572                if (GetUniformFieldValues(interp, nPoints, count, &p, pend, fieldObjPtr)
573                    != TCL_OK) {
574                    return TCL_ERROR;
575                }
576            } else if (isStructuredGrid) {
577                if (GetStructuredGridFieldValues(interp, nPoints, &p, pend, fieldObjPtr)
578                    != TCL_OK) {
579                    return TCL_ERROR;
580                }
581            } else {
582                if (GetCloudFieldValues(interp, nXYPoints, count[2], &p, pend, fieldObjPtr)
583                    != TCL_OK) {
584                    return TCL_ERROR;
585                }
586            }
587#ifdef notdef
588        } else {
589            fprintf(stderr, "unknown line (%.80s)\n", line);
590#endif
591        }
592    }
593
594    if (isUniform) {
595        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
596        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
597        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
598        Tcl_AppendToObj(objPtr, "DATASET STRUCTURED_POINTS\n", -1);
599        sprintf(mesg, "DIMENSIONS %d %d %d\n", count[0], count[1], count[2]);
600        Tcl_AppendToObj(objPtr, mesg, -1);
601        sprintf(mesg, "ORIGIN %g %g %g\n", origin[0], origin[1], origin[2]);
602        Tcl_AppendToObj(objPtr, mesg, -1);
603        sprintf(mesg, "SPACING %g %g %g\n", dx, dy, dz);
604        Tcl_AppendToObj(objPtr, mesg, -1);
605        sprintf(mesg, "POINT_DATA %d\n", nPoints);
606        Tcl_AppendToObj(objPtr, mesg, -1);
607        if (hasVectors) {
608            sprintf(mesg, "VECTORS %s double\n", name);
609            Tcl_AppendToObj(objPtr, mesg, -1);
610        } else {
611            sprintf(mesg, "SCALARS %s double 1\n", name);
612            Tcl_AppendToObj(objPtr, mesg, -1);
613            sprintf(mesg, "LOOKUP_TABLE default\n");
614            Tcl_AppendToObj(objPtr, mesg, -1);
615        }
616        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
617    } else if (isStructuredGrid) {
618#ifdef notdef
619        fprintf(stderr, "dv0 %g %g %g\n", dv0[0], dv0[1], dv0[2]);
620        fprintf(stderr, "dv1 %g %g %g\n", dv1[0], dv1[1], dv1[2]);
621        fprintf(stderr, "dv2 %g %g %g\n", dv2[0], dv2[1], dv2[2]);
622#endif
623        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
624        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
625        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
626        Tcl_AppendToObj(objPtr, "DATASET STRUCTURED_GRID\n", -1);
627        sprintf(mesg, "DIMENSIONS %d %d %d\n", count[0], count[1], count[2]);
628        Tcl_AppendToObj(objPtr, mesg, -1);
629        for (iz = 0; iz < count[2]; iz++) {
630            for (iy = 0; iy < count[1]; iy++) {
631                for (ix = 0; ix < count[0]; ix++) {
632                    double x, y, z;
633                    x = origin[0] + dv2[0] * iz + dv1[0] * iy + dv0[0] * ix;
634                    y = origin[1] + dv2[1] * iz + dv1[1] * iy + dv0[1] * ix;
635                    z = origin[2] + dv2[2] * iz + dv1[2] * iy + dv0[2] * ix;
636                    sprintf(mesg, "%g %g %g\n", x, y, z);
637                    Tcl_AppendToObj(pointsObjPtr, mesg, -1);
638                }
639            }
640        }
641        sprintf(mesg, "POINTS %d double\n", nPoints);
642        Tcl_AppendToObj(objPtr, mesg, -1);
643        Tcl_AppendObjToObj(objPtr, pointsObjPtr);
644        sprintf(mesg, "POINT_DATA %d\n", nPoints);
645        Tcl_AppendToObj(objPtr, mesg, -1);
646        sprintf(mesg, "SCALARS %s double 1\n", name);
647        Tcl_AppendToObj(objPtr, mesg, -1);
648        sprintf(mesg, "LOOKUP_TABLE default\n");
649        Tcl_AppendToObj(objPtr, mesg, -1);
650        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
651    } else {
652        /* Fill points.  Have to wait to do this since origin, delta can come after
653         * the point list in the file.
654         */
655        for (iz = 0; iz < count[2]; iz++) {
656            for (i = 0; i < nXYPoints; i++) {
657                sprintf(mesg, "%g %g %g\n", points[i*2], points[i*2+1], (origin[2] + dz * iz));
658                Tcl_AppendToObj(pointsObjPtr, mesg, -1);
659            }
660        }
661
662        objPtr = Tcl_NewStringObj("# vtk DataFile Version 2.0\n", -1);
663        Tcl_AppendToObj(objPtr, "Converted from DX file\n", -1);
664        Tcl_AppendToObj(objPtr, "ASCII\n", -1);
665        Tcl_AppendToObj(objPtr, "DATASET UNSTRUCTURED_GRID\n", -1);
666        sprintf(mesg, "POINTS %d double\n", nPoints);
667        Tcl_AppendToObj(objPtr, mesg, -1);
668        Tcl_AppendObjToObj(objPtr, pointsObjPtr);
669#ifdef DO_WEDGES
670        {
671            double xmin, xmax, ymin, ymax;
672            xmin = xmax = points[0];
673            ymin = ymax = points[1];
674            for (i = 1; i < nXYPoints; i++) {
675                double x, y;
676                x = points[i*2];
677                y = points[i*2+1];
678                if (x < xmin) xmin = x;
679                if (x > xmax) xmax = x;
680                if (y < ymin) ymin = y;
681                if (y > ymax) ymax = y;
682            }
683
684            char fpts[128];
685            sprintf(fpts, "/tmp/tmppts%d", getpid());
686            char fcells[128];
687            sprintf(fcells, "/tmp/tmpcells%d", getpid());
688
689            FILE *ftmp = fopen(fpts, "w");
690            // save corners of bounding box first, to work around meshing
691            // problems in voronoi utility
692            int numBoundaryPoints = 4;
693
694            // Bump out the bounds by an epsilon to avoid problem when
695            // corner points are already nodes
696            double XEPS = (xmax - xmin) / 10.0f;
697            double YEPS = (ymax - ymin) / 10.0f;
698
699            fprintf(ftmp, "%g %g\n", xmin - XEPS, ymin - YEPS);
700            fprintf(ftmp, "%g %g\n", xmax + XEPS, ymin - YEPS);
701            fprintf(ftmp, "%g %g\n", xmax + XEPS, ymax + YEPS);
702            fprintf(ftmp, "%g %g\n", xmin - XEPS, ymax + YEPS);
703
704            for (i = 0; i < nXYPoints; i++) {
705                fprintf(ftmp, "%g %g\n", points[i*2], points[i*2+1]);
706            }
707            fclose(ftmp);
708#ifdef CHECK_WINDINGS
709            double normal[3];
710            normal[0] = normal[1] = 0.0;
711            if (dx >= 0) {
712                normal[2] = 1.0;
713            } else {
714                normal[2] = -1.0;
715            }
716#endif
717            char cmdstr[512];
718            sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
719            if (system(cmdstr) == 0) {
720                int c0, c1, c2;
721                ftmp = fopen(fcells, "r");
722                while (!feof(ftmp)) {
723                    if (fscanf(ftmp, "%d %d %d", &c0, &c1, &c2) == 3) {
724                        c0 -= numBoundaryPoints;
725                        c1 -= numBoundaryPoints;
726                        c2 -= numBoundaryPoints;
727                        // skip boundary points we added
728                        if (c0 >= 0 && c1 >= 0 && c2 >= 0) {
729#ifdef CHECK_WINDINGS
730                            /* Winding of base triangle should point to
731                               outside of cell using right hand rule */
732                            double v1[3], v2[3], c[3];
733                            v1[0] = points[c1*2] - points[c0*2];
734                            v1[1] = points[c1*2+1] - points[c0*2+1];
735                            v1[2] = 0;
736                            v2[0] = points[c2*2] - points[c0*2];
737                            v2[1] = points[c2*2+1] - points[c0*2+1];
738                            v2[2] = 0;
739                            Cross(v1, v2, c);
740                            if (Dot(normal, c) > 0) {
741                                int tmp = c0;
742                                c0 = c2;
743                                c2 = tmp;
744                            }
745#endif
746                            for (iz = 0; iz < count[2] - 1; iz++) {
747                                int base = nXYPoints * iz;
748                                int top = base + nXYPoints;
749                                nCells++;
750                                sprintf(mesg, "%d %d %d %d %d %d %d\n", 6,
751                                        base + c0, base + c1, base + c2,
752                                        top + c0, top + c1, top + c2);
753                                Tcl_AppendToObj(cellsObjPtr, mesg, -1);
754                            }
755                        }
756                    } else {
757                        break;
758                    }
759                }
760                fclose(ftmp);
761             } else {
762                fprintf(stderr, "voronoi mesher failed");
763            }
764            unlink(fpts);
765            unlink(fcells);
766        }
767        sprintf(mesg, "CELLS %d %d\n", nCells, 7*nCells);
768        Tcl_AppendToObj(objPtr, mesg, -1);
769        Tcl_AppendObjToObj(objPtr, cellsObjPtr);
770        sprintf(mesg, "CELL_TYPES %d\n", nCells);
771        Tcl_AppendToObj(objPtr, mesg, -1);
772        sprintf(mesg, "%d\n", 13);
773        for (i = 0; i < nCells; i++) {
774            Tcl_AppendToObj(objPtr, mesg, -1);
775        }
776#endif
777        if (points != NULL) {
778            free(points);
779        }
780        sprintf(mesg, "POINT_DATA %d\n", nPoints);
781        Tcl_AppendToObj(objPtr, mesg, -1);
782        sprintf(mesg, "SCALARS %s double 1\n", name);
783        Tcl_AppendToObj(objPtr, mesg, -1);
784        sprintf(mesg, "LOOKUP_TABLE default\n");
785        Tcl_AppendToObj(objPtr, mesg, -1);
786        Tcl_AppendObjToObj(objPtr, fieldObjPtr);
787    }
788
789    Tcl_DecrRefCount(pointsObjPtr);
790    Tcl_DecrRefCount(cellsObjPtr);
791    Tcl_DecrRefCount(fieldObjPtr);
792    Tcl_SetObjResult(interp, objPtr);
793    return TCL_OK;
794}
795
796/*
797 * ------------------------------------------------------------------------
798 *  RpDxToVtk_Init --
799 *
800 *  Invoked when the Rappture GUI library is being initialized
801 *  to install the "DxToVtk" command into the interpreter.
802 *
803 *  Returns TCL_OK if successful, or TCL_ERROR (along with an error
804 *  message in the interp) if anything goes wrong.
805 * ------------------------------------------------------------------------
806 */
807int
808RpDxToVtk_Init(Tcl_Interp *interp)
809{
810    /* install the widget command */
811    Tcl_CreateObjCommand(interp, "Rappture::DxToVtk", DxToVtkCmd,
812        NULL, NULL);
813    return TCL_OK;
814}
Note: See TracBrowser for help on using the repository browser.