diff --git a/CMakeLists.txt b/CMakeLists.txt index 093b2bf..62d01c0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -48,13 +48,14 @@ target_sources(pfe PRIVATE src/dihedral_angles_filter.h src/external_points_filter.cc src/external_points_filter.h - src/kd_tree.cc src/kd_tree.h src/mesh_fit_filter.cc src/mesh_fit_filter.h src/point_tris_dist.cc - src/point_tris_dist.h) + src/point_tris_dist.h + src/project_external_points_on_surface.cc + src/project_external_points_on_surface.h) target_link_libraries(pfe PRIVATE ${VTK_COMPONENTS}) diff --git a/src/main.cc b/src/main.cc index 3e0c2fb..4ad8fb5 100644 --- a/src/main.cc +++ b/src/main.cc @@ -2,7 +2,7 @@ #include "aspect_ratio_filter.h" #include "dihedral_angles_filter.h" #include "external_points_filter.h" - +#include "project_external_points_on_surface.h" #include "mesh_fit_filter.h" #include "point_tris_dist.h" @@ -12,6 +12,7 @@ #include #include #include +#include #ifdef USE_VIEWER #include @@ -30,8 +31,8 @@ int main(int argc, char **argv) { - if (argc != 2) { - std::cerr << "Usage: " << argv[0] << " FILE.vtk" << std::endl; + if (argc != 3) { + std::cerr << "Usage: " << argv[0] << " tetmesh polydata" << std::endl; return EXIT_FAILURE; } @@ -95,49 +96,56 @@ int main(int argc, char **argv) { externalPointsFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); - vtkNew meshFitFilter; - meshFitFilter->SetInputConnection(externalPointsFilter->GetOutputPort()); + // vtkNew meshFitFilter; + // meshFitFilter->SetInputConnection(externalPointsFilter->GetOutputPort()); + + vtkNew pdReader; + pdReader->SetFileName(argv[2]); + + vtkNew pepouse; + pepouse->SetInputConnection(0, externalPointsFilter->GetOutputPort()); + pepouse->SetInputConnection(1, pdReader->GetOutputPort()); vtkNew writer; - writer->SetInputConnection(meshFitFilter->GetOutputPort()); + writer->SetInputConnection(pepouse->GetOutputPort()); writer->SetFileTypeToASCII(); writer->SetFileName("out.vtk"); writer->Write(); - vtkNew points; - points->InsertPoint(0, 0., 0., 0.); - points->InsertPoint(1, 1., 0., 0.); - points->InsertPoint(2, 1., 1., 0.); + // vtkNew points; + // points->InsertPoint(0, 0., 0., 0.); + // points->InsertPoint(1, 1., 0., 0.); + // points->InsertPoint(2, 1., 1., 0.); - vtkNew strips; - strips->InsertNextCell(3); - strips->InsertCellPoint(0); - strips->InsertCellPoint(1); - strips->InsertCellPoint(2); + // vtkNew strips; + // strips->InsertNextCell(3); + // strips->InsertCellPoint(0); + // strips->InsertCellPoint(1); + // strips->InsertCellPoint(2); - vtkNew profile; - profile->SetPoints(points); - profile->SetStrips(strips); + // vtkNew profile; + // profile->SetPoints(points); + // profile->SetStrips(strips); - vtkCell *cell = profile->GetCell(0); - // vtkTriangle *triangle = vtkTriangle::SafeDownCast(cell); - // vtkTriangle *triangle = (vtkTriangle *) cell; + // vtkCell *cell = profile->GetCell(0); + // // vtkTriangle *triangle = vtkTriangle::SafeDownCast(cell); + // // vtkTriangle *triangle = (vtkTriangle *) cell; - double direction[3];// = {0, 0, 0}; - double point[3] = {0, 1, 0}; + // double direction[3];// = {0, 0, 0}; + // double point[3] = {0, 1, 0}; - double d = pointTriangleDistance(point, cell, direction); - //double d = 5; + // double d = pointTriangleDistance(point, cell, direction); + // //double d = 5; - std::cout << "[" << point[0] - << ", " << point[1] - << ", " << point[2] - << "] (" << d << ")" - << "[" << direction[0] - << ", " << direction[1] - << ", " << direction[2] << "]\n"; + // std::cout << "[" << point[0] + // << ", " << point[1] + // << ", " << point[2] + // << "] (" << d << ")" + // << "[" << direction[0] + // << ", " << direction[1] + // << ", " << direction[2] << "]\n"; #ifdef USE_VIEWER /* Volume rendering properties */ diff --git a/src/project_external_points_on_surface.cc b/src/project_external_points_on_surface.cc new file mode 100644 index 0000000..1f29492 --- /dev/null +++ b/src/project_external_points_on_surface.cc @@ -0,0 +1,93 @@ +#include "point_tris_dist.h" + +#include "project_external_points_on_surface.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +vtkStandardNewMacro(ProjectExternalPointsOnSurface); + + +ProjectExternalPointsOnSurface::ProjectExternalPointsOnSurface() { + SetNumberOfInputPorts(2); +} + + +int ProjectExternalPointsOnSurface::FillInputPortInformation(int port, vtkInformation *info) { + if (port == 0) { + vtkUnstructuredGridAlgorithm::FillInputPortInformation(port, info); + } else if (port == 1) { + vtkNew input2; + input2->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData"); + info->Append(input2); + } + return 1; +} + + +vtkTypeBool ProjectExternalPointsOnSurface::RequestData( + vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) { + (void) request; + vtkUnstructuredGrid* us = vtkUnstructuredGrid::GetData(inputVector[0]); + vtkPolyData* pd = vtkPolyData::GetData(inputVector[1]); + vtkUnstructuredGrid* output = + vtkUnstructuredGrid::GetData(outputVector); + + output->CopyStructure(us); + output->GetPointData()->PassData(us->GetPointData()); + output->GetCellData()->PassData(us->GetCellData()); + + /* Generate cell links, these tell us which cell a point belongs to. */ + vtkNew links; + links->BuildLinks(pd); + + vtkNew kdTree; + kdTree->BuildLocatorFromPoints(pd->GetPoints()); + + vtkIdTypeArray *externalPoints = + vtkIdTypeArray::SafeDownCast(us->GetFieldData()->GetArray("external_points")); + for (vtkIdType i = 0; i < externalPoints->GetNumberOfTuples(); i++) { + double p[3]; + us->GetPoint(externalPoints->GetTuple1(i), p); + double dist2; + vtkIdType closest = kdTree->FindClosestPoint(p, dist2); + vtkIdType *cellIds = links->GetCells(closest); + vtkIdType nCells = links->GetNcells(closest); + double min_dist = std::numeric_limits::infinity(); + double min_direction[3]; + + /* For each tri the closest point belongs to. */ + for (vtkIdType j = 0; j < nCells; j++) { + vtkIdType cellId = cellIds[j]; + double direction[3]; + double dist = pointTriangleDistance(p, pd->GetCell(cellId), direction); + /* Find the closest one. */ + if (dist < min_dist) { + min_dist = dist; + min_direction[0] = direction[0]; + min_direction[1] = direction[1]; + min_direction[2] = direction[2]; + } + } + + vtkMath::MultiplyScalar(min_direction, min_dist); + vtkMath::Subtract(p, min_direction, p); + us->GetPoints()->SetPoint(externalPoints->GetTuple1(i), p); + us->Modified(); + } + + return true; +} diff --git a/src/project_external_points_on_surface.h b/src/project_external_points_on_surface.h new file mode 100644 index 0000000..c19077f --- /dev/null +++ b/src/project_external_points_on_surface.h @@ -0,0 +1,23 @@ +#ifndef PROJECT_EXTERNAL_POINTS_ON_SURFACE_H +#define PROJECT_EXTERNAL_POINTS_ON_SURFACE_H + +#include +#include + + +class ProjectExternalPointsOnSurface : public vtkUnstructuredGridAlgorithm { +public: + static ProjectExternalPointsOnSurface *New(); + vtkTypeMacro(ProjectExternalPointsOnSurface, vtkUnstructuredGridAlgorithm); + ProjectExternalPointsOnSurface(); + int FillInputPortInformation(int, vtkInformation *info) override; + vtkTypeBool RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; + +protected: + ~ProjectExternalPointsOnSurface() override = default; +}; + + +#endif