Paraview support for GenerateSDF.exe in InteractiveComputerGraphics/Discregrid

GenerateSDF.cpp

//
//https://github.com/InteractiveComputerGraphics/Discregrid


#include <Discregrid/All>
#include <Eigen/Dense>

#include "resource_path.hpp"

#include <string>
#include <iostream>
#include <array>

using namespace Eigen;

std::istream& operator>>(std::istream& is, std::array<unsigned int, 3>& data)  
{  
	is >> data[0] >> data[1] >> data[2];  
	return is;  
}  

std::istream& operator>>(std::istream& is, AlignedBox3d& data)  
{  
	is	>> data.min()[0] >> data.min()[1] >> data.min()[2]
		>> data.max()[0] >> data.max()[1] >> data.max()[2];  
	return is;  
}  

void createBoxMesh(int nx, int ny, int nz, Eigen::AlignedBox3d boundingbox,
	std::vector<Vector3d>& vertex, std::vector<Vector4i>& elements);
void writeSDF2VTK(std::string fname, std::vector<Vector3d>& vertex, std::vector<Vector4i>& elements, Eigen::VectorXd& sdf);

void createPlaneMesh(int nx, int ny, Eigen::AlignedBox3d boundingbox, std::string plane, double depth,
	std::vector<Vector3d>& vertex, std::vector<Vector4i>& elements);
void wrtePlaneSDF2VTK(std::string fname, const std::vector<Vector3d>& vertex, const std::vector<Vector4i>& elements, Eigen::VectorXd& sdf);


void SDFVolume_mesh(Discregrid::MeshDistance &md, std::array<unsigned int, 3> const& resolution,
				    Eigen::AlignedBox3d &boundingbox) {
	std::vector<Vector3d> vertex;
	std::vector<Vector4i> elements;

	createBoxMesh(resolution[0], resolution[1], resolution[2], 
		boundingbox, vertex, elements);

	Eigen::VectorXd sdf(vertex.size());

#pragma omp parallel for
		for (int l = 0; l < vertex.size(); ++l)
		{
			auto x = vertex[l];
			sdf[l] = md.signedDistanceCached(x);
		}

	writeSDF2VTK("sdfVolume_exact.vtk",vertex,elements,sdf);
}

void SDFVolume_cubic(Discregrid::CubicLagrangeDiscreteGrid& grid, std::array<unsigned int, 3> const& resolution,
	Eigen::AlignedBox3d& boundingbox) {
	std::vector<Vector3d> vertex;
	std::vector<Vector4i> elements;

	createBoxMesh(resolution[0], resolution[1], resolution[2],
		boundingbox, vertex, elements);

	Eigen::VectorXd sdf(vertex.size());

#pragma omp parallel for
		for (int l = 0; l < vertex.size(); ++l)
		{
			auto x = vertex[l];
			sdf[l] = grid.interpolate(0,x);
			if (sdf[l] == std::numeric_limits<double>::max())
			{
				sdf[l] = 0.0;
			}
		}

	writeSDF2VTK("sdfVolume_grid.vtk",vertex, elements, sdf);
}

void SDFPlane_mesh(Discregrid::MeshDistance& md, std::array<unsigned int, 3> const& resolution,
	Eigen::AlignedBox3d& boundingbox,std::string plane="xy", double depth=0.5) {
	std::vector<Vector3d> vertex;
	std::vector<Vector4i> elements;

	createPlaneMesh(resolution[0], resolution[1], boundingbox, plane, depth,
		vertex, elements);

	Eigen::VectorXd sdf(vertex.size());

#pragma omp parallel for
	for (int l = 0; l < vertex.size(); ++l)
	{
		auto x = vertex[l];
		sdf[l] = md.signedDistanceCached(x);
	}

	std::string fname = "SDFPlane_exact_" + std::string(plane) + std::to_string(depth) + ".vtk";
	wrtePlaneSDF2VTK(fname, vertex, elements, sdf);
}

#include <cxxopts/cxxopts.hpp>

int main(int argc, char* argv[])
{
	cxxopts::Options options(argv[0], "Generates a signed distance field from a closed two-manifold triangle mesh.");
	options.positional_help("[input OBJ file]");

	options.add_options()
	("h,help", "Prints this help text")
	("r,resolution", "Grid resolution", cxxopts::value<std::array<unsigned int, 3>>()->default_value("10 10 10"))
	("d,domain", "Domain extents (bounding box), format: \"minX minY minZ maxX maxY maxZ\"", cxxopts::value<AlignedBox3d>())
	("i,invert", "Invert SDF")
	("o,output", "Ouput file in cdf format", cxxopts::value<std::string>()->default_value(""))
	("input", "OBJ file containing input triangle mesh", cxxopts::value<std::vector<std::string>>())
	;

	try
	{
		options.parse_positional("input");
		auto result = options.parse(argc, argv);

		if (result.count("help"))
		{
			std::cout << options.help() << std::endl;
			std::cout << std::endl << std::endl << "Example: GenerateSDF -r \"50 50 50\" dragon.obj" << std::endl;
			exit(0);
		}
		if (!result.count("input"))
		{
			std::cout << "ERROR: No input mesh given." << std::endl;
			std::cout << options.help() << std::endl;
			std::cout << std::endl << std::endl << "Example: GenerateSDF -r \"50 50 50\" dragon.obj" << std::endl;
			exit(1);
		}
		auto resolution = result["r"].as<std::array<unsigned int, 3>>();
		auto filename = result["input"].as<std::vector<std::string>>().front();

		if (!std::ifstream(filename).good())
		{
			std::cerr << "ERROR: Input file does not exist!" << std::endl;
			exit(1);
		}

		std::cout << "Load mesh...";
		Discregrid::TriangleMesh mesh(filename);
		std::cout << "DONE" << std::endl;

		std::cout << "Set up data structures...";
		Discregrid::MeshDistance md(mesh);
		std::cout << "DONE" << std::endl;

		Eigen::AlignedBox3d domain;
		domain.setEmpty();
		if (result.count("d"))
		{
			domain = result["d"].as<Eigen::AlignedBox3d>();
		}
		if (domain.isEmpty())
		{
			for (auto const& x : mesh.vertices())
			{
				domain.extend(x);
			}
			domain.max() += 1.0e-3 * domain.diagonal().norm() * Vector3d::Ones();
			domain.min() -= 1.0e-3 * domain.diagonal().norm() * Vector3d::Ones();
		}

		Discregrid::CubicLagrangeDiscreteGrid sdf(domain, resolution);
		auto func = Discregrid::DiscreteGrid::ContinuousFunction{};
		if (result.count("invert"))
		{
			func = [&md](Vector3d const& xi) {return -1.0 * md.signedDistanceCached(xi); };
		}
		else
		{
			func = [&md](Vector3d const& xi) {return md.signedDistanceCached(xi); };
		}

		std::cout << "Save signed distance function into vtk..." << std::endl;
		//SDFVolume_mesh(md, resolution, domain);
		SDFPlane_mesh(md, resolution, domain, "xz", 0.5);
		SDFPlane_mesh(md, resolution, domain, "xy", 0.7);
		std::cout << "DONE" << std::endl;

		std::cout << "Generate discretization..." << std::endl;
		sdf.addFunction(func, true);
		std::cout << "DONE" << std::endl;

		std::cout << "Serialize discretization...";
		auto output_file = result["o"].as<std::string>();
		if (output_file == "")
		{
			output_file = filename;
			if (output_file.find(".") != std::string::npos)
			{
				auto lastindex = output_file.find_last_of(".");
				output_file = output_file.substr(0, lastindex);
			}
			output_file += ".cdf";
		}
		sdf.save(output_file);
		std::cout << "DONE" << std::endl;


		std::cout << "Save signed distance function into vtk..." << std::endl;
		SDFVolume_cubic(sdf, resolution, domain);
		std::cout << "DONE" << std::endl;
	}
	catch (cxxopts::OptionException const& e)
	{
		std::cout << "error parsing options: " << e.what() << std::endl;
		exit(1);
	}
	
	return 0;
}


//------------------------------------------------------------------------
//------------------------------Helper Funcs------------------------------
//------------------------------------------------------------------------
void createBoxMesh(int nx, int ny, int nz, Eigen::AlignedBox3d boundingbox,
	std::vector<Vector3d>& vertex, std::vector<Vector4i>& elements) {

	// Extract minimum and maximum coordinates
	const double x0 = boundingbox.min()[0];
	const double x1 = boundingbox.max()[0];
	const double y0 = boundingbox.min()[1];
	const double y1 = boundingbox.max()[1];
	const double z0 = boundingbox.min()[2];
	const double z1 = boundingbox.max()[2];

	const double a = x0;
	const double b = x1;
	const double c = y0;
	const double d = y1;
	const double e = z0;
	const double f = z1;

	assert(nx >= 1 && ny >= 1 && nz >= 1, "BoxMesh: number of vertices must be at least 1 in each dimension");

	// Create vertices
	Vector3d x;
	std::size_t nverts = 0;
	for (std::size_t iz = 0; iz <= nz; iz++)
	{
		x[2] = e + (static_cast<double>(iz)) * (f - e) / static_cast<double>(nz);
		for (std::size_t iy = 0; iy <= ny; iy++)
		{
			x[1] = c + (static_cast<double>(iy)) * (d - c) / static_cast<double>(ny);
			for (std::size_t ix = 0; ix <= nx; ix++)
			{
				x[0] = a + (static_cast<double>(ix)) * (b - a) / static_cast<double>(nx);

				Vector3d pos = x;
				vertex.push_back(pos);
				nverts++;
			}
		}
	}

	// Create tetrahedra
	std::size_t cell = 0;
	for (std::size_t iz = 0; iz < nz; iz++)
	{
		for (std::size_t iy = 0; iy < ny; iy++)
		{
			for (std::size_t ix = 0; ix < nx; ix++)
			{
				const std::size_t v0 = iz * (nx + 1) * (ny + 1) + iy * (nx + 1) + ix;
				const std::size_t v1 = v0 + 1;
				const std::size_t v2 = v0 + (nx + 1);
				const std::size_t v3 = v1 + (nx + 1);
				const std::size_t v4 = v0 + (nx + 1) * (ny + 1);
				const std::size_t v5 = v1 + (nx + 1) * (ny + 1);
				const std::size_t v6 = v2 + (nx + 1) * (ny + 1);
				const std::size_t v7 = v3 + (nx + 1) * (ny + 1);

				// Add cells
				// Note that v0 < v1 < v2 < v3 < vmid.
				elements.push_back(Vector4i(v0, v1, v3, v7));
				elements.push_back(Vector4i(v0, v1, v7, v5));
				elements.push_back(Vector4i(v0, v5, v7, v4));
				elements.push_back(Vector4i(v0, v3, v2, v7));
				elements.push_back(Vector4i(v0, v6, v4, v7));
				elements.push_back(Vector4i(v0, v2, v6, v7));
			}
		}
	}


	int NumOfVertex = vertex.size();
	int NumOfElements = elements.size();
	printf("#Box Mesh has %d vertex\n", NumOfVertex);
	printf("                    %d elements\n", NumOfElements);
	std::cout
		<< "                    bounding box " << boundingbox.min().transpose()
		<< ", " << boundingbox.max().transpose() << std::endl;
}


void writeSDF2VTK(std::string fname, std::vector<Vector3d>& vertex, std::vector<Vector4i>& elements, Eigen::VectorXd& sdf)
{
	int i;
	FILE* fp;
	char fileName[1024];

	sprintf(fileName, fname.c_str());
	printf("#GCPS: Write mesh to file %s...\n", fileName);

	fp = fopen(fileName, "w");
	fprintf(fp, "# vtk DataFile Version 4.1\n");
	fprintf(fp, "vtk output\n");
	fprintf(fp, "ASCII\n");
	fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");

	//Output vertices
	int NumVerts = vertex.size();
	fprintf(fp, "POINTS %d float\n", NumVerts);
	for (auto vi = 0; vi < NumVerts; ++vi) {
		fprintf(fp, "%lf %lf %lf\n", vertex[vi][0],
			vertex[vi][1],
			vertex[vi][2]);
	}
	fprintf(fp, "\n");

	//Output elements
	int NumEles = elements.size();
	fprintf(fp, "CELLS %d %d\n", NumEles, NumEles * 5);//5 value define a tet ele
	for (auto ei = 0; ei < NumEles; ++ei) {
		Vector4i ele = elements[ei];
		fprintf(fp, "4 %d %d %d %d\n", ele[0], ele[1], ele[2], ele[3]);
	}
	fprintf(fp, "\n");

	fprintf(fp, "CELL_TYPES %d\n", NumEles);
	for (auto ei = 0; ei < NumEles; ++ei) {
		if (ei % 10 == 0 && ei > 0) fprintf(fp, "\n");
		fprintf(fp, "%d ", 10);
	}
	fprintf(fp, "\n");
	fprintf(fp, "\n");

	//Output scalars/vector fields
	fprintf(fp, "POINT_DATA %d\n", NumVerts);
	fprintf(fp, "FIELD FieldData %d\n", 1);

	fprintf(fp, "sdf 1 %d float\n", NumVerts);
	for (auto vi = 0; vi < NumVerts; ++vi) {
		fprintf(fp, "%lf\n", sdf[vi]);
	}
	fprintf(fp, "\n");

	fclose(fp);
}


void createPlaneMesh(int nx, int ny, Eigen::AlignedBox3d boundingbox, std::string plane, double depth,
	std::vector<Vector3d>& vertex, std::vector<Vector4i>& elements)
{
	Vector3i dir;
	if (plane == "xy")
		dir = Vector3i(0, 1, 2);
	else if (plane == "xz")
		dir = Vector3i(0, 2, 1);
	else//yz
		dir = Vector3i(1, 2, 0);
	int _x = dir[0], _y = dir[1], _z = dir[2];

	// Extract minimum and maximum coordinates
	Vector3d boxsize = boundingbox.sizes();
	const double x0 = boundingbox.min()[_x];
	const double x1 = boundingbox.max()[_x];
	const double y0 = boundingbox.min()[_y];
	const double y1 = boundingbox.max()[_y];
	const double zp = boundingbox.min()[_z] + depth * boxsize[_z];

	const double a = x0;
	const double b = x1;
	const double c = y0;
	const double d = y1;

	// Create vertices:
	Vector3d x;
	std::size_t nverts = 0;

	x[2] = zp;
	for (std::size_t iy = 0; iy <= ny; iy++)
	{
		x[1] = c + ((static_cast<double>(iy)) * (d - c) / static_cast<double>(ny));
		for (std::size_t ix = 0; ix <= nx; ix++)
		{
			x[0] = a + ((static_cast<double>(ix)) * (b - a) / static_cast<double>(nx));

			Vector3d pos;
			pos[_x] = x[0]; pos[_y] = x[1]; pos[_z] = x[2];

			vertex.push_back(pos);
			nverts++;
		}
	}

	// Create triangles
	for (std::size_t iy = 0; iy < ny; iy++)
	{
		for (std::size_t ix = 0; ix < nx; ix++)
		{
			const std::size_t v0 = iy * (nx + 1) + ix;
			const std::size_t v1 = v0 + 1;
			const std::size_t v2 = v0 + (nx + 1);
			const std::size_t v3 = v1 + (nx + 1);
			std::vector<std::size_t> cell_data;

			//two triangles per lattice site
			elements.push_back(Vector4i(v0, v1, v2, -1));
			elements.push_back(Vector4i(v1, v2, v3, -1));
		}
	}


	int NumOfVertex = vertex.size();
	int NumOfElements = elements.size();
	printf("#GCPS: Plane Mesh has %d vertex\n", NumOfVertex);
	printf("                    %d elements\n", NumOfElements);
	printf("			   dir%s (%d,%d,%d)\n", plane.c_str(), dir[0], dir[1], dir[2]);
}

void wrtePlaneSDF2VTK(std::string fname, const std::vector<Vector3d>& vertex, const std::vector<Vector4i>& elements, Eigen::VectorXd& sdf)
{
	int i;
	FILE* fp;
	char fileName[1024];

	sprintf(fileName, fname.c_str());
	printf("#GCPS: Write sdf plane to file %s...\n", fileName);

	fp = fopen(fileName, "w");
	fprintf(fp, "# vtk DataFile Version 4.1\n");
	fprintf(fp, "vtk output\n");
	fprintf(fp, "ASCII\n");
	fprintf(fp, "DATASET POLYDATA\n");

	//Output vertices
	int NumVerts = vertex.size();
	fprintf(fp, "POINTS %d float\n", NumVerts);
	for (auto vi = 0; vi < NumVerts; ++vi) {
		fprintf(fp, "%lf %lf %lf\n", vertex[vi][0],
			vertex[vi][1],
			vertex[vi][2]);
	}
	fprintf(fp, "\n");

	//Output bdfacets
	int NumFacets = elements.size();
	fprintf(fp, "POLYGONS %d %d\n", NumFacets, NumFacets * 4);//4 value define a triangle facet
	for (auto fi = 0; fi < NumFacets; ++fi) {
		Vector4i face = elements[fi];
		fprintf(fp, "3 %d %d %d\n", face[0], face[1], face[2]);
	}
	fprintf(fp, "\n");


	//Output scalars/vector fields
	fprintf(fp, "POINT_DATA %d\n", NumVerts);
	fprintf(fp, "FIELD FieldData %d\n", 1);

	fprintf(fp, "sdf 1 %d float\n", NumVerts);
	for (auto vi = 0; vi < NumVerts; ++vi) {
		fprintf(fp, "%lf\n", sdf[vi]);
	}
	fprintf(fp, "\n");

	fclose(fp);
}