implement basic angle computation and analysis

CMakeLists.txt

@@ -6,6 +6,8 @@ add_subdirectory(external/vtk)
 
 add_executable(pfe)
 
+target_compile_features(pfe PRIVATE cxx_std_17)
+
 target_sources(pfe PRIVATE src/main.cc)
 target_link_libraries(pfe PRIVATE
   VTK::CommonCore

data/TetMesh2.vtk

@@ -0,0 +1,23 @@
+# vtk DataFile Version 3.0
+vtk output
+ASCII
+DATASET UNSTRUCTURED_GRID
+POINTS 4 float
+0 0 0
+1 0 0
+0 1 0
+0 0 1
+CELLS 1 5
+4 0 1 2 3
+CELL_TYPES 1
+10
+POINT_DATA 4
+SCALARS Indication float
+LOOKUP_TABLE default
+1
+1
+1
+1
+FIELD law 1
+cell_law 1 1 int
+0

src/main.cc

@@ -1,3 +1,6 @@
+#include "vtkCell.h"
+#include "vtkUnstructuredGrid.h"
+#include <algorithm>
 #include <vtkActor.h>
 #include <vtkCamera.h>
 #include <vtkCylinderSource.h>
@@ -16,8 +19,80 @@
 #include <vtkPolyDataReader.h>
 #include <vtkXMLPolyDataReader.h>
 #include <vtkPiecewiseFunction.h>
+#include <vtkCellIterator.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());
+}
+
+
+void cellAngles(vtkDataSet *dataSet, vtkIdList *idList, double *angles) {
+	// std::cout << "nb points: " << idList->GetNumberOfIds() << std::endl;
+	double a[3], b[3], c[3], d[3];
+	dataSet->GetPoint(idList->GetId(0), a);
+	dataSet->GetPoint(idList->GetId(1), b);
+	dataSet->GetPoint(idList->GetId(2), c);
+	dataSet->GetPoint(idList->GetId(3), d);
+	// std::cout << "ids " << idList->GetId(0)
+	// 		  << " " << idList->GetId(1)
+	// 		  << " " << idList->GetId(2)
+	// 		  << " " << idList->GetId(3) << std::endl;
+	// std::cout << "coords" << std::endl
+	// 		  << a[0] << ", " << a[1] << ", " << a[2] << std::endl
+	// 		  << b[0] << ", " << b[1] << ", " << b[2] << std::endl
+	// 		  << c[0] << ", " << c[1] << ", " << c[2] << std::endl
+	// 		  << d[0] << ", " << d[1] << ", " << d[2] << std::endl;
+	double ab[3], ac[3], ad[3],
+		ba[3], bc[3], bd[3],
+		ca[3], cb[3], cd[3],
+		da[3], db[3], dc[3];
+	vtkMath::Subtract(b, a, ab);
+	vtkMath::Subtract(c, a, ac);
+	vtkMath::Subtract(d, a, ad);
+	vtkMath::Subtract(a, b, ba);
+	vtkMath::Subtract(c, b, bc);
+	vtkMath::Subtract(d, b, bd);
+	vtkMath::Subtract(a, c, ca);
+	vtkMath::Subtract(b, c, cb);
+	vtkMath::Subtract(d, c, cd);
+	vtkMath::Subtract(a, d, da);
+	vtkMath::Subtract(b, d, db);
+	vtkMath::Subtract(c, d, dc);
+	angles[0] = vtkMath::AngleBetweenVectors(ab, ac);
+	angles[1] = vtkMath::AngleBetweenVectors(ac, ad);
+	angles[2] = vtkMath::AngleBetweenVectors(ad, ab);
+	angles[3] = vtkMath::AngleBetweenVectors(ba, bc);
+	angles[4] = vtkMath::AngleBetweenVectors(bc, bd);
+	angles[5] = vtkMath::AngleBetweenVectors(bd, ba);
+	angles[6] = vtkMath::AngleBetweenVectors(ca, cb);
+	angles[7] = vtkMath::AngleBetweenVectors(cb, cd);
+	angles[8] = vtkMath::AngleBetweenVectors(cd, ca);
+	angles[9] = vtkMath::AngleBetweenVectors(da, db);
+	angles[10] = vtkMath::AngleBetweenVectors(db, dc);
+	angles[11] = vtkMath::AngleBetweenVectors(dc, da);
+}
+
 
 int main(int argc, char **argv) {
 	if (argc != 2) {
@@ -33,6 +108,25 @@ int main(int argc, char **argv) {
 	// vtkNew<vtkXMLPolyDataReader> reader;
 	vtkNew<vtkUnstructuredGridReader> reader;
 	reader->SetFileName(argv[1]);
+	reader->Update();
+	vtkUnstructuredGrid *grid = reader->GetOutput();
+	auto it = grid->NewCellIterator();
+	std::vector<double> angles(grid->GetNumberOfCells() * 12);
+	size_t i = 0;
+	for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
+		if (it->GetCellType() != VTK_TETRA) continue;
+		// double angles[12];
+		cellAngles(grid, it->GetPointIds(), angles.data() + i);
+		i += 12;
+		// for (size_t i = 0; i < 12; i++) {
+		// 	std::cout << angles[i] << ", ";
+		// }
+		// std::cout << "\b\b \n";
+	}
+	std::cout << "avg: " << average(angles)
+			  << ", stddev: " << standard_deviation(angles)
+			  << ", min: " << *std::min_element(angles.begin(), angles.end())
+			  << ", max: " << *std::max_element(angles.begin(), angles.end()) << std::endl;
 
 	vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
 	volumeMapper->SetInputConnection(reader->GetOutputPort());
@@ -42,9 +136,8 @@ int main(int argc, char **argv) {
 	vtkNew<vtkVolume> volume;
 	volume->SetMapper(volumeMapper);
 	vtkNew<vtkPiecewiseFunction> transferFunction;
-	transferFunction->AddPoint(-1, 100);
-	transferFunction->AddPoint(0, .5);
-	// transferFunction->AddPoint(1, 1);
+	transferFunction->AddPoint(-1, 0);
+	transferFunction->AddPoint(1, 1);
 	volume->GetProperty()->SetScalarOpacity(transferFunction);
 	volume->GetProperty()->SetColor(transferFunction);
 	// vtkNew<vtkActor> actor;