source: trunk/packages/vizservers/vtkvis/LIC.cpp @ 3999

Last change on this file since 3999 was 3999, checked in by ldelgass, 6 years ago

Fix for building against VTK head

  • Property svn:eol-style set to native
File size: 15.9 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 <cassert>
9
10#include <vtkActor.h>
11#include <vtkProperty.h>
12#include <vtkPointData.h>
13#include <vtkCellData.h>
14#include <vtkCellDataToPointData.h>
15#include <vtkImageData.h>
16#include <vtkPolyData.h>
17#include <vtkDataSetMapper.h>
18#include <vtkPainterPolyDataMapper.h>
19#include <vtkDataSetSurfaceFilter.h>
20#include <vtkCutter.h>
21#include <vtkImageMask.h>
22#include <vtkImageCast.h>
23#include <vtkDelaunay2D.h>
24#include <vtkDelaunay3D.h>
25#include <vtkTransform.h>
26
27#include "LIC.h"
28#include "Trace.h"
29
30using namespace VtkVis;
31
32LIC::LIC() :
33    GraphicsObject(),
34    _sliceAxis(Z_AXIS),
35    _colorMap(NULL)
36{
37}
38
39LIC::~LIC()
40{
41#ifdef WANT_TRACE
42    if (_dataSet != NULL)
43        TRACE("Deleting LIC for %s", _dataSet->getName().c_str());
44    else
45        TRACE("Deleting LIC with NULL DataSet");
46#endif
47}
48
49void LIC::initProp()
50{
51    if (_prop == NULL) {
52        _prop = vtkSmartPointer<vtkActor>::New();
53        vtkProperty *property = getActor()->GetProperty();
54        property->SetOpacity(_opacity);
55        property->SetEdgeColor(_edgeColor[0], _edgeColor[1], _edgeColor[2]);
56        property->SetLineWidth(_edgeWidth);
57        property->EdgeVisibilityOff();
58        property->SetAmbient(.2);
59        if (!_lighting)
60            property->LightingOff();
61    }
62}
63
64void LIC::update()
65{
66    if (_dataSet == NULL)
67        return;
68
69    vtkDataSet *ds = _dataSet->getVtkDataSet();
70
71    vtkSmartPointer<vtkCellDataToPointData> cellToPtData;
72    if (ds->GetPointData() == NULL ||
73        ds->GetPointData()->GetVectors() == NULL) {
74         TRACE("No vector point data in dataset %s", _dataSet->getName().c_str());
75         if (ds->GetCellData() != NULL &&
76             ds->GetCellData()->GetVectors() != NULL) {
77             cellToPtData = 
78                 vtkSmartPointer<vtkCellDataToPointData>::New();
79#ifdef USE_VTK6
80             cellToPtData->SetInputData(ds);
81#else
82             cellToPtData->SetInput(ds);
83#endif
84             //cellToPtData->PassCellDataOn();
85             cellToPtData->Update();
86             ds = cellToPtData->GetOutput();
87         } else {
88             USER_ERROR("No vector field was found in the data set");
89             return;
90         }
91    }
92
93    vtkImageData *imageData = vtkImageData::SafeDownCast(ds);
94    if (imageData != NULL) {
95        if (_lic == NULL) {
96            _lic = vtkSmartPointer<vtkImageDataLIC2D>::New();
97        }
98        if (!_dataSet->is2D()) {
99            // 3D image/volume/uniform grid
100            if (_volumeSlicer == NULL)
101                _volumeSlicer = vtkSmartPointer<vtkExtractVOI>::New();
102            int dims[3];
103            imageData->GetDimensions(dims);
104            TRACE("Image data dimensions: %d %d %d", dims[0], dims[1], dims[2]);
105#ifdef USE_VTK6
106            _volumeSlicer->SetInputData(ds);
107#else
108            _volumeSlicer->SetInput(ds);
109#endif
110            _volumeSlicer->SetVOI(0, dims[0]-1, 0, dims[1]-1, (dims[2]-1)/2, (dims[2]-1)/2);
111            _volumeSlicer->SetSampleRate(1, 1, 1);
112            _lic->SetInputConnection(_volumeSlicer->GetOutputPort());
113        } else {
114            // 2D image/uniform grid
115#ifdef USE_VTK6
116            _lic->SetInputData(ds);
117#else
118            _lic->SetInput(ds);
119#endif
120        }
121        if (_mapper == NULL) {
122            _mapper = vtkSmartPointer<vtkDataSetMapper>::New();
123        }
124        _mapper->SetInputConnection(_lic->GetOutputPort());
125    } else if (vtkPolyData::SafeDownCast(ds) == NULL ||
126               _dataSet->isCloud() ||
127               _dataSet->is2D()) {
128        // structured/rectilinear/unstructured grid, cloud or 2D polydata
129        if (_lic == NULL) {
130            _lic = vtkSmartPointer<vtkImageDataLIC2D>::New();
131        }
132
133        // Need to convert to vtkImageData
134        double bounds[6];
135        ds->GetBounds(bounds);
136        double xSize = bounds[1] - bounds[0];
137        double ySize = bounds[3] - bounds[2];
138        double zSize = bounds[5] - bounds[4];
139        int minDir = 2;
140        if (xSize < ySize && xSize < zSize)
141            minDir = 0;
142        if (ySize < xSize && ySize < zSize)
143            minDir = 1;
144        if (_probeFilter == NULL)
145            _probeFilter = vtkSmartPointer<vtkProbeFilter>::New();
146
147        PrincipalPlane plane;
148        double offset;
149        if (!_dataSet->is2D(&plane, &offset)) {
150            // 3D Data: Sample a plane within the bounding box
151            vtkSmartPointer<vtkCutter> cutter = vtkSmartPointer<vtkCutter>::New();
152            if (_cutPlane == NULL) {
153                _cutPlane = vtkSmartPointer<vtkPlane>::New();
154            }
155            if (minDir == 0) {
156                _cutPlane->SetNormal(1, 0, 0);
157                _cutPlane->SetOrigin(bounds[0] + (bounds[1]-bounds[0])/2.,
158                                     0,
159                                     0);
160            } else if (minDir == 1) {
161                _cutPlane->SetNormal(0, 1, 0);
162                _cutPlane->SetOrigin(0,
163                                     bounds[2] + (bounds[3]-bounds[2])/2.,
164                                     0);
165            } else {
166                _cutPlane->SetNormal(0, 0, 1);
167                _cutPlane->SetOrigin(0,
168                                     0,
169                                     bounds[4] + (bounds[5]-bounds[4])/2.);
170            }
171
172            if (_dataSet->isCloud()) {
173                // 3D cloud -- Need to mesh it before we can resample
174                vtkSmartPointer<vtkDelaunay3D> mesher = vtkSmartPointer<vtkDelaunay3D>::New();
175#ifdef USE_VTK6
176                mesher->SetInputData(ds);
177#else
178                mesher->SetInput(ds);
179#endif
180                cutter->SetInputConnection(mesher->GetOutputPort());
181            } else {
182#ifdef USE_VTK6
183                cutter->SetInputData(ds);
184#else
185                cutter->SetInput(ds);
186#endif
187            }
188            cutter->SetCutFunction(_cutPlane);
189            _probeFilter->SetSourceConnection(cutter->GetOutputPort());
190        } else {
191            if (_dataSet->isCloud()) {
192                // 2D cloud -- Need to mesh it before we can resample
193                vtkSmartPointer<vtkDelaunay2D> mesher = vtkSmartPointer<vtkDelaunay2D>::New();
194                if (plane == PLANE_ZY) {
195                    vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
196                    trans->RotateWXYZ(90, 0, 1, 0);
197                    if (offset != 0.0) {
198                        trans->Translate(-offset, 0, 0);
199                    }
200                    mesher->SetTransform(trans);
201                } else if (plane == PLANE_XZ) {
202                    vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
203                    trans->RotateWXYZ(-90, 1, 0, 0);
204                    if (offset != 0.0) {
205                        trans->Translate(0, -offset, 0);
206                    }
207                    mesher->SetTransform(trans);
208                } else if (offset != 0.0) {
209                    // XY with Z offset
210                    vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
211                    trans->Translate(0, 0, -offset);
212                    mesher->SetTransform(trans);
213                }
214#ifdef USE_VTK6
215                mesher->SetInputData(ds);
216#else
217                mesher->SetInput(ds);
218#endif
219                _probeFilter->SetSourceConnection(mesher->GetOutputPort());
220            } else {
221#ifdef USE_VTK6
222                _probeFilter->SetSourceData(ds);
223#else
224                _probeFilter->SetSource(ds);
225#endif
226            }
227        }
228
229        vtkSmartPointer<vtkImageData> imageData = vtkSmartPointer<vtkImageData>::New();
230        int xdim, ydim, zdim;
231        double origin[3], spacing[3];
232        if (minDir == 0) {
233            xdim = 1;
234            ydim = 128;
235            zdim = 128;
236            origin[0] = bounds[0] + xSize/2.;
237            origin[1] = bounds[2];
238            origin[2] = bounds[4];
239            spacing[0] = 0;
240            spacing[1] = ySize/((double)(ydim-1));
241            spacing[2] = zSize/((double)(zdim-1));
242            _sliceAxis = X_AXIS;
243        } else if (minDir == 1) {
244            xdim = 128;
245            ydim = 1;
246            zdim = 128;
247            origin[0] = bounds[0];
248            origin[1] = bounds[2] + ySize/2.;
249            origin[2] = bounds[4];
250            spacing[0] = xSize/((double)(xdim-1));
251            spacing[1] = 0;
252            spacing[2] = zSize/((double)(zdim-1));
253            _sliceAxis = Y_AXIS;
254        } else {
255            xdim = 128;
256            ydim = 128;
257            zdim = 1;
258            origin[0] = bounds[0];
259            origin[1] = bounds[2];
260            origin[2] = bounds[4] + zSize/2.;
261            spacing[0] = xSize/((double)(xdim-1));
262            spacing[1] = ySize/((double)(ydim-1));
263            spacing[2] = 0;
264            _sliceAxis = Z_AXIS;
265        }
266        imageData->SetDimensions(xdim, ydim, zdim);
267        imageData->SetOrigin(origin);
268        imageData->SetSpacing(spacing);
269#ifdef USE_VTK6
270        _probeFilter->SetInputData(imageData);
271#else
272        _probeFilter->SetInput(imageData);
273#endif
274        _lic->SetInputConnection(_probeFilter->GetOutputPort());
275
276        if (_mapper == NULL) {
277            _mapper = vtkSmartPointer<vtkDataSetMapper>::New();
278            _mapper->SetColorModeToMapScalars();
279        }
280        _lic->Update();
281        vtkSmartPointer<vtkImageCast> cast = vtkSmartPointer<vtkImageCast>::New();
282        cast->SetInputConnection(_probeFilter->GetOutputPort());
283        cast->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS,
284                                     _probeFilter->GetValidPointMaskArrayName());
285        cast->SetOutputScalarTypeToUnsignedChar();
286        vtkSmartPointer<vtkImageMask> mask = vtkSmartPointer<vtkImageMask>::New();
287        mask->SetInputConnection(0, _lic->GetOutputPort());
288        mask->SetInputConnection(1, cast->GetOutputPort());
289        _mapper->SetInputConnection(mask->GetOutputPort());
290    } else {
291        // DataSet is a 3D PolyData with cells
292        vtkPolyData *pd = vtkPolyData::SafeDownCast(ds);
293        assert(pd != NULL);
294
295        // XXX: This mapper seems to conflict with offscreen rendering, so the main
296        // renderer needs to be set not to use FBOs for this to work
297        _mapper = vtkSmartPointer<vtkPainterPolyDataMapper>::New();
298        vtkPainterPolyDataMapper *ppdmapper = vtkPainterPolyDataMapper::SafeDownCast(_mapper);
299        if (_painter == NULL) {
300            _painter = vtkSmartPointer<vtkSurfaceLICPainter>::New();
301            _painter->SetDelegatePainter(ppdmapper->GetPainter());
302         }
303        ppdmapper->SetPainter(_painter);
304#ifdef USE_VTK6
305        ppdmapper->SetInputData(pd);
306#else
307        ppdmapper->SetInput(pd);
308#endif
309    }
310
311    setInterpolateBeforeMapping(true);
312
313    if (_lut == NULL) {
314        setColorMap(ColorMap::getGrayDefault());
315    }
316
317    initProp();
318    getActor()->SetMapper(_mapper);
319
320    _mapper->Update();
321}
322
323void LIC::setInterpolateBeforeMapping(bool state)
324{
325    if (_mapper != NULL) {
326        _mapper->SetInterpolateScalarsBeforeMapping((state ? 1 : 0));
327    }
328}
329
330/**
331 * \brief Select a 2D slice plane from a 3D DataSet
332 *
333 * \param[in] axis Axis of slice plane
334 * \param[in] ratio Position [0,1] of slice plane along axis
335 */
336void LIC::selectVolumeSlice(Axis axis, double ratio)
337{
338    if (_dataSet->is2D()) {
339        WARN("DataSet not 3D, returning");
340        return;
341    }
342
343    if (_volumeSlicer == NULL &&
344        _probeFilter == NULL) {
345        WARN("Called before update() or DataSet is not a volume");
346        return;
347    }
348
349    int dims[3];
350
351    vtkImageData *imageData = vtkImageData::SafeDownCast(_dataSet->getVtkDataSet());
352    if (imageData == NULL) {
353        if (_probeFilter != NULL) {
354            imageData = vtkImageData::SafeDownCast(_probeFilter->GetInput());
355            if (imageData == NULL) {
356                ERROR("Couldn't get probe filter input image");
357                return;
358            }
359        } else {
360            ERROR("Not a volume data set");
361            return;
362        }
363    }
364    imageData->GetDimensions(dims);
365
366    _sliceAxis = axis;
367
368    if (_probeFilter != NULL) {
369        vtkImageData *imageData = vtkImageData::SafeDownCast(_probeFilter->GetInput());
370        double bounds[6];
371        assert(vtkDataSet::SafeDownCast(_probeFilter->GetSource()) != NULL);
372        _dataSet->getBounds(bounds);
373        int dim = 128;
374
375        switch (axis) {
376        case X_AXIS:
377            imageData->SetDimensions(1, dim, dim);
378            imageData->SetOrigin(bounds[0] + (bounds[1]-bounds[0])*ratio, bounds[2], bounds[4]);
379            imageData->SetSpacing(0,
380                                  (bounds[3]-bounds[2])/((double)(dim-1)), 
381                                  (bounds[5]-bounds[4])/((double)(dim-1)));
382            if (_cutPlane != NULL) {
383                _cutPlane->SetNormal(1, 0, 0);
384                _cutPlane->SetOrigin(bounds[0] + (bounds[1]-bounds[0])*ratio, 0, 0);
385            }
386            break;
387        case Y_AXIS:
388            imageData->SetDimensions(dim, 1, dim);
389            imageData->SetOrigin(bounds[0], bounds[2] + (bounds[3]-bounds[2])*ratio, bounds[4]);
390            imageData->SetSpacing((bounds[1]-bounds[0])/((double)(dim-1)), 
391                                  0,
392                                  (bounds[5]-bounds[4])/((double)(dim-1)));
393            if (_cutPlane != NULL) {
394                _cutPlane->SetNormal(0, 1, 0);
395                _cutPlane->SetOrigin(0, bounds[2] + (bounds[3]-bounds[2])*ratio, 0);
396            }
397            break;
398        case Z_AXIS:
399            imageData->SetDimensions(dim, dim, 1);
400            imageData->SetOrigin(bounds[0], bounds[2], bounds[4] + (bounds[5]-bounds[4])*ratio);
401            imageData->SetSpacing((bounds[1]-bounds[0])/((double)(dim-1)), 
402                                  (bounds[3]-bounds[2])/((double)(dim-1)),
403                                  0);
404            if (_cutPlane != NULL) {
405                _cutPlane->SetNormal(0, 0, 1);
406                _cutPlane->SetOrigin(0, 0, bounds[4] + (bounds[5]-bounds[4])*ratio);
407            }
408            break;
409        default:
410            ERROR("Invalid Axis");
411            return;
412        }
413    } else {
414        int voi[6];
415
416        switch (axis) {
417        case X_AXIS:
418            voi[0] = voi[1] = (int)((dims[0]-1) * ratio);
419            voi[2] = 0;
420            voi[3] = dims[1]-1;
421            voi[4] = 0;
422            voi[5] = dims[2]-1;
423            break;
424        case Y_AXIS:
425            voi[2] = voi[3] = (int)((dims[1]-1) * ratio);
426            voi[0] = 0;
427            voi[1] = dims[0]-1;
428            voi[4] = 0;
429            voi[5] = dims[2]-1;
430            break;
431        case Z_AXIS:
432            voi[4] = voi[5] = (int)((dims[2]-1) * ratio);
433            voi[0] = 0;
434            voi[1] = dims[0]-1;
435            voi[2] = 0;
436            voi[3] = dims[1]-1;
437            break;
438        default:
439            ERROR("Invalid Axis");
440            return;
441        }
442
443        _volumeSlicer->SetVOI(voi);
444    }
445
446    if (_lic != NULL)
447        _lic->Update();
448
449    if (_mapper != NULL)
450        _mapper->Update();
451}
452
453/**
454 * \brief Called when the color map has been edited
455 */
456void LIC::updateColorMap()
457{
458    setColorMap(_colorMap);
459}
460
461/**
462 * \brief Associate a colormap lookup table with the DataSet
463 */
464void LIC::setColorMap(ColorMap *cmap)
465{
466    if (cmap == NULL)
467        return;
468
469    _colorMap = cmap;
470 
471    if (_lut == NULL) {
472        _lut = vtkSmartPointer<vtkLookupTable>::New();
473        if (_mapper != NULL) {
474            _mapper->UseLookupTableScalarRangeOff();
475            _mapper->SetLookupTable(_lut);
476        }
477    }
478
479    _lut->DeepCopy(cmap->getLookupTable());
480    _lut->SetRange(_dataRange);
481    _lut->Modified();
482}
483
484void LIC::updateRanges(Renderer *renderer)
485{
486    GraphicsObject::updateRanges(renderer);
487
488    if (_lut != NULL) {
489        _lut->SetRange(_dataRange);
490    }
491}
492
493/**
494 * \brief Set a group of world coordinate planes to clip rendering
495 *
496 * Passing NULL for planes will remove all cliping planes
497 */
498void LIC::setClippingPlanes(vtkPlaneCollection *planes)
499{
500    if (_mapper != NULL) {
501        _mapper->SetClippingPlanes(planes);
502    }
503}
Note: See TracBrowser for help on using the repository browser.