trimesh.cpp

来自「ＲＢＦ平台」· C++ 代码 · 共 1,522 行 · 第 1/3 页

#include "stdafx.h"
#include "TriMesh.h"
#include "meshio.h"
#include "bsphere.h"
#include "timestamp.h"
#include "lineqn.h"
#include <algorithm>
#include <utility>
#include <queue>
using std::pair;
using std::make_pair;
using std::priority_queue;
using std::find;


// Forward declarations
static void tstrip_build(TriMesh &mesh, int f, vector<signed char> &face_avail,
			             vector<int> &todo);

static void collect_tris_in_strips(vector<int> &tstrips);

static void rot_coord_sys(const vec &old_u, const vec &old_v, const vec &new_norm,
				 		  vec &new_u, vec &new_v);

static float cosmaxangle(const point &v1, const point &v2, const point &v3);
static float flip_benefit(const point &v1, const point &v2, const point &v3, const point &v4);
static float flip_benefit(const TriMesh *mesh, int f, int e);
int edge_flip(TriMesh *mesh, int f1, int e);

TriMesh::TriMesh()
{

}

TriMesh::~TriMesh()
{
	if(!vertices.empty()) vertices.clear();
	if(!faces.empty()) faces.clear();
	if(!tstrips.empty()) tstrips.clear();
	if(!colors.empty()) colors.clear();
	if(!confidences.empty()) confidences.clear();
	if(!flags.empty()) flags.clear();
	if(!normals.empty()) normals.clear();
	if(!pdir1.empty()) pdir1.clear();
	if(!pdir2.empty()) pdir2.clear();
	if(!curv1.empty()) curv1.clear();
	if(!curv2.empty()) curv2.clear();
	if(!dcurv.empty()) dcurv.clear();
	if(!cornerareas.empty()) cornerareas.clear();
	if(!pointareas.empty()) pointareas.clear();
	if(!neighbors.empty()) neighbors.clear();
	if(!adjacentfaces.empty()) adjacentfaces.clear();
	if(!across_edge.empty()) across_edge.clear();
}

// Read a TriMesh from a file.  Defined to use a helper function to make
// subclassing easier.
TriMesh *TriMesh::read(const char *filename)
{
	TriMesh *mesh = new TriMesh();
	CMeshIO meshio;

	if (meshio.read_helper(filename, mesh))
		return mesh;

	delete mesh;
	return NULL;
}

void TriMesh::write(TriMesh * trimesh, const char *filename, bool istrips)
{
	CMeshIO meshio;
	meshio.write(filename,trimesh);
}

// Debugging printout, controllable by a "verbose"ness parameter
int TriMesh::verbose = 1;

void TriMesh::set_verbose(int verbose_)
{
	verbose = verbose_;
}

int TriMesh::dprintf(const char *format, ...)
{
	if (!verbose)
		return 0;

	va_list ap;
	va_start(ap, format);
	int ret = vfprintf(stderr, format, ap);
	va_end(ap);
	fflush(stderr);
	return ret;
}

// Convert faces to tstrips
void TriMesh::need_tstrips()
{
	if (faces.empty() || !tstrips.empty())
		return;

	need_across_edge();

	dprintf("Building triangle strips... ");
	int nf = faces.size();

	vector<int> todo;
	vector<signed char> face_avail;
	face_avail.reserve(nf);
	for (int i = 0; i < nf; i++) {
		face_avail[i] = (across_edge[i][0] != -1) +
				(across_edge[i][1] != -1) +
				(across_edge[i][2] != -1);
		if (face_avail[i] == 1)
			todo.push_back(i);
	}

	tstrips.reserve(faces.size() * 2);

	int nstrips = 0;
	i = 0;
	while (i < faces.size()) {
		int next;
		if (todo.empty()) {
			next = i++;
		} else {
			next = todo.back();
			todo.pop_back();
		}
		if (face_avail[next] < 0)
			continue;
		tstrip_build(*this, next, face_avail, todo);
		nstrips++;
	}

	convert_strips(TSTRIP_LENGTH);

	dprintf("Done.\n  %d strips (Avg. length %.1f)\n",
		nstrips, (float) faces.size() / nstrips);
}


// Build a triangle strip starting with the given face
static void tstrip_build(TriMesh &mesh, int f, vector<signed char> &face_avail,
	vector<int> &todo)
{
	TriMesh::Face &v = mesh.faces[f];
	if (face_avail[f] == 0) {
		mesh.tstrips.push_back(v[0]);
		mesh.tstrips.push_back(v[1]);
		mesh.tstrips.push_back(v[2]);
		mesh.tstrips.push_back(-1);
		face_avail[f] = -1;
		return;
	}

	int score[3];
	for (int i = 0; i < 3; i++) {
		score[i] = 0;
		int ae = mesh.across_edge[f][i];
		if (ae == -1 || face_avail[ae] < 0)
			continue;
		score[i]++;
		int next_edge = mesh.faces[ae].indexof(v[(i+1)%3]);
		int nae = mesh.across_edge[ae][next_edge];
		if (nae == -1 || face_avail[nae] < 0)
			continue;
		score[i]++;
		if (face_avail[ae] == 2)
			score[i]++;
	}

	int best_score = __max(__max(score[0], score[1]), score[2]);
	int best = (score[0] == best_score) ? 0 :
		   (score[1] == best_score) ? 1 : 2;

	int vlast2 = v[ best     ];
	int vlast1 = v[(best+1)%3];
	int vnext  = v[(best+2)%3];
	int dir = 1;
	mesh.tstrips.push_back(vlast2);
	mesh.tstrips.push_back(vlast1);

	while (1) {
		mesh.tstrips.push_back(vnext);
		face_avail[f] = -1;
		for (int j = 0; j < 3; j++) {
			int ae = mesh.across_edge[f][j];
			if (ae == -1)
				continue;
			if (face_avail[ae] > 0)
				face_avail[ae]--;
			if (face_avail[ae] == 1)
				todo.push_back(ae);
		}

		f = mesh.across_edge[f][mesh.faces[f].indexof(vlast2)];
		if (f == -1 || face_avail[f] < 0)
			break;
		vlast2 = vlast1;
		vlast1 = vnext;
		vnext = mesh.faces[f][(mesh.faces[f].indexof(vlast2)+3+dir)%3];
		dir = -dir;
	}
	
	mesh.tstrips.push_back(-1);
}

// Convert triangle strips to faces
void TriMesh::need_faces()
{
	if (tstrips.empty() || !faces.empty())
		return;

	//dprintf("Unpacking triangle strips... ");
	int nfaces = 0;
	int i = 0;
	while (i < tstrips.size()) {
		nfaces += tstrips[i] - 2;
		i += tstrips[i] + 1;
	}

	faces.reserve(nfaces);

	int len = 0;
	bool flip = false;
	for (i = 0; i < tstrips.size(); i++) {
		if (len == 0) {
			len = tstrips[i] - 2;
			flip = false;
			i += 2;
			continue;
		}
		if (flip)
			faces.push_back(Face(tstrips[i-1],
					     tstrips[i-2],
					     tstrips[i]));
		else
			faces.push_back(Face(tstrips[i-2],
					     tstrips[i-1],
					     tstrips[i]));
		flip = !flip;
		len--;
	}
	//dprintf("Done.\n  %d triangles\n", nfaces);
}

// Convert between "length preceding strip" and "-1 following strip"
// representations
void TriMesh::convert_strips(tstrip_rep rep)
{
	if (tstrips.empty())
		return;
	if (rep == TSTRIP_TERM && tstrips.back() == -1)
		return;
	if (rep == TSTRIP_LENGTH && tstrips.back() != -1) {
		//collect_tris_in_strips(tstrips);
		return;
	}

	if (rep == TSTRIP_TERM) {
		int len = tstrips[0];
		for (int i = 1; i < tstrips.size(); i++) {
			if (len) {
				tstrips[i-1] = tstrips[i];
				len--;
			} else {
				tstrips[i-1] = -1;
				len = tstrips[i];
			}
		}
		tstrips.back() = -1;
	} else {
		int len = 0;
		for (int i = tstrips.size() - 2; i >= 0; i--) {
			if (tstrips[i] == -1) {
				tstrips[i+1] = len;
				len = 0;
			} else {
				tstrips[i+1] = tstrips[i];
				len++;
			}
		}
		tstrips[0] = len;
		//collect_tris_in_strips(tstrips);
	}
}

// Collect all single triangles to be at the end of the list of tstrips
static void collect_tris_in_strips(vector<int> &tstrips)
{
	if (tstrips.empty())
		return;
	vector<int> tris;

	int n = 0, offset = 0;
	bool have_tri = false, bad_strip = false;
	for (int i = 0; i < tstrips.size(); i++) {
		if (n == 0) {
			n = tstrips[i];
			bad_strip = (n < 3);
			have_tri = (n == 3);
			n++;
		}
		if (bad_strip) {
			offset++;
		} else if (have_tri) {
			tris.push_back(tstrips[i]);
			offset++;
		} else if (offset > 0) {
			tstrips[i-offset] = tstrips[i];
		}
		n--;
	}
	if (offset == 0)
		return;
	tstrips.erase(tstrips.end() - offset, tstrips.end());
	tstrips.insert(tstrips.end(), tris.begin(), tris.end());
}

// Compute per-vertex normals
void TriMesh::need_normals()
{
	if (normals.size() == vertices.size())
		return;
	need_faces();

	dprintf("Computing normals... ");

	normals.clear();
	normals.resize(vertices.size());

	int nf = faces.size(), nv = vertices.size();
	for (int i = 0; i < nf; i++) {
		const point &p0 = vertices[faces[i][0]];
		const point &p1 = vertices[faces[i][1]];
		const point &p2 = vertices[faces[i][2]];
		vec a = p0-p1, b = p1-p2, c = p2-p0;
		float l2a = len2(a), l2b = len2(b), l2c = len2(c);
		vec facenormal = a CROSS b;
		normals[faces[i][0]] += facenormal * (1.0f / (l2a * l2c));
		normals[faces[i][1]] += facenormal * (1.0f / (l2b * l2a));
		normals[faces[i][2]] += facenormal * (1.0f / (l2c * l2b));
	}

	for (i = 0; i < nv; i++)
		normalize(normals[i]);
}

// Find the direct neighbors of each vertex
void TriMesh::need_neighbors()
{
	if (!neighbors.empty())
		return;
	need_faces();

	//dprintf("Finding vertex neighbors... ");
	int nv = vertices.size(), nf = faces.size();

	vector<int> numneighbors(nv);
	for (int i = 0; i < nf; i++) {
		numneighbors[faces[i][0]]++;
		numneighbors[faces[i][1]]++;
		numneighbors[faces[i][2]]++;
	}

	neighbors.resize(nv);
	for (i = 0; i < nv; i++)
		neighbors[i].reserve(numneighbors[i]+2); // Slop for boundaries

	for (i = 0; i < nf; i++) {
		for (int j = 0; j < 3; j++) {
			vector<int> &me = neighbors[faces[i][j]];
			int n1 = faces[i][(j+1)%3];
			int n2 = faces[i][(j+2)%3];
			if (find(me.begin(), me.end(), n1) == me.end())
				me.push_back(n1);
			if (find(me.begin(), me.end(), n2) == me.end())
				me.push_back(n2);
		}
	}
}


// Find the faces touching each vertex
void TriMesh::need_adjacentfaces()
{
	if (!adjacentfaces.empty())
		return;
	need_faces();

	dprintf("Finding vertex to triangle maps... ");
	int nv = vertices.size(), nf = faces.size();

	vector<int> numadjacentfaces(nv);
	for (int i = 0; i < nf; i++) {
		numadjacentfaces[faces[i][0]]++;
		numadjacentfaces[faces[i][1]]++;
		numadjacentfaces[faces[i][2]]++;
	}

	adjacentfaces.resize(vertices.size());
	for (i = 0; i < nv; i++)
		adjacentfaces[i].reserve(numadjacentfaces[i]);

	for (i = 0; i < nf; i++) {
		for (int j = 0; j < 3; j++)
			adjacentfaces[faces[i][j]].push_back(i);
	}
}

// Find the face across each edge from each other face (-1 on boundary)
// If topology is bad, not necessarily what one would expect...
void TriMesh::need_across_edge()
{
	if (!across_edge.empty())
		return;
	need_adjacentfaces();

	dprintf("Finding across-edge maps... ");

	int nf = faces.size();
	across_edge.resize(nf, Face(-1,-1,-1));

	for (int i = 0; i < nf; i++) {
		for (int j = 0; j < 3; j++) {
			if (across_edge[i][j] != -1)
				continue;
			int v1 = faces[i][(j+1)%3];
			int v2 = faces[i][(j+2)%3];
			const vector<int> &a1 = adjacentfaces[v1];
			const vector<int> &a2 = adjacentfaces[v2];
			for (int k1 = 0; k1 < a1.size(); k1++) {
				int other = a1[k1];
				if (other == i)
					continue;
				std::vector<int>::const_iterator it =
					find(a2.begin(), a2.end(), other);
				if (it == a2.end())
					continue;
				int ind = (faces[other].indexof(v1)+1)%3;
				if (faces[other][(ind+1)%3] != v2)
					continue;
				across_edge[i][j] = other;
				across_edge[other][ind] = i;
				break;
			}
		}
	}
}


// Find axis-aligned bounding box of the vertices
void TriMesh::need_bbox()
{
	if (vertices.empty())
		return;

	dprintf("Computing bounding box... ");

	bbox.min = bbox.max = vertices[0];
	for (int i = 1; i < vertices.size(); i++) {
		if (vertices[i][0] < bbox.min[0])  bbox.min[0] = vertices[i][0];
		if (vertices[i][0] > bbox.max[0])  bbox.max[0] = vertices[i][0];
		if (vertices[i][1] < bbox.min[1])  bbox.min[1] = vertices[i][1];
		if (vertices[i][1] > bbox.max[1])  bbox.max[1] = vertices[i][1];
		if (vertices[i][2] < bbox.min[2])  bbox.min[2] = vertices[i][2];
		if (vertices[i][2] > bbox.max[2])  bbox.max[2] = vertices[i][2];
	}

	dprintf("Done.\n  x = %g .. %g, y = %g .. %g, z = %g .. %g\n",
		bbox.min[0], bbox.max[0], bbox.min[1],
		bbox.max[1], bbox.min[2], bbox.max[2]);
}


// Change this to #if 0 to enable a simpler (approximate) bsphere computation
// that does not use the Miniball code
#if 1
// Compute bounding sphere of the vertices.
void TriMesh::need_bsphere()
{
	if (vertices.empty() || bsphere.valid)
		return;

	dprintf("Computing bounding sphere... ");

	Miniball<3,float> mb;
	mb.check_in(vertices.begin(), vertices.end());
	mb.build();
	bsphere.center = mb.center();
	bsphere.r = sqrt(mb.squared_radius());
	bsphere.valid = true; 

	dprintf("Done.\n  center = (%g, %g, %g), radius = %g\n",
		bsphere.center[0], bsphere.center[1],
		bsphere.center[2], bsphere.r);
}
#else
// Find extreme vertex in a given direction
static int farthest_vertex_along(const TriMesh &t, const vec &dir)
{
	const vector<point> &v = t.vertices;
	int nv = v.size();

	int farthest = 0;
	float farthest_dot = v[0] DOT dir;