pfe/src/main.cc

#include "analysis/dihedral_angles_filter.h"
#include "analysis/max_distance_filter.h"
#include "fitting/remove_external_cells_filter.h"
#include "fitting/project_surface_points_on_poly.h"
#include "fitting/relaxation_filter.h"
#include "surface_points_filter.h"

#include <vtkGeometryFilter.h>
#include <vtkCellArrayIterator.h>
#include <vtkCellData.h>
#include <vtkDataSetReader.h>
#include <vtkOBJReader.h>
#include <vtkPolyData.h>
#include <vtkPolyDataReader.h>
#include <vtkSetGet.h>
#include <vtkUnstructuredGrid.h>
#include <vtkUnstructuredGridReader.h>
#include <vtkXMLPolyDataReader.h>
#include <vtkXMLUnstructuredGridReader.h>
#include <vtkXMLUnstructuredGridWriter.h>
#include <vtksys/SystemTools.hxx>
#include <string>

#ifdef USE_VIEWER
#include <vtkCamera.h>
#include <vtkNamedColors.h>
#include <vtkOpenGLProjectedTetrahedraMapper.h>
#include <vtkPiecewiseFunction.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkVolume.h>
#include <vtkVolumeMapper.h>
#include <vtkVolumeProperty.h>
#endif // USE_VIEWER

#include <array>


vtkAlgorithm *readerFromFileName(const char *fileName) {
	std::string extension =
		vtksys::SystemTools::GetFilenameLastExtension(fileName);
	if (extension == ".vtk") {
		auto reader = vtkDataSetReader::New();
		reader->SetFileName(fileName);
		return {reader};
	} else if (extension == ".vtu") {
		auto reader = vtkXMLUnstructuredGridReader::New();
		reader->SetFileName(fileName);
		return {reader};
	} else if (extension == ".vtp") {
		auto reader = vtkXMLPolyDataReader::New();
		reader->SetFileName(fileName);
		return {reader};
	} else if (extension == ".obj") {
		auto reader = vtkOBJReader::New();
		reader->SetFileName(fileName);
		return {reader};
	}
	throw std::runtime_error("Invalid file extension.");
}


static void usage(char **argv) {
	std::cerr << "Usage: " << argv[0] << " analyze|fit tetmesh polydata [radiusScale relaxationIterCount]" << std::endl;
}


int main(int argc, char **argv) {
	if (argc < 4) {
		usage(argv);
		return EXIT_FAILURE;
	}
	bool fit = false;
	if (std::string(argv[1]) == "fit") {
		fit = true;
	} else if (std::string(argv[1]) != "analyze") {
		usage(argv);
		return EXIT_FAILURE;
	}
	if ((fit && argc > 6) || (!fit && argc != 4)) {
		usage(argv);
		return EXIT_FAILURE;
	}

	double radiusScale = 2.;
	int iterCount = 1;
	if(argc > 4) {
		radiusScale = std::stod(argv[4]);
	}
	if (argc > 5) {
		iterCount = std::stoi(argv[5]);
	}

	auto tetMeshReader = readerFromFileName(argv[2]);
	auto polyMeshReader = readerFromFileName(argv[3]);

	vtkNew<RemoveExternalCellsFilter> removeExternalCellsFilter;
	removeExternalCellsFilter->SetInputConnection(0,
		tetMeshReader->GetOutputPort());
	removeExternalCellsFilter->SetInputConnection(1,
		polyMeshReader->GetOutputPort());

	vtkNew<SurfacePointsFilter> surfacePointsFilter;
	surfacePointsFilter->SetInputConnection(
		removeExternalCellsFilter->GetOutputPort());

	vtkNew<ProjectSurfacePointsOnPoly> project;
	project->SetRadiusScale(radiusScale);
	project->SetInputConnection(0, surfacePointsFilter->GetOutputPort());
	project->SetInputConnection(1, polyMeshReader->GetOutputPort());

	vtkNew<RelaxationFilter> relaxationFilter;
	relaxationFilter->SetIterCount(iterCount);
	relaxationFilter->SetInputConnection(project->GetOutputPort());

	vtkNew<DihedralAnglesFilter> dihedralAnglesFilter;
	if (fit) {
		dihedralAnglesFilter->SetInputConnection(
			relaxationFilter->GetOutputPort());
	} else {
		dihedralAnglesFilter->SetInputConnection(
			tetMeshReader->GetOutputPort());
	}

	vtkNew<vtkGeometryFilter> geometryFilter;
	geometryFilter->SetInputConnection(dihedralAnglesFilter->GetOutputPort());

	vtkNew<MaxDistanceFilter> maxDistanceFilter;
	maxDistanceFilter->SetInputConnection(0, geometryFilter->GetOutputPort());
	maxDistanceFilter->SetInputConnection(1, polyMeshReader->GetOutputPort());

	vtkNew<vtkXMLUnstructuredGridWriter> writer;
	writer->SetInputConnection(dihedralAnglesFilter->GetOutputPort());
	writer->SetDataModeToAscii();
	writer->SetFileName("out.vtu");
	writer->Write();

	maxDistanceFilter->Update();
	std::cerr << "Max distance: " << maxDistanceFilter->GetMaxDist() << "\n"
	          << "Max dist point id: " << maxDistanceFilter->GetMaxId()
	          << " in " << (maxDistanceFilter->GetMaxInput() == 0 ? "tetMesh\n" : "polyMesh\n")
	          << "Average min angle: " << dihedralAnglesFilter->GetAverageMinDegrees() << "\n"
	          << "Min min angle: " << dihedralAnglesFilter->GetMinMinDegrees() << "\n";

#ifdef USE_VIEWER
	/* Volume rendering properties */
	vtkNew<vtkOpenGLProjectedTetrahedraMapper> volumeMapper;
	volumeMapper->SetInputConnection(externalPointsFilter->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);

	/* 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();
#endif
	tetMeshReader->Delete();
	polyMeshReader->Delete();
	return EXIT_SUCCESS;
}