From 62e4cd16296c38f9b2a090e7faf26fd9fe8362ac Mon Sep 17 00:00:00 2001
From: papush! <contact@ccolin.fr>
Date: Mon, 14 Mar 2022 14:30:06 +0100
Subject: [PATCH] fix projection

---
 src/project_surface_points_on_poly.cc | 63 ++++++++++++++++++---------
 1 file changed, 42 insertions(+), 21 deletions(-)

diff --git a/src/project_surface_points_on_poly.cc b/src/project_surface_points_on_poly.cc
index 9e6ebca..c650dde 100644
--- a/src/project_surface_points_on_poly.cc
+++ b/src/project_surface_points_on_poly.cc
@@ -71,6 +71,19 @@ static void closestPolyMeshPoint(vtkPolyData *polyMesh,
 }
 
 
+void findPointsWithinRadius(double radius, vtkPointSet *pointSet,
+                            const double *point, vtkIdList *points) {
+	for (vtkIdType i = 0; i < pointSet->GetNumberOfPoints(); i++) {
+		double vec[3];
+		vtkMath::Subtract(point, pointSet->GetPoint(i), vec);
+		double dist = vtkMath::Norm(vec);
+		if (dist < radius) {
+			points->InsertNextId(i);
+		}
+	}
+}
+
+
 static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
                              vtkPolyData *polyMesh,
                              vtkIdType pointId,
@@ -87,14 +100,15 @@ static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
 
 	closestPolyMeshPoint(polyMesh, point, kdTree, links, direction, &distance);
 	vtkNew<vtkIdList> affectedPoints;
-	kdTree->FindPointsWithinRadius(distance * radiusScale, point, affectedPoints);
-	
+	findPointsWithinRadius(distance * radiusScale, tetMesh, point, affectedPoints);
+	// kdTree->FindPointsWithinRadius(distance * radiusScale, point, affectedPoints);
+
 	if(distance > 0) for (vtkIdType j = 0; j < affectedPoints->GetNumberOfIds(); j++) {
 		vtkIdType affectedPointId = affectedPoints->GetId(j);
 		double affectedPoint[3];
 		tetMesh->GetPoint(affectedPointId, affectedPoint);
 		double dist2 = vtkMath::Distance2BetweenPoints(affectedPoint, point);
-		if(dist2 > (distance * radiusScale) * (distance * radiusScale)) continue;
+		// if(dist2 > (distance * radiusScale) * (distance * radiusScale)) continue;
 
 		double motion[3] = {direction[0], direction[1], direction[2]};
 		double factor = 1. - std::sqrt(dist2) / (distance * radiusScale);
@@ -112,6 +126,29 @@ static void moveSurfacePoint(vtkUnstructuredGrid *tetMesh,
 }
 
 
+static void applyMotionVectors(vtkUnstructuredGrid *tetMesh,
+                               vtkIdTypeArray *surfacePoints,
+                               std::vector<double> &motionVectors,
+                               const std::vector<unsigned> &numberOfAffectors) {
+	for (vtkIdType i = 0; i < tetMesh->GetNumberOfPoints(); i++) {
+		bool skip = false;
+		for (vtkIdType j = 0; j < surfacePoints->GetNumberOfTuples(); j++) {
+			if (i == surfacePoints->GetValue(j)) skip = true;
+		}
+		if (skip) continue;
+		if(numberOfAffectors[i] == 0) continue;
+		double point[3];
+		tetMesh->GetPoint(i, point);
+		double *motion = motionVectors.data() + (i * 3);
+
+		vtkMath::MultiplyScalar(motion, 1. / numberOfAffectors[i]);
+		vtkMath::Subtract(point, motion, point);
+		tetMesh->GetPoints()->SetPoint(i, point);
+		tetMesh->Modified();
+	}
+}
+
+
 vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 	vtkInformation *request,
 	vtkInformationVector **inputVector,
@@ -144,24 +181,8 @@ vtkTypeBool ProjectSurfacePointsOnPoly::RequestData(
 		                 numberOfAffectors,
 		                 kdTree, links);
 	}
-	for (vtkIdType i = 0; i < tetMesh->GetNumberOfPoints(); i++) {
-		bool skip = false;
-		for (vtkIdType j = 0; j < surfacePoints->GetNumberOfTuples(); j++) {
-			if (i == surfacePoints->GetValue(j)) skip = true;
-		}
-		if (skip) continue;
-		
-		if(numberOfAffectors[i] != 0) {
-			double point[3];
-			tetMesh->GetPoint(i, point);
-			double *motion = motionVectors.data() + (i * 3);
-
-			vtkMath::MultiplyScalar(motion, 1. / numberOfAffectors[i]);
-			vtkMath::Subtract(point, motion, point);
-			tetMesh->GetPoints()->SetPoint(i, point);
-			tetMesh->Modified();
-		}
-	}
+	applyMotionVectors(tetMesh, surfacePoints,
+	                   motionVectors, numberOfAffectors);
 
 	return true;
 }