From b894ac2d95de75b181e4aff067d55c221ae5ae3f Mon Sep 17 00:00:00 2001 From: ccolin Date: Tue, 15 Feb 2022 12:55:12 +0100 Subject: [PATCH] implement basic angle computation and analysis --- CMakeLists.txt | 2 + src/main.cc | 99 ++++++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 98 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c0c2ed7..a2cc230 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,6 +6,8 @@ add_subdirectory(external/vtk) add_executable(pfe) +target_compile_features(pfe PRIVATE cxx_std_17) + target_sources(pfe PRIVATE src/main.cc) target_link_libraries(pfe PRIVATE VTK::CommonCore diff --git a/src/main.cc b/src/main.cc index 9b7af79..ecda49a 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,3 +1,6 @@ +#include "vtkCell.h" +#include "vtkUnstructuredGrid.h" +#include #include #include #include @@ -16,8 +19,80 @@ #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()); +} + + +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); +} + int main(int argc, char **argv) { if (argc != 2) { @@ -33,6 +108,25 @@ int main(int argc, char **argv) { // vtkNew reader; vtkNew reader; reader->SetFileName(argv[1]); + reader->Update(); + vtkUnstructuredGrid *grid = reader->GetOutput(); + auto it = grid->NewCellIterator(); + std::vector angles(grid->GetNumberOfCells() * 12); + size_t i = 0; + for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) { + if (it->GetCellType() != VTK_TETRA) continue; + // double angles[12]; + cellAngles(grid, it->GetPointIds(), angles.data() + i); + i += 12; + // for (size_t i = 0; i < 12; i++) { + // std::cout << angles[i] << ", "; + // } + // std::cout << "\b\b \n"; + } + std::cout << "avg: " << average(angles) + << ", stddev: " << standard_deviation(angles) + << ", min: " << *std::min_element(angles.begin(), angles.end()) + << ", max: " << *std::max_element(angles.begin(), angles.end()) << std::endl; vtkNew volumeMapper; volumeMapper->SetInputConnection(reader->GetOutputPort()); @@ -42,9 +136,8 @@ int main(int argc, char **argv) { vtkNew volume; volume->SetMapper(volumeMapper); vtkNew transferFunction; - transferFunction->AddPoint(-1, 100); - transferFunction->AddPoint(0, .5); - // transferFunction->AddPoint(1, 1); + transferFunction->AddPoint(-1, 0); + transferFunction->AddPoint(1, 1); volume->GetProperty()->SetScalarOpacity(transferFunction); volume->GetProperty()->SetColor(transferFunction); // vtkNew actor;