From 8e0c2317084d9381c28c2180bff71802baab5855 Mon Sep 17 00:00:00 2001 From: ccolin Date: Thu, 17 Feb 2022 14:02:20 +0100 Subject: [PATCH] refactor angle filter --- CMakeLists.txt | 6 +- src/angles_filter.cc | 87 ++++++++++++++++++++++++++ src/angles_filter.h | 25 ++++++++ src/main.cc | 144 ++++++++----------------------------------- 4 files changed, 141 insertions(+), 121 deletions(-) create mode 100644 src/angles_filter.cc create mode 100644 src/angles_filter.h diff --git a/CMakeLists.txt b/CMakeLists.txt index a2cc230..a31b39c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,7 +8,11 @@ add_executable(pfe) target_compile_features(pfe PRIVATE cxx_std_17) -target_sources(pfe PRIVATE src/main.cc) +target_sources(pfe PRIVATE + src/main.cc + src/angles_filter.cc + src/angles_filter.h) + target_link_libraries(pfe PRIVATE VTK::CommonCore VTK::ViewsCore diff --git a/src/angles_filter.cc b/src/angles_filter.cc new file mode 100644 index 0000000..5df391e --- /dev/null +++ b/src/angles_filter.cc @@ -0,0 +1,87 @@ +#include "angles_filter.h" + +#include +#include +#include +#include +#include + + +vtkStandardNewMacro(AnglesFilter); + + +vtkTypeBool AnglesFilter::RequestData( + vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) { + (void) request; + vtkUnstructuredGrid* input = + vtkUnstructuredGrid::GetData(inputVector[0]); + vtkUnstructuredGrid* output = + vtkUnstructuredGrid::GetData(outputVector); + output->CopyStructure(input); + output->GetPointData()->PassData(input->GetPointData()); + vtkCellData *cellData = output->GetCellData(); + cellData->PassData(input->GetCellData()); + vtkNew anglesArray; + anglesArray->SetName("angles"); + anglesArray->SetNumberOfComponents(12); + anglesArray->SetNumberOfTuples(input->GetNumberOfCells()); + double *anglesBase = anglesArray->GetPointer(0); + size_t i = 0; + auto it = input->NewCellIterator(); + for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) { + if (it->GetCellType() != VTK_TETRA) continue; + cellAngles(input, it->GetPointIds(), anglesBase + i); + i += 12; + } + cellData->AddArray((vtkAbstractArray *) anglesArray); + return true; +} + + +void cellAngles(vtkDataSet *dataSet, vtkIdList *idList, double *angles) { + // std::cout << "nb points: " << idList->GetNumberOfIds() << std::endl; + double a[3], b[3], c[3], d[3]; + dataSet->GetPoint(idList->GetId(0), a); + dataSet->GetPoint(idList->GetId(1), b); + dataSet->GetPoint(idList->GetId(2), c); + dataSet->GetPoint(idList->GetId(3), d); + // std::cout << "ids " << idList->GetId(0) + // << " " << idList->GetId(1) + // << " " << idList->GetId(2) + // << " " << idList->GetId(3) << std::endl; + // std::cout << "coords" << std::endl + // << a[0] << ", " << a[1] << ", " << a[2] << std::endl + // << b[0] << ", " << b[1] << ", " << b[2] << std::endl + // << c[0] << ", " << c[1] << ", " << c[2] << std::endl + // << d[0] << ", " << d[1] << ", " << d[2] << std::endl; + double ab[3], ac[3], ad[3], + ba[3], bc[3], bd[3], + ca[3], cb[3], cd[3], + da[3], db[3], dc[3]; + vtkMath::Subtract(b, a, ab); + vtkMath::Subtract(c, a, ac); + vtkMath::Subtract(d, a, ad); + vtkMath::Subtract(a, b, ba); + vtkMath::Subtract(c, b, bc); + vtkMath::Subtract(d, b, bd); + vtkMath::Subtract(a, c, ca); + vtkMath::Subtract(b, c, cb); + vtkMath::Subtract(d, c, cd); + vtkMath::Subtract(a, d, da); + vtkMath::Subtract(b, d, db); + vtkMath::Subtract(c, d, dc); + angles[0] = vtkMath::AngleBetweenVectors(ab, ac); + angles[1] = vtkMath::AngleBetweenVectors(ac, ad); + angles[2] = vtkMath::AngleBetweenVectors(ad, ab); + angles[3] = vtkMath::AngleBetweenVectors(ba, bc); + angles[4] = vtkMath::AngleBetweenVectors(bc, bd); + angles[5] = vtkMath::AngleBetweenVectors(bd, ba); + angles[6] = vtkMath::AngleBetweenVectors(ca, cb); + angles[7] = vtkMath::AngleBetweenVectors(cb, cd); + angles[8] = vtkMath::AngleBetweenVectors(cd, ca); + angles[9] = vtkMath::AngleBetweenVectors(da, db); + angles[10] = vtkMath::AngleBetweenVectors(db, dc); + angles[11] = vtkMath::AngleBetweenVectors(dc, da); +} diff --git a/src/angles_filter.h b/src/angles_filter.h new file mode 100644 index 0000000..f5c6956 --- /dev/null +++ b/src/angles_filter.h @@ -0,0 +1,25 @@ +#ifndef ANGLES_FILTER_H +#define ANGLES_FILTER_H + + +#include +#include + + +void cellAngles(vtkDataSet *dataSet, vtkIdList *idList, double *angles); + + +class AnglesFilter : public vtkUnstructuredGridAlgorithm { +public: + static AnglesFilter *New(); + vtkTypeMacro(AnglesFilter, vtkUnstructuredGridAlgorithm); + vtkTypeBool RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; + +protected: + ~AnglesFilter() override = default; +}; + + +#endif diff --git a/src/main.cc b/src/main.cc index 33ad8b5..a472af0 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,15 +1,9 @@ -#include -#include +#include "angles_filter.h" + #include #include -#include -#include #include -#include #include -#include -#include -#include #include #include #include @@ -18,11 +12,7 @@ #include #include #include -#include -#include #include -#include -#include #include #include @@ -30,111 +20,25 @@ #include -void cellAngles(vtkDataSet *dataSet, vtkIdList *idList, double *angles) { - // std::cout << "nb points: " << idList->GetNumberOfIds() << std::endl; - double a[3], b[3], c[3], d[3]; - dataSet->GetPoint(idList->GetId(0), a); - dataSet->GetPoint(idList->GetId(1), b); - dataSet->GetPoint(idList->GetId(2), c); - dataSet->GetPoint(idList->GetId(3), d); - // std::cout << "ids " << idList->GetId(0) - // << " " << idList->GetId(1) - // << " " << idList->GetId(2) - // << " " << idList->GetId(3) << std::endl; - // std::cout << "coords" << std::endl - // << a[0] << ", " << a[1] << ", " << a[2] << std::endl - // << b[0] << ", " << b[1] << ", " << b[2] << std::endl - // << c[0] << ", " << c[1] << ", " << c[2] << std::endl - // << d[0] << ", " << d[1] << ", " << d[2] << std::endl; - double ab[3], ac[3], ad[3], - ba[3], bc[3], bd[3], - ca[3], cb[3], cd[3], - da[3], db[3], dc[3]; - vtkMath::Subtract(b, a, ab); - vtkMath::Subtract(c, a, ac); - vtkMath::Subtract(d, a, ad); - vtkMath::Subtract(a, b, ba); - vtkMath::Subtract(c, b, bc); - vtkMath::Subtract(d, b, bd); - vtkMath::Subtract(a, c, ca); - vtkMath::Subtract(b, c, cb); - vtkMath::Subtract(d, c, cd); - vtkMath::Subtract(a, d, da); - vtkMath::Subtract(b, d, db); - vtkMath::Subtract(c, d, dc); - angles[0] = vtkMath::AngleBetweenVectors(ab, ac); - angles[1] = vtkMath::AngleBetweenVectors(ac, ad); - angles[2] = vtkMath::AngleBetweenVectors(ad, ab); - angles[3] = vtkMath::AngleBetweenVectors(ba, bc); - angles[4] = vtkMath::AngleBetweenVectors(bc, bd); - angles[5] = vtkMath::AngleBetweenVectors(bd, ba); - angles[6] = vtkMath::AngleBetweenVectors(ca, cb); - angles[7] = vtkMath::AngleBetweenVectors(cb, cd); - angles[8] = vtkMath::AngleBetweenVectors(cd, ca); - angles[9] = vtkMath::AngleBetweenVectors(da, db); - angles[10] = vtkMath::AngleBetweenVectors(db, dc); - angles[11] = vtkMath::AngleBetweenVectors(dc, da); -} +// template +// T average(const std::vector &data) { +// T avg = 0; +// for (const T &t : data) { +// avg += t; +// } +// return avg / data.size(); +// } -class AddAngleComputationAlgorithm : public vtkUnstructuredGridAlgorithm { -public: - static AddAngleComputationAlgorithm *New(); - vtkTypeMacro(AddAngleComputationAlgorithm, vtkUnstructuredGridAlgorithm); - - vtkTypeBool RequestData(vtkInformation *request, - vtkInformationVector **inputVector, - vtkInformationVector *outputVector) override { - (void) request; - vtkUnstructuredGrid* input = - vtkUnstructuredGrid::GetData(inputVector[0]); - vtkUnstructuredGrid* output = - vtkUnstructuredGrid::GetData(outputVector); - output->CopyStructure(input); - output->GetPointData()->PassData(input->GetPointData()); - vtkCellData *cellData = output->GetCellData(); - cellData->PassData(input->GetCellData()); - vtkNew anglesArray; - anglesArray->SetName("angles"); - anglesArray->SetNumberOfComponents(12); - anglesArray->SetNumberOfTuples(input->GetNumberOfCells()); - double *anglesBase = anglesArray->GetPointer(0); - size_t i = 0; - auto it = input->NewCellIterator(); - for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) { - if (it->GetCellType() != VTK_TETRA) continue; - cellAngles(input, it->GetPointIds(), anglesBase + i); - i += 12; - } - cellData->AddArray((vtkAbstractArray *) anglesArray); - return true; - } - -protected: - ~AddAngleComputationAlgorithm() override = default; -}; - -vtkStandardNewMacro(AddAngleComputationAlgorithm); - -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()); -} +// 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) { @@ -152,10 +56,10 @@ int main(int argc, char **argv) { vtkNew reader; reader->SetFileName(argv[1]); // reader->Update(); - vtkNew addAngles; - addAngles->SetInputConnection(reader->GetOutputPort()); - addAngles->Update(); - vtkUnstructuredGrid *grid = addAngles->GetOutput(); + vtkNew anglesFilter; + anglesFilter->SetInputConnection(reader->GetOutputPort()); + anglesFilter->Update(); + vtkUnstructuredGrid *grid = anglesFilter->GetOutput(); auto *angles = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("angles")); for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) { std::cout << "cell " << i << ": "; @@ -170,7 +74,7 @@ int main(int argc, char **argv) { // << ", max: " << *std::max_element(angles.begin(), angles.end()) << std::endl; vtkNew volumeMapper; - volumeMapper->SetInputConnection(addAngles->GetOutputPort()); + volumeMapper->SetInputConnection(anglesFilter->GetOutputPort()); // vtkNew mapper; // mapper->SetInputConnection(reader->GetOutputPort());