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