1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkConnectivityFilter.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 "vtkConnectivityFilter.h"
16
17 #include "vtkCell.h"
18 #include "vtkCellData.h"
19 #include "vtkDataSet.h"
20 #include "vtkFloatArray.h"
21 #include "vtkIdList.h"
22 #include "vtkInformation.h"
23 #include "vtkInformationVector.h"
24 #include "vtkMath.h"
25 #include "vtkObjectFactory.h"
26 #include "vtkPointData.h"
27 #include "vtkPointData.h"
28 #include "vtkPoints.h"
29 #include "vtkUnstructuredGrid.h"
30 #include "vtkIdTypeArray.h"
31
32 vtkStandardNewMacro(vtkConnectivityFilter);
33
34 // Construct with default extraction mode to extract largest regions.
vtkConnectivityFilter()35 vtkConnectivityFilter::vtkConnectivityFilter()
36 {
37 this->RegionSizes = vtkIdTypeArray::New();
38 this->ExtractionMode = VTK_EXTRACT_LARGEST_REGION;
39 this->ColorRegions = 0;
40
41 this->ScalarConnectivity = 0;
42 this->ScalarRange[0] = 0.0;
43 this->ScalarRange[1] = 1.0;
44
45 this->ClosestPoint[0] = this->ClosestPoint[1] = this->ClosestPoint[2] = 0.0;
46
47 this->CellScalars = vtkFloatArray::New();
48 this->CellScalars->Allocate(8);
49
50 this->NeighborCellPointIds = vtkIdList::New();
51 this->NeighborCellPointIds->Allocate(8);
52
53 this->Seeds = vtkIdList::New();
54 this->SpecifiedRegionIds = vtkIdList::New();
55
56 this->NewScalars = 0;
57 this->NewCellScalars = 0;
58
59 this->OutputPointsPrecision = vtkAlgorithm::DEFAULT_PRECISION;
60 }
61
~vtkConnectivityFilter()62 vtkConnectivityFilter::~vtkConnectivityFilter()
63 {
64 this->RegionSizes->Delete();
65 this->CellScalars->Delete();
66 this->NeighborCellPointIds->Delete();
67 this->Seeds->Delete();
68 this->SpecifiedRegionIds->Delete();
69 }
70
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)71 int vtkConnectivityFilter::RequestData(
72 vtkInformation *vtkNotUsed(request),
73 vtkInformationVector **inputVector,
74 vtkInformationVector *outputVector)
75 {
76 // get the info objects
77 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
78 vtkInformation *outInfo = outputVector->GetInformationObject(0);
79
80 // get the input and output
81 vtkDataSet *input = vtkDataSet::SafeDownCast(
82 inInfo->Get(vtkDataObject::DATA_OBJECT()));
83 vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
84 outInfo->Get(vtkDataObject::DATA_OBJECT()));
85
86 vtkIdType numPts, numCells, cellId, newCellId, i, j, pt;
87 vtkPoints *newPts;
88 int id;
89 int maxCellsInRegion;
90 int largestRegionId = 0;
91 vtkPointData *pd=input->GetPointData(), *outputPD=output->GetPointData();
92 vtkCellData *cd=input->GetCellData(), *outputCD=output->GetCellData();
93
94 vtkDebugMacro(<<"Executing connectivity filter.");
95
96 // Check input/allocate storage
97 //
98 numCells=input->GetNumberOfCells();
99 if ( (numPts=input->GetNumberOfPoints()) < 1 || numCells < 1 )
100 {
101 vtkDebugMacro(<<"No data to connect!");
102 return 1;
103 }
104 output->Allocate(numCells,numCells);
105
106 // See whether to consider scalar connectivity
107 //
108 this->InScalars = input->GetPointData()->GetScalars();
109 if ( !this->ScalarConnectivity )
110 {
111 this->InScalars = NULL;
112 }
113 else
114 {
115 if ( this->ScalarRange[1] < this->ScalarRange[0] )
116 {
117 this->ScalarRange[1] = this->ScalarRange[0];
118 }
119 }
120
121 // Initialize. Keep track of points and cells visited.
122 //
123 this->RegionSizes->Reset();
124 this->Visited = new vtkIdType[numCells];
125 for ( i=0; i < numCells; i++ )
126 {
127 this->Visited[i] = -1;
128 }
129 this->PointMap = new vtkIdType[numPts];
130 for ( i=0; i < numPts; i++ )
131 {
132 this->PointMap[i] = -1;
133 }
134
135 this->NewScalars = vtkIdTypeArray::New();
136 this->NewScalars->SetName("RegionId");
137 this->NewScalars->SetNumberOfTuples(numPts);
138
139 this->NewCellScalars = vtkIdTypeArray::New();
140 this->NewCellScalars->SetName("RegionId");
141 this->NewCellScalars->SetNumberOfTuples(numCells);
142
143 newPts = vtkPoints::New();
144
145 // Set the desired precision for the points in the output.
146 if(this->OutputPointsPrecision == vtkAlgorithm::DEFAULT_PRECISION)
147 {
148 vtkPointSet *inputPointSet = vtkPointSet::SafeDownCast(input);
149 if(inputPointSet)
150 {
151 newPts->SetDataType(inputPointSet->GetPoints()->GetDataType());
152 }
153 else
154 {
155 newPts->SetDataType(VTK_FLOAT);
156 }
157 }
158 else if(this->OutputPointsPrecision == vtkAlgorithm::SINGLE_PRECISION)
159 {
160 newPts->SetDataType(VTK_FLOAT);
161 }
162 else if(this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
163 {
164 newPts->SetDataType(VTK_DOUBLE);
165 }
166
167 newPts->Allocate(numPts);
168
169 // Traverse all cells marking those visited. Each new search
170 // starts a new connected region. Connected region grows
171 // using a connected wave propagation.
172 //
173 this->Wave = vtkIdList::New();
174 this->Wave->Allocate(numPts/4+1,numPts);
175 this->Wave2 = vtkIdList::New();
176 this->Wave2->Allocate(numPts/4+1,numPts);
177
178 this->PointNumber = 0;
179 this->RegionNumber = 0;
180 maxCellsInRegion = 0;
181
182 this->CellIds = vtkIdList::New();
183 this->CellIds->Allocate(8, VTK_CELL_SIZE);
184 this->PointIds = vtkIdList::New();
185 this->PointIds->Allocate(8, VTK_CELL_SIZE);
186
187 if ( this->ExtractionMode != VTK_EXTRACT_POINT_SEEDED_REGIONS &&
188 this->ExtractionMode != VTK_EXTRACT_CELL_SEEDED_REGIONS &&
189 this->ExtractionMode != VTK_EXTRACT_CLOSEST_POINT_REGION )
190 { //visit all cells marking with region number
191 for (cellId=0; cellId < numCells; cellId++)
192 {
193 if ( cellId && !(cellId % 5000) )
194 {
195 this->UpdateProgress (0.1 + 0.8*cellId/numCells);
196 }
197
198 if ( this->Visited[cellId] < 0 )
199 {
200 this->NumCellsInRegion = 0;
201 this->Wave->InsertNextId(cellId);
202 this->TraverseAndMark (input);
203
204 if ( this->NumCellsInRegion > maxCellsInRegion )
205 {
206 maxCellsInRegion = this->NumCellsInRegion;
207 largestRegionId = this->RegionNumber;
208 }
209
210 this->RegionSizes->InsertValue(this->RegionNumber++,
211 this->NumCellsInRegion);
212 this->Wave->Reset();
213 this->Wave2->Reset();
214 }
215 }
216 }
217 else // regions have been seeded, everything considered in same region
218 {
219 this->NumCellsInRegion = 0;
220
221 if ( this->ExtractionMode == VTK_EXTRACT_POINT_SEEDED_REGIONS )
222 {
223 for (i=0; i < this->Seeds->GetNumberOfIds(); i++)
224 {
225 pt = this->Seeds->GetId(i);
226 if ( pt >= 0 )
227 {
228 input->GetPointCells(pt,this->CellIds);
229 for (j=0; j < this->CellIds->GetNumberOfIds(); j++)
230 {
231 this->Wave->InsertNextId(this->CellIds->GetId(j));
232 }
233 }
234 }
235 }
236 else if ( this->ExtractionMode == VTK_EXTRACT_CELL_SEEDED_REGIONS )
237 {
238 for (i=0; i < this->Seeds->GetNumberOfIds(); i++)
239 {
240 cellId = this->Seeds->GetId(i);
241 if ( cellId >= 0 )
242 {
243 this->Wave->InsertNextId(cellId);
244 }
245 }
246 }
247 else if ( this->ExtractionMode == VTK_EXTRACT_CLOSEST_POINT_REGION )
248 {//loop over points, find closest one
249 double minDist2, dist2, x[3];
250 vtkIdType minId = 0;
251 for (minDist2=VTK_DOUBLE_MAX, i=0; i<numPts; i++)
252 {
253 input->GetPoint(i,x);
254 dist2 = vtkMath::Distance2BetweenPoints(x,this->ClosestPoint);
255 if ( dist2 < minDist2 )
256 {
257 minId = i;
258 minDist2 = dist2;
259 }
260 }
261 input->GetPointCells(minId,this->CellIds);
262 for (j=0; j < this->CellIds->GetNumberOfIds(); j++)
263 {
264 this->Wave->InsertNextId(this->CellIds->GetId(j));
265 }
266 }
267 this->UpdateProgress (0.5);
268
269 //mark all seeded regions
270 this->TraverseAndMark (input);
271 this->RegionSizes->InsertValue(this->RegionNumber,this->NumCellsInRegion);
272 this->UpdateProgress (0.9);
273 }
274
275 vtkDebugMacro (<<"Extracted " << this->RegionNumber << " region(s)");
276 this->Wave->Delete();
277 this->Wave2->Delete();
278
279 // Now that points and cells have been marked, traverse these lists pulling
280 // everything that has been visited.
281 //
282 //Pass through point data that has been visited
283 outputPD->CopyAllocate(pd);
284 outputCD->CopyAllocate(cd);
285
286 for (i=0; i < numPts; i++)
287 {
288 if ( this->PointMap[i] > -1 )
289 {
290 newPts->InsertPoint(this->PointMap[i],input->GetPoint(i));
291 outputPD->CopyData(pd,i,this->PointMap[i]);
292 }
293 }
294
295 // if coloring regions; send down new scalar data
296 if ( this->ColorRegions )
297 {
298 int idx = outputPD->AddArray(this->NewScalars);
299 outputPD->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
300 idx = outputCD->AddArray(this->NewCellScalars);
301 outputCD->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
302 }
303 this->NewScalars->Delete();
304 this->NewCellScalars->Delete();
305
306 output->SetPoints(newPts);
307 newPts->Delete();
308
309 // Create output cells
310 //
311 if ( this->ExtractionMode == VTK_EXTRACT_POINT_SEEDED_REGIONS ||
312 this->ExtractionMode == VTK_EXTRACT_CELL_SEEDED_REGIONS ||
313 this->ExtractionMode == VTK_EXTRACT_CLOSEST_POINT_REGION ||
314 this->ExtractionMode == VTK_EXTRACT_ALL_REGIONS)
315 { // extract any cell that's been visited
316 for (cellId=0; cellId < numCells; cellId++)
317 {
318 if ( this->Visited[cellId] >= 0 )
319 {
320 // special handling for polyhedron cells
321 if (vtkUnstructuredGrid::SafeDownCast(input) &&
322 input->GetCellType(cellId) == VTK_POLYHEDRON)
323 {
324 vtkUnstructuredGrid::SafeDownCast(input)->
325 GetFaceStream(cellId, this->PointIds);
326 vtkUnstructuredGrid::ConvertFaceStreamPointIds(this->PointIds,
327 this->PointMap);
328 }
329 else
330 {
331 input->GetCellPoints(cellId, this->PointIds);
332 for (i=0; i < this->PointIds->GetNumberOfIds(); i++)
333 {
334 id = this->PointMap[this->PointIds->GetId(i)];
335 this->PointIds->InsertId(i,id);
336 }
337 }
338 newCellId = output->InsertNextCell(input->GetCellType(cellId),
339 this->PointIds);
340 outputCD->CopyData(cd,cellId,newCellId);
341 }
342 }
343 }
344 else if ( this->ExtractionMode == VTK_EXTRACT_SPECIFIED_REGIONS )
345 {
346 for (cellId=0; cellId < numCells; cellId++)
347 {
348 int inReg, regionId;
349 if ( (regionId=this->Visited[cellId]) >= 0 )
350 {
351 for (inReg=0,i=0; i<this->SpecifiedRegionIds->GetNumberOfIds(); i++)
352 {
353 if ( regionId == this->SpecifiedRegionIds->GetId(i) )
354 {
355 inReg = 1;
356 break;
357 }
358 }
359 if ( inReg )
360 {
361 // special handling for polyhedron cells
362 if (vtkUnstructuredGrid::SafeDownCast(input) &&
363 input->GetCellType(cellId) == VTK_POLYHEDRON)
364 {
365 vtkUnstructuredGrid::SafeDownCast(input)->
366 GetFaceStream(cellId, this->PointIds);
367 vtkUnstructuredGrid::ConvertFaceStreamPointIds(this->PointIds,
368 this->PointMap);
369 }
370 else
371 {
372 input->GetCellPoints(cellId, this->PointIds);
373 for (i=0; i < this->PointIds->GetNumberOfIds(); i++)
374 {
375 id = this->PointMap[this->PointIds->GetId(i)];
376 this->PointIds->InsertId(i,id);
377 }
378 }
379 newCellId = output->InsertNextCell(input->GetCellType(cellId),
380 this->PointIds);
381 outputCD->CopyData(cd,cellId,newCellId);
382 }
383 }
384 }
385 }
386 else //extract largest region
387 {
388 for (cellId=0; cellId < numCells; cellId++)
389 {
390 if ( this->Visited[cellId] == largestRegionId )
391 {
392 // special handling for polyhedron cells
393 if (vtkUnstructuredGrid::SafeDownCast(input) &&
394 input->GetCellType(cellId) == VTK_POLYHEDRON)
395 {
396 vtkUnstructuredGrid::SafeDownCast(input)->
397 GetFaceStream(cellId, this->PointIds);
398 vtkUnstructuredGrid::ConvertFaceStreamPointIds(this->PointIds,
399 this->PointMap);
400 }
401 else
402 {
403 input->GetCellPoints(cellId, this->PointIds);
404 for (i=0; i < this->PointIds->GetNumberOfIds(); i++)
405 {
406 id = this->PointMap[this->PointIds->GetId(i)];
407 this->PointIds->InsertId(i,id);
408 }
409 }
410 newCellId = output->InsertNextCell(input->GetCellType(cellId),
411 this->PointIds);
412 outputCD->CopyData(cd,cellId,newCellId);
413 }
414 }
415 }
416
417 delete [] this->Visited;
418 delete [] this->PointMap;
419 this->PointIds->Delete();
420 this->CellIds->Delete();
421 output->Squeeze();
422 vtkDataArray* outScalars = 0;
423 if (this->ColorRegions && (outScalars=output->GetPointData()->GetScalars()))
424 {
425 outScalars->Resize(output->GetNumberOfPoints());
426 }
427
428 int num = this->GetNumberOfExtractedRegions();
429 int count = 0;
430
431 for (int ii = 0; ii < num; ii++)
432 {
433 count += this->RegionSizes->GetValue (ii);
434 }
435 vtkDebugMacro (<< "Total # of cells accounted for: " << count);
436 vtkDebugMacro (<< "Extracted " << output->GetNumberOfCells() << " cells");
437
438 return 1;
439 }
440
441
442 // Mark current cell as visited and assign region number. Note:
443 // traversal occurs across shared vertices.
444 //
TraverseAndMark(vtkDataSet * input)445 void vtkConnectivityFilter::TraverseAndMark (vtkDataSet *input)
446 {
447 vtkIdType i, j, k, cellId, numIds, ptId, numPts, numCells;
448 vtkIdList *tmpWave;
449
450 while ( (numIds=this->Wave->GetNumberOfIds()) > 0 )
451 {
452 for ( i=0; i < numIds; i++ )
453 {
454 cellId = this->Wave->GetId(i);
455 if ( this->Visited[cellId] < 0 )
456 {
457 this->NewCellScalars->SetValue(cellId, this->RegionNumber);
458 this->Visited[cellId] = this->RegionNumber;
459 this->NumCellsInRegion++;
460 input->GetCellPoints(cellId, this->PointIds);
461
462 numPts = this->PointIds->GetNumberOfIds();
463 for (j=0; j < numPts; j++)
464 {
465 if ( this->PointMap[ptId=this->PointIds->GetId(j)] < 0 )
466 {
467 this->PointMap[ptId] = this->PointNumber++;
468 this->NewScalars->SetValue(this->PointMap[ptId],
469 this->RegionNumber);
470 }
471
472 input->GetPointCells(ptId,this->CellIds);
473
474 // check connectivity criterion (geometric + scalar)
475 numCells = this->CellIds->GetNumberOfIds();
476 for (k=0; k < numCells; k++)
477 {
478 cellId = this->CellIds->GetId(k);
479 if ( this->InScalars )
480 {
481 int numScalars, ii;
482 double s, range[2];
483
484 input->GetCellPoints(cellId, this->NeighborCellPointIds);
485 numScalars = this->NeighborCellPointIds->GetNumberOfIds();
486 this->CellScalars->SetNumberOfComponents(this->InScalars->GetNumberOfComponents());
487 this->CellScalars->SetNumberOfTuples(numScalars);
488 this->InScalars->GetTuples(this->NeighborCellPointIds,
489 this->CellScalars);
490 range[0] = VTK_DOUBLE_MAX; range[1] = -VTK_DOUBLE_MAX;
491 for (ii=0; ii < numScalars; ii++)
492 {
493 s = this->CellScalars->GetComponent(ii,0);
494 if ( s < range[0] )
495 {
496 range[0] = s;
497 }
498 if ( s > range[1] )
499 {
500 range[1] = s;
501 }
502 }
503 if ( range[1] >= this->ScalarRange[0] &&
504 range[0] <= this->ScalarRange[1] )
505 {
506 this->Wave2->InsertNextId(cellId);
507 }
508 }
509 else
510 {
511 this->Wave2->InsertNextId(cellId);
512 }
513 }//for all cells using this point
514 }//for all points of this cell
515 }//if cell not yet visited
516 }//for all cells in this wave
517
518 tmpWave = this->Wave;
519 this->Wave = this->Wave2;
520 this->Wave2 = tmpWave;
521 tmpWave->Reset();
522 } //while wave is not empty
523
524 return;
525 }
526
527 // Obtain the number of connected regions.
GetNumberOfExtractedRegions()528 int vtkConnectivityFilter::GetNumberOfExtractedRegions()
529 {
530 return this->RegionSizes->GetMaxId() + 1;
531 }
532
533 // Initialize list of point ids/cell ids used to seed regions.
InitializeSeedList()534 void vtkConnectivityFilter::InitializeSeedList()
535 {
536 this->Modified();
537 this->Seeds->Reset();
538 }
539
540 // Add a seed id (point or cell id). Note: ids are 0-offset.
AddSeed(vtkIdType id)541 void vtkConnectivityFilter::AddSeed(vtkIdType id)
542 {
543 this->Modified();
544 this->Seeds->InsertNextId(id);
545 }
546
547 // Delete a seed id (point or cell id). Note: ids are 0-offset.
DeleteSeed(vtkIdType id)548 void vtkConnectivityFilter::DeleteSeed(vtkIdType id)
549 {
550 this->Modified();
551 this->Seeds->DeleteId(id);
552 }
553
554 // Initialize list of region ids to extract.
InitializeSpecifiedRegionList()555 void vtkConnectivityFilter::InitializeSpecifiedRegionList()
556 {
557 this->Modified();
558 this->SpecifiedRegionIds->Reset();
559 }
560
561 // Add a region id to extract. Note: ids are 0-offset.
AddSpecifiedRegion(int id)562 void vtkConnectivityFilter::AddSpecifiedRegion(int id)
563 {
564 this->Modified();
565 this->SpecifiedRegionIds->InsertNextId(id);
566 }
567
568 // Delete a region id to extract. Note: ids are 0-offset.
DeleteSpecifiedRegion(int id)569 void vtkConnectivityFilter::DeleteSpecifiedRegion(int id)
570 {
571 this->Modified();
572 this->SpecifiedRegionIds->DeleteId(id);
573 }
574
FillInputPortInformation(int,vtkInformation * info)575 int vtkConnectivityFilter::FillInputPortInformation(int, vtkInformation *info)
576 {
577 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
578 return 1;
579 }
580
PrintSelf(ostream & os,vtkIndent indent)581 void vtkConnectivityFilter::PrintSelf(ostream& os, vtkIndent indent)
582 {
583 this->Superclass::PrintSelf(os,indent);
584
585 os << indent << "Extraction Mode: ";
586 os << this->GetExtractionModeAsString() << "\n";
587
588 os << indent << "Closest Point: (" << this->ClosestPoint[0] << ", "
589 << this->ClosestPoint[1] << ", " << this->ClosestPoint[2] << ")\n";
590
591 os << indent << "Color Regions: " << (this->ColorRegions ? "On\n" : "Off\n");
592
593 os << indent << "Scalar Connectivity: "
594 << (this->ScalarConnectivity ? "On\n" : "Off\n");
595
596 double *range = this->GetScalarRange();
597 os << indent << "Scalar Range: (" << range[0] << ", " << range[1] << ")\n";
598 os << indent << "Output Points Precision: " << this->OutputPointsPrecision
599 << "\n";
600 }
601
602