source: branches/1.3/packages/vizservers/vtkvis/LIC.cpp @ 4530

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

Merge r4024,4049,4051-4052,4060 from trunk

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