From 621f24dd8baf3eff15f0b023877adc66f4782c47 Mon Sep 17 00:00:00 2001 From: papush! Date: Mon, 7 Mar 2022 23:09:59 +0100 Subject: [PATCH] start implementing RemoveExternalCellsFilter --- CMakeLists.txt | 8 +- src/main.cc | 126 ++++------------------- src/remove_external_cells_filter.cc | 152 ++++++++++++++++++++++++++++ src/remove_external_cells_filter.h | 29 ++++++ 4 files changed, 209 insertions(+), 106 deletions(-) create mode 100644 src/remove_external_cells_filter.cc create mode 100644 src/remove_external_cells_filter.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 9b3462a..1280cb9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,7 +9,9 @@ set(VTK_COMPONENTS VTK::CommonColor VTK::CommonDataModel VTK::IOLegacy - VTK::IOXML) + VTK::IOGeometry + VTK::IOXML + VTK::FiltersModeling) set(ENABLE_VIEWER OFF CACHE BOOL "Enable the 3D viewer, depends on Qt.") if(ENABLE_VIEWER) list(APPEND VTK_COMPONENTS @@ -55,7 +57,9 @@ target_sources(pfe PRIVATE src/point_tris_dist.cc src/point_tris_dist.h src/project_surface_points_on_poly.cc - src/project_surface_points_on_poly.h) + src/project_surface_points_on_poly.h + src/remove_external_cells_filter.cc + src/remove_external_cells_filter.h) target_link_libraries(pfe PRIVATE ${VTK_COMPONENTS}) diff --git a/src/main.cc b/src/main.cc index 9f731fc..5e2107b 100644 --- a/src/main.cc +++ b/src/main.cc @@ -5,14 +5,16 @@ #include "project_surface_points_on_poly.h" #include "mesh_fit_filter.h" #include "point_tris_dist.h" +#include "remove_external_cells_filter.h" #include #include #include -#include +#include #include #include #include +#include #ifdef USE_VIEWER #include @@ -36,117 +38,33 @@ int main(int argc, char **argv) { return EXIT_FAILURE; } - /* Reader */ - vtkNew reader; - reader->SetFileName(argv[1]); + vtkNew tetMeshReader; + tetMeshReader->SetFileName(argv[1]); - - /* Angles filter */ - vtkNew anglesFilter; - anglesFilter->SetInputConnection(reader->GetOutputPort()); - - /* Display angles in console. */ - // anglesFilter->Update(); - // vtkUnstructuredGrid *grid = anglesFilter->GetOutput(); - // auto *angles = grid->GetCellData()->GetArray("angles"); - // for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) { - // std::cerr << "cell " << i << ": "; - // for (ssize_t j = 0; j < angles->GetNumberOfComponents(); j++) { - // std::cerr << angles->GetTuple(i)[j] << ", "; - // } - // std::cerr << "\b\b \n"; - // } - - /* Aspect ratio */ - vtkNew aspectRatioFilter; - aspectRatioFilter->SetInputConnection(anglesFilter->GetOutputPort()); - - /* Display aspect ratios in console. */ - // aspectRatioFilter->Update(); - // grid = aspectRatioFilter->GetOutput(); - // auto *aspectRatios = grid->GetCellData()->GetArray("aspect_ratio"); - // for (ssize_t i = 0; i < aspectRatios->GetNumberOfTuples(); i++) { - // std::cerr << "cell " << i << ": " - // << aspectRatios->GetTuple1(i) << "\n"; - // } - - /* Dihedral angles filter */ vtkNew dihedralAnglesFilter; - dihedralAnglesFilter->SetInputConnection(aspectRatioFilter->GetOutputPort()); + dihedralAnglesFilter->SetInputConnection(tetMeshReader->GetOutputPort()); - /* Display dihedral angles in console. */ - // dihedralAnglesFilter->Update(); - // grid = dihedralAnglesFilter->GetOutput(); - // auto dihedralAngles = grid->GetCellData()->GetArray("dihedral_angles"); - // for (vtkIdType i = 0; i < grid->GetNumberOfCells(); i++) { - // std::cerr << "dihedral angles "; - // for (vtkIdType j = 0; j < 6; j++) { - // std::cerr << dihedralAngles->GetTuple(i)[j] << ", "; - // } - // std::cerr << "\b\b\n"; - // } + // vtkNew polyMeshReader; + // polyMeshReader->SetFileName(argv[2]); + vtkNew polyMeshReader; + polyMeshReader->SetFileName(argv[2]); - /* External faces filter. */ - // vtkNew externalFacesFilter; - // externalFacesFilter->DebugOn(); - // externalFacesFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); + vtkNew removeExternalCellsFilter; + removeExternalCellsFilter->SetInputConnection(0, + dihedralAnglesFilter->GetOutputPort()); + removeExternalCellsFilter->SetInputConnection(1, + polyMeshReader->GetOutputPort()); - /* Surface points filter. */ - vtkNew surfacePointsFilter; - surfacePointsFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); + // vtkNew project; + // project->SetInputConnection(0, removeExternalCellsFilter->GetOutputPort()); + // project->SetInputConnection(1, pdReader->GetOutputPort()); - - // vtkNew meshFitFilter; - // meshFitFilter->SetInputConnection(surfacePointsFilter->GetOutputPort()); - - vtkNew pdReader; - pdReader->SetFileName(argv[2]); - - vtkNew project; - project->SetInputConnection(0, surfacePointsFilter->GetOutputPort()); - project->SetInputConnection(1, pdReader->GetOutputPort()); - - - vtkNew writer; - writer->SetInputConnection(project->GetOutputPort()); - writer->SetFileTypeToASCII(); - writer->SetFileName("out.vtk"); + vtkNew writer; + writer->SetInputConnection(removeExternalCellsFilter->GetOutputPort()); + writer->SetDataModeToAscii(); + writer->SetFileName("out.vtu"); writer->Write(); - - // 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 profile; - // profile->SetPoints(points); - // profile->SetStrips(strips); - - // 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 d = pointTriangleDistance(point, cell, direction); - // //double d = 5; - - // std::cout << "[" << point[0] - // << ", " << point[1] - // << ", " << point[2] - // << "] (" << d << ")" - // << "[" << direction[0] - // << ", " << direction[1] - // << ", " << direction[2] << "]\n"; - #ifdef USE_VIEWER /* Volume rendering properties */ vtkNew volumeMapper; diff --git a/src/remove_external_cells_filter.cc b/src/remove_external_cells_filter.cc new file mode 100644 index 0000000..c16f2e0 --- /dev/null +++ b/src/remove_external_cells_filter.cc @@ -0,0 +1,152 @@ +#include "remove_external_cells_filter.h" +#include "point_tris_dist.h" + +#include +#include +#include +#include +#include +#include +#include + + +vtkStandardNewMacro(RemoveExternalCellsFilter); + + +RemoveExternalCellsFilter::RemoveExternalCellsFilter() { + SetNumberOfInputPorts(2); +} + + +int RemoveExternalCellsFilter::FillInputPortInformation( + int port, vtkInformation *info) { + if (port == 0) { + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid"); + } else if (port == 1) { + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData"); + } + return 1; +} + + +// static bool isSurfaceCell(vtkUnstructuredGrid *tetMesh, vtkIdType cellId) { +// vtkNew facePoints; +// vtkNew sharedCells; +// for (vtkIdType faceId = 0; faceId < 4; faceId++) { +// facePoints->InsertNextId((faceId + 0) % 4); +// facePoints->InsertNextId((faceId + 1) % 4); +// facePoints->InsertNextId((faceId + 2) % 4); +// tetMesh->GetCellNeighbors(cellId, facePoints, sharedCells); +// if (sharedCells->GetNumberOfIds() == 1) { // Surface cell. +// return true; +// } +// facePoints->Reset(); +// sharedCells->Reset(); +// } +// return false; +// } + + +// static void closestTri(vtkPolyData *polyMesh, +// vtkKdTree *kdTree, +// vtkStaticCellLinks *polyMeshLinks, +// double *point, +// double *closestTriDist, +// double *closestTriDirection) { +// vtkIdType closestPoint; +// { +// double dummy; +// closestPoint = kdTree->FindClosestPoint(point, dummy); +// } +// vtkIdType *triIds = polyMeshLinks->GetCells(closestPoint); +// vtkIdType nTris = polyMeshLinks->GetNcells(closestPoint); + +// *closestTriDist = std::numeric_limits::infinity(); +// for (vtkIdType j = 0; j < nTris; j++) { +// vtkIdType cellId = triIds[j]; +// double direction[3]; +// double dist = pointTriangleDistance(point, polyMesh->GetCell(cellId), +// direction); +// if (dist < *closestTriDist) { +// *closestTriDist = dist; +// closestTriDirection[0] = direction[0]; +// closestTriDirection[1] = direction[1]; +// closestTriDirection[2] = direction[2]; +// } +// } +// } + + +// static bool isExternalCell(vtkUnstructuredGrid *tetMesh, +// vtkPolyData *polyMesh, +// vtkStaticCellLinks *polyMeshLinks, +// vtkKdTree *kdTree, +// vtkIdType cellId) { +// vtkNew pointIds; +// tetMesh->GetCellPoints(cellId, pointIds); +// double a[3], b[3], c[3]; +// tetMesh->GetPoint(pointIds->GetId(0), a); +// tetMesh->GetPoint(pointIds->GetId(1), b); +// tetMesh->GetPoint(pointIds->GetId(2), c); + +// double aDist, bDist, cDist; +// double aDirection[3], bDirection[3], cDirection[3]; +// closestTri(polyMesh, kdTree, polyMeshLinks, a, +// &aDist, aDirection); +// closestTri(polyMesh, kdTree, polyMeshLinks, b, +// &bDist, bDirection); +// closestTri(polyMesh, kdTree, polyMeshLinks, c, +// &cDist, cDirection); +// if (aDist > 0 && bDist > 0 && cDist > 0) { +// return true; +// } +// return false; +// } + + +vtkTypeBool RemoveExternalCellsFilter::RequestData( + vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) { + (void) request; + vtkUnstructuredGrid* tetMesh = vtkUnstructuredGrid::GetData(inputVector[0]); + vtkPolyData *polyMesh = vtkPolyData::GetData(inputVector[1]); + vtkUnstructuredGrid* output = vtkUnstructuredGrid::GetData(outputVector); + // output->CopyStructure(tetMesh); + // output->GetPointData()->PassData(tetMesh->GetPointData()); + // output->GetCellData()->PassData(tetMesh->GetCellData()); + + vtkNew selectEnclosedPoints; + selectEnclosedPoints->SetInputDataObject(0, tetMesh); + selectEnclosedPoints->SetInputDataObject(1, polyMesh); + selectEnclosedPoints->Update(); + unsigned char *enclosedPoints = (unsigned char *) + selectEnclosedPoints->GetOutput()->GetPointData() + ->GetArray("SelectedPoints")->GetVoidPointer(0); + vtkNew outPoints; + vtkNew outPointIds; + vtkIdType outPointId = 0; + outPointIds->Allocate(3); + output->Allocate(); + auto it = tetMesh->NewCellIterator(); + for (it->InitTraversal(); + !it->IsDoneWithTraversal(); + it->GoToNextCell()) { + vtkIdList *pointIds = it->GetPointIds(); + if (!enclosedPoints[pointIds->GetId(0)] + && !enclosedPoints[pointIds->GetId(1)] + && !enclosedPoints[pointIds->GetId(2)] + && !enclosedPoints[pointIds->GetId(3)]) { + outPointIds->SetId(0, outPointId + 0); + outPointIds->SetId(1, outPointId + 1); + outPointIds->SetId(2, outPointId + 2); + outPointId += 3; + outPoints->InsertPoints(outPointIds, pointIds, + tetMesh->GetPoints()); + output->InsertNextCell(it->GetCellType(), outPointIds); + } + } + output->SetPoints(outPoints); + + return true; +} diff --git a/src/remove_external_cells_filter.h b/src/remove_external_cells_filter.h new file mode 100644 index 0000000..8b76764 --- /dev/null +++ b/src/remove_external_cells_filter.h @@ -0,0 +1,29 @@ +#ifndef REMOVE_EXTERNAL_CELLS_FILTER_H +#define REMOVE_EXTERNAL_CELLS_FILTER_H + + +#include +#include +#include + + +/* + Removes cells from the UnstructuredGrid passed as first input that + lay outside the PolyData passed as second input. +*/ +class RemoveExternalCellsFilter : public vtkUnstructuredGridAlgorithm { +public: + static RemoveExternalCellsFilter *New(); + vtkTypeMacro(RemoveExternalCellsFilter, vtkUnstructuredGridAlgorithm); + RemoveExternalCellsFilter(); + int FillInputPortInformation(int port, vtkInformation *info) override; + vtkTypeBool RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; + +protected: + ~RemoveExternalCellsFilter() override = default; +}; + + +#endif