diff --git a/CMakeLists.txt b/CMakeLists.txt index eafd6db..3e03e9e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,6 +12,7 @@ set(VTK_COMPONENTS VTK::IOGeometry VTK::IOXML VTK::FiltersModeling + VTK::FiltersGeometry VTK::vtksys) set(ENABLE_VIEWER OFF CACHE BOOL "Enable the 3D viewer, depends on Qt.") if(ENABLE_VIEWER) @@ -62,7 +63,11 @@ target_sources(pfe PRIVATE src/remove_external_cells_filter.cc src/remove_external_cells_filter.h src/relaxation_filter.cc - src/relaxation_filter.h) + src/relaxation_filter.h + src/max_distance_filter.cc + src/max_distance_filter.h + src/closest_polymesh_point.cc + src/closest_polymesh_point.h) target_link_libraries(pfe PRIVATE ${VTK_COMPONENTS}) diff --git a/data/tetrahedron.obj b/data/tetrahedron.obj new file mode 100644 index 0000000..5ac91e3 --- /dev/null +++ b/data/tetrahedron.obj @@ -0,0 +1,12 @@ +# Blender v2.79 (sub 0) OBJ File: '' +# www.blender.org +o Solid +v 0.000000 0.769688 0.000000 +v 0.000000 0.000000 -0.769688 +v 0.769688 0.000000 0.000000 +v 0.000000 0.000000 0.000000 +s off +f 1 2 3 +f 1 3 4 +f 1 4 2 +f 2 4 3 diff --git a/src/closest_polymesh_point.cc b/src/closest_polymesh_point.cc new file mode 100644 index 0000000..7053093 --- /dev/null +++ b/src/closest_polymesh_point.cc @@ -0,0 +1,60 @@ +#include "closest_polymesh_point.h" +#include "point_tris_dist.h" +#include + + +vtkIdType findClosestPoint(const double *point, vtkPointSet *pointSet) { + vtkIdType minPoint = 0; + double minDist2 = std::numeric_limits::infinity(); + for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) { + double point2[3]; + pointSet->GetPoint(i, point2); + double dist2 = vtkMath::Distance2BetweenPoints(point, point2); + if (dist2 < minDist2) { + minDist2 = dist2; + minPoint = i; + } + } + return minPoint; +} + + +void closestPolyMeshPoint(vtkPolyData *polyMesh, + const double *point, + vtkKdTree *kdTree, + vtkStaticCellLinks *links, + double *toClosestPoint, + double *distanceToClosestPoint) { + vtkIdType closestPolyMeshPoint; + { + double dummy; + // closestPolyMeshPoint = kdTree->FindClosestPoint((double *) point, dummy); + closestPolyMeshPoint = findClosestPoint(point, polyMesh); + } + vtkIdType *cellIds = links->GetCells(closestPolyMeshPoint); + vtkIdType nCells = links->GetNcells(closestPolyMeshPoint); + if (nCells == 0) { + *distanceToClosestPoint = -1; + return; + } + double minDist = std::numeric_limits::infinity(); + + /* For each tri the closest point belongs to. */ + for (vtkIdType i = 0; i < nCells; i++) { + vtkIdType cellId = cellIds[i]; + double direction[3]; + double dist = pointTriangleDistance( + (double *) point, polyMesh->GetCell(cellId), direction); + /* Find the closest one. */ + if (dist < minDist) { + minDist = dist; + toClosestPoint[0] = direction[0]; + toClosestPoint[1] = direction[1]; + toClosestPoint[2] = direction[2]; + } + } + vtkMath::MultiplyScalar(toClosestPoint, minDist); + *distanceToClosestPoint = minDist; +} + + diff --git a/src/closest_polymesh_point.h b/src/closest_polymesh_point.h new file mode 100644 index 0000000..fef4b32 --- /dev/null +++ b/src/closest_polymesh_point.h @@ -0,0 +1,17 @@ +#ifndef CLOSEST_POLYMESH_POINT_H +#define CLOSEST_POLYMESH_POINT_H + +#include +#include +#include + + +void closestPolyMeshPoint(vtkPolyData *polyMesh, + const double *point, + vtkKdTree *kdTree, + vtkStaticCellLinks *links, + double *toClosestPoint, + double *distanceToClosestPoint); + + +#endif \ No newline at end of file diff --git a/src/dihedral_angles_filter.cc b/src/dihedral_angles_filter.cc index 60581ee..9040ad7 100644 --- a/src/dihedral_angles_filter.cc +++ b/src/dihedral_angles_filter.cc @@ -1,5 +1,7 @@ #include "dihedral_angles_filter.h" +#include +#include #include #include #include @@ -36,6 +38,9 @@ vtkTypeBool DihedralAnglesFilter::RequestData( minDihedralAngleArray->SetNumberOfTuples(input->GetNumberOfCells()); double *minDihedralAngleBase = minDihedralAngleArray->GetPointer(0); + AverageMinDegrees = 0; + MinMinDegrees = std::numeric_limits::infinity(); + size_t i = 0; vtkCellIterator *it = input->NewCellIterator(); for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) { @@ -62,9 +67,14 @@ vtkTypeBool DihedralAnglesFilter::RequestData( minDihedralAngleBase[i / 6] = *std::min_element(dihedralAnglesBase + i, dihedralAnglesBase + i + 6); + AverageMinDegrees += minDihedralAngleBase[i / 6]; + MinMinDegrees = std::min(MinMinDegrees, minDihedralAngleBase[i / 6]); i += 6; } + AverageMinDegrees /= input->GetNumberOfCells(); + AverageMinDegrees = AverageMinDegrees * 180. / 3.141592653589793; + MinMinDegrees = MinMinDegrees * 180. / 3.141592653589793; cellData->AddArray((vtkAbstractArray *) dihedralAnglesArray); cellData->AddArray((vtkAbstractArray *) minDihedralAngleArray); return true; diff --git a/src/dihedral_angles_filter.h b/src/dihedral_angles_filter.h index 2a441fb..ad46c81 100644 --- a/src/dihedral_angles_filter.h +++ b/src/dihedral_angles_filter.h @@ -12,9 +12,13 @@ public: vtkTypeBool RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) override; + vtkGetMacro(AverageMinDegrees, double); + vtkGetMacro(MinMinDegrees, double); protected: ~DihedralAnglesFilter() override = default; + double AverageMinDegrees; + double MinMinDegrees; }; diff --git a/src/main.cc b/src/main.cc index cd357b8..efbebcc 100644 --- a/src/main.cc +++ b/src/main.cc @@ -3,7 +3,9 @@ #include "surface_points_filter.h" #include "project_surface_points_on_poly.h" #include "relaxation_filter.h" +#include "max_distance_filter.h" +#include #include #include #include @@ -98,12 +100,24 @@ int main(int argc, char **argv) { dihedralAnglesFilter->SetInputConnection( relaxationFilter->GetOutputPort()); + vtkNew geometryFilter; + geometryFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); + + vtkNew maxDistanceFilter; + maxDistanceFilter->SetInputConnection(0, geometryFilter->GetOutputPort()); + maxDistanceFilter->SetInputConnection(1, polyMeshReader->GetOutputPort()); + vtkNew writer; writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); writer->SetDataModeToAscii(); writer->SetFileName("out.vtu"); writer->Write(); + maxDistanceFilter->Update(); + std::cerr << "Max distance: " << maxDistanceFilter->GetMaxDist() << "\n" + << "Average min angle: " << dihedralAnglesFilter->GetAverageMinDegrees() << "\n" + << "Min min angle: " << dihedralAnglesFilter->GetMinMinDegrees() << "\n"; + #ifdef USE_VIEWER /* Volume rendering properties */ vtkNew volumeMapper; diff --git a/src/max_distance_filter.cc b/src/max_distance_filter.cc new file mode 100644 index 0000000..6732a98 --- /dev/null +++ b/src/max_distance_filter.cc @@ -0,0 +1,63 @@ +#include "max_distance_filter.h" +#include "closest_polymesh_point.h" + +#include +#include +#include +#include + + +vtkStandardNewMacro(MaxDistanceFilter); + + +MaxDistanceFilter::MaxDistanceFilter() { + SetNumberOfInputPorts(2); + SetNumberOfOutputPorts(0); +} + + +int MaxDistanceFilter::FillInputPortInformation(int port, + vtkInformation *info) { + if (port == 0) { + vtkPolyDataAlgorithm::FillInputPortInformation(port, info); + } else if (port == 1) { + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData"); + } + return 1; +} + + +vtkTypeBool MaxDistanceFilter::RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) { + (void) request; + (void) outputVector; + vtkPolyData* input1 = vtkPolyData::GetData(inputVector[0]); + vtkPolyData* input2 = vtkPolyData::GetData(inputVector[1]); + vtkNew tree1; + tree1->BuildLocatorFromPoints(input1->GetPoints()); + vtkNew tree2; + tree2->BuildLocatorFromPoints(input2->GetPoints()); + vtkNew links1; + links1->BuildLinks(input1); + vtkNew links2; + links2->BuildLinks(input2); + MaxDist = 0; + for (vtkIdType i = 0; i < input1->GetNumberOfPoints(); i++) { + double point[3]; + input1->GetPoint(i, point); + double vec[3]; + double dist; + closestPolyMeshPoint(input2, point, tree2, links2, vec, &dist); + MaxDist = std::max(dist, MaxDist); + } + for (vtkIdType i = 0; i < input2->GetNumberOfPoints(); i++) { + double point[3]; + input2->GetPoint(i, point); + double vec[3]; + double dist; + closestPolyMeshPoint(input1, point, tree1, links1, vec, &dist); + MaxDist = std::max(dist, MaxDist); + } + return true; +} diff --git a/src/max_distance_filter.h b/src/max_distance_filter.h new file mode 100644 index 0000000..6ecf0e8 --- /dev/null +++ b/src/max_distance_filter.h @@ -0,0 +1,23 @@ +#ifndef MAX_DISTANCE_FILTER_H +#define MAX_DISTANCE_FILTER_H + +#include + + +class MaxDistanceFilter : public vtkPolyDataAlgorithm { +public: + static MaxDistanceFilter *New(); + vtkTypeMacro(MaxDistanceFilter, vtkPolyDataAlgorithm); + MaxDistanceFilter(); + int FillInputPortInformation(int, vtkInformation *info) override; + vtkTypeBool RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; + vtkGetMacro(MaxDist, double); + +protected: + double MaxDist; +}; + + +#endif \ No newline at end of file diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc index 0e9891f..5ae8850 100644 --- a/src/project_surface_points_on_poly.cc +++ b/src/project_surface_points_on_poly.cc @@ -1,6 +1,5 @@ -#include "point_tris_dist.h" - #include "project_surface_points_on_poly.h" +#include "closest_polymesh_point.h" #include #include @@ -36,40 +35,6 @@ int ProjectSurfacePointsOnPoly::FillInputPortInformation(int port, vtkInformatio } -static void closestPolyMeshPoint(vtkPolyData *polyMesh, - const double *point, - vtkKdTree *kdTree, - vtkStaticCellLinks *links, - double *toClosestPoint, - double *distanceToClosestPoint) { - vtkIdType closestPolyMeshPoint; - { - double dummy; - closestPolyMeshPoint = kdTree->FindClosestPoint((double *) point, dummy); - } - vtkIdType *cellIds = links->GetCells(closestPolyMeshPoint); - vtkIdType nCells = links->GetNcells(closestPolyMeshPoint); - double minDist = std::numeric_limits::infinity(); - - /* For each tri the closest point belongs to. */ - for (vtkIdType i = 0; i < nCells; i++) { - vtkIdType cellId = cellIds[i]; - double direction[3]; - double dist = pointTriangleDistance( - (double *) point, polyMesh->GetCell(cellId), direction); - /* Find the closest one. */ - if (dist < minDist) { - minDist = dist; - toClosestPoint[0] = direction[0]; - toClosestPoint[1] = direction[1]; - toClosestPoint[2] = direction[2]; - } - } - vtkMath::MultiplyScalar(toClosestPoint, minDist); - *distanceToClosestPoint = minDist; -} - - void findPointsWithinRadius(double radius, vtkPointSet *pointSet, const double *point, vtkIdList *points) { for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) {