1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkMarchingContourFilter.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkMarchingContourFilter.h"
16
17 #include "vtkCell.h"
18 #include "vtkContourFilter.h"
19 #include "vtkContourValues.h"
20 #include "vtkImageMarchingCubes.h"
21 #include "vtkIncrementalPointLocator.h"
22 #include "vtkInformation.h"
23 #include "vtkInformationVector.h"
24 #include "vtkMarchingCubes.h"
25 #include "vtkMarchingSquares.h"
26 #include "vtkMergePoints.h"
27 #include "vtkObjectFactory.h"
28 #include "vtkPointData.h"
29 #include "vtkPolyData.h"
30 #include "vtkScalarTree.h"
31 #include "vtkStreamingDemandDrivenPipeline.h"
32 #include "vtkStructuredPoints.h"
33 #include "vtkTrivialProducer.h"
34
35 #include <cmath>
36
37 vtkStandardNewMacro(vtkMarchingContourFilter);
38
39 // Construct object with initial range (0,1) and single contour value
40 // of 0.0.
vtkMarchingContourFilter()41 vtkMarchingContourFilter::vtkMarchingContourFilter()
42 {
43 this->ContourValues = vtkContourValues::New();
44
45 this->ComputeNormals = 1;
46 this->ComputeGradients = 0;
47 this->ComputeScalars = 1;
48
49 this->Locator = nullptr;
50
51 this->UseScalarTree = 0;
52 this->ScalarTree = nullptr;
53 }
54
~vtkMarchingContourFilter()55 vtkMarchingContourFilter::~vtkMarchingContourFilter()
56 {
57 this->ContourValues->Delete();
58 if (this->Locator)
59 {
60 this->Locator->UnRegister(this);
61 this->Locator = nullptr;
62 }
63 if (this->ScalarTree)
64 {
65 this->ScalarTree->Delete();
66 }
67 }
68
69 // Overload standard modified time function. If contour values are modified,
70 // then this object is modified as well.
GetMTime()71 vtkMTimeType vtkMarchingContourFilter::GetMTime()
72 {
73 vtkMTimeType mTime = this->Superclass::GetMTime();
74 vtkMTimeType time;
75
76 if (this->ContourValues)
77 {
78 time = this->ContourValues->GetMTime();
79 mTime = (time > mTime ? time : mTime);
80 }
81 if (this->Locator)
82 {
83 time = this->Locator->GetMTime();
84 mTime = (time > mTime ? time : mTime);
85 }
86
87 return mTime;
88 }
89
90 //
91 // General contouring filter. Handles arbitrary input.
92 //
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)93 int vtkMarchingContourFilter::RequestData(vtkInformation* vtkNotUsed(request),
94 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
95 {
96 // get the info objects
97 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
98 vtkInformation* outInfo = outputVector->GetInformationObject(0);
99
100 // get the input and output
101 vtkDataSet* input = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
102 vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
103
104 vtkDataArray* inScalars;
105 vtkIdType numCells;
106
107 vtkDebugMacro(<< "Executing marching contour filter");
108
109 numCells = input->GetNumberOfCells();
110 inScalars = input->GetPointData()->GetScalars();
111 if (!inScalars || numCells < 1)
112 {
113 vtkErrorMacro(<< "No data to contour");
114 return 1;
115 }
116
117 // If structured points, use more efficient algorithms
118 if ((input->GetDataObjectType() == VTK_STRUCTURED_POINTS))
119 {
120 if (inScalars->GetDataType() != VTK_BIT)
121 {
122 int dim = input->GetCell(0)->GetCellDimension();
123
124 if (input->GetCell(0)->GetCellDimension() >= 2)
125 {
126 vtkDebugMacro(<< "Structured Points");
127 this->StructuredPointsContour(dim, input, output);
128 return 1;
129 }
130 }
131 }
132
133 if ((input->GetDataObjectType() == VTK_IMAGE_DATA))
134 {
135 if (inScalars->GetDataType() != VTK_BIT)
136 {
137 int dim = input->GetCell(0)->GetCellDimension();
138
139 if (input->GetCell(0)->GetCellDimension() >= 2)
140 {
141 vtkDebugMacro(<< "Image");
142 this->ImageContour(dim, input, output);
143 return 1;
144 }
145 }
146 }
147
148 vtkDebugMacro(<< "Unoptimized");
149 this->DataSetContour(input, output);
150
151 return 1;
152 }
153
StructuredPointsContour(int dim,vtkDataSet * input,vtkPolyData * thisOutput)154 void vtkMarchingContourFilter::StructuredPointsContour(
155 int dim, vtkDataSet* input, vtkPolyData* thisOutput)
156 {
157 vtkPolyData* output;
158 vtkIdType numContours = this->ContourValues->GetNumberOfContours();
159 double* values = this->ContourValues->GetValues();
160
161 if (dim == 2) // marching squares
162 {
163 vtkMarchingSquares* msquares;
164 int i;
165
166 msquares = vtkMarchingSquares::New();
167 msquares->SetInputData((vtkImageData*)input);
168 msquares->SetDebug(this->Debug);
169 msquares->SetNumberOfContours(numContours);
170 for (i = 0; i < numContours; i++)
171 {
172 msquares->SetValue(i, values[i]);
173 }
174
175 msquares->Update();
176 output = msquares->GetOutput();
177 output->Register(this);
178 msquares->Delete();
179 }
180
181 else // marching cubes
182 {
183 vtkMarchingCubes* mcubes;
184 int i;
185
186 mcubes = vtkMarchingCubes::New();
187 mcubes->SetInputData((vtkImageData*)input);
188 mcubes->SetComputeNormals(this->ComputeNormals);
189 mcubes->SetComputeGradients(this->ComputeGradients);
190 mcubes->SetComputeScalars(this->ComputeScalars);
191 mcubes->SetDebug(this->Debug);
192 mcubes->SetNumberOfContours(numContours);
193 for (i = 0; i < numContours; i++)
194 {
195 mcubes->SetValue(i, values[i]);
196 }
197
198 mcubes->Update();
199 output = mcubes->GetOutput();
200 output->Register(this);
201 mcubes->Delete();
202 }
203
204 thisOutput->CopyStructure(output);
205 thisOutput->GetPointData()->ShallowCopy(output->GetPointData());
206 output->UnRegister(this);
207 }
208
DataSetContour(vtkDataSet * input,vtkPolyData * output)209 void vtkMarchingContourFilter::DataSetContour(vtkDataSet* input, vtkPolyData* output)
210 {
211 vtkIdType numContours = this->ContourValues->GetNumberOfContours();
212 double* values = this->ContourValues->GetValues();
213
214 vtkContourFilter* contour = vtkContourFilter::New();
215 contour->SetInputData((vtkImageData*)input);
216 contour->SetComputeNormals(this->ComputeNormals);
217 contour->SetComputeGradients(this->ComputeGradients);
218 contour->SetComputeScalars(this->ComputeScalars);
219 contour->SetDebug(this->Debug);
220 contour->SetNumberOfContours(numContours);
221 for (int i = 0; i < numContours; i++)
222 {
223 contour->SetValue(i, values[i]);
224 }
225
226 contour->Update();
227 output->ShallowCopy(contour->GetOutput());
228 this->SetOutput(output);
229 contour->Delete();
230 }
231
ImageContour(int dim,vtkDataSet * input,vtkPolyData * output)232 void vtkMarchingContourFilter::ImageContour(int dim, vtkDataSet* input, vtkPolyData* output)
233 {
234 vtkIdType numContours = this->ContourValues->GetNumberOfContours();
235 double* values = this->ContourValues->GetValues();
236 vtkPolyData* contourOutput;
237
238 vtkNew<vtkTrivialProducer> producer;
239 producer->SetOutput(input);
240 // Explicitly update the update extent in the trivial producer to prevent
241 // an error when downstream algorithms request a different extent.
242 producer->UpdateWholeExtent();
243
244 if (dim == 2) // marching squares
245 {
246 vtkMarchingSquares* msquares;
247 int i;
248
249 msquares = vtkMarchingSquares::New();
250 msquares->SetInputConnection(producer->GetOutputPort());
251 msquares->SetDebug(this->Debug);
252 msquares->SetNumberOfContours(numContours);
253 for (i = 0; i < numContours; i++)
254 {
255 msquares->SetValue(i, values[i]);
256 }
257
258 contourOutput = msquares->GetOutput();
259 msquares->Update();
260 output->ShallowCopy(contourOutput);
261 msquares->Delete();
262 }
263
264 else // image marching cubes
265 {
266 vtkImageMarchingCubes* mcubes;
267 int i;
268
269 mcubes = vtkImageMarchingCubes::New();
270
271 mcubes->SetInputConnection(producer->GetOutputPort());
272 mcubes->SetComputeNormals(this->ComputeNormals);
273 mcubes->SetComputeGradients(this->ComputeGradients);
274 mcubes->SetComputeScalars(this->ComputeScalars);
275 mcubes->SetDebug(this->Debug);
276 mcubes->SetNumberOfContours(numContours);
277 for (i = 0; i < numContours; i++)
278 {
279 mcubes->SetValue(i, values[i]);
280 }
281
282 contourOutput = mcubes->GetOutput();
283 mcubes->Update();
284 output->ShallowCopy(contourOutput);
285 mcubes->Delete();
286 }
287 }
288
289 // Specify a spatial locator for merging points. By default,
290 // an instance of vtkMergePoints is used.
SetLocator(vtkIncrementalPointLocator * locator)291 void vtkMarchingContourFilter::SetLocator(vtkIncrementalPointLocator* locator)
292 {
293 if (this->Locator == locator)
294 {
295 return;
296 }
297 if (this->Locator)
298 {
299 this->Locator->UnRegister(this);
300 this->Locator = nullptr;
301 }
302 if (locator)
303 {
304 locator->Register(this);
305 }
306 this->Locator = locator;
307 this->Modified();
308 }
309
CreateDefaultLocator()310 void vtkMarchingContourFilter::CreateDefaultLocator()
311 {
312 if (this->Locator == nullptr)
313 {
314 this->Locator = vtkMergePoints::New();
315 }
316 }
317
FillInputPortInformation(int,vtkInformation * info)318 int vtkMarchingContourFilter::FillInputPortInformation(int, vtkInformation* info)
319 {
320 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
321 return 1;
322 }
323
PrintSelf(ostream & os,vtkIndent indent)324 void vtkMarchingContourFilter::PrintSelf(ostream& os, vtkIndent indent)
325 {
326 this->Superclass::PrintSelf(os, indent);
327
328 os << indent << "Compute Gradients: " << (this->ComputeGradients ? "On\n" : "Off\n");
329 os << indent << "Compute Normals: " << (this->ComputeNormals ? "On\n" : "Off\n");
330 os << indent << "Compute Scalars: " << (this->ComputeScalars ? "On\n" : "Off\n");
331 os << indent << "Use Scalar Tree: " << (this->UseScalarTree ? "On\n" : "Off\n");
332
333 this->ContourValues->PrintSelf(os, indent.GetNextIndent());
334
335 if (this->Locator)
336 {
337 os << indent << "Locator: " << this->Locator << "\n";
338 }
339 else
340 {
341 os << indent << "Locator: (none)\n";
342 }
343 }
344