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

Last change on this file since 4515 was 4515, checked in by gah, 10 years ago

add new utilities (squeezer, tweener), update RpDxToVtk?.c, experimental RpDicomToVtk?.c not used yet
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.