From e789941169da5cfbdb58a3ee538329faebc3631e Mon Sep 17 00:00:00 2001
From: papush! <contact@ccolin.fr>
Date: Thu, 10 Mar 2022 16:49:49 +0100
Subject: [PATCH] improve dihedral angles filter

---
 src/dihedral_angles_filter.cc | 12 ++++++++++++
 src/main.cc                   | 10 +++++-----
 2 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/src/dihedral_angles_filter.cc b/src/dihedral_angles_filter.cc
index 73a08c2..60581ee 100644
--- a/src/dihedral_angles_filter.cc
+++ b/src/dihedral_angles_filter.cc
@@ -24,11 +24,18 @@ vtkTypeBool DihedralAnglesFilter::RequestData(
 	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);
+
+	vtkNew<vtkDoubleArray> minDihedralAngleArray;
+	minDihedralAngleArray->SetName("min_dihedral_angle");
+	minDihedralAngleArray->SetNumberOfTuples(input->GetNumberOfCells());
+	double *minDihedralAngleBase = minDihedralAngleArray->GetPointer(0);
+
 	size_t i = 0;
 	vtkCellIterator *it = input->NewCellIterator();
 	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
@@ -52,8 +59,13 @@ vtkTypeBool DihedralAnglesFilter::RequestData(
 		dihedralAnglesBase[i+4] = vtkMath::AngleBetweenVectors(nbdc, nadc);
 		dihedralAnglesBase[i+5] = vtkMath::AngleBetweenVectors(nbdc, nabd);
 
+		minDihedralAngleBase[i / 6] =
+			*std::min_element(dihedralAnglesBase + i,
+			                  dihedralAnglesBase + i + 6);
+
 		i += 6;
 	}
 	cellData->AddArray((vtkAbstractArray *) dihedralAnglesArray);
+	cellData->AddArray((vtkAbstractArray *) minDihedralAngleArray);
 	return true;
 }
diff --git a/src/main.cc b/src/main.cc
index 09833d7..6d54bcb 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -65,12 +65,9 @@ int main(int argc, char **argv) {
 	auto tetMeshReader = readerFromFileName(argv[1]);
 	auto polyMeshReader = readerFromFileName(argv[2]);
 
-	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
-	dihedralAnglesFilter->SetInputConnection(tetMeshReader->GetOutputPort());
-
 	vtkNew<RemoveExternalCellsFilter> removeExternalCellsFilter;
 	removeExternalCellsFilter->SetInputConnection(0,
-		dihedralAnglesFilter->GetOutputPort());
+		tetMeshReader->GetOutputPort());
 	removeExternalCellsFilter->SetInputConnection(1,
 		polyMeshReader->GetOutputPort());
 
@@ -82,8 +79,11 @@ int main(int argc, char **argv) {
 	project->SetInputConnection(0, surfacePointsFilter->GetOutputPort());
 	project->SetInputConnection(1, polyMeshReader->GetOutputPort());
 
+	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
+	dihedralAnglesFilter->SetInputConnection(project->GetOutputPort());
+
 	vtkNew<vtkXMLUnstructuredGridWriter> writer;
-	writer->SetInputConnection(project->GetOutputPort());
+	writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
 	writer->SetDataModeToAscii();
 	writer->SetFileName("out.vtu");
 	writer->Write();