add external points filter

CMakeLists.txt

@@ -15,7 +15,9 @@ target_sources(pfe PRIVATE
   src/aspect_ratio_filter.cc
   src/aspect_ratio_filter.h
   src/dihedral_angles_filter.cc
-  src/dihedral_angles_filter.h)
+  src/dihedral_angles_filter.h
+  src/external_points_filter.cc
+  src/external_points_filter.h)
 
 target_link_libraries(pfe PRIVATE
   VTK::CommonCore

data/TetMesh2.vtk

@@ -18,6 +18,3 @@ LOOKUP_TABLE default
 1
 1
 1
-FIELD law 1
-cell_law 1 1 int
-0

src/external_points_filter.cc

@@ -0,0 +1,93 @@
+#include "external_points_filter.h"
+
+#include <vtkIdList.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkIdTypeArray.h>
+#include <vtkCellIterator.h>
+#include <vtkTriangle.h>
+#include <vtkStaticCellLinks.h>
+#include <vtkIntArray.h>
+
+
+
+/* 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(ExternalPointsFilter);
+
+
+vtkTypeBool ExternalPointsFilter::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* input =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkUnstructuredGrid* output =
+		vtkUnstructuredGrid::GetData(outputVector);
+
+	/* Copy point and cell data unmodified. */
+	output->CopyStructure(input);
+	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("external_points");
+	externalPoints->SetNumberOfComponents(1);
+	vtkNew<vtkIntArray> isExternal;
+	isExternal->SetName("bl");
+	isExternal->SetNumberOfComponents(1);
+	isExternal->SetNumberOfTuples(input->GetNumberOfPoints());
+	isExternal->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
+	auto *it = input->NewCellIterator();
+	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
+		vtkIdList *pointIds = 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);
+				}
+			}
+		}
+	}
+
+	output->GetPointData()->SetScalars(isExternal);
+	output->GetFieldData()->AddArray(externalPoints);
+	return true;
+}

src/external_points_filter.h

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

src/main.cc

@@ -1,6 +1,7 @@
 #include "angles_filter.h"
 #include "aspect_ratio_filter.h"
 #include "dihedral_angles_filter.h"
+#include "external_points_filter.h"
 
 #include <vtkCellData.h>
 #include <vtkUnstructuredGrid.h>
@@ -14,32 +15,10 @@
 #include <vtkVolume.h>
 #include <vtkOpenGLProjectedTetrahedraMapper.h>
 #include <vtkUnstructuredGridReader.h>
+#include <vtkUnstructuredGridWriter.h>
 #include <vtkPiecewiseFunction.h>
 
 #include <array>
-#include <vector>
-#include <algorithm>
-
-
-// template <typename T>
-// T average(const std::vector<T> &data) {
-// 	T avg = 0;
-// 	for (const T &t : data) {
-// 		avg += t;
-// 	}
-// 	return avg / data.size();
-// }
-
-
-// template <typename T>
-// T standard_deviation(const std::vector<T> &data) {
-// 	T avg = average(data);
-// 	T stddev = 0;
-// 	for (const T &t : data) {
-// 		stddev += (t - avg) * (t - avg);
-// 	}
-// 	return std::sqrt(stddev / data.size());
-// }
 
 
 int main(int argc, char **argv) {
@@ -58,74 +37,89 @@ int main(int argc, char **argv) {
 	anglesFilter->SetInputConnection(reader->GetOutputPort());
 
 	/* Display angles in console. */
-	anglesFilter->Update();
+	// anglesFilter->Update();
-	vtkUnstructuredGrid *grid = anglesFilter->GetOutput();
+	// vtkUnstructuredGrid *grid = anglesFilter->GetOutput();
-	auto *angles = grid->GetCellData()->GetArray("angles");
+	// auto *angles = grid->GetCellData()->GetArray("angles");
-	for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) {
+	// for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) {
-		std::cout << "cell " << i << ": ";
+	// 	std::cerr << "cell " << i << ": ";
-		for (ssize_t j = 0; j < angles->GetNumberOfComponents(); j++) {
+	// 	for (ssize_t j = 0; j < angles->GetNumberOfComponents(); j++) {
-			std::cout << angles->GetTuple(i)[j] << ", ";
+	// 		std::cerr << angles->GetTuple(i)[j] << ", ";
-		}
+	// 	}
-		std::cout << "\b\b \n";
+	// 	std::cerr << "\b\b \n";
-	}
+	// }
 
 	/* Aspect ratio */
 	vtkNew<AspectRatioFilter> aspectRatioFilter;
 	aspectRatioFilter->SetInputConnection(anglesFilter->GetOutputPort());
 
 	/* Display aspect ratios in console. */
-	aspectRatioFilter->Update();
+	// aspectRatioFilter->Update();
-	grid = aspectRatioFilter->GetOutput();
+	// grid = aspectRatioFilter->GetOutput();
-	auto *aspectRatios = grid->GetCellData()->GetArray("aspect_ratio");
+	// auto *aspectRatios = grid->GetCellData()->GetArray("aspect_ratio");
-	for (ssize_t i = 0; i < aspectRatios->GetNumberOfTuples(); i++) {
+	// for (ssize_t i = 0; i < aspectRatios->GetNumberOfTuples(); i++) {
-		std::cout << "cell " << i << ": "
+	// 	std::cerr << "cell " << i << ": "
-				  << aspectRatios->GetTuple1(i) << "\n";
+	// 			  << aspectRatios->GetTuple1(i) << "\n";
-	}
+	// }
 
 	/* Dihedral angles filter */
 	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
 	dihedralAnglesFilter->SetInputConnection(aspectRatioFilter->GetOutputPort());
 
 	/* Display dihedral angles in console. */
-	dihedralAnglesFilter->Update();
+	// dihedralAnglesFilter->Update();
-	grid = dihedralAnglesFilter->GetOutput();
+	// grid = dihedralAnglesFilter->GetOutput();
-	auto dihedralAngles = grid->GetCellData()->GetArray("dihedral_angles");
+	// auto dihedralAngles = grid->GetCellData()->GetArray("dihedral_angles");
-	for (vtkIdType i = 0; i < grid->GetNumberOfCells(); i++) {
+	// for (vtkIdType i = 0; i < grid->GetNumberOfCells(); i++) {
-		std::cout << "dihedral angles ";
+	// 	std::cerr << "dihedral angles ";
-		for (vtkIdType j = 0; j < 6; j++) {
+	// 	for (vtkIdType j = 0; j < 6; j++) {
-			std::cout << dihedralAngles->GetTuple(i)[j] << ", ";
+	// 		std::cerr << dihedralAngles->GetTuple(i)[j] << ", ";
-		}
+	// 	}
-		std::cout << "\b\b\n";
+	// 	std::cerr << "\b\b\n";
-	}
+	// }
 
-	/* Volume rendering properties */
+	/* External faces filter. */
-	vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
+	// vtkNew<vtkmExternalFaces> externalFacesFilter;
-	volumeMapper->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
+	// externalFacesFilter->DebugOn();
-	vtkNew<vtkVolume> volume;
+	// externalFacesFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
-	volume->SetMapper(volumeMapper);
-	vtkNew<vtkPiecewiseFunction> transferFunction;
-	transferFunction->AddPoint(-1, 0);
-	transferFunction->AddPoint(1, 1);
-	volume->GetProperty()->SetScalarOpacity(transferFunction);
-	volume->GetProperty()->SetColor(transferFunction);
 
-	/* Renderer */
+	/* External points filter. */
-	vtkNew<vtkNamedColors> colors;
+	vtkNew<ExternalPointsFilter> externalPointsFilter;
-	std::array<unsigned char, 4> bkg{{26, 51, 102, 255}};
+	externalPointsFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
-	colors->SetColor("BkgColor", bkg.data());
+
-	vtkNew<vtkRenderer> renderer;
+	vtkNew<vtkUnstructuredGridWriter> writer;
-	renderer->AddVolume(volume);
+	writer->SetInputConnection(externalPointsFilter->GetOutputPort());
-	renderer->SetBackground(colors->GetColor3d("BkgColor").GetData());
+	writer->SetFileTypeToASCII();
-	renderer->ResetCamera();
+	writer->SetFileName("out.vtk");
-	renderer->GetActiveCamera()->Zoom(1.5);
+	writer->Write();
-	vtkNew<vtkRenderWindow> renderWindow;
+
-	renderWindow->SetSize(300, 300);
+	// /* Volume rendering properties */
-	renderWindow->AddRenderer(renderer);
+	// vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
-	renderWindow->SetWindowName("PFE");
+	// volumeMapper->SetInputConnection(externalPointsFilter->GetOutputPort());
-	vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
+	// vtkNew<vtkVolume> volume;
-	renderWindowInteractor->SetRenderWindow(renderWindow);
+	// volume->SetMapper(volumeMapper);
-	renderWindow->Render();
+	// vtkNew<vtkPiecewiseFunction> transferFunction;
-	renderWindowInteractor->Start();
+	// transferFunction->AddPoint(-1, 0);
+	// transferFunction->AddPoint(1, 1);
+	// volume->GetProperty()->SetScalarOpacity(transferFunction);
+	// volume->GetProperty()->SetColor(transferFunction);
+
+	// /* Renderer */
+	// vtkNew<vtkNamedColors> colors;
+	// std::array<unsigned char, 4> bkg{{26, 51, 102, 255}};
+	// colors->SetColor("BkgColor", bkg.data());
+	// vtkNew<vtkRenderer> renderer;
+	// renderer->AddVolume(volume);
+	// renderer->SetBackground(colors->GetColor3d("BkgColor").GetData());
+	// renderer->ResetCamera();
+	// renderer->GetActiveCamera()->Zoom(1.5);
+	// vtkNew<vtkRenderWindow> renderWindow;
+	// renderWindow->SetSize(300, 300);
+	// renderWindow->AddRenderer(renderer);
+	// renderWindow->SetWindowName("PFE");
+	// vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
+	// renderWindowInteractor->SetRenderWindow(renderWindow);
+	// renderWindow->Render();
+	// renderWindowInteractor->Start();
 
 	return EXIT_SUCCESS;
 }