source: vtkvis/trunk/Streamlines.cpp @ 5073

Last change on this file since 5073 was 5073, checked in by ldelgass, 5 years ago

merge r5072 from release branch

  • Property svn:eol-style set to native
File size: 56.3 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2/*
3 * Copyright (C) 2004-2012  HUBzero Foundation, LLC
4 *
5 * Author: Leif Delgass <ldelgass@purdue.edu>
6 */
7
8#include <cstdlib>
9#include <ctime>
10#include <cfloat>
11#include <cmath>
12#include <cstring>
13
14#include <vtkMath.h>
15#include <vtkActor.h>
16#include <vtkProperty.h>
17#include <vtkPoints.h>
18#include <vtkCellArray.h>
19#include <vtkPolyLine.h>
20#include <vtkRegularPolygonSource.h>
21#include <vtkPointData.h>
22#include <vtkCellData.h>
23#include <vtkCellDataToPointData.h>
24#include <vtkDelaunay2D.h>
25#include <vtkDelaunay3D.h>
26#include <vtkUnstructuredGrid.h>
27#include <vtkPolygon.h>
28#include <vtkPolyData.h>
29#include <vtkTubeFilter.h>
30#include <vtkRibbonFilter.h>
31#include <vtkTransform.h>
32#include <vtkTransformPolyDataFilter.h>
33#include <vtkVertexGlyphFilter.h>
34#include <vtkMaskPoints.h>
35
36#include "Streamlines.h"
37#include "Renderer.h"
38#include "Trace.h"
39
40using namespace VtkVis;
41
42Streamlines::Streamlines() :
43    GraphicsObject(),
44    _lineType(LINES),
45    _colorMap(NULL),
46    _colorMode(COLOR_BY_VECTOR_MAGNITUDE),
47    _colorFieldType(DataSet::POINT_DATA),
48    _seedVisible(true),
49    _renderer(NULL),
50    _dataScale(1)
51{
52    _faceCulling = true;
53    _color[0] = 1.0f;
54    _color[1] = 1.0f;
55    _color[2] = 1.0f;
56    _seedColor[0] = 1.0f;
57    _seedColor[1] = 1.0f;
58    _seedColor[2] = 1.0f;
59    _colorFieldRange[0] = DBL_MAX;
60    _colorFieldRange[1] = -DBL_MAX;
61    vtkMath::RandomSeed((int)time(NULL));
62    srand((unsigned int)time(NULL));
63}
64
65Streamlines::~Streamlines()
66{
67#ifdef WANT_TRACE
68    if (_dataSet != NULL)
69        TRACE("Deleting Streamlines for %s", _dataSet->getName().c_str());
70    else
71        TRACE("Deleting Streamlines with NULL DataSet");
72#endif
73}
74
75void Streamlines::setDataSet(DataSet *dataSet,
76                             Renderer *renderer)
77{
78    if (_dataSet != dataSet) {
79        _dataSet = dataSet;
80
81        _renderer = renderer;
82
83        if (renderer->getUseCumulativeRange()) {
84            renderer->getCumulativeDataRange(_dataRange,
85                                             _dataSet->getActiveScalarsName(),
86                                             1);
87            renderer->getCumulativeDataRange(_vectorMagnitudeRange,
88                                             _dataSet->getActiveVectorsName(),
89                                             3);
90            for (int i = 0; i < 3; i++) {
91                renderer->getCumulativeDataRange(_vectorComponentRange[i],
92                                                 _dataSet->getActiveVectorsName(),
93                                                 3, i);
94            }
95        } else {
96            _dataSet->getScalarRange(_dataRange);
97            _dataSet->getVectorRange(_vectorMagnitudeRange);
98            for (int i = 0; i < 3; i++) {
99                _dataSet->getVectorRange(_vectorComponentRange[i], i);
100            }
101        }
102
103        update();
104    }
105}
106
107/**
108 * \brief Create and initialize a VTK Prop to render Streamlines
109 */
110void Streamlines::initProp()
111{
112    if (_linesActor == NULL) {
113        _linesActor = vtkSmartPointer<vtkActor>::New();
114        _linesActor->GetProperty()->SetColor(_color[0], _color[1], _color[2]);
115        _linesActor->GetProperty()->SetEdgeColor(_edgeColor[0], _edgeColor[1], _edgeColor[2]);
116        _linesActor->GetProperty()->SetLineWidth(_edgeWidth);
117        _linesActor->GetProperty()->SetOpacity(_opacity);
118        _linesActor->GetProperty()->SetAmbient(.2);
119        if (!_lighting)
120            _linesActor->GetProperty()->LightingOff();
121        switch (_lineType) {
122        case LINES:
123            setCulling(_linesActor->GetProperty(), false);
124            _linesActor->GetProperty()->SetRepresentationToWireframe();
125            _linesActor->GetProperty()->EdgeVisibilityOff();
126            break;
127        case TUBES:
128            if (_faceCulling && _opacity == 1.0)
129                setCulling(_linesActor->GetProperty(), true);
130            _linesActor->GetProperty()->SetRepresentationToSurface();
131            _linesActor->GetProperty()->EdgeVisibilityOff();
132            break;
133        case RIBBONS:
134            setCulling(_linesActor->GetProperty(), false);
135            _linesActor->GetProperty()->SetRepresentationToSurface();
136            _linesActor->GetProperty()->EdgeVisibilityOff();
137            break;
138        default:
139            ;
140        }
141    }
142    if (_seedActor == NULL) {
143        _seedActor = vtkSmartPointer<vtkActor>::New();
144        _seedActor->GetProperty()->SetColor(_seedColor[0], _seedColor[1], _seedColor[2]);
145        _seedActor->GetProperty()->SetEdgeColor(_seedColor[0], _seedColor[1], _seedColor[2]);
146        _seedActor->GetProperty()->SetLineWidth(1);
147        _seedActor->GetProperty()->SetPointSize(2);
148        _seedActor->GetProperty()->SetOpacity(_opacity);
149        _seedActor->GetProperty()->SetRepresentationToWireframe();
150        _seedActor->GetProperty()->LightingOff();
151        setSeedVisibility(_seedVisible);
152    }
153    if (_prop == NULL) {
154        _prop = vtkSmartPointer<vtkAssembly>::New();
155        getAssembly()->AddPart(_linesActor);
156        getAssembly()->AddPart(_seedActor);
157    }
158}
159
160/**
161 * \brief Get a pseudo-random number in range [min,max]
162 */
163double Streamlines::getRandomNum(double min, double max)
164{
165#if 1
166    return vtkMath::Random(min, max);
167#else
168    int r = rand();
169    return (min + ((double)r / RAND_MAX) * (max - min));
170#endif
171}
172
173/**
174 * \brief Get a random 3D point within an AABB
175 *
176 * \param[out] pt The random point
177 * \param[in] bounds The bounds of the AABB
178 */
179void Streamlines::getRandomPoint(double pt[3], const double bounds[6])
180{
181    pt[0] = getRandomNum(bounds[0], bounds[1]);
182    pt[1] = getRandomNum(bounds[2], bounds[3]);
183    pt[2] = getRandomNum(bounds[4], bounds[5]);
184}
185
186/**
187 * \brief Get a random point within a triangle (including edges)
188 *
189 * \param[out] pt The random point
190 * \param[in] v1 Triangle vertex 1
191 * \param[in] v2 Triangle vertex 2
192 * \param[in] v3 Triangle vertex 3
193 */
194void Streamlines::getRandomPointInTriangle(double pt[3],
195                                           const double v0[3],
196                                           const double v1[3],
197                                           const double v2[3])
198{
199    // Choose random barycentric coordinates
200    double bary[3];
201    bary[0] = getRandomNum(0, 1);
202    bary[1] = getRandomNum(0, 1);
203    if (bary[0] + bary[1] > 1.0) {
204        bary[0] = 1.0 - bary[0];
205        bary[1] = 1.0 - bary[1];
206    }
207    bary[2] = 1.0 - bary[0] - bary[1];
208
209    //TRACE("bary %g %g %g", bary[0], bary[1], bary[2]);
210    // Convert to cartesian coords
211    for (int i = 0; i < 3; i++) {
212        pt[i] = v0[i] * bary[0] + v1[i] * bary[1] + v2[i] * bary[2];
213    }
214}
215
216void Streamlines::getRandomPointInTetrahedron(double pt[3],
217                                              const double v0[3],
218                                              const double v1[3],
219                                              const double v2[3],
220                                              const double v3[3])
221{
222    // Choose random barycentric coordinates
223    double bary[4];
224    bary[0] = getRandomNum(0, 1);
225    bary[1] = getRandomNum(0, 1);
226    bary[2] = getRandomNum(0, 1);
227    if (bary[0] + bary[1] > 1.0) {
228        bary[0] = 1.0 - bary[0];
229        bary[1] = 1.0 - bary[1];
230    }
231    if (bary[1] + bary[2] > 1.0) {
232        double tmp = bary[2];
233        bary[2] = 1.0 - bary[0] - bary[1];
234        bary[1] = 1.0 - tmp;
235    } else if (bary[0] + bary[1] + bary[2] > 1.0) {
236        double tmp = bary[2];
237        bary[2] = bary[0] + bary[1] + bary[2] - 1.0;
238        bary[0] = 1.0 - bary[1] - tmp;
239    }
240    bary[3] = 1.0 - bary[0] - bary[1] - bary[2];
241    //TRACE("bary %g %g %g %g", bary[0], bary[1], bary[2], bary[3]);
242    // Convert to cartesian coords
243    for (int i = 0; i < 3; i++) {
244#if 0
245        pt[i] = (v0[i] - v3[i]) * bary[0] +
246                (v1[i] - v3[i]) * bary[1] +
247                (v2[i] - v3[i]) * bary[2] + v3[i];
248#else
249        pt[i] = v0[i] * bary[0] + v1[i] * bary[1] +
250                v2[i] * bary[2] + v3[i] * bary[3];
251#endif
252    }
253}
254
255/**
256 * \brief Get a random point on a line segment (including endpoints)
257 */
258void Streamlines::getRandomPointOnLineSegment(double pt[3],
259                                              const double endpt[3],
260                                              const double endpt2[3])
261{
262    double ratio = getRandomNum(0, 1);
263    pt[0] = endpt[0] + ratio * (endpt2[0] - endpt[0]);
264    pt[1] = endpt[1] + ratio * (endpt2[1] - endpt[1]);
265    pt[2] = endpt[2] + ratio * (endpt2[2] - endpt[2]);
266}
267
268/**
269 * \brief Get a random point within a vtkDataSet's mesh
270 *
271 * Note: This currently doesn't always give a uniform distribution
272 * of points in space and can generate points outside the mesh for
273 * unusual cell types
274 */
275void Streamlines::getRandomCellPt(double pt[3], vtkDataSet *ds)
276{
277    int numCells = (int)ds->GetNumberOfCells();
278    // XXX: Not uniform distribution (shouldn't use mod, and assumes
279    // all cells are equal area/volume)
280    int cell = rand() % numCells;
281    int type = ds->GetCellType(cell);
282    switch (type) {
283    case VTK_VERTEX: {
284        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
285        ds->GetCellPoints(cell, ptIds);
286        assert(ptIds->GetNumberOfIds() == 1);
287        ds->GetPoint(ptIds->GetId(0), pt);
288    }
289        break;
290    case VTK_POLY_VERTEX: {
291        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
292        ds->GetCellPoints(cell, ptIds);
293        assert(ptIds->GetNumberOfIds() >= 1);
294        int id = rand() % ptIds->GetNumberOfIds();
295        ds->GetPoint(ptIds->GetId(id), pt);
296    }
297        break;
298    case VTK_LINE: {
299        double v[2][3];
300        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
301        ds->GetCellPoints(cell, ptIds);
302        assert(ptIds->GetNumberOfIds() == 2);
303        for (int i = 0; i < 2; i++) {
304            ds->GetPoint(ptIds->GetId(i), v[i]);
305        }
306        getRandomPointOnLineSegment(pt, v[0], v[1]);
307    }
308        break;
309    case VTK_POLY_LINE: {
310        double v[2][3];
311        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
312        ds->GetCellPoints(cell, ptIds);
313        assert(ptIds->GetNumberOfIds() >= 2);
314        int id = rand() % (ptIds->GetNumberOfIds()-1);
315        for (int i = 0; i < 2; i++) {
316            ds->GetPoint(ptIds->GetId(id+i), v[i]);
317        }
318        getRandomPointOnLineSegment(pt, v[0], v[1]);
319    }
320        break;
321    case VTK_TRIANGLE: {
322        double v[3][3];
323        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
324        ds->GetCellPoints(cell, ptIds);
325        assert(ptIds->GetNumberOfIds() == 3);
326        for (int i = 0; i < 3; i++) {
327            ds->GetPoint(ptIds->GetId(i), v[i]);
328        }
329        getRandomPointInTriangle(pt, v[0], v[1], v[2]);
330    }
331        break;
332    case VTK_TRIANGLE_STRIP: {
333        double v[3][3];
334        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
335        ds->GetCellPoints(cell, ptIds);
336        assert(ptIds->GetNumberOfIds() >= 3);
337        int id = rand() % (ptIds->GetNumberOfIds()-2);
338        for (int i = 0; i < 3; i++) {
339            ds->GetPoint(ptIds->GetId(id+i), v[i]);
340        }
341        getRandomPointInTriangle(pt, v[0], v[1], v[2]);
342    }
343        break;
344    case VTK_POLYGON: {
345        vtkPolygon *poly = vtkPolygon::SafeDownCast(ds->GetCell(cell));
346        assert (poly != NULL);
347        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
348        poly->Triangulate(ptIds);
349        assert(ptIds->GetNumberOfIds() >= 3 && ptIds->GetNumberOfIds() % 3 == 0);
350        int tri = rand() % (ptIds->GetNumberOfIds()/3);
351        double v[3][3];
352        for (int i = 0; i < 3; i++) {
353            ds->GetPoint(ptIds->GetId(i + tri * 3), v[i]);
354        }
355        getRandomPointInTriangle(pt, v[0], v[1], v[2]);
356    }
357        break;
358    case VTK_QUAD: {
359        double v[4][3];
360        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
361        ds->GetCellPoints(cell, ptIds);
362        assert(ptIds->GetNumberOfIds() == 4);
363        for (int i = 0; i < 4; i++) {
364            ds->GetPoint(ptIds->GetId(i), v[i]);
365        }
366        int tri = rand() & 0x1;
367        if (tri) {
368            getRandomPointInTriangle(pt, v[0], v[1], v[2]);
369        } else {
370            getRandomPointInTriangle(pt, v[0], v[2], v[3]);
371        }
372    }
373        break;
374    case VTK_TETRA: {
375        double v[4][3];
376        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
377        ds->GetCellPoints(cell, ptIds);
378        assert(ptIds->GetNumberOfIds() == 4);
379        for (int i = 0; i < 4; i++) {
380            ds->GetPoint(ptIds->GetId(i), v[i]);
381        }
382        getRandomPointInTetrahedron(pt, v[0], v[1], v[2], v[3]);
383    }
384        break;
385    case VTK_WEDGE: {
386        double v[6][3];
387        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
388        ds->GetCellPoints(cell, ptIds);
389        assert(ptIds->GetNumberOfIds() == 6);
390        for (int i = 0; i < 6; i++) {
391            ds->GetPoint(ptIds->GetId(i), v[i]);
392        }
393        double vv[3][3];
394        getRandomPointOnLineSegment(vv[0], v[0], v[3]);
395        getRandomPointOnLineSegment(vv[1], v[1], v[4]);
396        getRandomPointOnLineSegment(vv[2], v[2], v[5]);
397        getRandomPointInTriangle(pt, vv[0], vv[1], vv[2]);
398    }
399        break;
400    case VTK_PYRAMID: {
401        double v[5][3];
402        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
403        ds->GetCellPoints(cell, ptIds);
404        assert(ptIds->GetNumberOfIds() == 5);
405        for (int i = 0; i < 5; i++) {
406            ds->GetPoint(ptIds->GetId(i), v[i]);
407        }
408        int tetra = rand() & 0x1;
409        if (tetra) {
410            getRandomPointInTetrahedron(pt, v[0], v[1], v[2], v[4]);
411        } else {
412            getRandomPointInTetrahedron(pt, v[0], v[2], v[3], v[4]);
413        }
414    }
415        break;
416    case VTK_PIXEL:
417    case VTK_VOXEL: {
418        double bounds[6];
419        ds->GetCellBounds(cell, bounds);
420        getRandomPoint(pt, bounds);
421    }
422        break;
423    default: {
424        vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
425        vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
426        ds->GetCell(cell)->Triangulate(0, ptIds, pts);
427        if (ptIds->GetNumberOfIds() % 4 == 0) {
428            int tetra = rand() % (ptIds->GetNumberOfIds()/4);
429            double v[4][3];
430            for (int i = 0; i < 4; i++) {
431                ds->GetPoint(ptIds->GetId(i + tetra * 4), v[i]);
432            }
433            getRandomPointInTetrahedron(pt, v[0], v[1], v[2], v[4]);
434        } else {
435            assert(ptIds->GetNumberOfIds() % 3 == 0);
436            int tri = rand() % (ptIds->GetNumberOfIds()/3);
437            double v[3][3];
438            for (int i = 0; i < 3; i++) {
439                ds->GetPoint(ptIds->GetId(i + tri * 3), v[i]);
440            }
441            getRandomPointInTriangle(pt, v[0], v[1], v[2]);
442        }
443    }
444    }
445}
446
447/**
448 * \brief Internal method to set up pipeline after a state change
449 */
450void Streamlines::update()
451{
452    if (_dataSet == NULL) {
453        return;
454    }
455
456    vtkDataSet *ds = _dataSet->getVtkDataSet();
457
458    double bounds[6];
459    _dataSet->getBounds(bounds);
460    double xLen = bounds[1] - bounds[0];
461    double yLen = bounds[3] - bounds[2];
462    double zLen = bounds[5] - bounds[4];
463    double maxBound = 0.0;
464    if (xLen > maxBound) {
465        maxBound = xLen;
466    }
467    if (yLen > maxBound) {
468        maxBound = yLen;
469    }
470    if (zLen > maxBound) {
471        maxBound = zLen;
472    }
473
474    double cellSizeRange[2];
475    double avgSize;
476    _dataSet->getCellSizeRange(cellSizeRange, &avgSize);
477    _dataScale = avgSize / 8.;
478
479    vtkSmartPointer<vtkCellDataToPointData> cellToPtData;
480
481    if (ds->GetPointData() == NULL ||
482        ds->GetPointData()->GetVectors() == NULL) {
483        TRACE("No vector point data found in DataSet %s", _dataSet->getName().c_str());
484        if (ds->GetCellData() != NULL &&
485            ds->GetCellData()->GetVectors() != NULL) {
486            cellToPtData =
487                vtkSmartPointer<vtkCellDataToPointData>::New();
488#ifdef USE_VTK6
489            cellToPtData->SetInputData(ds);
490#else
491            cellToPtData->SetInput(ds);
492#endif
493            //cellToPtData->PassCellDataOn();
494            cellToPtData->Update();
495            ds = cellToPtData->GetOutput();
496        } else {
497            USER_ERROR("No vector field was found in the data set.");
498            return;
499        }
500    }
501
502    if (_streamTracer == NULL) {
503        _streamTracer = vtkSmartPointer<vtkStreamTracer>::New();
504    }
505
506    if (_dataSet->isCloud()) {
507        // DataSet is a point cloud
508        PrincipalPlane plane;
509        double offset;
510        vtkDataSet *mesherOutput = NULL;
511        if (_dataSet->is2D(&plane, &offset)) {
512            // Generate a 2D PolyData
513            _mesher = vtkSmartPointer<vtkDelaunay2D>::New();
514            if (plane == PLANE_ZY) {
515                vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
516                trans->RotateWXYZ(90, 0, 1, 0);
517                if (offset != 0.0) {
518                    trans->Translate(-offset, 0, 0);
519                }
520                vtkDelaunay2D::SafeDownCast(_mesher)->SetTransform(trans);
521            } else if (plane == PLANE_XZ) {
522                vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
523                trans->RotateWXYZ(-90, 1, 0, 0);
524                if (offset != 0.0) {
525                    trans->Translate(0, -offset, 0);
526                }
527                vtkDelaunay2D::SafeDownCast(_mesher)->SetTransform(trans);
528            } else if (offset != 0.0) {
529                // XY with Z offset
530                vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
531                trans->Translate(0, 0, -offset);
532                vtkDelaunay2D::SafeDownCast(_mesher)->SetTransform(trans);
533            }
534#ifdef USE_VTK6
535            vtkDelaunay2D::SafeDownCast(_mesher)->SetInputData(ds);
536#else
537            _mesher->SetInput(ds);
538#endif
539            _streamTracer->SetInputConnection(_mesher->GetOutputPort());
540            _mesher->Update();
541            mesherOutput = vtkDelaunay2D::SafeDownCast(_mesher)->GetOutput();
542        } else {
543            // Generate a 3D unstructured grid
544            _mesher = vtkSmartPointer<vtkDelaunay3D>::New();
545#ifdef USE_VTK6
546            vtkDelaunay3D::SafeDownCast(_mesher)->SetInputData(ds);
547#else
548            _mesher->SetInput(ds);
549#endif
550            _streamTracer->SetInputConnection(_mesher->GetOutputPort());
551            _mesher->Update();
552            mesherOutput = vtkDelaunay3D::SafeDownCast(_mesher)->GetOutput();
553        }
554        DataSet::getCellSizeRange(mesherOutput, cellSizeRange, &avgSize);
555        _dataScale = avgSize / 8.;
556     } else {
557        // DataSet is NOT a cloud (has cells)
558#ifdef USE_VTK6
559        _streamTracer->SetInputData(ds);
560#else
561        _streamTracer->SetInput(ds);
562#endif
563    }
564
565    _streamTracer->SetMaximumPropagation(xLen + yLen + zLen);
566    _streamTracer->SetIntegratorTypeToRungeKutta45();
567
568    TRACE("Setting streamlines max length to %g",
569          _streamTracer->GetMaximumPropagation());
570
571    if (_pdMapper == NULL) {
572        _pdMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
573        _pdMapper->SetResolveCoincidentTopologyToPolygonOffset();
574        _pdMapper->ScalarVisibilityOn();
575    }
576    if (_seedMapper == NULL) {
577        _seedMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
578        _seedMapper->SetResolveCoincidentTopologyToPolygonOffset();
579        _seedMapper->ScalarVisibilityOff();
580    }
581
582    // Set up seed source object
583    setSeedToFilledMesh(200);
584
585    switch (_lineType) {
586    case LINES: {
587        _streamTracer->SetComputeVorticity(false);
588        _pdMapper->SetInputConnection(_streamTracer->GetOutputPort());
589    }
590        break;
591    case TUBES: {
592        _streamTracer->SetComputeVorticity(true);
593        _lineFilter = vtkSmartPointer<vtkTubeFilter>::New();
594        vtkTubeFilter *tubeFilter = vtkTubeFilter::SafeDownCast(_lineFilter);
595        tubeFilter->SetNumberOfSides(5);
596        _lineFilter->SetInputConnection(_streamTracer->GetOutputPort());
597        _pdMapper->SetInputConnection(_lineFilter->GetOutputPort());
598    }
599        break;
600    case RIBBONS: {
601        _streamTracer->SetComputeVorticity(true);
602        _lineFilter = vtkSmartPointer<vtkRibbonFilter>::New();
603        _lineFilter->SetInputConnection(_streamTracer->GetOutputPort());
604        _pdMapper->SetInputConnection(_lineFilter->GetOutputPort());
605    }
606        break;
607    default:
608        ERROR("Unknown LineType: %d", _lineType);
609    }
610
611#if 0 && defined(WANT_TRACE)
612    _streamTracer->Update();
613    vtkPolyData *pd = _streamTracer->GetOutput();
614    DataSet::print(pd);
615    TRACE("Points: %d Verts: %d Lines: %d Polys: %d Strips: %d",
616          pd->GetNumberOfPoints(),
617          pd->GetNumberOfVerts(),
618          pd->GetNumberOfLines(),
619          pd->GetNumberOfPolys(),
620          pd->GetNumberOfStrips());
621#if 0
622    vtkCellArray *arr = pd->GetLines();
623    arr->InitTraversal();
624    vtkIdType npts, *pts;
625    arr->GetNextCell(npts, pts);
626    for (int i = 0; i < npts; i++) {
627        TRACE("Pt: %d", pts[i]);
628    }
629#endif
630#endif
631
632    initProp();
633
634    _seedActor->SetMapper(_seedMapper);
635
636    if (_lut == NULL) {
637        setColorMap(ColorMap::getDefault());
638    }
639
640    setColorMode(_colorMode);
641
642    _linesActor->SetMapper(_pdMapper);
643    _pdMapper->Update();
644    _seedMapper->Update();
645#ifdef WANT_TRACE
646    double *b = getBounds();
647    TRACE("bounds: %g %g %g %g %g %g", b[0], b[1], b[2], b[3], b[4], b[5]);
648#endif
649}
650
651void Streamlines::setNumberOfSeedPoints(int numPoints)
652{
653    switch(_seedType) {
654    case DATASET_FILLED_MESH:
655        setSeedToFilledMesh(numPoints);
656        break;
657    case FILLED_MESH:
658        if (_seedMesh != NULL) {
659            setSeedToFilledMesh(_seedMesh, numPoints);
660        } else {
661            ERROR("NULL _seedMesh");
662        }
663        break;
664    case DATASET_MESH_POINTS:
665        setSeedToMeshPoints(numPoints);
666        break;
667    case MESH_POINTS:
668        setSeedToMeshPoints(_seedMesh, numPoints);
669        break;
670    default:
671        ERROR("Can't set number of points for seed type %d", _seedType);
672        break;
673    }
674    _pdMapper->Update();
675    _seedMapper->Update();
676}
677
678/**
679 * \brief Use points of the DataSet associated with this
680 * Streamlines as seeds
681 */
682void Streamlines::setSeedToMeshPoints(int maxPoints)
683{
684    _seedType = DATASET_MESH_POINTS;
685    setSeedToMeshPoints(_dataSet->getVtkDataSet(), maxPoints);
686}
687
688/**
689 * \brief Use seed points randomly distributed within the cells
690 * of the DataSet associated with this Streamlines
691 *
692 * Note: The current implementation doesn't give a uniform
693 * distribution of points, and points outside the mesh bounds
694 * may be generated
695 *
696 * \param[in] numPoints Number of random seed points to generate
697 */
698void Streamlines::setSeedToFilledMesh(int numPoints)
699{
700    _seedType = DATASET_FILLED_MESH;
701    setSeedToFilledMesh(_dataSet->getVtkDataSet(), numPoints);
702}
703
704/**
705 * \brief Use points of a supplied vtkDataSet as seeds
706 *
707 * \param[in] seed vtkDataSet with points to use as seeds
708 * \param[in] maxPoints Maximum number of points to be used as seeds
709 */
710void Streamlines::setSeedToMeshPoints(vtkDataSet *seed, int maxPoints)
711{
712    if (seed != _dataSet->getVtkDataSet()) {
713        _seedType = MESH_POINTS;
714        _seedMesh = seed;
715    } else {
716        _seedType = DATASET_MESH_POINTS;
717        _seedMesh = NULL;
718    }
719
720    if (_streamTracer == NULL)
721        return;
722
723#ifndef USE_VTK6
724    vtkSmartPointer<vtkDataSet> oldSeed;
725    if (_streamTracer->GetSource() != NULL) {
726        oldSeed = _streamTracer->GetSource();
727    }
728#endif
729
730    if (maxPoints > 0 && seed->GetNumberOfPoints() > maxPoints) {
731        TRACE("Seed points: %d", maxPoints);
732        vtkSmartPointer<vtkMaskPoints> mask = vtkSmartPointer<vtkMaskPoints>::New();
733#ifdef USE_VTK6
734        mask->SetInputData(seed);
735#else
736        mask->SetInput(seed);
737#endif
738        mask->SetMaximumNumberOfPoints(maxPoints);
739        mask->RandomModeOn();
740        mask->GenerateVerticesOff();
741        _streamTracer->SetSourceConnection(mask->GetOutputPort());
742
743        vtkSmartPointer<vtkVertexGlyphFilter> vertFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
744        vertFilter->SetInputConnection(mask->GetOutputPort());
745        _seedMapper->SetInputConnection(vertFilter->GetOutputPort());
746    } else {
747        TRACE("Seed points: %d", seed->GetNumberOfPoints());
748
749#ifdef USE_VTK6
750        _streamTracer->SetSourceData(seed);
751#else
752        _streamTracer->SetSource(seed);
753#endif
754        if (vtkPolyData::SafeDownCast(seed) != NULL) {
755#ifdef USE_VTK6
756            _seedMapper->SetInputData(vtkPolyData::SafeDownCast(seed));
757#else
758            _seedMapper->SetInput(vtkPolyData::SafeDownCast(seed));
759#endif
760        } else {
761            vtkSmartPointer<vtkVertexGlyphFilter> vertFilter = vtkSmartPointer<vtkVertexGlyphFilter>::New();
762#ifdef USE_VTK6
763            vertFilter->SetInputData(seed);
764#else
765            vertFilter->SetInput(seed);
766#endif
767            _seedMapper->SetInputConnection(vertFilter->GetOutputPort());
768        }
769    }
770
771#ifndef USE_VTK6
772    if (oldSeed != NULL) {
773        oldSeed->SetPipelineInformation(NULL);
774    }
775#endif
776}
777
778/**
779 * \brief Use seed points randomly distributed within the cells
780 * of a supplied vtkDataSet
781 *
782 * Note: The current implementation doesn't give a uniform
783 * distribution of points, and points outside the mesh bounds
784 * may be generated
785 *
786 * \param[in] ds vtkDataSet containing cells
787 * \param[in] numPoints Number of random seed points to generate
788 */
789void Streamlines::setSeedToFilledMesh(vtkDataSet *ds, int numPoints)
790{
791    if (ds != _dataSet->getVtkDataSet()) {
792        _seedMesh = ds;
793        _seedType = FILLED_MESH;
794#ifdef DEBUG
795        DataSet::print(ds);
796#endif
797    } else {
798        _seedMesh = NULL;
799        if (_mesher != NULL) {
800            if (vtkDelaunay2D::SafeDownCast(_mesher) != NULL) {
801                ds = vtkDelaunay2D::SafeDownCast(_mesher)->GetOutput();
802            } else {
803                ds = vtkDelaunay3D::SafeDownCast(_mesher)->GetOutput();
804            }
805        }
806    }
807    if (_streamTracer != NULL) {
808        // Set up seed source object
809        vtkSmartPointer<vtkPolyData> seed = vtkSmartPointer<vtkPolyData>::New();
810        vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
811        vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
812
813        if (ds->GetNumberOfCells() < 1) {
814            ERROR("No cells in mesh");
815        }
816
817        for (int i = 0; i < numPoints; i++) {
818            double pt[3];
819            getRandomCellPt(pt, ds);
820            //TRACE("Seed pt: %g %g %g", pt[0], pt[1], pt[2]);
821            pts->InsertNextPoint(pt);
822            cells->InsertNextCell(1);
823            cells->InsertCellPoint(i);
824        }
825
826        seed->SetPoints(pts);
827        seed->SetVerts(cells);
828
829        TRACE("Seed points: %d", seed->GetNumberOfPoints());
830        vtkSmartPointer<vtkDataSet> oldSeed;
831        if (_streamTracer->GetSource() != NULL) {
832            oldSeed = _streamTracer->GetSource();
833        }
834
835#ifdef USE_VTK6
836        _streamTracer->SetSourceData(seed);
837        _seedMapper->SetInputData(seed);
838#else
839        _streamTracer->SetSource(seed);
840
841        if (oldSeed != NULL) {
842            oldSeed->SetPipelineInformation(NULL);
843        }
844        _seedMapper->SetInput(seed);
845#endif
846    }
847}
848
849/**
850 * \brief Use seed points along a line
851 *
852 * \param[in] start Starting point of rake line
853 * \param[in] end End point of rake line
854 * \param[in] numPoints Number of points along line to generate
855 */
856void Streamlines::setSeedToRake(double start[3], double end[3], int numPoints)
857{
858    if (numPoints < 2)
859        return;
860    _seedType = RAKE;
861    _seedMesh = NULL;
862    if (_streamTracer != NULL) {
863        // Set up seed source object
864        vtkSmartPointer<vtkPolyData> seed = vtkSmartPointer<vtkPolyData>::New();
865        vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
866        vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
867        vtkSmartPointer<vtkPolyLine> polyline = vtkSmartPointer<vtkPolyLine>::New();
868
869        double dir[3];
870        for (int i = 0; i < 3; i++) {
871            dir[i] = end[i] - start[i];
872        }
873
874        polyline->GetPointIds()->SetNumberOfIds(numPoints);
875        for (int i = 0; i < numPoints; i++) {
876            double pt[3];
877            for (int ii = 0; ii < 3; ii++) {
878                pt[ii] = start[ii] + dir[ii] * ((double)i / (numPoints-1));
879            }
880            //TRACE("Seed pt: %g %g %g", pt[0], pt[1], pt[2]);
881            pts->InsertNextPoint(pt);
882            polyline->GetPointIds()->SetId(i, i);
883        }
884
885        cells->InsertNextCell(polyline);
886        seed->SetPoints(pts);
887        seed->SetLines(cells);
888
889        TRACE("Seed points: %d", seed->GetNumberOfPoints());
890        vtkSmartPointer<vtkDataSet> oldSeed;
891        if (_streamTracer->GetSource() != NULL) {
892            oldSeed = _streamTracer->GetSource();
893        }
894#ifdef USE_VTK6
895        _streamTracer->SetSourceData(seed);
896        _seedMapper->SetInputData(seed);
897#else
898        _streamTracer->SetSource(seed);
899        if (oldSeed != NULL) {
900            oldSeed->SetPipelineInformation(NULL);
901        }
902
903        _seedMapper->SetInput(seed);
904#endif
905    }
906}
907
908/**
909 * \brief Create seed points inside a disk with an optional hole
910 *
911 * \param[in] center Center point of disk
912 * \param[in] normal Normal vector to orient disk
913 * \param[in] radius Radius of disk
914 * \param[in] innerRadius Radius of hole at center of disk
915 * \param[in] numPoints Number of random points to generate
916 */
917void Streamlines::setSeedToDisk(double center[3],
918                                double normal[3],
919                                double radius,
920                                double innerRadius,
921                                int numPoints)
922{
923    _seedType = DISK;
924    _seedMesh = NULL;
925    if (_streamTracer != NULL) {
926        // Set up seed source object
927        vtkSmartPointer<vtkPolyData> seed = vtkSmartPointer<vtkPolyData>::New();
928        vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
929        vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
930
931        // The following code is based on vtkRegularPolygonSource::RequestData
932
933        double px[3];
934        double py[3];
935        double axis[3] = {1., 0., 0.};
936
937        if (vtkMath::Normalize(normal) == 0.0) {
938            normal[0] = 0.0;
939            normal[1] = 0.0;
940            normal[2] = 1.0;
941        }
942
943        // Find axis in plane (orthogonal to normal)
944        bool done = false;
945        vtkMath::Cross(normal, axis, px);
946        if (vtkMath::Normalize(px) > 1.0e-3) {
947            done = true;
948        }
949        if (!done) {
950            axis[0] = 0.0;
951            axis[1] = 1.0;
952            axis[2] = 0.0;
953            vtkMath::Cross(normal, axis, px);
954            if (vtkMath::Normalize(px) > 1.0e-3) {
955                done = true;
956            }
957        }
958        if (!done) {
959            axis[0] = 0.0;
960            axis[1] = 0.0;
961            axis[2] = 1.0;
962            vtkMath::Cross(normal, axis, px);
963            vtkMath::Normalize(px);
964        }
965        // Create third orthogonal basis vector
966        vtkMath::Cross(px, normal, py);
967
968        double minSquared = (innerRadius*innerRadius)/(radius*radius);
969        for (int j = 0; j < numPoints; j++) {
970            // Get random sweep angle and radius
971#ifdef USE_VTK6
972            double angle = getRandomNum(0, 2.0 * vtkMath::Pi());
973#else
974            double angle = getRandomNum(0, 2.0 * vtkMath::DoublePi());
975#endif
976            // Need sqrt to get uniform distribution
977            double r = sqrt(getRandomNum(minSquared, 1)) * radius;
978            double pt[3];
979            for (int i = 0; i < 3; i++) {
980                pt[i] = center[i] + r * (px[i] * cos(angle) + py[i] * sin(angle));
981            }
982            //TRACE("Seed pt: %g %g %g", pt[0], pt[1], pt[2]);
983            pts->InsertNextPoint(pt);
984            cells->InsertNextCell(1);
985            cells->InsertCellPoint(j);
986        }
987
988        seed->SetPoints(pts);
989        seed->SetVerts(cells);
990
991        TRACE("Seed points: %d", seed->GetNumberOfPoints());
992        vtkSmartPointer<vtkDataSet> oldSeed;
993        if (_streamTracer->GetSource() != NULL) {
994            oldSeed = _streamTracer->GetSource();
995        }
996#ifdef USE_VTK6
997        _streamTracer->SetSourceData(seed);
998        _seedMapper->SetInputData(seed);
999#else
1000        _streamTracer->SetSource(seed);
1001        if (oldSeed != NULL) {
1002            oldSeed->SetPipelineInformation(NULL);
1003        }
1004
1005        _seedMapper->SetInput(seed);
1006#endif
1007    }
1008}
1009
1010/**
1011 * \brief Use seed points from an n-sided polygon
1012 *
1013 * \param[in] center Center point of polygon
1014 * \param[in] normal Normal vector to orient polygon
1015 * \param[in] angle Angle in degrees to rotate about normal
1016 * \param[in] radius Radius of circumscribing circle
1017 * \param[in] numSides Number of polygon sides (and points) to generate
1018 */
1019void Streamlines::setSeedToPolygon(double center[3],
1020                                   double normal[3],
1021                                   double angle,
1022                                   double radius,
1023                                   int numSides)
1024{
1025    _seedType = POLYGON;
1026    _seedMesh = NULL;
1027    if (_streamTracer != NULL) {
1028        // Set up seed source object
1029        vtkSmartPointer<vtkRegularPolygonSource> seed = vtkSmartPointer<vtkRegularPolygonSource>::New();
1030
1031        seed->SetCenter(center);
1032        seed->SetNormal(normal);
1033        seed->SetRadius(radius);
1034        seed->SetNumberOfSides(numSides);
1035        seed->GeneratePolygonOn();
1036
1037        if (angle != 0.0) {
1038            vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
1039            trans->RotateWXYZ(angle, normal);
1040            vtkSmartPointer<vtkTransformPolyDataFilter> transFilt =
1041                vtkSmartPointer<vtkTransformPolyDataFilter>::New();
1042            transFilt->SetInputConnection(seed->GetOutputPort());
1043            transFilt->SetTransform(trans);
1044        }
1045
1046        TRACE("Seed points: %d", numSides);
1047        vtkSmartPointer<vtkDataSet> oldSeed;
1048        if (_streamTracer->GetSource() != NULL) {
1049            oldSeed = _streamTracer->GetSource();
1050        }
1051
1052        if (angle != 0.0) {
1053            vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
1054            trans->Translate(+center[0], +center[1], +center[2]);
1055            trans->RotateWXYZ(angle, normal);
1056            trans->Translate(-center[0], -center[1], -center[2]);
1057            vtkSmartPointer<vtkTransformPolyDataFilter> transFilt =
1058                vtkSmartPointer<vtkTransformPolyDataFilter>::New();
1059            transFilt->SetInputConnection(seed->GetOutputPort());
1060            transFilt->SetTransform(trans);
1061            _streamTracer->SetSourceConnection(transFilt->GetOutputPort());
1062            _seedMapper->SetInputConnection(transFilt->GetOutputPort());
1063        } else {
1064            _streamTracer->SetSourceConnection(seed->GetOutputPort());
1065            _seedMapper->SetInputConnection(seed->GetOutputPort());
1066        }
1067
1068#ifndef USE_VTK6
1069        if (oldSeed != NULL) {
1070            oldSeed->SetPipelineInformation(NULL);
1071        }
1072#endif
1073    }
1074}
1075
1076/**
1077 * \brief Use seed points from an n-sided polygon
1078 *
1079 * \param[in] center Center point of polygon
1080 * \param[in] normal Normal vector to orient polygon
1081 * \param[in] angle Angle in degrees to rotate about normal
1082 * \param[in] radius Radius of circumscribing circle
1083 * \param[in] numSides Number of polygon sides (and points) to generate
1084 * \param[in] numPoints Number of random points to generate
1085 */
1086void Streamlines::setSeedToFilledPolygon(double center[3],
1087                                         double normal[3],
1088                                         double angle,
1089                                         double radius,
1090                                         int numSides,
1091                                         int numPoints)
1092{
1093    _seedType = FILLED_POLYGON;
1094    _seedMesh = NULL;
1095    if (_streamTracer != NULL) {
1096         // Set up seed source object
1097        vtkSmartPointer<vtkPolyData> seed = vtkSmartPointer<vtkPolyData>::New();
1098        vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
1099        vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
1100
1101        // The following code is based on vtkRegularPolygonSource::RequestData
1102
1103        double px[3];
1104        double py[3];
1105        double axis[3] = {1., 0., 0.};
1106
1107        if (vtkMath::Normalize(normal) == 0.0) {
1108            normal[0] = 0.0;
1109            normal[1] = 0.0;
1110            normal[2] = 1.0;
1111        }
1112
1113        // Find axis in plane (orthogonal to normal)
1114        bool done = false;
1115        vtkMath::Cross(normal, axis, px);
1116        if (vtkMath::Normalize(px) > 1.0e-3) {
1117            done = true;
1118        }
1119        if (!done) {
1120            axis[0] = 0.0;
1121            axis[1] = 1.0;
1122            axis[2] = 0.0;
1123            vtkMath::Cross(normal, axis, px);
1124            if (vtkMath::Normalize(px) > 1.0e-3) {
1125                done = true;
1126            }
1127        }
1128        if (!done) {
1129            axis[0] = 0.0;
1130            axis[1] = 0.0;
1131            axis[2] = 1.0;
1132            vtkMath::Cross(normal, axis, px);
1133            vtkMath::Normalize(px);
1134        }
1135        // Create third orthogonal basis vector
1136        vtkMath::Cross(px, normal, py);
1137
1138        double verts[numSides][3];
1139#ifdef USE_VTK6
1140        double sliceTheta = 2.0 * vtkMath::Pi() / (double)numSides;
1141#else
1142        double sliceTheta = 2.0 * vtkMath::DoublePi() / (double)numSides;
1143#endif
1144        angle = vtkMath::RadiansFromDegrees(angle);
1145        for (int j = 0; j < numSides; j++) {
1146            for (int i = 0; i < 3; i++) {
1147                double theta = sliceTheta * (double)j - angle;
1148                verts[j][i] = center[i] + radius * (px[i] * cos(theta) +
1149                                                    py[i] * sin(theta));
1150            }
1151            //TRACE("Vert %d: %g %g %g", j, verts[j][0], verts[j][1], verts[j][2]);
1152        }
1153
1154        // Note: this gives a uniform distribution because the polygon is regular and
1155        // the triangular sections have equal area
1156        if (numSides == 3) {
1157            for (int j = 0; j < numPoints; j++) {
1158                double pt[3];
1159                getRandomPointInTriangle(pt, verts[0], verts[1], verts[2]);
1160                //TRACE("Seed pt: %g %g %g", pt[0], pt[1], pt[2]);
1161                pts->InsertNextPoint(pt);
1162                cells->InsertNextCell(1);
1163                cells->InsertCellPoint(j);
1164            }
1165        } else {
1166            for (int j = 0; j < numPoints; j++) {
1167                // Get random triangle section
1168                int tri = rand() % numSides;
1169                double pt[3];
1170                getRandomPointInTriangle(pt, center, verts[tri], verts[(tri+1) % numSides]);
1171                //TRACE("Seed pt: %g %g %g", pt[0], pt[1], pt[2]);
1172                pts->InsertNextPoint(pt);
1173                cells->InsertNextCell(1);
1174                cells->InsertCellPoint(j);
1175            }
1176        }
1177
1178        seed->SetPoints(pts);
1179        seed->SetVerts(cells);
1180
1181        TRACE("Seed points: %d", seed->GetNumberOfPoints());
1182        vtkSmartPointer<vtkDataSet> oldSeed;
1183        if (_streamTracer->GetSource() != NULL) {
1184            oldSeed = _streamTracer->GetSource();
1185        }
1186
1187#ifdef USE_VTK6
1188        _streamTracer->SetSourceData(seed);
1189        _seedMapper->SetInputData(seed);
1190#else
1191        _streamTracer->SetSource(seed);
1192        if (oldSeed != NULL) {
1193            oldSeed->SetPipelineInformation(NULL);
1194        }
1195
1196        _seedMapper->SetInput(seed);
1197#endif
1198    }
1199}
1200
1201/**
1202 * \brief Set the integration method used
1203 */
1204void Streamlines::setIntegrator(IntegratorType integrator)
1205{
1206    if (_streamTracer != NULL) {
1207        switch (integrator) {
1208        case RUNGE_KUTTA2:
1209            _streamTracer->SetIntegratorTypeToRungeKutta2();
1210            break;
1211        case RUNGE_KUTTA4:
1212            _streamTracer->SetIntegratorTypeToRungeKutta4();
1213            break;
1214        case RUNGE_KUTTA45:
1215            _streamTracer->SetIntegratorTypeToRungeKutta45();
1216            break;
1217        default:
1218            ;
1219        }
1220    }
1221}
1222
1223/**
1224 * \brief Set the direction of integration
1225 */
1226void Streamlines::setIntegrationDirection(IntegrationDirection dir)
1227{
1228    if (_streamTracer != NULL) {
1229        switch (dir) {
1230        case FORWARD:
1231            _streamTracer->SetIntegrationDirectionToForward();
1232            break;
1233        case BACKWARD:
1234            _streamTracer->SetIntegrationDirectionToBackward();
1235            break;
1236        case BOTH:
1237            _streamTracer->SetIntegrationDirectionToBoth();
1238            break;
1239        default:
1240            ;
1241        }
1242    }
1243}
1244
1245/**
1246 * \brief Set the step size units.  Length units are world
1247 * coordinates, and cell units means steps are from cell to
1248 * cell.  Default is cell units.
1249 *
1250 * Note: calling this function will not convert existing
1251 * initial, minimum or maximum step value settings to the
1252 * new units, so this function should be called before
1253 * setting step values.
1254 */
1255void Streamlines::setIntegrationStepUnit(StepUnit unit)
1256{
1257    if (_streamTracer != NULL) {
1258        switch (unit) {
1259        case LENGTH_UNIT:
1260            _streamTracer->SetIntegrationStepUnit(vtkStreamTracer::LENGTH_UNIT);
1261            break;
1262        case CELL_LENGTH_UNIT:
1263            _streamTracer->SetIntegrationStepUnit(vtkStreamTracer::CELL_LENGTH_UNIT);
1264            break;
1265        default:
1266            ;
1267        }
1268    }
1269}
1270
1271/**
1272 * \brief Set initial step size for adaptive step integrator in
1273 * step units (see setIntegrationStepUnit).  For non-adaptive
1274 * integrators, this is the fixed step size.
1275 */
1276void Streamlines::setInitialIntegrationStep(double step)
1277{
1278    if (_streamTracer != NULL) {
1279        _streamTracer->SetInitialIntegrationStep(step);
1280    }
1281}
1282
1283/**
1284 * \brief Set minimum step for adaptive step integrator in
1285 * step units (see setIntegrationStepUnit)
1286 */
1287void Streamlines::setMinimumIntegrationStep(double step)
1288{
1289    if (_streamTracer != NULL) {
1290        _streamTracer->SetMinimumIntegrationStep(step);
1291    }
1292}
1293
1294/**
1295 * \brief Set maximum step for adaptive step integrator in
1296 * step units (see setIntegrationStepUnit)
1297 */
1298void Streamlines::setMaximumIntegrationStep(double step)
1299{
1300    if (_streamTracer != NULL) {
1301        _streamTracer->SetMaximumIntegrationStep(step);
1302    }
1303}
1304
1305/**
1306 * \brief Set maximum error tolerance
1307 */
1308void Streamlines::setMaximumError(double error)
1309{
1310    if (_streamTracer != NULL) {
1311        _streamTracer->SetMaximumError(error);
1312    }
1313}
1314
1315/**
1316 * \brief Set maximum length of stream lines in world
1317 * coordinates
1318 */
1319void Streamlines::setMaxPropagation(double length)
1320{
1321    if (_streamTracer != NULL) {
1322        _streamTracer->SetMaximumPropagation(length);
1323    }
1324}
1325
1326/**
1327 * \brief Set maximum number of integration steps
1328 */
1329void Streamlines::setMaxNumberOfSteps(int steps)
1330{
1331    if (_streamTracer != NULL) {
1332        _streamTracer->SetMaximumNumberOfSteps(steps);
1333    }
1334}
1335
1336/**
1337 * \brief Set the minimum speed before integration stops
1338 */
1339void Streamlines::setTerminalSpeed(double speed)
1340{
1341    if (_streamTracer != NULL) {
1342        _streamTracer->SetTerminalSpeed(speed);
1343    }
1344}
1345
1346/**
1347 * \brief Set streamline type to polylines
1348 */
1349void Streamlines::setLineTypeToLines()
1350{
1351    _lineType = LINES;
1352    if (_streamTracer != NULL &&
1353        _pdMapper != NULL) {
1354        _streamTracer->SetComputeVorticity(false);
1355        _pdMapper->SetInputConnection(_streamTracer->GetOutputPort());
1356        _lineFilter = NULL;
1357        setCulling(_linesActor->GetProperty(), false);
1358        _linesActor->GetProperty()->SetRepresentationToWireframe();
1359        _linesActor->GetProperty()->LightingOff();
1360    }
1361}
1362
1363/**
1364 * \brief Set streamline type to 3D tubes
1365 *
1366 * \param[in] numSides Number of sides (>=3) for tubes
1367 * \param[in] radius World coordinate minimum tube radius
1368 */
1369void Streamlines::setLineTypeToTubes(int numSides, double radius)
1370{
1371    _lineType = TUBES;
1372    if (_streamTracer != NULL) {
1373        _streamTracer->SetComputeVorticity(true);
1374        if (vtkTubeFilter::SafeDownCast(_lineFilter) == NULL) {
1375            _lineFilter = vtkSmartPointer<vtkTubeFilter>::New();
1376            _lineFilter->SetInputConnection(_streamTracer->GetOutputPort());
1377        }
1378        vtkTubeFilter *tubeFilter = vtkTubeFilter::SafeDownCast(_lineFilter);
1379        if (numSides < 3)
1380            numSides = 3;
1381        tubeFilter->SetNumberOfSides(numSides);
1382        tubeFilter->SetRadius(_dataScale * radius);
1383        _pdMapper->SetInputConnection(_lineFilter->GetOutputPort());
1384        if (_faceCulling && _opacity == 1.0)
1385            setCulling(_linesActor->GetProperty(), true);
1386        _linesActor->GetProperty()->SetRepresentationToSurface();
1387        _linesActor->GetProperty()->LightingOn();
1388     }
1389}
1390
1391/**
1392 * \brief Set streamline type to 3D ribbons
1393 *
1394 * \param[in] width Minimum half-width of ribbons
1395 * \param[in] angle Default ribbon angle in degrees from normal
1396 */
1397void Streamlines::setLineTypeToRibbons(double width, double angle)
1398{
1399    _lineType = RIBBONS;
1400    if (_streamTracer != NULL) {
1401        _streamTracer->SetComputeVorticity(true);
1402        if (vtkRibbonFilter::SafeDownCast(_lineFilter) == NULL) {
1403            _lineFilter = vtkSmartPointer<vtkRibbonFilter>::New();
1404            _lineFilter->SetInputConnection(_streamTracer->GetOutputPort());
1405        }
1406        vtkRibbonFilter *ribbonFilter = vtkRibbonFilter::SafeDownCast(_lineFilter);
1407        ribbonFilter->SetWidth(_dataScale * width);
1408        ribbonFilter->SetAngle(angle);
1409        ribbonFilter->UseDefaultNormalOn();
1410        _pdMapper->SetInputConnection(_lineFilter->GetOutputPort());
1411        setCulling(_linesActor->GetProperty(), false);
1412        _linesActor->GetProperty()->SetRepresentationToSurface();
1413        _linesActor->GetProperty()->LightingOn();
1414    }
1415}
1416
1417void Streamlines::updateRanges(Renderer *renderer)
1418{
1419    if (_dataSet == NULL) {
1420        ERROR("called before setDataSet");
1421        return;
1422    }
1423
1424    if (renderer->getUseCumulativeRange()) {
1425        renderer->getCumulativeDataRange(_dataRange,
1426                                         _dataSet->getActiveScalarsName(),
1427                                         1);
1428        renderer->getCumulativeDataRange(_vectorMagnitudeRange,
1429                                         _dataSet->getActiveVectorsName(),
1430                                         3);
1431        for (int i = 0; i < 3; i++) {
1432            renderer->getCumulativeDataRange(_vectorComponentRange[i],
1433                                             _dataSet->getActiveVectorsName(),
1434                                             3, i);
1435        }
1436    } else {
1437        _dataSet->getScalarRange(_dataRange);
1438        _dataSet->getVectorRange(_vectorMagnitudeRange);
1439        for (int i = 0; i < 3; i++) {
1440            _dataSet->getVectorRange(_vectorComponentRange[i], i);
1441        }
1442    }
1443
1444    // Need to update color map ranges and/or active vector field
1445    double *rangePtr = _colorFieldRange;
1446    if (_colorFieldRange[0] > _colorFieldRange[1]) {
1447        rangePtr = NULL;
1448    }
1449    setColorMode(_colorMode, _colorFieldType, _colorFieldName.c_str(), rangePtr);
1450}
1451
1452void Streamlines::setColorMode(ColorMode mode)
1453{
1454    _colorMode = mode;
1455    if (_dataSet == NULL)
1456        return;
1457
1458    switch (mode) {
1459    case COLOR_BY_SCALAR:
1460        setColorMode(mode,
1461                     _dataSet->getActiveScalarsType(),
1462                     _dataSet->getActiveScalarsName());
1463        break;
1464    case COLOR_BY_VECTOR_MAGNITUDE:
1465        setColorMode(mode,
1466                     _dataSet->getActiveVectorsType(),
1467                     _dataSet->getActiveVectorsName());
1468        break;
1469    case COLOR_BY_VECTOR_X:
1470        setColorMode(mode,
1471                     _dataSet->getActiveVectorsType(),
1472                     _dataSet->getActiveVectorsName());
1473        break;
1474    case COLOR_BY_VECTOR_Y:
1475        setColorMode(mode,
1476                     _dataSet->getActiveVectorsType(),
1477                     _dataSet->getActiveVectorsName());
1478        break;
1479    case COLOR_BY_VECTOR_Z:
1480        setColorMode(mode,
1481                     _dataSet->getActiveVectorsType(),
1482                     _dataSet->getActiveVectorsName());
1483        break;
1484    case COLOR_CONSTANT:
1485    default:
1486        setColorMode(mode, DataSet::POINT_DATA, NULL, NULL);
1487        break;
1488    }
1489}
1490
1491void Streamlines::setColorMode(ColorMode mode,
1492                               const char *name, double range[2])
1493{
1494    if (_dataSet == NULL)
1495        return;
1496    DataSet::DataAttributeType type = DataSet::POINT_DATA;
1497    int numComponents = 1;
1498    if (name != NULL && strlen(name) > 0 &&
1499        !_dataSet->getFieldInfo(name, &type, &numComponents)) {
1500        ERROR("Field not found: %s", name);
1501        return;
1502    }
1503    setColorMode(mode, type, name, range);
1504}
1505
1506void Streamlines::setColorMode(ColorMode mode, DataSet::DataAttributeType type,
1507                               const char *name, double range[2])
1508{
1509    _colorMode = mode;
1510    _colorFieldType = type;
1511    if (name == NULL)
1512        _colorFieldName.clear();
1513    else
1514        _colorFieldName = name;
1515    if (range == NULL) {
1516        _colorFieldRange[0] = DBL_MAX;
1517        _colorFieldRange[1] = -DBL_MAX;
1518    } else {
1519        memcpy(_colorFieldRange, range, sizeof(double)*2);
1520    }
1521
1522    if (_dataSet == NULL || _pdMapper == NULL)
1523        return;
1524
1525    switch (type) {
1526    case DataSet::POINT_DATA:
1527        _pdMapper->SetScalarModeToUsePointFieldData();
1528        break;
1529    case DataSet::CELL_DATA:
1530        _pdMapper->SetScalarModeToUseCellFieldData();
1531        break;
1532    default:
1533        ERROR("Unsupported DataAttributeType: %d", type);
1534        return;
1535    }
1536
1537    if (name != NULL && strlen(name) > 0) {
1538        _pdMapper->SelectColorArray(name);
1539    } else {
1540        _pdMapper->SetScalarModeToDefault();
1541    }
1542
1543    if (_lut != NULL) {
1544        if (range != NULL) {
1545            _lut->SetRange(range);
1546        } else if (name != NULL && strlen(name) > 0) {
1547            double r[2];
1548            int comp = -1;
1549            if (mode == COLOR_BY_VECTOR_X)
1550                comp = 0;
1551            else if (mode == COLOR_BY_VECTOR_Y)
1552                comp = 1;
1553            else if (mode == COLOR_BY_VECTOR_Z)
1554                comp = 2;
1555
1556            if (_renderer->getUseCumulativeRange()) {
1557                int numComponents;
1558                if  (!_dataSet->getFieldInfo(name, type, &numComponents)) {
1559                    ERROR("Field not found: %s, type: %d", name, type);
1560                    return;
1561                } else if (numComponents < comp+1) {
1562                    ERROR("Request for component %d in field with %d components",
1563                          comp, numComponents);
1564                    return;
1565                }
1566                _renderer->getCumulativeDataRange(r, name, type, numComponents, comp);
1567            } else {
1568                _dataSet->getDataRange(r, name, type, comp);
1569            }
1570            _lut->SetRange(r);
1571        }
1572    }
1573
1574    switch (mode) {
1575    case COLOR_BY_SCALAR:
1576        _pdMapper->ScalarVisibilityOn();
1577        break;
1578    case COLOR_BY_VECTOR_MAGNITUDE:
1579        _pdMapper->ScalarVisibilityOn();
1580        if (_lut != NULL) {
1581            _lut->SetVectorModeToMagnitude();
1582        }
1583        break;
1584    case COLOR_BY_VECTOR_X:
1585        _pdMapper->ScalarVisibilityOn();
1586        if (_lut != NULL) {
1587            _lut->SetVectorModeToComponent();
1588            _lut->SetVectorComponent(0);
1589        }
1590        break;
1591    case COLOR_BY_VECTOR_Y:
1592        _pdMapper->ScalarVisibilityOn();
1593        if (_lut != NULL) {
1594            _lut->SetVectorModeToComponent();
1595            _lut->SetVectorComponent(1);
1596        }
1597        break;
1598    case COLOR_BY_VECTOR_Z:
1599        _pdMapper->ScalarVisibilityOn();
1600        if (_lut != NULL) {
1601            _lut->SetVectorModeToComponent();
1602            _lut->SetVectorComponent(2);
1603        }
1604        break;
1605    case COLOR_CONSTANT:
1606    default:
1607        _pdMapper->ScalarVisibilityOff();
1608        break;
1609    }
1610}
1611
1612/**
1613 * \brief Called when the color map has been edited
1614 */
1615void Streamlines::updateColorMap()
1616{
1617    setColorMap(_colorMap);
1618}
1619
1620/**
1621 * \brief Associate a colormap lookup table with the DataSet
1622 */
1623void Streamlines::setColorMap(ColorMap *cmap)
1624{
1625    if (cmap == NULL)
1626        return;
1627
1628    _colorMap = cmap;
1629 
1630    if (_lut == NULL) {
1631        _lut = vtkSmartPointer<vtkLookupTable>::New();
1632        if (_pdMapper != NULL) {
1633            _pdMapper->UseLookupTableScalarRangeOn();
1634            _pdMapper->SetLookupTable(_lut);
1635        }
1636        _lut->DeepCopy(cmap->getLookupTable());
1637        switch (_colorMode) {
1638        case COLOR_CONSTANT:
1639        case COLOR_BY_SCALAR:
1640            _lut->SetRange(_dataRange);
1641            break;
1642        case COLOR_BY_VECTOR_MAGNITUDE:
1643            _lut->SetRange(_vectorMagnitudeRange);
1644            break;
1645        case COLOR_BY_VECTOR_X:
1646            _lut->SetRange(_vectorComponentRange[0]);
1647            break;
1648        case COLOR_BY_VECTOR_Y:
1649            _lut->SetRange(_vectorComponentRange[1]);
1650            break;
1651        case COLOR_BY_VECTOR_Z:
1652            _lut->SetRange(_vectorComponentRange[2]);
1653            break;
1654        default:
1655            break;
1656        }
1657    } else {
1658        double range[2];
1659        _lut->GetTableRange(range);
1660        _lut->DeepCopy(cmap->getLookupTable());
1661        _lut->SetRange(range);
1662        _lut->Modified();
1663    }
1664
1665    switch (_colorMode) {
1666    case COLOR_BY_VECTOR_MAGNITUDE:
1667        _lut->SetVectorModeToMagnitude();
1668        break;
1669    case COLOR_BY_VECTOR_X:
1670        _lut->SetVectorModeToComponent();
1671        _lut->SetVectorComponent(0);
1672        break;
1673    case COLOR_BY_VECTOR_Y:
1674        _lut->SetVectorModeToComponent();
1675        _lut->SetVectorComponent(1);
1676        break;
1677    case COLOR_BY_VECTOR_Z:
1678        _lut->SetVectorModeToComponent();
1679        _lut->SetVectorComponent(2);
1680        break;
1681    default:
1682         break;
1683    }
1684}
1685
1686/**
1687 * \brief Turn on/off lighting of this object
1688 */
1689void Streamlines::setLighting(bool state)
1690{
1691    _lighting = state;
1692    if (_linesActor != NULL)
1693        _linesActor->GetProperty()->SetLighting((state ? 1 : 0));
1694}
1695
1696/**
1697 * \brief Set opacity of this object
1698 */
1699void Streamlines::setOpacity(double opacity)
1700{
1701    _opacity = opacity;
1702    if (_linesActor != NULL) {
1703        _linesActor->GetProperty()->SetOpacity(_opacity);
1704        if (_opacity < 1.0)
1705            setCulling(_linesActor->GetProperty(), false);
1706        else if (_faceCulling && _lineType == TUBES)
1707            setCulling(_linesActor->GetProperty(), true);
1708    }
1709    if (_seedActor != NULL) {
1710        _seedActor->GetProperty()->SetOpacity(_opacity);
1711    }
1712}
1713
1714/**
1715 * \brief Turn on/off rendering of this Streamlines
1716 */
1717void Streamlines::setVisibility(bool state)
1718{
1719    if (_linesActor != NULL) {
1720        _linesActor->SetVisibility((state ? 1 : 0));
1721    }
1722    if (_seedActor != NULL) {
1723        if (!state ||
1724            (state && _seedVisible)) {
1725            _seedActor->SetVisibility((state ? 1 : 0));
1726        }
1727    }
1728}
1729
1730/**
1731 * \brief Turn on/off rendering of the seed geometry
1732 */
1733void Streamlines::setSeedVisibility(bool state)
1734{
1735    _seedVisible = state;
1736    if (_seedActor != NULL) {
1737        _seedActor->SetVisibility((state ? 1 : 0));
1738    }
1739}
1740
1741/**
1742 * \brief Get visibility state of the Streamlines
1743 *
1744 * \return Are the Streamlines visible?
1745 */
1746bool Streamlines::getVisibility()
1747{
1748    if (_linesActor == NULL) {
1749        return false;
1750    } else {
1751        return (_linesActor->GetVisibility() != 0);
1752    }
1753}
1754
1755/**
1756 * \brief Turn on/off rendering of edges
1757 */
1758void Streamlines::setEdgeVisibility(bool state)
1759{
1760    if (_linesActor != NULL) {
1761        _linesActor->GetProperty()->SetEdgeVisibility((state ? 1 : 0));
1762    }
1763}
1764
1765/**
1766 * \brief Set RGB color of stream lines
1767 */
1768void Streamlines::setColor(float color[3])
1769{
1770    _color[0] = color[0];
1771    _color[1] = color[1];
1772    _color[2] = color[2];
1773    if (_linesActor != NULL)
1774        _linesActor->GetProperty()->SetColor(_color[0], _color[1], _color[2]);
1775}
1776
1777/**
1778 * \brief Set RGB color of stream line edges
1779 */
1780void Streamlines::setEdgeColor(float color[3])
1781{
1782    _edgeColor[0] = color[0];
1783    _edgeColor[1] = color[1];
1784    _edgeColor[2] = color[2];
1785    if (_linesActor != NULL)
1786        _linesActor->GetProperty()->SetEdgeColor(_edgeColor[0], _edgeColor[1], _edgeColor[2]);
1787}
1788
1789/**
1790 * \brief Set RGB color of seed geometry
1791 */
1792void Streamlines::setSeedColor(float color[3])
1793{
1794    _seedColor[0] = color[0];
1795    _seedColor[1] = color[1];
1796    _seedColor[2] = color[2];
1797    if (_seedActor != NULL) {
1798        _seedActor->GetProperty()->SetColor(_seedColor[0], _seedColor[1], _seedColor[2]);
1799        _seedActor->GetProperty()->SetEdgeColor(_seedColor[0], _seedColor[1], _seedColor[2]);
1800    }
1801}
1802
1803/**
1804 * \brief Set point size of seed geometry (may be a no-op)
1805 */
1806void Streamlines::setSeedPointSize(float size)
1807{
1808    if (_seedActor != NULL) {
1809        _seedActor->GetProperty()->SetPointSize(size);
1810    }
1811}
1812
1813/**
1814 * \brief Set pixel width of stream lines (may be a no-op)
1815 */
1816void Streamlines::setEdgeWidth(float edgeWidth)
1817{
1818    _edgeWidth = edgeWidth;
1819    if (_linesActor != NULL)
1820        _linesActor->GetProperty()->SetLineWidth(_edgeWidth);
1821}
1822
1823/**
1824 * \brief Set a group of world coordinate planes to clip rendering
1825 *
1826 * Passing NULL for planes will remove all cliping planes
1827 */
1828void Streamlines::setClippingPlanes(vtkPlaneCollection *planes)
1829{
1830    if (_pdMapper != NULL) {
1831        _pdMapper->SetClippingPlanes(planes);
1832    }
1833    if (_seedMapper != NULL) {
1834        _seedMapper->SetClippingPlanes(planes);
1835    }
1836}
Note: See TracBrowser for help on using the repository browser.