source: nanovis/trunk/DataSetResample.cpp @ 5704

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

Add VTK reader/resampler to nanovis

File size: 10.2 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2/*
3 * Copyright (C) 2004-2013  HUBzero Foundation, LLC
4 *
5 */
6
7#include <iostream>
8#include <float.h>
9
10#include <vtkSmartPointer.h>
11#include <vtkDataSet.h>
12#include <vtkProbeFilter.h>
13#include <vtkCutter.h>
14#include <vtkImageData.h>
15#include <vtkImageResample.h>
16#include <vtkPolyData.h>
17#include <vtkUnstructuredGrid.h>
18#include <vtkDelaunay2D.h>
19#include <vtkDelaunay3D.h>
20#include <vtkTransform.h>
21#include <vtkGaussianSplatter.h>
22#include <vtkExtractVOI.h>
23#include <vtkShepardMethod.h>
24
25#include "DataSetResample.h"
26#include "Trace.h"
27
28enum PrincipalPlane {
29    PLANE_XY,
30    PLANE_ZY,
31    PLANE_XZ
32};
33
34static bool is2D(vtkDataSet *dataSet, PrincipalPlane *plane, double *offset)
35{
36    double bounds[6];
37    dataSet->GetBounds(bounds);
38    if (bounds[4] == bounds[5]) {
39        // Zmin = Zmax, XY plane
40        if (plane != NULL) {
41            *plane = PLANE_XY;
42        }
43        if (offset != NULL)
44            *offset = bounds[4];
45        return true;
46    } else if (bounds[0] == bounds[1]) {
47        // Xmin = Xmax, ZY plane
48        if (plane != NULL) {
49            *plane = PLANE_ZY;
50        }
51        if (offset != NULL)
52            *offset = bounds[0];
53        return true;
54    } else if (bounds[2] == bounds[3]) {
55        // Ymin = Ymax, XZ plane
56        if (plane != NULL) {
57            *plane = PLANE_XZ;
58         }
59        if (offset != NULL)
60            *offset = bounds[2];
61        return true;
62    }
63    return false;
64}
65
66/**
67 * \brief Determines if mesh is a point cloud (has no cells)
68 */
69static bool isCloud(vtkDataSet *dataSet)
70{
71    vtkPolyData *pd = vtkPolyData::SafeDownCast(dataSet);
72    if (pd) {
73        // If PolyData, is a cloud if there are no cells other than vertices
74        if (pd->GetNumberOfLines() == 0 &&
75            pd->GetNumberOfPolys() == 0 &&
76            pd->GetNumberOfStrips() == 0) {
77            return true;
78        } else {
79            // Has cells
80            return false;
81        }
82    } else {
83        return (dataSet->GetNumberOfCells() == 0);
84    }
85}
86
87vtkSmartPointer<vtkImageData>
88nv::resampleVTKDataSet(vtkDataSet *dataSetIn, int maxDim, nv::CloudStyle cloudStyle)
89{
90    double bounds[6];
91    dataSetIn->GetBounds(bounds);
92
93    int xdim, ydim, zdim;
94    double xSpacing, ySpacing, zSpacing;
95
96    vtkSmartPointer<vtkImageData> dataSetOut;
97    if (vtkImageData::SafeDownCast(dataSetIn) != NULL) {
98        // We have a uniform grid, use vtkImageResample
99        vtkImageData *imgIn = vtkImageData::SafeDownCast(dataSetIn);
100
101        TRACE("Image input Dims: %d,%d,%d Spacing: %g,%g,%g",
102              imgIn->GetDimensions()[0],
103              imgIn->GetDimensions()[1],
104              imgIn->GetDimensions()[2],
105              imgIn->GetSpacing()[0],
106              imgIn->GetSpacing()[1],
107              imgIn->GetSpacing()[2]);
108
109        // Input DataSet is vtkImageData
110        xdim = (imgIn->GetDimensions()[0] > maxDim) ? maxDim : imgIn->GetDimensions()[0];
111        ydim = (imgIn->GetDimensions()[1] > maxDim) ? maxDim : imgIn->GetDimensions()[1];
112        zdim = (imgIn->GetDimensions()[2] > maxDim) ? maxDim : imgIn->GetDimensions()[2];
113
114        xSpacing = (bounds[1]-bounds[0])/((double)(xdim-1));
115        ySpacing = (bounds[3]-bounds[2])/((double)(ydim-1));
116        zSpacing = (bounds[5]-bounds[4])/((double)(zdim-1));
117
118        vtkSmartPointer<vtkImageResample> filter = vtkSmartPointer<vtkImageResample>::New();
119        filter->SetInputData(dataSetIn);
120        filter->SetAxisOutputSpacing(0, xSpacing);
121        filter->SetAxisOutputSpacing(1, ySpacing);
122        filter->SetAxisOutputSpacing(2, zSpacing);
123        filter->Update();
124        dataSetOut = filter->GetOutput();
125    } else {
126        double xLen = bounds[1] - bounds[0];
127        double yLen = bounds[3] - bounds[2];
128        double zLen = bounds[5] - bounds[4];
129
130        xdim = (xLen == 0.0 ? 1 : maxDim);
131        ydim = (yLen == 0.0 ? 1 : maxDim);
132        zdim = (zLen == 0.0 ? 1 : maxDim);
133
134        xSpacing = (xLen == 0.0 ? 0.0 : xLen/((double)(xdim-1)));
135        ySpacing = (yLen == 0.0 ? 0.0 : yLen/((double)(ydim-1)));
136        zSpacing = (zLen == 0.0 ? 0.0 : zLen/((double)(zdim-1)));
137
138        vtkSmartPointer<vtkImageData> imageData = vtkSmartPointer<vtkImageData>::New();
139        imageData->SetDimensions(xdim, ydim, zdim);
140        imageData->SetOrigin(bounds[0], bounds[2], bounds[4]);
141        imageData->SetSpacing(xSpacing, ySpacing, zSpacing); 
142
143        vtkSmartPointer<vtkProbeFilter> filter = vtkSmartPointer<vtkProbeFilter>::New();
144        filter->SetInputData(imageData);
145
146        if (isCloud(dataSetIn)) {
147            // Need to mesh or splat cloud
148            if (cloudStyle == CLOUD_MESH) {
149                PrincipalPlane plane;
150                double offset;
151                if (is2D(dataSetIn, &plane, &offset)) {
152                    vtkSmartPointer<vtkDelaunay2D> mesher = vtkSmartPointer<vtkDelaunay2D>::New();
153                    if (plane == PLANE_ZY) {
154                        vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
155                        trans->RotateWXYZ(90, 0, 1, 0);
156                        if (offset != 0.0) {
157                            trans->Translate(-offset, 0, 0);
158                        }
159                        mesher->SetTransform(trans);
160                    } else if (plane == PLANE_XZ) {
161                        vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
162                        trans->RotateWXYZ(-90, 1, 0, 0);
163                        if (offset != 0.0) {
164                            trans->Translate(0, -offset, 0);
165                        }
166                        mesher->SetTransform(trans);
167                    } else if (offset != 0.0) {
168                        // XY with Z offset
169                        vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
170                        trans->Translate(0, 0, -offset);
171                        mesher->SetTransform(trans);
172                    }
173                    mesher->SetInputData(dataSetIn);
174                    filter->SetSourceConnection(mesher->GetOutputPort());
175                    //filter->SpatialMatchOn();
176                    filter->Update();
177                    dataSetOut = filter->GetImageDataOutput();
178                } else {
179                    vtkSmartPointer<vtkDelaunay3D> mesher = vtkSmartPointer<vtkDelaunay3D>::New();
180                    mesher->SetInputData(dataSetIn);
181                    filter->SetSourceConnection(mesher->GetOutputPort());
182                    //filter->SpatialMatchOn();
183                    filter->Update();
184                    dataSetOut = filter->GetImageDataOutput();
185                }
186            } else if (cloudStyle == CLOUD_SPLAT) {
187                // CLOUD_SPLAT
188                vtkSmartPointer<vtkGaussianSplatter> splatter = vtkSmartPointer<vtkGaussianSplatter>::New();
189                splatter->SetInputData(dataSetIn);
190                PrincipalPlane plane;
191                double offset;
192                int dims[3];
193                dims[0] = xdim;
194                dims[1] = ydim;
195                dims[2] = zdim;
196                if (is2D(dataSetIn, &plane, &offset)) {
197                    vtkSmartPointer<vtkExtractVOI> volumeSlicer = vtkSmartPointer<vtkExtractVOI>::New();
198                    if (plane == PLANE_ZY) {
199                        dims[0] = 3;
200                        volumeSlicer->SetVOI(1, 1, 0, dims[1]-1, 0, dims[1]-1);
201                    } else if (plane == PLANE_XZ) {
202                        dims[1] = 3;
203                        volumeSlicer->SetVOI(0, dims[0]-1, 1, 1, 0, dims[2]-1);
204                    } else {
205                        dims[2] = 3;
206                        volumeSlicer->SetVOI(0, dims[0]-1, 0, dims[1]-1, 1, 1);
207                    }
208                    splatter->SetSampleDimensions(dims);
209                    volumeSlicer->SetInputConnection(splatter->GetOutputPort());
210                    volumeSlicer->SetSampleRate(1, 1, 1);
211                    volumeSlicer->Update();
212                    dataSetOut = volumeSlicer->GetOutput();
213                } else {
214                    splatter->SetSampleDimensions(dims);
215                    splatter->Update();
216                    dataSetOut = splatter->GetOutput();
217                }
218            } else {
219                // CLOUD_SHEAPRD_METHOD
220                vtkSmartPointer<vtkShepardMethod> splatter = vtkSmartPointer<vtkShepardMethod>::New();
221                splatter->SetInputData(dataSetIn);
222                PrincipalPlane plane;
223                double offset;
224                int dims[3];
225                dims[0] = xdim;
226                dims[1] = ydim;
227                dims[2] = zdim;
228                if (is2D(dataSetIn, &plane, &offset)) {
229                    vtkSmartPointer<vtkExtractVOI> volumeSlicer = vtkSmartPointer<vtkExtractVOI>::New();
230                    if (plane == PLANE_ZY) {
231                        dims[0] = 3;
232                        volumeSlicer->SetVOI(1, 1, 0, dims[1]-1, 0, dims[1]-1);
233                    } else if (plane == PLANE_XZ) {
234                        dims[1] = 3;
235                        volumeSlicer->SetVOI(0, dims[0]-1, 1, 1, 0, dims[2]-1);
236                    } else {
237                        dims[2] = 3;
238                        volumeSlicer->SetVOI(0, dims[0]-1, 0, dims[1]-1, 1, 1);
239                    }
240                    splatter->SetSampleDimensions(dims);
241                    volumeSlicer->SetInputConnection(splatter->GetOutputPort());
242                    volumeSlicer->SetSampleRate(1, 1, 1);
243                    volumeSlicer->Update();
244                    dataSetOut = volumeSlicer->GetOutput();
245                } else {
246                    splatter->SetSampleDimensions(dims);
247                    splatter->Update();
248                    dataSetOut = splatter->GetOutput();
249                }
250            }
251        } else {
252            filter->SetSourceData(dataSetIn);
253            //filter->SpatialMatchOn();
254            filter->Update();
255            dataSetOut = filter->GetImageDataOutput();
256        }
257    }
258
259    TRACE("Image output Dims: %d,%d,%d Spacing: %g,%g,%g",
260          dataSetOut->GetDimensions()[0],
261          dataSetOut->GetDimensions()[1],
262          dataSetOut->GetDimensions()[2],
263          dataSetOut->GetSpacing()[0],
264          dataSetOut->GetSpacing()[1],
265          dataSetOut->GetSpacing()[2]);
266
267    return dataSetOut;
268}
Note: See TracBrowser for help on using the repository browser.