mod_geo-tp/external/MeshReconstruction/lib/Cube.cpp

#include "Cube.h"
#include <cmath>

using namespace MeshReconstruction;

namespace
{
	// Cube has 8 vertices. Each vertex can have positive or negative sign.
	// 2^8 = 256 possible configurations mapped to intersected edges in each case.
	// The 12 edges are numbered as 1, 2, 4, ..., 2048 and are stored as a 12-bit bitstring for each configuration.
	const int signConfigToIntersectedEdges[256] = {
		0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
		0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
		0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
		0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
		0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
		0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
		0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
		0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
		0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
		0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
		0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
		0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
		0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
		0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
		0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
		0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
		0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
		0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
		0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
		0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
		0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
		0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
		0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
		0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
		0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
		0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
		0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
		0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
		0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
		0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
		0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
		0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 };

	struct Edge
	{
		int edgeFlag : 12; // flag: 1, 2, 4, ... 2048
		int vert0; // 0-7
		int vert1; // 0-7
	};

	const Edge edges[12] =
	{
		{ 1, 0, 1 }, // edge 0
		{ 2, 1, 2 }, // edge 1
		{ 4, 2, 3 }, // ...
		{ 8, 3, 0 },
		{ 16, 4, 5 },
		{ 32, 5, 6 },
		{ 64, 6, 7 },
		{ 128, 7, 4 },
		{ 256, 0, 4 },
		{ 512, 1, 5 },
		{ 1024, 2, 6 },
		{ 2048, 3, 7 } // edge 11
	};
}

Vec3 Cube::LerpVertex(double isoLevel, int i1, int i2) const
{
	auto const Eps = 1e-5;
	auto const v1 = sdf[i1];
	auto const v2 = sdf[i2];
	auto const& p1 = pos[i1];
	auto const& p2 = pos[i2];

	if (abs(isoLevel - v1) < Eps) return p1;
	if (abs(isoLevel - v2) < Eps) return p2;
	if (abs(v1 - v2) < Eps) return p1;

	auto mu = (isoLevel - v1) / (v2 - v1);
	return p1 + (p2 - p1)*mu;
}

Cube::Cube(Rect3 const& space, Fun3s const& sdf)
{
	auto mx = space.min.x;
	auto my = space.min.y;
	auto mz = space.min.z;

	auto sx = space.size.x;
	auto sy = space.size.y;
	auto sz = space.size.z;

	pos[0] = space.min;
	pos[1] = { mx + sx, my, mz };
	pos[2] = { mx + sx, my, mz + sz };
	pos[3] = { mx, my, mz + sz };
	pos[4] = { mx, my + sy, mz };
	pos[5] = { mx + sx, my + sy, mz };
	pos[6] = { mx + sx, my + sy, mz + sz };
	pos[7] = { mx, my + sy, mz + sz };

	for (auto i = 0; i < 8; ++i)
	{
		auto sd = sdf(pos[i]);
		if (sd == 0) sd += 1e-6;
		this->sdf[i] = sd;
	}
}

int Cube::SignConfig(double isoLevel) const
{
	auto edgeIndex = 0;

	for (auto i = 0; i < 8; ++i)
	{
		if (sdf[i] < isoLevel)
		{
			edgeIndex |= 1 << i;
		}
	}

	return edgeIndex;
}

IntersectInfo Cube::Intersect(double iso) const
{
	// idea:
	// from signs at 8 corners of cube a sign configuration (256 possible ones) is computed
	// this configuration can be used to index into a table that tells which of the 12 edges are intersected
	// find vertices adjacent to edges and interpolate cut vertex and store it in IntersectionInfo object

	IntersectInfo intersect;
	intersect.signConfig = SignConfig(iso);

	for (auto e = 0; e<12; ++e)
	{
		if (signConfigToIntersectedEdges[intersect.signConfig] & edges[e].edgeFlag)
		{
			auto v0 = edges[e].vert0;
			auto v1 = edges[e].vert1;
			auto vert = LerpVertex(iso, v0, v1);
			intersect.edgeVertIndices[e] = vert;
		}
	}

	return intersect;
}