add external point mapping

CMakeLists.txt

@@ -48,13 +48,14 @@ target_sources(pfe PRIVATE
   src/dihedral_angles_filter.h
   src/external_points_filter.cc
   src/external_points_filter.h
-  
   src/kd_tree.cc
   src/kd_tree.h
   src/mesh_fit_filter.cc
   src/mesh_fit_filter.h
   src/point_tris_dist.cc
-  src/point_tris_dist.h)
+  src/point_tris_dist.h
+  src/project_external_points_on_surface.cc
+  src/project_external_points_on_surface.h)
 
 target_link_libraries(pfe PRIVATE ${VTK_COMPONENTS})
 

src/main.cc

@@ -2,7 +2,7 @@
 #include "aspect_ratio_filter.h"
 #include "dihedral_angles_filter.h"
 #include "external_points_filter.h"
-
+#include "project_external_points_on_surface.h"
 #include "mesh_fit_filter.h"
 #include "point_tris_dist.h"
 
@@ -12,6 +12,7 @@
 #include <vtkUnstructuredGridWriter.h>
 #include <vtkPolyData.h>
 #include <vtkCellArrayIterator.h>
+#include <vtkXMLPolyDataReader.h>
 
 #ifdef USE_VIEWER
 #include <vtkNamedColors.h>
@@ -30,8 +31,8 @@
 
 
 int main(int argc, char **argv) {
-	if (argc != 2) {
-		std::cerr << "Usage: " << argv[0] << " FILE.vtk" << std::endl;
+	if (argc != 3) {
+		std::cerr << "Usage: " << argv[0] << " tetmesh polydata" << std::endl;
 		return EXIT_FAILURE;
 	}
 
@@ -95,49 +96,56 @@ int main(int argc, char **argv) {
 	externalPointsFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
 
 
-	vtkNew<MeshFitFilter> meshFitFilter;
-	meshFitFilter->SetInputConnection(externalPointsFilter->GetOutputPort());
+	// vtkNew<MeshFitFilter> meshFitFilter;
+	// meshFitFilter->SetInputConnection(externalPointsFilter->GetOutputPort());
+
+	vtkNew<vtkXMLPolyDataReader> pdReader;
+	pdReader->SetFileName(argv[2]);
+
+	vtkNew<ProjectExternalPointsOnSurface> pepouse;
+	pepouse->SetInputConnection(0, externalPointsFilter->GetOutputPort());
+	pepouse->SetInputConnection(1, pdReader->GetOutputPort());
 
 
 	vtkNew<vtkUnstructuredGridWriter> writer;
-	writer->SetInputConnection(meshFitFilter->GetOutputPort());
+	writer->SetInputConnection(pepouse->GetOutputPort());
 	writer->SetFileTypeToASCII();
 	writer->SetFileName("out.vtk");
 	writer->Write();
 
 
-	vtkNew<vtkPoints> points;
-	points->InsertPoint(0, 0., 0., 0.);
-	points->InsertPoint(1, 1., 0., 0.);
-	points->InsertPoint(2, 1., 1., 0.);
+	// vtkNew<vtkPoints> points;
+	// points->InsertPoint(0, 0., 0., 0.);
+	// points->InsertPoint(1, 1., 0., 0.);
+	// points->InsertPoint(2, 1., 1., 0.);
 
-	vtkNew<vtkCellArray> strips;
-	strips->InsertNextCell(3);
-	strips->InsertCellPoint(0);
-	strips->InsertCellPoint(1);
-	strips->InsertCellPoint(2);
+	// vtkNew<vtkCellArray> strips;
+	// strips->InsertNextCell(3);
+	// strips->InsertCellPoint(0);
+	// strips->InsertCellPoint(1);
+	// strips->InsertCellPoint(2);
 
-	vtkNew<vtkPolyData> profile;
-	profile->SetPoints(points);
-	profile->SetStrips(strips);
+	// vtkNew<vtkPolyData> profile;
+	// profile->SetPoints(points);
+	// profile->SetStrips(strips);
 
-	vtkCell *cell = profile->GetCell(0);
-	// vtkTriangle *triangle = vtkTriangle::SafeDownCast(cell);
-	// vtkTriangle *triangle = (vtkTriangle *) cell;
+	// vtkCell *cell = profile->GetCell(0);
+	// // vtkTriangle *triangle = vtkTriangle::SafeDownCast(cell);
+	// // vtkTriangle *triangle = (vtkTriangle *) cell;
 
-	double direction[3];// = {0, 0, 0};
-	double point[3] = {0, 1, 0};
+	// double direction[3];// = {0, 0, 0};
+	// double point[3] = {0, 1, 0};
 
-	double d = pointTriangleDistance(point, cell, direction);
-	//double d = 5;
+	// double d = pointTriangleDistance(point, cell, direction);
+	// //double d = 5;
 
-	std::cout << "[" << point[0]
-			  << ", " << point[1]
-			  << ", " << point[2]
-			  << "] (" << d << ")"
-			  <<  "[" << direction[0]
-			  << ", " << direction[1]
-			  << ", " << direction[2] << "]\n";
+	// std::cout << "[" << point[0]
+	// 		  << ", " << point[1]
+	// 		  << ", " << point[2]
+	// 		  << "] (" << d << ")"
+	// 		  <<  "[" << direction[0]
+	// 		  << ", " << direction[1]
+	// 		  << ", " << direction[2] << "]\n";
 
 #ifdef USE_VIEWER
 	/* Volume rendering properties */

src/project_external_points_on_surface.cc

@@ -0,0 +1,93 @@
+#include "point_tris_dist.h"
+
+#include "project_external_points_on_surface.h"
+
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkDoubleArray.h>
+#include <vtkCellIterator.h>
+#include <vtkInformation.h>
+#include <vtkInformationVector.h>
+#include <vtkPolyData.h>
+#include <vtkKdTree.h>
+#include <vtkStaticCellLinks.h>
+
+#include <limits>
+
+
+vtkStandardNewMacro(ProjectExternalPointsOnSurface);
+
+
+ProjectExternalPointsOnSurface::ProjectExternalPointsOnSurface() {
+	SetNumberOfInputPorts(2);
+}
+
+
+int ProjectExternalPointsOnSurface::FillInputPortInformation(int port, vtkInformation *info) {
+	if (port == 0) {
+		vtkUnstructuredGridAlgorithm::FillInputPortInformation(port, info);
+	} else if (port == 1) {
+		vtkNew<vtkInformation> input2;
+		input2->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData");
+		info->Append(input2);
+	}
+	return 1;
+}
+
+
+vtkTypeBool ProjectExternalPointsOnSurface::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* us = vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkPolyData* pd = vtkPolyData::GetData(inputVector[1]);
+	vtkUnstructuredGrid* output =
+		vtkUnstructuredGrid::GetData(outputVector);
+
+	output->CopyStructure(us);
+	output->GetPointData()->PassData(us->GetPointData());
+	output->GetCellData()->PassData(us->GetCellData());
+
+	/* Generate cell links, these tell us which cell a point belongs to. */
+	vtkNew<vtkStaticCellLinks> links;
+	links->BuildLinks(pd);
+
+	vtkNew<vtkKdTree> kdTree;
+	kdTree->BuildLocatorFromPoints(pd->GetPoints());
+
+	vtkIdTypeArray *externalPoints =
+		vtkIdTypeArray::SafeDownCast(us->GetFieldData()->GetArray("external_points"));
+	for (vtkIdType i = 0; i < externalPoints->GetNumberOfTuples(); i++) {
+		double p[3];
+		us->GetPoint(externalPoints->GetTuple1(i), p);
+		double dist2;
+		vtkIdType closest = kdTree->FindClosestPoint(p, dist2);
+		vtkIdType *cellIds = links->GetCells(closest);
+		vtkIdType nCells = links->GetNcells(closest);
+		double min_dist = std::numeric_limits<double>::infinity();
+		double min_direction[3];
+
+		/* For each tri the closest point belongs to. */
+		for (vtkIdType j = 0; j < nCells; j++) {
+			vtkIdType cellId = cellIds[j];
+			double direction[3];
+			double dist = pointTriangleDistance(p, pd->GetCell(cellId), direction);
+			/* Find the closest one. */
+			if (dist < min_dist) {
+				min_dist = dist;
+				min_direction[0] = direction[0];
+				min_direction[1] = direction[1];
+				min_direction[2] = direction[2];
+			}
+		}
+
+		vtkMath::MultiplyScalar(min_direction, min_dist);
+		vtkMath::Subtract(p, min_direction, p);
+		us->GetPoints()->SetPoint(externalPoints->GetTuple1(i), p);
+		us->Modified();
+	}
+
+	return true;
+}

src/project_external_points_on_surface.h

@@ -0,0 +1,23 @@
+#ifndef PROJECT_EXTERNAL_POINTS_ON_SURFACE_H
+#define PROJECT_EXTERNAL_POINTS_ON_SURFACE_H
+
+#include <vtkUnstructuredGridAlgorithm.h>
+#include <vtkIdList.h>
+
+
+class ProjectExternalPointsOnSurface : public vtkUnstructuredGridAlgorithm {
+public:
+	static ProjectExternalPointsOnSurface *New();
+	vtkTypeMacro(ProjectExternalPointsOnSurface, vtkUnstructuredGridAlgorithm);
+	ProjectExternalPointsOnSurface();
+	int FillInputPortInformation(int, vtkInformation *info) override;
+	vtkTypeBool RequestData(vtkInformation *request,
+							vtkInformationVector **inputVector,
+							vtkInformationVector *outputVector) override;
+
+protected:
+	~ProjectExternalPointsOnSurface() override = default;
+};
+
+
+#endif