From 357068d1b81abf697b209af2c7de642d4ee25a68 Mon Sep 17 00:00:00 2001 From: ccolin Date: Wed, 23 Feb 2022 13:01:50 +0100 Subject: [PATCH] add external points filter --- CMakeLists.txt | 4 +- data/TetMesh2.vtk | 3 - src/external_points_filter.cc | 93 +++++++++++++++++++++ src/external_points_filter.h | 22 +++++ src/main.cc | 148 ++++++++++++++++------------------ 5 files changed, 189 insertions(+), 81 deletions(-) create mode 100644 src/external_points_filter.cc create mode 100644 src/external_points_filter.h diff --git a/CMakeLists.txt b/CMakeLists.txt index aed3ad2..f34243f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,7 +15,9 @@ target_sources(pfe PRIVATE src/aspect_ratio_filter.cc src/aspect_ratio_filter.h src/dihedral_angles_filter.cc - src/dihedral_angles_filter.h) + src/dihedral_angles_filter.h + src/external_points_filter.cc + src/external_points_filter.h) target_link_libraries(pfe PRIVATE VTK::CommonCore diff --git a/data/TetMesh2.vtk b/data/TetMesh2.vtk index 0c45140..7c9c98b 100644 --- a/data/TetMesh2.vtk +++ b/data/TetMesh2.vtk @@ -18,6 +18,3 @@ LOOKUP_TABLE default 1 1 1 -FIELD law 1 -cell_law 1 1 int -0 diff --git a/src/external_points_filter.cc b/src/external_points_filter.cc new file mode 100644 index 0000000..d54836d --- /dev/null +++ b/src/external_points_filter.cc @@ -0,0 +1,93 @@ +#include "external_points_filter.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +/* Remove from ids_a the elements that are not also in ids_b. */ +static void keep_only_dupes(std::vector &ids_a, + vtkIdType n_ids_b, const vtkIdType *ids_b) { + for (auto it = ids_a.begin(); it != ids_a.end();) { + bool duplicate = false; + for (vtkIdType i = 0; i < n_ids_b; i++) { + if (*it == ids_b[i]) { + duplicate = true; + break; + } + } + if (!duplicate) { + it = ids_a.erase(it); + } else { + ++it; + } + } +} + + +vtkStandardNewMacro(ExternalPointsFilter); + + +vtkTypeBool ExternalPointsFilter::RequestData( + vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) { + (void) request; + vtkUnstructuredGrid* input = + vtkUnstructuredGrid::GetData(inputVector[0]); + vtkUnstructuredGrid* output = + vtkUnstructuredGrid::GetData(outputVector); + + /* Copy point and cell data unmodified. */ + output->CopyStructure(input); + output->GetPointData()->PassData(input->GetPointData()); + output->GetCellData()->PassData(input->GetCellData()); + + /* Create an array to store the IDs of external points. */ + vtkNew externalPoints; + externalPoints->SetName("external_points"); + externalPoints->SetNumberOfComponents(1); + vtkNew isExternal; + isExternal->SetName("bl"); + isExternal->SetNumberOfComponents(1); + isExternal->SetNumberOfTuples(input->GetNumberOfPoints()); + isExternal->FillValue(0); + + /* Generate cell links, these tell us which cell a point belongs to. */ + vtkNew links; + links->BuildLinks(input); + + // TODO: try to use vtkDataSet::GetCellNeighbors instead + auto *it = input->NewCellIterator(); + for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) { + vtkIdList *pointIds = it->GetPointIds(); + for (vtkIdType faceId = 0; faceId < 4; faceId++) { + vtkIdType firstPoint = pointIds->GetId(faceId); + vtkIdType n_cells = links->GetNcells(firstPoint); + vtkIdType *cells = links->GetCells(firstPoint); + std::vector intersection(cells, cells + n_cells); + for (vtkIdType i = 1; i < 3; i++) { + vtkIdType pointId = pointIds->GetId((faceId + i) % 4); + keep_only_dupes(intersection, + links->GetNcells(pointId), links->GetCells(pointId)); + } + if (intersection.size() == 1) { // The face only belongs to one cell. + for (vtkIdType i = 0; i < 3; i++) { + externalPoints->InsertNextValue(pointIds->GetId((faceId + i) % 4)); + isExternal->SetValue(pointIds->GetId((faceId + i) % 4), 1); + } + } + } + } + + output->GetPointData()->SetScalars(isExternal); + output->GetFieldData()->AddArray(externalPoints); + return true; +} diff --git a/src/external_points_filter.h b/src/external_points_filter.h new file mode 100644 index 0000000..9c6e5b7 --- /dev/null +++ b/src/external_points_filter.h @@ -0,0 +1,22 @@ +#ifndef EXTERNAL_POINTS_FILTER_H +#define EXTERNAL_POINTS_FILTER_H + + +#include +#include + + +class ExternalPointsFilter : public vtkUnstructuredGridAlgorithm { +public: + static ExternalPointsFilter *New(); + vtkTypeMacro(ExternalPointsFilter, vtkUnstructuredGridAlgorithm); + vtkTypeBool RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; + +protected: + ~ExternalPointsFilter() override = default; +}; + + +#endif diff --git a/src/main.cc b/src/main.cc index bcf8a47..857e8d0 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,6 +1,7 @@ #include "angles_filter.h" #include "aspect_ratio_filter.h" #include "dihedral_angles_filter.h" +#include "external_points_filter.h" #include #include @@ -14,32 +15,10 @@ #include #include #include +#include #include #include -#include -#include - - -// template -// T average(const std::vector &data) { -// T avg = 0; -// for (const T &t : data) { -// avg += t; -// } -// return avg / data.size(); -// } - - -// template -// T standard_deviation(const std::vector &data) { -// T avg = average(data); -// T stddev = 0; -// for (const T &t : data) { -// stddev += (t - avg) * (t - avg); -// } -// return std::sqrt(stddev / data.size()); -// } int main(int argc, char **argv) { @@ -58,74 +37,89 @@ int main(int argc, char **argv) { 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::cout << "cell " << i << ": "; - for (ssize_t j = 0; j < angles->GetNumberOfComponents(); j++) { - std::cout << angles->GetTuple(i)[j] << ", "; - } - std::cout << "\b\b \n"; - } + // 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::cout << "cell " << i << ": " - << aspectRatios->GetTuple1(i) << "\n"; - } + // 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()); /* 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::cout << "dihedral angles "; - for (vtkIdType j = 0; j < 6; j++) { - std::cout << dihedralAngles->GetTuple(i)[j] << ", "; - } - std::cout << "\b\b\n"; - } + // 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"; + // } - /* Volume rendering properties */ - vtkNew volumeMapper; - volumeMapper->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); - vtkNew volume; - volume->SetMapper(volumeMapper); - vtkNew transferFunction; - transferFunction->AddPoint(-1, 0); - transferFunction->AddPoint(1, 1); - volume->GetProperty()->SetScalarOpacity(transferFunction); - volume->GetProperty()->SetColor(transferFunction); + /* External faces filter. */ + // vtkNew externalFacesFilter; + // externalFacesFilter->DebugOn(); + // externalFacesFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); - /* Renderer */ - vtkNew colors; - std::array bkg{{26, 51, 102, 255}}; - colors->SetColor("BkgColor", bkg.data()); - vtkNew renderer; - renderer->AddVolume(volume); - renderer->SetBackground(colors->GetColor3d("BkgColor").GetData()); - renderer->ResetCamera(); - renderer->GetActiveCamera()->Zoom(1.5); - vtkNew renderWindow; - renderWindow->SetSize(300, 300); - renderWindow->AddRenderer(renderer); - renderWindow->SetWindowName("PFE"); - vtkNew renderWindowInteractor; - renderWindowInteractor->SetRenderWindow(renderWindow); - renderWindow->Render(); - renderWindowInteractor->Start(); + /* External points filter. */ + vtkNew externalPointsFilter; + externalPointsFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); + + vtkNew writer; + writer->SetInputConnection(externalPointsFilter->GetOutputPort()); + writer->SetFileTypeToASCII(); + writer->SetFileName("out.vtk"); + writer->Write(); + + // /* Volume rendering properties */ + // vtkNew volumeMapper; + // volumeMapper->SetInputConnection(externalPointsFilter->GetOutputPort()); + // vtkNew volume; + // volume->SetMapper(volumeMapper); + // vtkNew transferFunction; + // transferFunction->AddPoint(-1, 0); + // transferFunction->AddPoint(1, 1); + // volume->GetProperty()->SetScalarOpacity(transferFunction); + // volume->GetProperty()->SetColor(transferFunction); + + // /* Renderer */ + // vtkNew colors; + // std::array bkg{{26, 51, 102, 255}}; + // colors->SetColor("BkgColor", bkg.data()); + // vtkNew renderer; + // renderer->AddVolume(volume); + // renderer->SetBackground(colors->GetColor3d("BkgColor").GetData()); + // renderer->ResetCamera(); + // renderer->GetActiveCamera()->Zoom(1.5); + // vtkNew renderWindow; + // renderWindow->SetSize(300, 300); + // renderWindow->AddRenderer(renderer); + // renderWindow->SetWindowName("PFE"); + // vtkNew renderWindowInteractor; + // renderWindowInteractor->SetRenderWindow(renderWindow); + // renderWindow->Render(); + // renderWindowInteractor->Start(); return EXIT_SUCCESS; }