pfe/src/main.cc

#include <vtkCell.h>
#include <vtkPointData.h>
#include <vtkCellData.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 <vtkUnstructuredGridAlgorithm.h>
#include <vtkDoubleArray.h>

#include <array>
#include <vector>
#include <algorithm>


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);
}


class AddAngleComputationAlgorithm : public vtkUnstructuredGridAlgorithm {
public:
	static AddAngleComputationAlgorithm *New();
	vtkTypeMacro(AddAngleComputationAlgorithm, vtkUnstructuredGridAlgorithm);

	vtkTypeBool RequestData(vtkInformation *request,
							vtkInformationVector **inputVector,
							vtkInformationVector *outputVector) override {
		(void) request;
		vtkUnstructuredGrid* input =
			vtkUnstructuredGrid::GetData(inputVector[0]);
		vtkUnstructuredGrid* output =
			vtkUnstructuredGrid::GetData(outputVector);
		output->CopyStructure(input);
		output->GetPointData()->PassData(input->GetPointData());
		vtkCellData *cellData = output->GetCellData();
		cellData->PassData(input->GetCellData());
		vtkNew<vtkDoubleArray> anglesArray;
		anglesArray->SetName("angles");
		anglesArray->SetNumberOfComponents(12);
		anglesArray->SetNumberOfTuples(input->GetNumberOfCells());
		double *anglesBase = anglesArray->GetPointer(0);
		size_t i = 0;
		auto it = input->NewCellIterator();
		for (it->InitTraversal(); !it->IsDoneWithTraversal(); it->GoToNextCell()) {
			if (it->GetCellType() != VTK_TETRA) continue;
			cellAngles(input, it->GetPointIds(), anglesBase + i);
			i += 12;
		}
		cellData->AddArray((vtkAbstractArray *) anglesArray);
		return true;
	}

protected:
	~AddAngleComputationAlgorithm() override = default;
};

vtkStandardNewMacro(AddAngleComputationAlgorithm);

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) {
	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();
	vtkNew<AddAngleComputationAlgorithm> addAngles;
	addAngles->SetInputConnection(reader->GetOutputPort());
	addAngles->Update();
	vtkUnstructuredGrid *grid = addAngles->GetOutput();
	auto *angles = vtkDoubleArray::SafeDownCast(grid->GetCellData()->GetArray("angles"));
	for (ssize_t i = 0; i < angles->GetNumberOfTuples(); i++) {
		std::cout << "cell " << i << ": ";
		for (ssize_t j = 0; j < angles->GetNumberOfComponents(); j++) {
			std::cout << angles->GetTuple(i)[j] << ", ";
		}
		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(addAngles->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;
}