1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2013-2015 Imperial College London
5 * Copyright 2013-2015 Andreas Schuh
6 *
7 * Licensed under the Apache License, Version 2.0 (the "License");
8 * you may not use this file except in compliance with the License.
9 * You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing, software
14 * distributed under the License is distributed on an "AS IS" BASIS,
15 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 * See the License for the specific language governing permissions and
17 * limitations under the License.
18 */
19
20 #include "mirtk/Common.h"
21 #include "mirtk/Options.h"
22
23 #include "mirtk/IOConfig.h"
24
25 #include "mirtk/BaseImage.h"
26 #include "mirtk/ImageReader.h"
27
28 #ifdef HAVE_MIRTK_Transformation
29 # include "mirtk/Transformation.h"
30 # include "mirtk/AffineTransformation.h"
31 #endif // HAVE_MIRTK_Transformation
32
33 #ifdef HAVE_MIRTK_PointSet
34 # include "mirtk/EdgeTable.h"
35 # include "mirtk/Triangle.h"
36 # include "mirtk/Polyhedron.h"
37 # include "mirtk/PointSetIO.h"
38 # include "mirtk/PointSetUtils.h"
39 # include "mirtk/DataStatistics.h"
40
41 # include "mirtk/Vtk.h"
42 # include "mirtk/VtkMath.h"
43
44 # include "vtkSmartPointer.h"
45 # include "vtkPolyData.h"
46 # include "vtkDataArray.h"
47 # include "vtkCellArray.h"
48 # include "vtkPointData.h"
49 # include "vtkCellData.h"
50 # include "vtkGenericCell.h"
51 # include "vtkOctreePointLocator.h"
52 # include "vtkUnsignedCharArray.h"
53 # include "vtkFloatArray.h"
54 # include "vtkDataSetSurfaceFilter.h"
55 # include "vtkCleanPolyData.h"
56 # include "vtkPolyDataConnectivityFilter.h"
57 # include "vtkIntersectionPolyDataFilter.h"
58 #endif // HAVE_MIRTK_PointSet
59
60 using namespace mirtk;
61
62
63
64 // =============================================================================
65 // Help
66 // =============================================================================
67
68 // -----------------------------------------------------------------------------
PrintHelp(const char * name)69 void PrintHelp(const char *name)
70 {
71 cout << "\n";
72 cout << "Usage: " << name << " <file> [options]\n";
73 cout << "\n";
74 cout << "Description:\n";
75 cout << " Prints some useful information about the given input file, which\n";
76 cout << " can be an image, transformation (requires MIRTK Transformation module),\n";
77 cout << " or point set (requires MIRTK Point Set module).\n";
78 cout << "\n";
79 cout << "Arguments:\n";
80 cout << " input Input image";
81 #ifdef HAVE_MIRTK_Transformation
82 cout << " or transformation";
83 #endif
84 #ifdef HAVE_MIRTK_PointSet
85 cout << " or point set";
86 #endif
87 cout << " file.\n";
88 cout << "\n";
89 cout << "Image options:\n";
90 // Note: No two first named options can be identical, otherwise Sphinx
91 // complains that the -attributes point set option below is a duplicate.
92 // (cf. help-rst.py output).
93 cout << " -a -attributes Print attributes of image.\n";
94 #ifdef HAVE_MIRTK_Transformation
95 cout << "\n";
96 cout << "Transformation options:\n";
97 cout << " -attributes Print attributes of transformation. (default)\n";
98 cout << " -type-name Print name of transformation type.\n";
99 cout << " -type-id Print enum of transformation type.\n";
100 #endif // HAVE_MIRTK_Transformation
101 #ifdef HAVE_MIRTK_PointSet
102 cout << "\n";
103 cout << "Point set options:\n";
104 cout << " -point <id>... Print info of points with the given zero-based indices.\n";
105 cout << " -cell-types Report cell types and how many of each. (default: off)\n";
106 cout << " -bounds Report point set bounds and center point. (default: off)\n";
107 cout << " -edgelength Report edge length statistics. (default: off)\n";
108 cout << " -self-intersections Check for self-intersections. (default: off)\n";
109 cout << " -area Display surface area information\n";
110 cout << " -output-surface <file> Write surface mesh to specified file. (default: none)\n";
111 #endif // HAVE_MIRTK_PointSet
112 PrintStandardOptions(cout);
113 cout << "\n";
114 cout.flush();
115 }
116
117 // =============================================================================
118 // Point set auxiliaries
119 // =============================================================================
120 #ifdef HAVE_MIRTK_PointSet
121
122 // -----------------------------------------------------------------------------
123 /// Convert input dataset to polygonal data, i.e., boundary surface
ConvertToPolyData(vtkDataSet * dataset)124 vtkSmartPointer<vtkPolyData> ConvertToPolyData(vtkDataSet *dataset)
125 {
126 vtkSmartPointer<vtkPolyData> polydata = vtkPolyData::SafeDownCast(dataset);
127 if (polydata) return polydata;
128 vtkSmartPointer<vtkDataSetSurfaceFilter> surface;
129 surface = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
130 SetVTKInput(surface, dataset);
131 surface->PassThroughCellIdsOff();
132 surface->PassThroughPointIdsOff();
133 surface->Update();
134 return surface->GetOutput();
135 }
136
137 // -----------------------------------------------------------------------------
138 /// Clean polydata
CleanPolyData(vtkSmartPointer<vtkPolyData> polydata)139 vtkSmartPointer<vtkPolyData> CleanPolyData(vtkSmartPointer<vtkPolyData> polydata)
140 {
141 vtkSmartPointer<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
142 cleaner->SetAbsoluteTolerance(.0);
143 cleaner->PointMergingOn();
144 cleaner->ConvertLinesToPointsOn();
145 cleaner->ConvertPolysToLinesOn();
146 cleaner->ConvertStripsToPolysOn();
147 SetVTKInput(cleaner, polydata);
148 cleaner->Update();
149 return cleaner->GetOutput();
150 }
151
152 // -----------------------------------------------------------------------------
153 /// Get center point and radius of bounding sphere
GetBoundingSphereRadius(vtkGenericCell * cell,double c[3])154 inline double GetBoundingSphereRadius(vtkGenericCell *cell, double c[3])
155 {
156 double p[3];
157 // Get center of bounding sphere
158 c[0] = .0, c[1] = .0, c[2] = .0;
159 double *weights = new double[cell->GetNumberOfPoints()];
160 int subId = cell->GetParametricCenter(p);
161 cell->EvaluateLocation(subId, p, c, weights);
162 delete[] weights;
163 // Get radius of bounding sphere
164 double r = .0;
165 vtkPoints *points = cell->GetPoints();
166 for (vtkIdType i = 0; i < points->GetNumberOfPoints(); ++i) {
167 points->GetPoint(i, p);
168 r = max(r, vtkMath::Distance2BetweenPoints(c, p));
169 }
170 return sqrt(r);
171 }
172
173 // -----------------------------------------------------------------------------
174 /// Auxiliary functor which counts the number of self-intersections of a triangulated surface mesh
175 struct CountTriangleTriangleIntersections
176 {
177 vtkPolyData *_DataSet;
178 vtkAbstractPointLocator *_PointLocator;
179 vtkDataArray *_Mask;
180 int _NumberOfIntersections;
181
CountTriangleTriangleIntersectionsCountTriangleTriangleIntersections182 CountTriangleTriangleIntersections() : _NumberOfIntersections(0) {}
183
CountTriangleTriangleIntersectionsCountTriangleTriangleIntersections184 CountTriangleTriangleIntersections(const CountTriangleTriangleIntersections &other, split)
185 :
186 _DataSet(other._DataSet),
187 _PointLocator(other._PointLocator),
188 _Mask(other._Mask),
189 _NumberOfIntersections(0)
190 {}
191
joinCountTriangleTriangleIntersections192 void join(const CountTriangleTriangleIntersections &other)
193 {
194 _NumberOfIntersections += other._NumberOfIntersections;
195 }
196
operator ()CountTriangleTriangleIntersections197 void operator ()(const blocked_range<vtkIdType> &re)
198 {
199 vtkSmartPointer<vtkGenericCell> cell1 = vtkSmartPointer<vtkGenericCell>::New();
200 vtkSmartPointer<vtkGenericCell> cell2 = vtkSmartPointer<vtkGenericCell>::New();
201 vtkSmartPointer<vtkIdList> ptIds = vtkSmartPointer<vtkIdList>::New();
202 vtkSmartPointer<vtkIdList> cellIds = vtkSmartPointer<vtkIdList>::New();
203
204 vtkIdType cellId2;
205 vtkNew<vtkIdList> ptCellIds;
206 double r1, c1[3], v1[3], v2[3], v3[3], u1[3], u2[3], u3[3];
207
208 for (vtkIdType cellId1 = re.begin(); cellId1 != re.end(); ++cellId1) {
209 _DataSet->GetCell(cellId1, cell1);
210 if (cell1->GetNumberOfPoints() != 3) continue;
211 _DataSet->GetPoint(cell1->GetPointId(0), v1);
212 _DataSet->GetPoint(cell1->GetPointId(1), v2);
213 _DataSet->GetPoint(cell1->GetPointId(2), v3);
214 r1 = GetBoundingSphereRadius(cell1, c1);
215 _PointLocator->FindPointsWithinRadius(3.0 * r1, c1, ptIds);
216 cellIds->Reset();
217 for (vtkIdType i = 0; i < ptIds->GetNumberOfIds(); ++i) {
218 _DataSet->GetPointCells(ptIds->GetId(i), ptCellIds.GetPointer());
219 for (vtkIdType j = 0; j < ptCellIds->GetNumberOfIds(); ++j) {
220 cellIds->InsertUniqueId(ptCellIds->GetId(j));
221 }
222 }
223 for (vtkIdType i = 0; i < cellIds->GetNumberOfIds(); ++i) {
224 cellId2 = cellIds->GetId(i);
225 if (cellId2 <= cellId1) continue; // unordered pairs of cells
226 _DataSet->GetCell(cellId2, cell2);
227 if (cell2->GetNumberOfPoints() != 3) continue;
228 if (cell1->GetPointIds()->IsId(cell2->GetPointId(0)) >= 0) continue;
229 if (cell1->GetPointIds()->IsId(cell2->GetPointId(1)) >= 0) continue;
230 if (cell1->GetPointIds()->IsId(cell2->GetPointId(2)) >= 0) continue;
231 _DataSet->GetPoint(cell2->GetPointId(0), u1);
232 _DataSet->GetPoint(cell2->GetPointId(1), u2);
233 _DataSet->GetPoint(cell2->GetPointId(2), u3);
234 if (Triangle::TriangleTriangleIntersection(v1, v2, v3, u1, u2, u3) != 0) {
235 _Mask->SetComponent(cellId1, 0, 1.0);
236 _Mask->SetComponent(cellId2, 0, 1.0);
237 ++_NumberOfIntersections;
238 }
239 }
240 }
241 }
242 };
243
244 // -----------------------------------------------------------------------------
245 /// Count number of self-intersections of triangulated surface mesh
NumberOfTriangleTriangleIntersections(vtkPolyData * polydata,const char * array_name=nullptr)246 int NumberOfTriangleTriangleIntersections(vtkPolyData *polydata, const char *array_name = nullptr)
247 {
248 vtkSmartPointer<vtkDataArray> mask = vtkSmartPointer<vtkUnsignedCharArray>::New();
249 mask->SetNumberOfComponents(1);
250 mask->SetNumberOfTuples(polydata->GetNumberOfCells());
251 mask->FillComponent(0, 0.);
252 vtkNew<vtkOctreePointLocator> octree;
253 octree->SetDataSet(polydata);
254 octree->BuildLocator();
255 CountTriangleTriangleIntersections count;
256 count._DataSet = polydata;
257 count._PointLocator = octree.GetPointer();
258 count._Mask = mask;
259 parallel_reduce(blocked_range<vtkIdType>(0, polydata->GetNumberOfCells()), count);
260 if (array_name) {
261 mask->SetName(array_name);
262 polydata->GetCellData()->AddArray(mask);
263 }
264 return count._NumberOfIntersections;
265 }
266
267 // -----------------------------------------------------------------------------
268 /// Auxiliary functor which counts the number of vertices inside the polyhedron
269 struct CountVerticesInsidePolyhedron
270 {
271 Polyhedron *_Polyhedron;
272 vtkDataArray *_Mask;
273 int _Num;
274
CountVerticesInsidePolyhedronCountVerticesInsidePolyhedron275 CountVerticesInsidePolyhedron() : _Num(0) {}
CountVerticesInsidePolyhedronCountVerticesInsidePolyhedron276 CountVerticesInsidePolyhedron(const CountVerticesInsidePolyhedron &other, split)
277 :
278 _Polyhedron(other._Polyhedron), _Mask(other._Mask), _Num(0)
279 {}
280
joinCountVerticesInsidePolyhedron281 void join(const CountVerticesInsidePolyhedron &other)
282 {
283 _Num += other._Num;
284 }
285
operator ()CountVerticesInsidePolyhedron286 void operator ()(const blocked_range<int> &re)
287 {
288 double p[3];
289 for (int ptId = re.begin(); ptId != re.end(); ++ptId) {
290 _Polyhedron->GetPoint(ptId, p);
291 if (_Polyhedron->IsInside(p)) {
292 _Mask->SetComponent(ptId, 0, 1);
293 ++_Num;
294 } else {
295 _Mask->SetComponent(ptId, 0, 0);
296 }
297 }
298 }
299 };
300
301 // -----------------------------------------------------------------------------
302 /// Point-inside-polyhedron test
NumberOfVerticesInsidePolyhedron(vtkPolyData * polydata)303 int NumberOfVerticesInsidePolyhedron(vtkPolyData *polydata)
304 {
305 vtkSmartPointer<vtkDataArray> mask = vtkSmartPointer<vtkUnsignedCharArray>::New();
306 mask->SetName("VertexInsidePolyhedron");
307 mask->SetNumberOfComponents(1);
308 mask->SetNumberOfTuples(polydata->GetNumberOfPoints());
309 polydata->GetPointData()->AddArray(mask);
310 Polyhedron polyhedron(polydata);
311 CountVerticesInsidePolyhedron count;
312 count._Polyhedron = &polyhedron;
313 count._Mask = mask;
314 parallel_reduce(blocked_range<int>(0, polyhedron.NumberOfPoints()), count);
315 return count._Num;
316 }
317
318 // -----------------------------------------------------------------------------
319 /// Compute volume of cells
320 struct ComputeCellVolumes
321 {
322 vtkPointSet *_PointSet;
323 vtkDataArray *_Volume;
324
operator ()ComputeCellVolumes325 void operator ()(const blocked_range<vtkIdType> &re) const
326 {
327 double volume;
328 vtkSmartPointer<vtkGenericCell> cell = vtkSmartPointer<vtkGenericCell>::New();
329 for (vtkIdType cellId = re.begin(); cellId != re.end(); ++cellId) {
330 _PointSet->GetCell(cellId, cell);
331 volume = ComputeVolume(cell);
332 _Volume->SetComponent(cellId, 0, volume);
333 }
334 }
335
RunComputeCellVolumes336 static vtkSmartPointer<vtkDataArray> Run(vtkPointSet *pointset)
337 {
338 vtkSmartPointer<vtkDataArray> volume;
339 volume = vtkSmartPointer<vtkFloatArray>::New();
340 volume->SetName("Volume");
341 volume->SetNumberOfComponents(1);
342 volume->SetNumberOfTuples(pointset->GetNumberOfCells());
343 pointset->GetCellData()->AddArray(volume);
344 ComputeCellVolumes body;
345 body._PointSet = pointset;
346 body._Volume = volume;
347 parallel_for(blocked_range<vtkIdType>(0, pointset->GetNumberOfCells()), body);
348 return volume;
349 }
350 };
351
352 // -----------------------------------------------------------------------------
DataObjectTypeString(int type)353 const char *DataObjectTypeString(int type)
354 {
355 switch (type) {
356 case VTK_IMAGE_DATA: return "Image";
357 case VTK_POLY_DATA: return "Polydata";
358 case VTK_UNSTRUCTURED_GRID: return "Unstructured grid";
359 case VTK_STRUCTURED_GRID: return "Structured grid";
360 case VTK_STRUCTURED_POINTS: return "Structured points";
361 default: return "unknown";
362 }
363 }
364
365 // -----------------------------------------------------------------------------
CellTypeString(int type)366 const char *CellTypeString(int type)
367 {
368 switch (type) {
369 case VTK_EMPTY_CELL: return "empty cells";
370 case VTK_VERTEX: return "vertices";
371 case VTK_POLY_VERTEX: return "poly-vertices";
372 case VTK_LINE: return "lines";
373 case VTK_TRIANGLE: return "triangles";
374 case VTK_TRIANGLE_STRIP: return "triangle strips";
375 case VTK_POLYGON: return "polygons";
376 case VTK_QUAD: return "quads";
377 case VTK_TETRA: return "tetrahedra";
378 case VTK_VOXEL: return "voxels";
379 case VTK_HEXAHEDRON: return "hexadrons";
380 case VTK_HEXAGONAL_PRISM: return "hexagonal prisms";
381 case VTK_PYRAMID: return "pyramids";
382 default: return "unknown cells";
383 }
384 }
385
386 #endif // HAVE_MIRTK_PointSet
387 // =============================================================================
388 // Main
389 // =============================================================================
390
391 // -----------------------------------------------------------------------------
main(int argc,char * argv[])392 int main(int argc, char *argv[])
393 {
394 EXPECTS_POSARGS(1);
395
396 bool image_info = false, dof_info = false, pointset_info = false;
397
398 // Read input file
399 const char *fname = POSARG(1);
400 const string fext = Extension(fname, EXT_LastWithoutGz);
401
402 UniquePtr<ImageReader> image_reader;
403 #ifdef HAVE_MIRTK_Transformation
404 UniquePtr<Transformation> dof;
405 if (Transformation::CheckHeader(fname)) {
406 dof.reset(Transformation::New(fname));
407 dof_info = true;
408 }
409 else
410 #endif // HAVE_MIRTK_Transformation
411 {
412 InitializeIOLibrary();
413 image_reader.reset(ImageReader::TryNew(fname));
414 image_info = (image_reader != nullptr);
415 }
416 #ifdef HAVE_MIRTK_PointSet
417 vtkSmartPointer<vtkPointSet> pointset;
418 vtkSmartPointer<vtkPolyData> polydata;
419 if (!image_info && !dof_info &&
420 // FIXME: Silence errors of VTK readers instead
421 fext != ".nii" && fext != ".hdr" && fext != ".img" && fext != ".png" &&
422 fext != ".gipl" && fext != ".pgm") {
423 const bool exit_on_failure = false;
424 pointset = ReadPointSet(POSARG(1), exit_on_failure);
425 if (pointset) {
426 polydata = ConvertToPolyData(pointset);
427 polydata->BuildLinks();
428 pointset_info = true;
429 }
430 }
431 #endif // HAVE_MIRTK_PointSet
432
433 if (!image_info && !dof_info && !pointset_info) {
434 FatalError("Unknown input file type or file format not supported!");
435 }
436
437 //////////////////////////////////////////////////////////////////////////////
438 // Image info
439 if (image_info) {
440
441 bool attributes = false;
442
443 for (ALL_OPTIONS) {
444 if (OPTION("-a") || OPTION("-attributes") || OPTION("-attr")) attributes = true;
445 else HANDLE_COMMON_OR_UNKNOWN_OPTION();
446 }
447
448 UniquePtr<BaseImage> image(image_reader->Run());
449 if (!attributes) {
450 cout << "Information from ImageReader::Print\n\n";
451 image_reader->Print();
452 cout << "\nInformation from BaseImage::Print\n\n";
453 }
454 image->Print();
455 if (!attributes) {
456 double min, max;
457 image->GetMinMaxAsDouble(&min, &max);
458 cout << "\nMinimum intensity: " << min << "\n";
459 cout << "Maximum intensity: " << max << "\n";
460 cout << "\nImage to world matrix\n";
461 image->GetImageToWorldMatrix().Print();
462 cout << "\nWorld to image matrix\n";
463 image->GetWorldToImageMatrix().Print();
464 if (!image->GetAffineMatrix().IsIdentity()) {
465 cout << "Including affine transformation\n";
466 #ifdef HAVE_MIRTK_Transformation
467 AffineTransformation transformation;
468 transformation.PutMatrix(image->GetAffineMatrix());
469 transformation.Print(Indent(1));
470 #else // HAVE_MIRTK_Transformation
471 image->GetAffineMatrix().Print();
472 #endif // HAVE_MIRTK_Transformation
473 }
474 }
475
476 } // image_info
477
478 //////////////////////////////////////////////////////////////////////////////
479 // Transformation info
480 #ifdef HAVE_MIRTK_Transformation
481 if (dof_info) {
482
483 bool attributes = false;
484 bool type_name = false;
485 bool type_id = false;
486
487 for (ALL_OPTIONS) {
488 if (OPTION("-attributes") || OPTION("-attr") || OPTION("-a")) attributes = true;
489 else if (OPTION("-type-name") || OPTION("-type")) type_name = true;
490 else if (OPTION("-type-id") || OPTION("-id")) type_id = true;
491 else HANDLE_COMMON_OR_UNKNOWN_OPTION();
492 }
493 if (!type_name && !type_id && !attributes) {
494 attributes = true;
495 }
496
497 if (type_name) {
498 cout << dof->NameOfClass();
499 }
500 if (type_id) {
501 if (type_name) cout << " (";
502 cout << dof->TypeOfClass();
503 if (type_name) cout << ")";
504 }
505 if (attributes) {
506 if (type_name || type_id) cout << "\n";
507 dof->Print();
508 }
509
510 } // dof_info
511 #endif // HAVE_MIRTK_Transformation
512
513 //////////////////////////////////////////////////////////////////////////////
514 // Point set info
515 #ifdef HAVE_MIRTK_PointSet
516 if (pointset_info) {
517
518 using namespace mirtk::data::statistic;
519
520 const char *self_intersections_name = nullptr;
521 bool report_bounds = false;
522 bool report_cell_types = false;
523 bool report_cell_volumes = false;
524 bool report_edge_lengths = false;
525 bool report_surface_area = false;
526 bool print_summary = false;
527 bool print_surface_summary = false;
528 bool list_pointdata_arrays = false;
529 bool list_celldata_arrays = false;
530
531 const char *output_pointset_name = NULL;
532 const char *output_surface_name = NULL;
533
534 Array<int> point_indices;
535 int max_point_index = -1;
536
537 for (ALL_OPTIONS) {
538 if (OPTION("-self-intersections")) {
539 if (HAS_ARGUMENT) self_intersections_name = ARGUMENT;
540 else self_intersections_name = "SelfIntersectionsMask";
541 }
542 else if (OPTION("-edgelength")) report_edge_lengths = true;
543 else if (OPTION("-bounds")) report_bounds = true;
544 else if (OPTION("-cell-types")) report_cell_types = true;
545 else if (OPTION("-pointdata")) list_pointdata_arrays = true;
546 else if (OPTION("-celldata")) list_celldata_arrays = true;
547 else if (OPTION("-data") || OPTION("-arrays") || OPTION("-data-arrays")) {
548 list_pointdata_arrays = list_celldata_arrays = true;
549 }
550 else if (OPTION("-summary")) print_summary = true;
551 else if (OPTION("-surface-graph") || OPTION("-surface-mesh") || OPTION("-surface-summary") || OPTION("-surface")) {
552 print_surface_summary = true;
553 }
554 else if (OPTION("-area") || OPTION("-surface-area")) {
555 report_surface_area = true;
556 }
557 else if (OPTION("-point")) {
558 int i;
559 do {
560 PARSE_ARGUMENT(i);
561 if (i < 0) FatalError("Invalid -point ID, must be non-negative: " << i);
562 point_indices.push_back(i);
563 if (i > max_point_index) max_point_index = i;
564 } while (HAS_ARGUMENT);
565 }
566 else if (OPTION("-vol") || OPTION("-volume") || OPTION("-volumes")) {
567 report_cell_volumes = true;
568 }
569 else if (OPTION("-o") || OPTION("-out") || OPTION("-output")) {
570 output_pointset_name = ARGUMENT;
571 }
572 else if (OPTION("-s") || OPTION("-output-surface")) {
573 output_surface_name = ARGUMENT;
574 }
575 else HANDLE_COMMON_OR_UNKNOWN_OPTION();
576 }
577
578 // Read input dataset
579 vtkSmartPointer<vtkPointSet> pointset = ReadPointSet(POSARG(1));
580 vtkSmartPointer<vtkPolyData> polydata = ConvertToPolyData(pointset);
581 polydata->BuildLinks();
582
583 // What info to print by default
584 if (!self_intersections_name &&
585 !report_bounds &&
586 !report_cell_types &&
587 !report_cell_volumes &&
588 !report_edge_lengths &&
589 !report_surface_area &&
590 !print_summary &&
591 !print_surface_summary &&
592 !list_pointdata_arrays &&
593 !list_celldata_arrays &&
594 point_indices.empty()) {
595 print_summary = true;
596 print_surface_summary = (polydata->GetNumberOfPolys() != 0);
597 list_pointdata_arrays = true;
598 list_celldata_arrays = true;
599 }
600
601 if (max_point_index >= pointset->GetNumberOfPoints()) {
602 FatalError("Invalid -point ID, point set has only " << pointset->GetNumberOfPoints() << " points!");
603 }
604
605 // Point set type and size
606 if (print_summary) {
607 cout << DataObjectTypeString(pointset->GetDataObjectType()) << ":" << endl;
608 cout << " No. of points: " << pointset->GetNumberOfPoints() << endl;
609 cout << " No. of edges: " << NumberOfEdges(pointset) << endl;
610 cout << " No. of cells: " << pointset->GetNumberOfCells() << endl;
611 }
612
613 // Cell types
614 if (report_cell_types && pointset->GetNumberOfCells() > 0) {
615 OrderedMap<int, int> cellTypeHist;
616 OrderedMap<int, int>::iterator it;
617 for (vtkIdType cellId = 0; cellId < pointset->GetNumberOfCells(); ++cellId) {
618 int type = pointset->GetCellType(cellId);
619 it = cellTypeHist.find(type);
620 if (it == cellTypeHist.end()) cellTypeHist[type] = 1;
621 else ++it->second;
622 }
623 for (it = cellTypeHist.begin(); it != cellTypeHist.end(); ++it) {
624 string s = CellTypeString(it->first);
625 cout << " No. of " << setw(18) << left << (s + ": ") << it->second << endl;
626 }
627 cout << " Maximum cell size: " << pointset->GetMaxCellSize() << endl;
628 }
629
630 // Bounding box
631 if (report_bounds) {
632 double bounds[6];
633 pointset->GetBounds(bounds);
634 cout << " Bounding box: [" << bounds[0] << ", " << bounds[1]
635 << "] x [" << bounds[2] << ", " << bounds[3]
636 << "] x [" << bounds[4] << ", " << bounds[5] << "]" << endl;
637 cout << " Center point: (" << .5 * (bounds[1] + bounds[0]) << ", "
638 << .5 * (bounds[3] + bounds[2]) << ", "
639 << .5 * (bounds[5] + bounds[4]) << ")" << endl;
640 }
641
642 if (pointset->GetNumberOfCells() > 0) {
643
644 // Compute cell volumes
645 if (report_cell_volumes && !IsSurfaceMesh(pointset)) {
646 vtkSmartPointer<vtkDataArray> volume = ComputeCellVolumes::Run(pointset);
647 double min_volume, max_volume, avg_volume, std_volume;
648 Extrema::Calculate(min_volume, max_volume, volume);
649 NormalDistribution::Calculate(avg_volume, std_volume, volume);
650 cout << "\n";
651 cout << " Average cell volume: " << avg_volume << "\n";
652 cout << " Cell volume StDev: " << std_volume << "\n";
653 cout << " Minimum cell volume: " << min_volume << "\n";
654 cout << " Maximum cell volume: " << max_volume << "\n";
655 cout.flush();
656 }
657
658 // Discard degenerate cells and unused points
659 polydata = CleanPolyData(polydata);
660 polydata->SetVerts(NULL);
661 polydata->SetLines(NULL);
662 polydata = CleanPolyData(polydata);
663 polydata->BuildLinks();
664
665 // Compute edge table
666 EdgeTable edgeTable(polydata);
667
668 // Euler characteristic / Genus
669 if (print_surface_summary) {
670
671 int npoints, nedges, nfaces, nbounds, ncomps, euler;
672 double genus = Genus(polydata, edgeTable, &npoints, &nedges, &nfaces,
673 &nbounds, &ncomps, &euler);
674 cout << "\n";
675 cout << "Surface mesh:\n";
676 cout << " V " << npoints << "\n";
677 cout << " E " << nedges << "\n";
678 cout << " F " << nfaces << "\n";
679 cout << " B " << nbounds << "\n";
680 cout << " C " << ncomps << "\n";
681 cout << "\n";
682 cout << " Euler characteristic / Genus (V - E + F = 2C - 2g - B)\n";
683 cout << " Euler: " << euler << "\n";
684 cout << " Genus: " << genus << "\n";
685 cout.flush();
686
687 if (nbounds > 0 && output_surface_name) {
688 UnorderedSet<int> boundaryPtIds = BoundaryPoints(polydata);
689 vtkSmartPointer<vtkDataArray> boundaryMask = NewVTKDataArray(VTK_UNSIGNED_CHAR);
690 boundaryMask->SetName("BoundaryMask");
691 boundaryMask->SetNumberOfComponents(1);
692 boundaryMask->SetNumberOfTuples(polydata->GetNumberOfPoints());
693 boundaryMask->FillComponent(0, .0);
694 for (auto it = boundaryPtIds.begin(); it != boundaryPtIds.end(); ++it) {
695 boundaryMask->SetComponent(*it, 0, 1.0);
696 }
697 polydata->GetPointData()->AddArray(boundaryMask);
698 }
699 }
700
701 if (report_edge_lengths) {
702 double min_length, max_length, avg_length, std_length;
703 EdgeLengthNormalDistribution(polydata->GetPoints(), edgeTable, avg_length, std_length);
704 GetMinMaxEdgeLength(polydata->GetPoints(), edgeTable, min_length, max_length);
705 cout << "\n";
706 cout << " Average edge length: " << avg_length << "\n";
707 cout << " Edge length StDev: " << std_length << "\n";
708 cout << " Minimum edge length: " << min_length << "\n";
709 cout << " Maximum edge length: " << max_length << "\n";
710 cout.flush();
711 }
712
713 if (report_surface_area) {
714 double sum_area = Area(polydata, true);
715 vtkDataArray *area = polydata->GetCellData()->GetArray("Area");
716 double min_area, max_area, avg_area, std_area;
717 Extrema::Calculate(min_area, max_area, area);
718 NormalDistribution::Calculate(avg_area, std_area, area);
719 double sum_volume = Volume(polydata);
720 cout << "\n";
721 cout << " Average cell area: " << avg_area << "\n";
722 cout << " Cell area StDev: " << std_area << "\n";
723 cout << " Minimum cell area: " << min_area << "\n";
724 cout << " Maximum cell area: " << max_area << "\n";
725 cout << " Total surface area: " << sum_area << "\n";
726 cout << " Total volume: " << sum_volume << "\n";
727 cout.flush();
728 }
729
730 // Self-intersections
731 if (self_intersections_name) {
732 cout << endl;
733 cout << " No. of triangle/triangle intersections: ", cout.flush();
734 cout << NumberOfTriangleTriangleIntersections(polydata, self_intersections_name) << endl;
735 if (verbose > 2) {
736 int count = 0;
737 vtkDataArray *tritri = polydata->GetCellData()->GetArray(self_intersections_name);
738 for (vtkIdType cellId = 0; cellId < polydata->GetNumberOfCells(); ++cellId) {
739 if (tritri->GetComponent(cellId, 0) != .0) {
740 ++count;
741 cout << " " << right << setw(5) << count;
742 if (count == 1) cout << "st";
743 else if (count == 2) cout << "nd";
744 else if (count == 3) cout << "rd";
745 else cout << "th";
746 cout << " intersected cell: ID = " << cellId;
747 double area = ComputeArea(polydata->GetCell(cellId));
748 cout << ", area = " << area << "\n";
749 }
750 }
751 cout.flush();
752 }
753 // cout << " No. of vertices inside polyhedron: ", cout.flush();
754 // cout << NumberOfVerticesInsidePolyhedron(polydata) << endl;
755 }
756 }
757
758 // Points
759 if (!point_indices.empty()) {
760 double p[3];
761 const int w = static_cast<int>(ToString(max_point_index).length());
762 for (size_t i = 0; i < point_indices.size(); ++i) {
763 if (point_indices[i] < 0 || point_indices[i] >= pointset->GetNumberOfPoints()) {
764 cout << "Invalid point index: " << point_indices[i] << "\n";
765 } else {
766 pointset->GetPoint(point_indices[i], p);
767 cout << "Point " << ToString(point_indices[i], w)
768 << ": [" << p[0] << ", " << p[1] << ", " << p[2] << "]\n";
769 }
770 }
771 }
772
773 // Attributes
774 if (list_pointdata_arrays) {
775 vtkPointData *pd = pointset->GetPointData();
776 if (pd->GetNumberOfArrays() > 0) {
777 cout << "\nPoint data attributes:\n";
778 for (int i = 0; i < pd->GetNumberOfArrays(); ++i) {
779 vtkDataArray *array = pd->GetArray(i);
780 printf(" %2d %-24s (dim: %2d, type: %-7s kind: %s)\n", i,
781 array->GetName(),
782 array->GetNumberOfComponents(),
783 (VtkDataTypeString(array->GetDataType()) + ",").c_str(),
784 VtkAttributeTypeString(pd->IsArrayAnAttribute(i)).c_str());
785 }
786 }
787 }
788 if (list_celldata_arrays) {
789 vtkCellData *cd = pointset->GetCellData();
790 if (cd->GetNumberOfArrays() > 0) {
791 cout << "\nCell data attributes:\n";
792 for (int i = 0; i < cd->GetNumberOfArrays(); ++i) {
793 vtkDataArray *array = cd->GetArray(i);
794 printf(" %2d %-24s (dim: %2d, type: %-7s kind: %s)\n", i,
795 array->GetName(),
796 array->GetNumberOfComponents(),
797 (VtkDataTypeString(array->GetDataType()) + ",").c_str(),
798 VtkAttributeTypeString(cd->IsArrayAnAttribute(i)).c_str());
799 }
800 }
801 }
802 cout.flush();
803
804 // Write output point set (possibly with additional attributes)
805 if (output_pointset_name) WritePointSet(output_pointset_name, pointset);
806
807 // Write output surface mesh (possibly with additional attributes)
808 if (output_surface_name && polydata->GetNumberOfCells() > 0) {
809 WritePointSet(output_surface_name, polydata);
810 }
811
812 } // dof_info
813 #endif // HAVE_MIRTK_PointSet
814
815 return 0;
816 }
817