:)

CMakeLists.txt

@@ -12,6 +12,7 @@ set(VTK_COMPONENTS
   VTK::IOGeometry
   VTK::IOXML
   VTK::FiltersModeling
+  VTK::FiltersPoints
   VTK::vtksys)
 set(ENABLE_VIEWER OFF CACHE BOOL "Enable the 3D viewer, depends on Qt.")
 if(ENABLE_VIEWER)

src/main.cc

@@ -15,6 +15,7 @@
 #include <vtkXMLUnstructuredGridReader.h>
 #include <vtkXMLUnstructuredGridWriter.h>
 #include <vtksys/SystemTools.hxx>
+#include <vtkPointSmoothingFilter.h>
 
 #ifdef USE_VIEWER
 #include <vtkCamera.h>
@@ -82,8 +83,11 @@ int main(int argc, char **argv) {
 	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
 	dihedralAnglesFilter->SetInputConnection(project->GetOutputPort());
 
+	vtkNew<vtkPointSmoothingFilter> pointSmoothingFilter;
+	pointSmoothingFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPoort());
+
 	vtkNew<vtkXMLUnstructuredGridWriter> writer;
-	writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+	writer->SetInputConnection(pointSmoothingFilter->GetOutputPort());
 	writer->SetDataModeToAscii();
 	writer->SetFileName("out.vtu");
 	writer->Write();

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++) {
+	applyMotionVectors(tetMesh, surfacePoints,
-		bool skip = false;
+	                   motionVectors, numberOfAffectors);
-		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();
-		}
-	}
 
 	return true;
 }