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 | |
---|
28 | enum PrincipalPlane { |
---|
29 | PLANE_XY, |
---|
30 | PLANE_ZY, |
---|
31 | PLANE_XZ |
---|
32 | }; |
---|
33 | |
---|
34 | static 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 | */ |
---|
69 | static 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 | |
---|
87 | vtkSmartPointer<vtkImageData> |
---|
88 | nv::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 | } |
---|