From d31d9629f014e5ef954d5dcfabc190a4b39d0bf9 Mon Sep 17 00:00:00 2001
From: papush! <contact@ccolin.fr>
Date: Tue, 15 Mar 2022 14:55:36 +0100
Subject: [PATCH] add relaxation (not very relaxing)

---
 src/main.cc                           |  7 +++-
 src/project_surface_points_on_poly.cc |  6 ++--
 src/relaxation_filter.cc              | 50 +++++++++++++++++++++++++--
 3 files changed, 57 insertions(+), 6 deletions(-)

diff --git a/src/main.cc b/src/main.cc
index 6d54bcb..f6ef659 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -2,6 +2,7 @@
 #include "remove_external_cells_filter.h"
 #include "surface_points_filter.h"
 #include "project_surface_points_on_poly.h"
+#include "relaxation_filter.h"
 
 #include <vtkCellArrayIterator.h>
 #include <vtkCellData.h>
@@ -79,8 +80,12 @@ int main(int argc, char **argv) {
 	project->SetInputConnection(0, surfacePointsFilter->GetOutputPort());
 	project->SetInputConnection(1, polyMeshReader->GetOutputPort());
 
+	vtkNew<RelaxationFilter> relaxationFilter;
+	relaxationFilter->SetInputConnection(project->GetOutputPort());
+
 	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
-	dihedralAnglesFilter->SetInputConnection(project->GetOutputPort());
+	dihedralAnglesFilter->SetInputConnection(
+		relaxationFilter->GetOutputPort());
 
 	vtkNew<vtkXMLUnstructuredGridWriter> writer;
 	writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc
index c650dde..6502ce3 100644
--- a/src/project_surface_points_on_poly.cc
+++ b/src/project_surface_points_on_poly.cc
@@ -154,10 +154,10 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 	vtkInformationVector **inputVector,
 	vtkInformationVector *outputVector) {
 	(void) request;
-	vtkUnstructuredGrid* tetMesh = vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkUnstructuredGrid* tetMesh =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
 	vtkPolyData* polyMesh = vtkPolyData::GetData(inputVector[1]);
-	vtkUnstructuredGrid* output =
-		vtkUnstructuredGrid::GetData(outputVector);
+	vtkUnstructuredGrid* output = vtkUnstructuredGrid::GetData(outputVector);
 
 	output->CopyStructure(tetMesh);
 	output->GetPointData()->PassData(tetMesh->GetPointData());
diff --git a/src/relaxation_filter.cc b/src/relaxation_filter.cc
index d70b0a1..13c2011 100644
--- a/src/relaxation_filter.cc
+++ b/src/relaxation_filter.cc
@@ -6,13 +6,59 @@
 #include <vtkDoubleArray.h>
 #include <vtkCellIterator.h>
 
+#include <set>
+
 vtkStandardNewMacro(RelaxationFilter);
 
 vtkTypeBool RelaxationFilter::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());
+	output->GetCellData()->PassData(input->GetCellData());
+	vtkNew<vtkPoints> newPoints;
+	newPoints->DeepCopy(output->GetPoints());
+	vtkNew<vtkIdList> neighborCells;
+	vtkNew<vtkIdList> cellPoints;
+	vtkIntArray *isSurface = vtkIntArray::SafeDownCast(
+		output->GetPointData()->GetArray("isSurface"));
+	std::set<vtkIdType> neighbors;
+	output->BuildLinks();
+	double avg[3];
+	for (vtkIdType i = 0; i < output->GetNumberOfPoints(); i++) {
+		output->GetPointCells(i, neighborCells);
+		if (neighborCells->GetNumberOfIds() != 0) {
+			for (vtkIdType j = 0; j < neighborCells->GetNumberOfIds(); j++) {
+				vtkIdType cellId = neighborCells->GetId(j);
+				output->GetCellPoints(cellId, cellPoints);
+				for (vtkIdType k = 0; k < cellPoints->GetNumberOfIds(); k++) {
+					vtkIdType neighborId = cellPoints->GetId(k);
+					if (isSurface->GetValue(neighborId)) {
+						neighbors.insert(neighborId);
+					}
+				}
+				cellPoints->Reset();
+			}
+			neighbors.erase(neighbors.find(i));
+			if (neighbors.size() != 0) {
+				avg[0] = 0; avg[1] = 0; avg[2] = 0;
+				for (const vtkIdType &j : neighbors) {
+					double point[3];
+					output->GetPoints()->GetPoint(j, point);
+					vtkMath::Add(point, avg, avg);
+				}
+				vtkMath::MultiplyScalar(avg, 1. / neighbors.size());
+				newPoints->SetPoint(i, avg);
+			}
+		}
+		neighborCells->Reset();
+		neighbors.clear();
+	}
+	output->SetPoints(newPoints);
 	return true;
 }