diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d2cc21..aed3ad2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,7 +13,9 @@ target_sources(pfe PRIVATE src/angles_filter.cc src/angles_filter.h src/aspect_ratio_filter.cc - src/aspect_ratio_filter.h) + src/aspect_ratio_filter.h + src/dihedral_angles_filter.cc + src/dihedral_angles_filter.h) target_link_libraries(pfe PRIVATE VTK::CommonCore diff --git a/src/dihedral_angles_filter.cc b/src/dihedral_angles_filter.cc new file mode 100644 index 0000000..73a08c2 --- /dev/null +++ b/src/dihedral_angles_filter.cc @@ -0,0 +1,59 @@ +#include "dihedral_angles_filter.h" + +#include +#include +#include +#include +#include +#include + + +vtkStandardNewMacro(DihedralAnglesFilter); + + +vtkTypeBool DihedralAnglesFilter::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 dihedralAnglesArray; + dihedralAnglesArray->SetName("dihedral_angles"); + dihedralAnglesArray->SetNumberOfComponents(6); + dihedralAnglesArray->SetNumberOfTuples(input->GetNumberOfCells()); + double *dihedralAnglesBase = dihedralAnglesArray->GetPointer(0); + size_t i = 0; + vtkCellIterator *it = input->NewCellIterator(); + for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) { + double a[3], b[3], c[3], d[3]; + vtkPoints *points = it->GetPoints(); + points->GetPoint(0, a); + points->GetPoint(1, b); + points->GetPoint(2, c); + points->GetPoint(3, d); + + double nacb[3], nadc[3], nabd[3], nbdc[3]; + vtkTriangle::ComputeNormal(a, c, b, nacb); + vtkTriangle::ComputeNormal(a, d, c, nadc); + vtkTriangle::ComputeNormal(a, b, d, nabd); + vtkTriangle::ComputeNormal(b, d, c, nbdc); + + dihedralAnglesBase[i+0] = vtkMath::AngleBetweenVectors(nacb, nadc); + dihedralAnglesBase[i+1] = vtkMath::AngleBetweenVectors(nadc, nabd); + dihedralAnglesBase[i+2] = vtkMath::AngleBetweenVectors(nabd, nacb); + dihedralAnglesBase[i+3] = vtkMath::AngleBetweenVectors(nbdc, nacb); + dihedralAnglesBase[i+4] = vtkMath::AngleBetweenVectors(nbdc, nadc); + dihedralAnglesBase[i+5] = vtkMath::AngleBetweenVectors(nbdc, nabd); + + i += 6; + } + cellData->AddArray((vtkAbstractArray *) dihedralAnglesArray); + return true; +} diff --git a/src/dihedral_angles_filter.h b/src/dihedral_angles_filter.h new file mode 100644 index 0000000..2a441fb --- /dev/null +++ b/src/dihedral_angles_filter.h @@ -0,0 +1,21 @@ +#ifndef DIHEDRAL_ANGLES_FILTER_H +#define DIHEDRAL_ANGLES_FILTER_H + +#include +#include + + +class DihedralAnglesFilter : public vtkUnstructuredGridAlgorithm { +public: + static DihedralAnglesFilter *New(); + vtkTypeMacro(DihedralAnglesFilter, vtkUnstructuredGridAlgorithm); + vtkTypeBool RequestData(vtkInformation *request, + vtkInformationVector **inputVector, + vtkInformationVector *outputVector) override; + +protected: + ~DihedralAnglesFilter() override = default; +}; + + +#endif \ No newline at end of file diff --git a/src/main.cc b/src/main.cc index f6c1720..bcf8a47 100644 --- a/src/main.cc +++ b/src/main.cc @@ -1,5 +1,6 @@ #include "angles_filter.h" #include "aspect_ratio_filter.h" +#include "dihedral_angles_filter.h" #include #include @@ -52,7 +53,7 @@ int main(int argc, char **argv) { reader->SetFileName(argv[1]); - /* Angle filter */ + /* Angles filter */ vtkNew anglesFilter; anglesFilter->SetInputConnection(reader->GetOutputPort()); @@ -81,9 +82,25 @@ int main(int argc, char **argv) { << 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"; + } + /* Volume rendering properties */ vtkNew volumeMapper; - volumeMapper->SetInputConnection(aspectRatioFilter->GetOutputPort()); + volumeMapper->SetInputConnection(dihedralAnglesFilter->GetOutputPort()); vtkNew volume; volume->SetMapper(volumeMapper); vtkNew transferFunction;