trimesh.cpp

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

	for (int i = 1; i < nv; i++) {
		float my_dot = v[i] DOT dir;
		if (my_dot > farthest_dot)
			farthest = i, farthest_dot = my_dot;
	}
	return farthest;
}
// Approximate bounding sphere code based on an algorithm by Ritter
void TriMesh::need_bsphere()
{
	if (vertices.empty())
		return;

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

	point best_min, best_max;
	vector<vec> dirs;
	dirs.push_back(vec(1,0,0));
	dirs.push_back(vec(0,1,0));
	dirs.push_back(vec(0,0,1));
	dirs.push_back(vec(1,1,1));
	dirs.push_back(vec(1,-1,1));
	dirs.push_back(vec(1,-1,-1));
	dirs.push_back(vec(1,1,-1));
	for (int i = 0; i < dirs.size(); i++) {
		point p1 = vertices[farthest_vertex_along(*this, -dirs[i])];
		point p2 = vertices[farthest_vertex_along(*this,  dirs[i])];
		if (dist2(p1, p2) > dist2(best_min, best_max)) {
			best_min = p1;
			best_max = p2;
		}
	}
	bsphere.center = 0.5f * (best_min + best_max);
	bsphere.r = 0.5f * dist(best_min, best_max);
	float r2 = sqr(bsphere.r);

	// Expand bsphere to contain all points
	for (int i = 0; i < vertices.size(); i++) {
		float d2 = dist2(vertices[i], bsphere.center);
		if (d2 <= r2)
			continue;
		float d = sqrt(d2);
		bsphere.r = 0.5f * (bsphere.r + d);
		r2 = sqr(bsphere.r);
		bsphere.center -= vertices[i];
		bsphere.center *= bsphere.r / d;
		bsphere.center += vertices[i];
	}

	dprintf("Done.\n  center = (%g, %g, %g), radius = %g\n",
		bsphere.center[0], bsphere.center[1],
		bsphere.center[2], bsphere.r);
}
#endif

// Remove the indicated vertices from the TriMesh.
//
// Faces are renumbered to reflect the new numbering of vertices, and any
// faces that included a vertex to be removed will also be removed.
//
// Any per-vertex properties are renumbered along with the vertices.
// Any triangle strips, neighbor lists, or adjacency lists are deleted.
void TriMesh::remove_vertices(TriMesh *mesh, const vector<bool> &toremove)
{
	int numvertices = mesh->vertices.size();

	// Build a table that tells how the vertices will be remapped
	if (!numvertices)
		return;

	TriMesh::dprintf("Removing vertices... ");
	vector<int> remap_table(numvertices);
	int next = 0;
	int i;
	for (i = 0; i < numvertices; i++) {
		if (toremove[i])
			remap_table[i] = -1;
		else
			remap_table[i] = next++;
	}

	// Nothing to delete?
	if (next == numvertices) {
		TriMesh::dprintf("None removed.\n");
		return;
	}

	// Nuke connectivity lists, as well as anything we can re-compute
	bool have_col = !mesh->colors.empty();
	bool have_conf = !mesh->confidences.empty();
	bool have_flags = !mesh->flags.empty();
	bool have_normals = !mesh->normals.empty();
	bool have_pdir1 = !mesh->pdir1.empty();
	bool have_pdir2 = !mesh->pdir2.empty();
	bool have_curv1 = !mesh->curv1.empty();
	bool have_curv2 = !mesh->curv2.empty();
	bool have_dcurv = !mesh->dcurv.empty();
	mesh->cornerareas.clear(); mesh->pointareas.clear();
	bool had_tstrips = !mesh->tstrips.empty();
	bool had_faces = !mesh->faces.empty();
	mesh->need_faces();
	mesh->tstrips.clear();
	mesh->adjacentfaces.clear();
	mesh->neighbors.clear();
	mesh->across_edge.clear();

	// Remap the vertices and per-vertex properties
#define REMAP(property) mesh->property[remap_table[i]] = mesh->property[i]
	for (i = 0; i < numvertices; i++) {
		if (remap_table[i] == -1 || remap_table[i] == i)
			continue;
		REMAP(vertices);
		if (have_col) REMAP(colors);
		if (have_conf) REMAP(confidences);
		if (have_flags) REMAP(flags);
		if (have_normals) REMAP(normals);
		if (have_pdir1) REMAP(pdir1);
		if (have_pdir2) REMAP(pdir2);
		if (have_curv1) REMAP(curv1);
		if (have_curv2) REMAP(curv2);
		if (have_dcurv) REMAP(dcurv);
	}
	TriMesh::dprintf("%d vertices removed... Done.\n", numvertices - next);

#define ERASE(property) mesh->property.erase(mesh->property.begin() + next, \
					     mesh->property.end())
	ERASE(vertices);
	if (have_col) ERASE(colors);
	if (have_conf) ERASE(confidences);
	if (have_flags) ERASE(flags);
	if (have_normals) ERASE(normals);
	if (have_pdir1) ERASE(pdir1);
	if (have_pdir2) ERASE(pdir2);
	if (have_curv1) ERASE(curv1);
	if (have_curv2) ERASE(curv2);
	if (have_dcurv) ERASE(dcurv);

	// Renumber faces
	int nextface = 0;
	for (i = 0; i < mesh->faces.size(); i++) {
		int n0 = (mesh->faces[nextface][0] = remap_table[mesh->faces[i][0]]);
		int n1 = (mesh->faces[nextface][1] = remap_table[mesh->faces[i][1]]);
		int n2 = (mesh->faces[nextface][2] = remap_table[mesh->faces[i][2]]);
		if ((n0 != -1) && (n1 != -1) && (n2 != -1))
			nextface++;
	}
	mesh->faces.erase(mesh->faces.begin() + nextface, mesh->faces.end());
	if (had_tstrips)
		mesh->need_tstrips();
	if (!had_faces)
		mesh->faces.clear();

	mesh->bsphere.valid = false;
}

// Remove vertices that aren't referenced by any face
void TriMesh::remove_unused_vertices(TriMesh *mesh)
{
	int numvertices = mesh->vertices.size();
	if (!numvertices)
		return;

	bool had_faces = !mesh->faces.empty();
	mesh->need_faces();
	vector<bool> unused(numvertices, true);
	for (int i = 0; i < mesh->faces.size(); i++) {
		unused[mesh->faces[i][0]] = false;
		unused[mesh->faces[i][1]] = false;
		unused[mesh->faces[i][2]] = false;
	}
	remove_vertices(mesh, unused);
	if (!had_faces)
		mesh->faces.clear();
}

// Remove faces as indicated by toremove.  Should probably be
// followed by a call to remove_unused_vertices()
void TriMesh::remove_faces(TriMesh *mesh, const vector<bool> &toremove)
{
	bool had_tstrips = !mesh->tstrips.empty();
	bool had_faces = !mesh->faces.empty();
	mesh->need_faces();
	int numfaces = mesh->faces.size();
	if (!numfaces)
		return;

	mesh->tstrips.clear();
	mesh->adjacentfaces.clear();
	mesh->neighbors.clear();
	mesh->across_edge.clear();
	mesh->cornerareas.clear();
	mesh->pointareas.clear();

	TriMesh::dprintf("Removing faces... ");
	int next = 0;
	for (int i = 0; i < numfaces; i++) {
		if (toremove[i])
			continue;
		mesh->faces[next++] = mesh->faces[i];
	}
	if (next == numfaces) {
		TriMesh::dprintf("None removed.\n");
		return;
	}

	mesh->faces.erase(mesh->faces.begin() + next, mesh->faces.end());
	TriMesh::dprintf("%d faces removed... Done.\n", numfaces - next);

	if (had_tstrips)
		mesh->need_tstrips();
	if (!had_faces)
		mesh->faces.clear();

	mesh->bsphere.valid = false;
}

// Remove long, skinny faces.  Should probably be followed by a
// call to remove_unused_vertices()
void TriMesh::remove_sliver_faces(TriMesh *mesh)
{
	bool had_faces = !mesh->faces.empty();
	mesh->need_faces();
	int numfaces = mesh->faces.size();

	const float l2thresh = sqr(4.0f * mesh->feature_size());
	const float cos2thresh = 0.85f;
	vector<bool> toremove(numfaces, false);
	for (int i = 0; i < numfaces; i++) {
		const point &v0 = mesh->vertices[mesh->faces[i][0]];
		const point &v1 = mesh->vertices[mesh->faces[i][1]];
		const point &v2 = mesh->vertices[mesh->faces[i][2]];
		float d01 = dist2(v0, v1);
		float d12 = dist2(v1, v2);
		float d20 = dist2(v2, v0);
		if (d01 < l2thresh && d12 < l2thresh && d20 < l2thresh)
			continue;
		// c2 is square of cosine of smallest angle
		float m = __min(__min(d01,d12),d20);
		float c2 = sqr(d01+d12+d20-2.0f*m) * m/(4.0f*d01*d12*d20);
		if (c2 < cos2thresh)
			continue;
		toremove[i] = true;
	}
	remove_faces(mesh, toremove);
	if (!had_faces)
		mesh->faces.clear();

	mesh->bsphere.valid = false;
}

// A characteristic "feature size" for the mesh.  Computed as an approximation
// to the median edge length
float TriMesh::feature_size()
{
	need_faces();
	if (faces.empty())
		return 0.0f;

	int nf = faces.size();
	int nsamp = __min(nf / 2, 333);

	vector<float> samples;
	samples.reserve(nsamp * 3);

	for (int i = 0; i < nsamp; i++) {
		// Quick 'n dirty portable random number generator
		static unsigned randq = 0;
		randq = unsigned(1664525) * randq + unsigned(1013904223);

		int ind = randq % nf;
		const point &p0 = vertices[faces[ind][0]];
		const point &p1 = vertices[faces[ind][1]];
		const point &p2 = vertices[faces[ind][2]];
		samples.push_back(dist2(p0,p1));
		samples.push_back(dist2(p1,p2));
		samples.push_back(dist2(p2,p0));
	}
	std::nth_element(samples.begin(),
		    samples.begin() + samples.size()/2,
		    samples.end());
	return sqrt(samples[samples.size()/2]);
}

//----------------------------Mesh smoothing-----------------------------//
// One iteration of umbrella-operator smoothing
void TriMesh::umbrella(TriMesh *mesh, float stepsize)
{
	mesh->need_neighbors();
	mesh->need_adjacentfaces();
	int nv = mesh->vertices.size();
	vector<vec> disp(nv);
	for (int i = 0; i < nv; i++) {
		if (mesh->is_bdy(i)) {
			// Change to #if 1 to smooth boundaries.
			// This way, we leave boundaries alone.
#if 0
			int nn = mesh->neighbors[i].size();
			int nnused = 0;
			if (!nn)
				continue;
			for (int j = 0; j < nn; j++) {
				if (!mesh->is_bdy(mesh->neighbors[i][j]))
					continue;
				disp[i] += mesh->vertices[mesh->neighbors[i][j]];
				nnused++;
			}
			disp[i] /= nnused;
			disp[i] -= mesh->vertices[i];
#else
			disp[i].clear();
#endif
		} else {
			int nn = mesh->neighbors[i].size();
			if (!nn)
				continue;
			for (int j = 0; j < nn; j++)
				disp[i] += mesh->vertices[mesh->neighbors[i][j]];
			disp[i] /= nn;
			disp[i] -= mesh->vertices[i];
		}
	}
	for (i = 0; i < nv; i++)
		mesh->vertices[i] += stepsize * disp[i];

	mesh->bsphere.valid = false;
}

// Several iterations of Taubin lambda/mu
void TriMesh::lmsmooth(TriMesh *mesh, int niters)
{
	mesh->need_neighbors();
	mesh->need_adjacentfaces();
	//riMesh::dprintf("Smoothing mesh... ");
	for (int i = 0; i < niters; i++) {
		umbrella(mesh, 0.330f);
		umbrella(mesh, -0.331f);
	}

	mesh->bsphere.valid = false;
}

//-------------------------Curvature computation------------------------//
// Rotate a coordinate system to be perpendicular to the given normal
static void rot_coord_sys(const vec &old_u, const vec &old_v,
			  const vec &new_norm,
			  vec &new_u, vec &new_v)
{
	new_u = old_u;
	new_v = old_v;
	vec old_norm = old_u CROSS old_v;
	float ndot = old_norm DOT new_norm;
	if (unlikely(ndot <= -1.0f)) {
		new_u = -new_u;
		new_v = -new_v;
		return;
	}
	vec perp_old = new_norm - ndot * old_norm;
	vec dperp = 1.0f / (1 + ndot) * (old_norm + new_norm);

	new_u -= dperp * (new_u DOT perp_old);
	new_v -= dperp * (new_v DOT perp_old);
}


// Reproject a curvature tensor from the basis spanned by old_u and old_v
// (which are assumed to be unit-length and perpendicular) to the
// new_u, new_v basis.
void proj_curv(const vec &old_u, const vec &old_v,
			   float old_ku, float old_kuv, float old_kv,
					   const vec &new_u, const vec &new_v,
					   float &new_ku, float &new_kuv, float &new_kv)
{
	vec r_new_u, r_new_v;
	rot_coord_sys(new_u, new_v, old_u CROSS old_v, r_new_u, r_new_v);

	float u1 = r_new_u DOT old_u;
	float v1 = r_new_u DOT old_v;
	float u2 = r_new_v DOT old_u;
	float v2 = r_new_v DOT old_v;
	new_ku  = old_ku * u1*u1 + old_kuv * (2.0f  * u1*v1) + old_kv * v1*v1;
	new_kuv = old_ku * u1*u2 + old_kuv * (u1*v2 + u2*v1) + old_kv * v1*v2;
	new_kv  = old_ku * u2*u2 + old_kuv * (2.0f  * u2*v2) + old_kv * v2*v2;
}


// Like the above, but for dcurv
void proj_dcurv(const vec &old_u, const vec &old_v,
						const Vec<4> old_dcurv,
						const vec &new_u, const vec &new_v,
						Vec<4> &new_dcurv)
{
	vec r_new_u, r_new_v;
	rot_coord_sys(new_u, new_v, old_u CROSS old_v, r_new_u, r_new_v);

	float u1 = r_new_u DOT old_u;
	float v1 = r_new_u DOT old_v;
	float u2 = r_new_v DOT old_u;
	float v2 = r_new_v DOT old_v;

	new_dcurv[0] = old_dcurv[0]*u1*u1*u1 +
		       old_dcurv[1]*3.0f*u1*u1*v1 +
		       old_dcurv[2]*3.0f*u1*v1*v1 +
		       old_dcurv[3]*v1*v1*v1;
	new_dcurv[1] = old_dcurv[0]*u1*u1*u2 +
		       old_dcurv[1]*(u1*u1*v2 + 2.0f*u2*u1*v1) +
		       old_dcurv[2]*(u2*v1*v1 + 2.0f*u1*v1*v2) +
		       old_dcurv[3]*v1*v1*v2;
	new_dcurv[2] = old_dcurv[0]*u1*u2*u2 +
		       old_dcurv[1]*(u2*u2*v1 + 2.0f*u1*u2*v2) +
		       old_dcurv[2]*(u1*v2*v2 + 2.0f*u2*v2*v1) +
		       old_dcurv[3]*v1*v2*v2;
	new_dcurv[3] = old_dcurv[0]*u2*u2*u2 +
		       old_dcurv[1]*3.0f*u2*u2*v2 +
		       old_dcurv[2]*3.0f*u2*v2*v2 +
		       old_dcurv[3]*v2*v2*v2;
}


// Given a curvature tensor, find principal directions and curvatures
// Makes sure that pdir1 and pdir2 are perpendicular to normal
void diagonalize_curv(const vec &old_u, const vec &old_v,
							  float ku, float kuv, float kv,
							  const vec &new_norm,
							  vec &pdir1, vec &pdir2, float &k1, float &k2)
{
	vec r_old_u, r_old_v;
	rot_coord_sys(old_u, old_v, new_norm, r_old_u, r_old_v);

	float c = 1, s = 0, tt = 0;
	if (likely(kuv != 0.0f)) {
		// Jacobi rotation to diagonalize
		float h = 0.5f * (kv - ku) / kuv;
		tt = (h < 0.0f) ?
			1.0f / (h - sqrt(1.0f + h*h)) :
			1.0f / (h + sqrt(1.0f + h*h));
		c = 1.0f / sqrt(1.0f + tt*tt);
		s = tt * c;
	}

	k1 = ku - tt * kuv;
	k2 = kv + tt * kuv;

	if (fabs(k1) >= fabs(k2)) {
		pdir1 = c*r_old_u - s*r_old_v;
	} else {
		swap(k1, k2);
		pdir1 = s*r_old_u + c*r_old_v;
	}
	pdir2 = new_norm CROSS pdir1;
}

// Compute principal curvatures and directions.
void TriMesh::need_curvatures()
{
	int i,j;
	if (curv1.size() == vertices.size())
		return;
	need_faces();
	need_normals();
	need_pointareas();

	//dprintf("Computing curvatures... ");

	// Resize the arrays we'll be using
	int nv = vertices.size(), nf = faces.size();
	curv1.clear(); curv1.resize(nv); curv2.clear(); curv2.resize(nv);
	pdir1.clear(); pdir1.resize(nv); pdir2.clear(); pdir2.resize(nv);
	vector<float> curv12(nv);

	// Set up an initial coordinate system per vertex
	for (i = 0; i < nf; i++) {
		pdir1[faces[i][0]] = vertices[faces[i][1]] -
				     vertices[faces[i][0]];
		pdir1[faces[i][1]] = vertices[faces[i][2]] -
				     vertices[faces[i][1]];
		pdir1[faces[i][2]] = vertices[faces[i][0]] -
				     vertices[faces[i][2]];
	}
	for (i = 0; i < nv; i++) {
		pdir1[i] = pdir1[i] CROSS normals[i];
		normalize(pdir1[i]);
		pdir2[i] = normals[i] CROSS pdir1[i];
	}

	
	LinEqn<float,3> leqn;
	// Compute curvature per-face
	for (i = 0; i < nf; i++) {
		// Edges
		vec e[3] = { vertices[faces[i][2]] - vertices[faces[i][1]],
			     vertices[faces[i][0]] - vertices[faces[i][2]],
			     vertices[faces[i][1]] - vertices[faces[i][0]] };

		// N-T-B coordinate system per face
		vec t = e[0];
		normalize(t);
		vec n = e[0] CROSS e[1];
		vec b = n CROSS t;
		normalize(b);

		float w[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
		float m[3] = {0,0,0};

		// Estimate curvature based on variation of normals
		// along edges
		for (j = 0; j < 3; j++) {