add aspect ratio computation

CMakeLists.txt

@@ -11,7 +11,9 @@ target_compile_features(pfe PRIVATE cxx_std_17)
 target_sources(pfe PRIVATE
   src/main.cc
   src/angles_filter.cc
-  src/angles_filter.h)
+  src/angles_filter.h
+  src/aspect_ratio_filter.cc
+  src/aspect_ratio_filter.h)
 
 target_link_libraries(pfe PRIVATE
   VTK::CommonCore

src/aspect_ratio_filter.cc

@@ -0,0 +1,69 @@
+#include "aspect_ratio_filter.h"
+
+#include <vtkUnstructuredGrid.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkDoubleArray.h>
+#include <vtkCellIterator.h>
+
+
+vtkStandardNewMacro(AspectRatioFilter);
+
+
+int AspectRatioFilter::FillInputPortInformation(int port, vtkInformation *info) {
+	vtkUnstructuredGridAlgorithm::FillInputPortInformation(port, info);
+	// Ensure there is an "angles" array in the DataCell.
+	return 1;
+}
+
+
+vtkTypeBool AspectRatioFilter::RequestData(
+	vtkInformation *request,
+	vtkInformationVector **inputVector,
+	vtkInformationVector *outputVector) {
+	(void) request;
+	vtkUnstructuredGrid* input =
+		vtkUnstructuredGrid::GetData(inputVector[0]);
+	vtkUnstructuredGrid* output =
+		vtkUnstructuredGrid::GetData(outputVector);
+
+	/* Pass the input structure and point data unmodified. */
+	output->CopyStructure(input);
+	output->GetPointData()->PassData(input->GetPointData());
+
+	/* Pass the existing cell data unmodified. */
+	vtkCellData *cellData = output->GetCellData();
+	cellData->PassData(input->GetCellData());
+
+	/* Then add an array to it. */
+	vtkNew<vtkDoubleArray> aspectRatios;
+	aspectRatios->SetName("aspect_ratio");
+	aspectRatios->SetNumberOfComponents(1);
+	aspectRatios->SetNumberOfTuples(input->GetNumberOfCells());
+
+	auto *angles = vtkDoubleArray::SafeDownCast(input->GetCellData()->GetArray("angles"));
+
+	double cellAngles[12];
+	double triangleRatios[4];
+	for (vtkIdType cellId = 0; cellId < input->GetNumberOfCells(); cellId++) {
+		if (input->GetCellType(cellId) != VTK_TETRA) continue;
+		angles->GetTuple(cellId, cellAngles);
+		for (vtkIdType angleId = 0;
+			 angleId < angles->GetNumberOfComponents();
+			 angleId += 3) {
+			double minAngle = std::min(
+				std::min(cellAngles[angleId], cellAngles[angleId + 1]),
+				cellAngles[angleId + 2]);
+			double maxAngle = std::max(
+				std::max(cellAngles[angleId], cellAngles[angleId + 1]),
+				cellAngles[angleId + 2]);
+			triangleRatios[angleId / 3] = minAngle / maxAngle;
+		}
+		double worst = std::min(
+			std::min(triangleRatios[0], triangleRatios[1]),
+			std::min(triangleRatios[2], triangleRatios[3]));
+		aspectRatios->SetTuple1(cellId, worst);
+	}
+	cellData->AddArray((vtkAbstractArray *) aspectRatios);
+	return true;
+}

src/aspect_ratio_filter.h

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

src/main.cc

@@ -1,4 +1,5 @@
 #include "angles_filter.h"
+#include "aspect_ratio_filter.h"
 
 #include <vtkCellData.h>
 #include <vtkUnstructuredGrid.h>
@@ -47,17 +48,16 @@ int main(int argc, char **argv) {
 		return EXIT_FAILURE;
 	}
 
-	vtkNew<vtkNamedColors> colors;
-
-	std::array<unsigned char, 4> bkg{{26, 51, 102, 255}};
-	colors->SetColor("BkgColor", bkg.data());
-
-	// vtkNew<vtkXMLPolyDataReader> reader;
+	/* Reader */
 	vtkNew<vtkUnstructuredGridReader> reader;
 	reader->SetFileName(argv[1]);
-	// reader->Update();
+
+
+	/* Angle filter */
 	vtkNew<AnglesFilter> anglesFilter;
 	anglesFilter->SetInputConnection(reader->GetOutputPort());
+
+	/* Display angles in console. */
 	anglesFilter->Update();
 	vtkUnstructuredGrid *grid = anglesFilter->GetOutput();
 	auto *angles = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("angles"));
@@ -68,16 +68,23 @@ int main(int argc, char **argv) {
 		}
 		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;
 
+	/* Aspect ratio */
+	vtkNew<AspectRatioFilter> aspectRatioFilter;
+	aspectRatioFilter->SetInputConnection(anglesFilter->GetOutputPort());
+
+	/* Display aspect ratios in console. */
+	aspectRatioFilter->Update();
+	grid = aspectRatioFilter->GetOutput();
+	auto *aspectRatios = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("aspect_ratio"));
+	for (ssize_t i = 0; i < aspectRatios->GetNumberOfTuples(); i++) {
+		std::cout << "cell " << i << ": "
+				  << aspectRatios->GetTuple1(i) << "\n";
+	}
+
+	/* Volume rendering properties */
 	vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
-	volumeMapper->SetInputConnection(anglesFilter->GetOutputPort());
-	// vtkNew<vtkPolyDataMapper> mapper;
-	// mapper->SetInputConnection(reader->GetOutputPort());
-
+	volumeMapper->SetInputConnection(aspectRatioFilter->GetOutputPort());
 	vtkNew<vtkVolume> volume;
 	volume->SetMapper(volumeMapper);
 	vtkNew<vtkPiecewiseFunction> transferFunction;
@@ -85,16 +92,12 @@ int main(int argc, char **argv) {
 	transferFunction->AddPoint(1, 1);
 	volume->GetProperty()->SetScalarOpacity(transferFunction);
 	volume->GetProperty()->SetColor(transferFunction);
-	// vtkNew<vtkActor> actor;
-	// actor->SetMapper(mapper);
-
-	// actor->GetProperty()->SetColor(
-	// 	colors->GetColor4d("Tomato").GetData());
-	// actor->RotateX(30.0);
-	// actor->RotateY(-45.0);
 
+	/* Renderer */
+	vtkNew<vtkNamedColors> colors;
+	std::array<unsigned char, 4> bkg{{26, 51, 102, 255}};
+	colors->SetColor("BkgColor", bkg.data());
 	vtkNew<vtkRenderer> renderer;
-	// renderer->AddActor(actor);
 	renderer->AddVolume(volume);
 	renderer->SetBackground(colors->GetColor3d("BkgColor").GetData());
 	renderer->ResetCamera();