pfe/src/main.cc

#include "vtkCell.h"
#include "vtkUnstructuredGrid.h"
#include <algorithm>
#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkCylinderSource.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkVolumeProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkVolumeMapper.h>
#include <vtkVolume.h>
#include <vtkOpenGLProjectedTetrahedraMapper.h>
#include <vtkUnstructuredGridReader.h>
#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) {
		std::cerr << "Usage: " << argv[0] << " FILE.vtk" << std::endl;
		return EXIT_FAILURE;
	}

	vtkNew<vtkNamedColors> colors;

	std::array<unsigned char, 4> bkg{{26, 51, 102, 255}};
	colors->SetColor("BkgColor", bkg.data());

	// 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());
	// vtkNew<vtkPolyDataMapper> mapper;
	// mapper->SetInputConnection(reader->GetOutputPort());

	vtkNew<vtkVolume> volume;
	volume->SetMapper(volumeMapper);
	vtkNew<vtkPiecewiseFunction> transferFunction;
	transferFunction->AddPoint(-1, 0);
	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);

	vtkNew<vtkRenderer> renderer;
	// renderer->AddActor(actor);
	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;
}