1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkAreaContourSpectrumFilter.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 "vtkAreaContourSpectrumFilter.h"
16
17 #include "vtkDataSetAttributes.h"
18 #include "vtkEdgeListIterator.h"
19 #include "vtkIdList.h"
20 #include "vtkInformation.h"
21 #include "vtkInformationVector.h"
22 #include "vtkObjectFactory.h"
23 #include "vtkPointData.h"
24 #include "vtkPolyData.h"
25 #include "vtkReebGraph.h"
26 #include "vtkTable.h"
27 #include "vtkTriangle.h"
28 #include "vtkVariantArray.h"
29
30 vtkStandardNewMacro(vtkAreaContourSpectrumFilter);
31
32 //------------------------------------------------------------------------------
vtkAreaContourSpectrumFilter()33 vtkAreaContourSpectrumFilter::vtkAreaContourSpectrumFilter()
34 {
35 this->SetNumberOfInputPorts(2);
36 this->ArcId = 0;
37 this->FieldId = 0;
38 this->NumberOfSamples = 100;
39 }
40
41 //------------------------------------------------------------------------------
42 vtkAreaContourSpectrumFilter::~vtkAreaContourSpectrumFilter() = default;
43
44 //------------------------------------------------------------------------------
FillInputPortInformation(int portNumber,vtkInformation * info)45 int vtkAreaContourSpectrumFilter::FillInputPortInformation(int portNumber, vtkInformation* info)
46 {
47 switch (portNumber)
48 {
49 case 0:
50 info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
51 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
52 break;
53 case 1:
54 info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
55 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkReebGraph");
56 break;
57 }
58 return 1;
59 }
60
61 //------------------------------------------------------------------------------
FillOutputPortInformation(int vtkNotUsed (portNumber),vtkInformation * info)62 int vtkAreaContourSpectrumFilter::FillOutputPortInformation(
63 int vtkNotUsed(portNumber), vtkInformation* info)
64 {
65
66 info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkTable");
67
68 return 1;
69 }
70
71 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)72 void vtkAreaContourSpectrumFilter::PrintSelf(ostream& os, vtkIndent indent)
73 {
74 this->Superclass::PrintSelf(os, indent);
75 os << indent << "Arc Id: " << this->ArcId << "\n";
76 os << indent << "Field Id: " << this->FieldId << "\n";
77 os << indent << "Number of Samples: " << this->NumberOfSamples << "\n";
78 }
79
80 //------------------------------------------------------------------------------
GetOutput()81 vtkTable* vtkAreaContourSpectrumFilter::GetOutput()
82 {
83 return vtkTable::SafeDownCast(this->GetOutputDataObject(0));
84 }
85
86 //------------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)87 int vtkAreaContourSpectrumFilter::RequestData(vtkInformation* vtkNotUsed(request),
88 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
89 {
90
91 vtkInformation *inInfoMesh = inputVector[0]->GetInformationObject(0),
92 *inInfoGraph = inputVector[1]->GetInformationObject(0);
93
94 if ((!inInfoMesh) || (!inInfoGraph))
95 {
96 return 0;
97 }
98 vtkPolyData* inputMesh = vtkPolyData::SafeDownCast(inInfoMesh->Get(vtkPolyData::DATA_OBJECT()));
99
100 vtkReebGraph* inputGraph =
101 vtkReebGraph::SafeDownCast(inInfoGraph->Get(vtkReebGraph::DATA_OBJECT()));
102
103 if ((inputMesh) && (inputGraph))
104 {
105 vtkInformation* outInfo = outputVector->GetInformationObject(0);
106 vtkTable* output = vtkTable::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
107
108 if (output)
109 {
110
111 // Retrieve the arc given by ArcId
112 vtkVariantArray* edgeInfo = vtkArrayDownCast<vtkVariantArray>(
113 inputGraph->GetEdgeData()->GetAbstractArray("Vertex Ids"));
114 // Invalid Reeb graph (no information associated to the edges)
115 if (!edgeInfo)
116 return 0;
117
118 // Retrieve the information to get the critical vertices Ids
119 vtkDataArray* criticalPointIds =
120 vtkArrayDownCast<vtkDataArray>(inputGraph->GetVertexData()->GetAbstractArray("Vertex Ids"));
121 // Invalid Reeb graph (no information associated to the vertices)
122 if (!criticalPointIds)
123 return 0;
124
125 vtkAbstractArray* vertexList = edgeInfo->GetPointer(ArcId)->ToArray();
126
127 // the arc defined by ArcId does not exist (out of bound?)
128 if (!vertexList)
129 return 0;
130
131 vtkDataArray* scalarField = inputMesh->GetPointData()->GetArray(FieldId);
132 if (!scalarField)
133 return 0;
134
135 // parse the input vertex list (region in which the connectivity of the
136 // level sets does not change) and compute the area signature
137 double cumulativeArea = 0;
138 std::vector<double> scalarValues, areaSignature;
139 std::vector<int> vertexIds;
140 std::vector<bool> visitedTriangles;
141
142 visitedTriangles.resize(inputMesh->GetNumberOfCells());
143 vertexIds.resize(vertexList->GetNumberOfTuples() + 2);
144 scalarValues.resize(vertexIds.size());
145 areaSignature.resize(vertexIds.size());
146
147 // include the critical points in the computation
148 // - iterates through the edges of the Reeb graph until we found the arc
149 // we're looking for
150 // - retrieve the Source and Target of the edge
151 // - pick the corresponding mesh vertex Ids in the VertexData.
152 std::pair<int, int> criticalPoints;
153 vtkEdgeListIterator* eIt = vtkEdgeListIterator::New();
154 inputGraph->GetEdges(eIt);
155
156 do
157 {
158 vtkEdgeType e = eIt->Next();
159 if (e.Id == ArcId)
160 {
161 if ((criticalPointIds->GetTuple(e.Source)) && (criticalPointIds->GetTuple(e.Target)))
162 {
163 criticalPoints.first = (int)*(criticalPointIds->GetTuple(e.Source));
164 criticalPoints.second = (int)*(criticalPointIds->GetTuple(e.Target));
165 }
166 else
167 {
168 // invalid Reeb graph
169 return 0;
170 }
171 }
172 } while (eIt->HasNext());
173
174 eIt->Delete();
175
176 vertexIds[0] = criticalPoints.first;
177 vertexIds[vertexIds.size() - 1] = criticalPoints.second;
178 // NB: the vertices of vertexList are already in sorted order of function
179 // value.
180 for (int i = 0; i < vertexList->GetNumberOfTuples(); i++)
181 vertexIds[i + 1] = vertexList->GetVariantValue(i).ToInt();
182
183 // mark all the input triangles as non visited.
184 for (unsigned int i = 0; i < visitedTriangles.size(); i++)
185 visitedTriangles[i] = false;
186
187 // now do the parsing
188 double min = scalarField->GetComponent(vertexIds[0], 0),
189 max = scalarField->GetComponent(vertexIds[vertexIds.size() - 1], 0);
190 for (unsigned int i = 0; i < vertexIds.size(); i++)
191 {
192 scalarValues[i] = scalarField->GetComponent(vertexIds[i], 0);
193
194 vtkIdList* starTriangleList = vtkIdList::New();
195 inputMesh->GetPointCells(vertexIds[i], starTriangleList);
196
197 for (int j = 0; j < starTriangleList->GetNumberOfIds(); j++)
198 {
199 vtkIdType tId = starTriangleList->GetId(j);
200 if (!visitedTriangles[tId])
201 {
202 vtkTriangle* t = vtkTriangle::SafeDownCast(inputMesh->GetCell(tId));
203
204 if ((scalarField->GetComponent(t->GetPointIds()->GetId(0), 0) <= scalarValues[i]) &&
205 (scalarField->GetComponent(t->GetPointIds()->GetId(1), 0) <= scalarValues[i]) &&
206 (scalarField->GetComponent(t->GetPointIds()->GetId(2), 0) <= scalarValues[i]) &&
207 (scalarField->GetComponent(t->GetPointIds()->GetId(0), 0) >= min) &&
208 (scalarField->GetComponent(t->GetPointIds()->GetId(1), 0) >= min) &&
209 (scalarField->GetComponent(t->GetPointIds()->GetId(2), 0) >= min))
210 {
211 // make sure the triangle is strictly in the covered function
212 // span.
213 cumulativeArea += t->ComputeArea();
214 visitedTriangles[tId] = true;
215 }
216 }
217 }
218 areaSignature[i] = cumulativeArea;
219 starTriangleList->Delete();
220 }
221
222 // now adjust to the desired sampling
223 std::vector<std::pair<int, double>> samples(NumberOfSamples);
224 unsigned int pos = 0;
225 for (int i = 0; i < NumberOfSamples; i++)
226 {
227 samples[i].first = 0;
228 samples[i].second = 0;
229 double temp = min + (i + 1.0) * ((max - min) / ((double)NumberOfSamples));
230 while ((pos < scalarValues.size()) && (scalarValues[pos] < temp))
231 {
232 samples[i].first++;
233 samples[i].second += areaSignature[pos];
234 pos++;
235 }
236 if (samples[i].first)
237 {
238 samples[i].second /= samples[i].first;
239 }
240 }
241
242 // no value at the start? put 0
243 if (!samples[0].first)
244 {
245 samples[0].first = 1;
246 samples[0].second = 0;
247 }
248 // no value at the end? put the cumulative area
249 if (!samples[samples.size() - 1].first)
250 {
251 samples[samples.size() - 1].first = 1;
252 samples[samples.size() - 1].second = cumulativeArea;
253 }
254
255 // fill out the blanks
256 int lastSample = 0;
257 for (int i = 0; i < NumberOfSamples; i++)
258 {
259 if (!samples[i].first)
260 {
261 // not enough vertices in the region for the number of desired
262 // samples. we have to interpolate.
263
264 // first, search for the next valid sample
265 int nextSample = i;
266 for (; nextSample < NumberOfSamples; nextSample++)
267 {
268 if (samples[nextSample].first)
269 break;
270 }
271
272 // next interpolate
273 samples[i].second = samples[lastSample].second +
274 (i - lastSample) * (samples[nextSample].second - samples[lastSample].second) /
275 (nextSample - lastSample);
276 }
277 else
278 lastSample = i;
279 }
280
281 // now prepare the output
282 vtkVariantArray* outputSignature = vtkVariantArray::New();
283 outputSignature->SetNumberOfTuples(static_cast<vtkIdType>(samples.size()));
284 for (unsigned int i = 0; i < samples.size(); i++)
285 {
286 outputSignature->SetValue(i, samples[i].second);
287 }
288 output->Initialize();
289 output->AddColumn(outputSignature);
290 outputSignature->Delete();
291 }
292
293 return 1;
294 }
295 return 0;
296 }
297