add dihedral angles filter

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

src/dihedral_angles_filter.cc

@@ -0,0 +1,59 @@
+#include "dihedral_angles_filter.h"
+
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkDoubleArray.h>
+#include <vtkCellIterator.h>
+#include <vtkTriangle.h>
+
+
+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<vtkDoubleArray> 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;
+}

src/dihedral_angles_filter.h

@@ -0,0 +1,21 @@
+#ifndef DIHEDRAL_ANGLES_FILTER_H
+#define DIHEDRAL_ANGLES_FILTER_H
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+
+
+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

src/main.cc

@@ -1,5 +1,6 @@
 #include "angles_filter.h"
 #include "aspect_ratio_filter.h"
+#include "dihedral_angles_filter.h"
 
 #include <vtkCellData.h>
 #include <vtkUnstructuredGrid.h>
@@ -52,7 +53,7 @@ int main(int argc, char **argv) {
 	reader->SetFileName(argv[1]);
 
 
-	/* Angle filter */
+	/* Angles filter */
 	vtkNew<AnglesFilter> 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;
+	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<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
-	volumeMapper->SetInputConnection(aspectRatioFilter->GetOutputPort());
+	volumeMapper->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
 	vtkNew<vtkVolume> volume;
 	volume->SetMapper(volumeMapper);
 	vtkNew<vtkPiecewiseFunction> transferFunction;