From d0381b755207a484b98ea074de05ec8859ddf603 Mon Sep 17 00:00:00 2001
From: papush! <contact@ccolin.fr>
Date: Tue, 8 Mar 2022 14:59:15 +0100
Subject: [PATCH] use vtkDataSet::GetCellNeighbors to find surface points

---
 src/surface_points_filter.cc | 77 ++++++++++++++----------------------
 1 file changed, 30 insertions(+), 47 deletions(-)

diff --git a/src/surface_points_filter.cc b/src/surface_points_filter.cc
index e87f625..3b240e3 100644
--- a/src/surface_points_filter.cc
+++ b/src/surface_points_filter.cc
@@ -12,26 +12,6 @@
 
 
 
-/* Remove from ids_a the elements that are not also in ids_b. */
-static void keep_only_dupes(std::vector<vtkIdType> &ids_a,
-                            vtkIdType n_ids_b, const vtkIdType *ids_b) {
-	for (auto it = ids_a.begin(); it != ids_a.end();) {
-		bool duplicate = false;
-		for (vtkIdType i = 0; i < n_ids_b; i++) {
-			if (*it == ids_b[i]) {
-				duplicate = true;
-				break;
-			}
-		}
-		if (!duplicate) {
-			it = ids_a.erase(it);
-		} else {
-			++it;
-		}
-	}
-}
-
-
 vtkStandardNewMacro(SurfacePointsFilter);
 
 
@@ -50,44 +30,47 @@ vtkTypeBool SurfacePointsFilter::RequestData(
 	output->GetPointData()->PassData(input->GetPointData());
 	output->GetCellData()->PassData(input->GetCellData());
 
-	/* Create an array to store the IDs of external points. */
-	vtkNew<vtkIdTypeArray> externalPoints;
-	externalPoints->SetName("surface_points");
-	externalPoints->SetNumberOfComponents(1);
-	vtkNew<vtkIntArray> isExternal;
-	isExternal->SetName("bl");
-	isExternal->SetNumberOfComponents(1);
-	isExternal->SetNumberOfTuples(input->GetNumberOfPoints());
-	isExternal->FillValue(0);
+	/* Create an array to store the IDs of surface points. */
+	vtkNew<vtkIdTypeArray> surfacePoints;
+	surfacePoints->SetName("surface_points");
+	surfacePoints->SetNumberOfComponents(1);
+	vtkNew<vtkIntArray> isSurface;
+	isSurface->SetName("isSurface");
+	isSurface->SetNumberOfComponents(1);
+	isSurface->SetNumberOfTuples(input->GetNumberOfPoints());
+	isSurface->FillValue(0);
 
 	/* Generate cell links, these tell us which cell a point belongs to. */
 	vtkNew<vtkStaticCellLinks> links;
 	links->BuildLinks(input);
 
-	// TODO: try to use vtkDataSet::GetCellNeighbors instead
+	vtkNew<vtkIdList> neighborCells;
+	vtkNew<vtkIdList> facePoints;
 	auto *it = input->NewCellIterator();
 	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
-		vtkIdList *pointIds = it->GetPointIds();
+		vtkIdList *cellPointIds = it->GetPointIds();
 		for (vtkIdType faceId = 0; faceId < 4; faceId++) {
-			vtkIdType firstPoint = pointIds->GetId(faceId);
-			vtkIdType n_cells = links->GetNcells(firstPoint);
-			vtkIdType *cells = links->GetCells(firstPoint);
-			std::vector<vtkIdType> intersection(cells, cells + n_cells);
-			for (vtkIdType i = 1; i < 3; i++) {
-				vtkIdType pointId = pointIds->GetId((faceId + i) % 4);
-				keep_only_dupes(intersection,
-				                links->GetNcells(pointId), links->GetCells(pointId));
-			}
-			if (intersection.size() == 1) { // The face only belongs to one cell.
-				for (vtkIdType i = 0; i < 3; i++) {
-					externalPoints->InsertNextValue(pointIds->GetId((faceId + i) % 4));
-					isExternal->SetValue(pointIds->GetId((faceId + i) % 4), 1);
-				}
+			vtkIdType idA = cellPointIds->GetId((faceId + 0) % 4);
+			vtkIdType idB = cellPointIds->GetId((faceId + 1) % 4);
+			vtkIdType idC = cellPointIds->GetId((faceId + 2) % 4);
+			facePoints->InsertNextId(idA);
+			facePoints->InsertNextId(idB);
+			facePoints->InsertNextId(idC);
+			input->GetCellNeighbors(it->GetCellId(), facePoints,
+			                        neighborCells);
+			facePoints->Reset();
+			if (neighborCells->GetNumberOfIds() == 1) {
+				surfacePoints->InsertNextValue(idA);
+				isSurface->SetValue(idA, 1);
+				surfacePoints->InsertNextValue(idB);
+				isSurface->SetValue(idB, 1);
+				surfacePoints->InsertNextValue(idC);
+				isSurface->SetValue(idC, 1);
 			}
 		}
 	}
 
-	output->GetPointData()->SetScalars(isExternal);
-	output->GetFieldData()->AddArray(externalPoints);
+	output->GetPointData()->SetScalars(isSurface);
+	output->GetFieldData()->AddArray(surfacePoints);
 	return true;
 }