C++ Manifold类代码示例

/*
* Checks if an edge is a boundary and if not, checks it for being locally Delaunay.
* If it is not locally delaunay it is flipped and all its neighbor edges are checked, ad nauseam.
*/
void recursiveDelaunayFlip(Manifold &m, Walker w, bool isAffected) {
	if ( !boundary(m, w.halfedge()) ) {
		// Check if the current halfedge is locally Delaunay using the inCircle function
		Vec3d p1 = m.pos(w.opp().vertex());
		Vec3d p2 = m.pos(w.vertex());
		Vec3d p3 = m.pos(w.next().vertex());
		Vec3d p4 = m.pos(w.opp().next().vertex()); // This seems to return erroneus values every time
		if (isAffected == true) {
			cout << "Affected quadrillateral:" << endl;
			cout << "p1: " << p1 << endl;
			cout << "p2: " << p2 << endl;
			cout << "p3: " << p3 << endl;
			cout << "p4: " << p4 << endl;
		}
		if ( inCircle( p1, p3, p2, p4 ) || inCircle(p1, p2, p4, p3) ) {
			// Since either point was in a triangle circumcircle, flip the edge
			m.flip_edge(w.halfedge());
			cout << "Edge to be flipped: " << p1 << ", " << p2 << ". Other vertices: " << p3 << ", " << p4 << endl;
			// Recursively check all the edges that share a neighbour with the flipped edge
			recursiveDelaunayFlip(m, w.next(), true);
			recursiveDelaunayFlip(m, w.prev(), true);
			recursiveDelaunayFlip(m, w.opp().next(), true);
			recursiveDelaunayFlip(m, w.opp().prev(), true);
		}
	}
}

int WireframeRenderer::maximum_face_valency(const Manifold& m)
{
    int max_val = 0;
    for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f)
        max_val = max(max_val, no_edges(m, *f));
    return max_val;
}

AmbientOcclusionRenderer::AmbientOcclusionRenderer(const Manifold& m, bool smooth, VertexAttributeVector<double>& field, double max_val):
    SimpleShaderRenderer(vss,fss)
{
    GLint old_prog;
    glGetIntegerv(GL_CURRENT_PROGRAM, &old_prog);
    glUseProgram(prog);

    GLuint scalar_attrib = glGetAttribLocation(prog, "scalar");
    glUniform1fARB(glGetUniformLocationARB(prog, "scalar_max"), max_val);

    glNewList(display_list,GL_COMPILE);

    for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f) {

        if(!smooth)
            glNormal3dv(normal(m, *f).get());
        if(no_edges(m, *f)== 3)
            glBegin(GL_TRIANGLES);
        else
            glBegin(GL_POLYGON);

        for(Walker w = m.walker(*f); !w.full_circle(); w = w.circulate_face_ccw())
        {
            Vec3d n(normal(m, w.vertex()));
            if(smooth)
                glNormal3dv(n.get());
            glVertexAttrib1d(scalar_attrib, field[w.vertex()]);
            glVertex3dv(m.pos(w.vertex()).get());
        }
        glEnd();
    }
    glEndList();
    glUseProgram(old_prog);

}

Manifold::GeodesicPtr Compound::PointOfReference::getFinalGeodesic(Vector3d vector) {
	/*Manifold* space = pointOfReference->getPosition()->getSpace();
	return space->getGeodesic(pointOfReference, vector);*/
	//std::cout << "getFinalGeodesic(vector)" << std::endl;
	assert(vector == vector);
	//This will need to be made more sophisticated when Compound is made to support more than one space.
	Manifold* space = pointOfReference->getPosition()->getSpace();
	Manifold::GeodesicPtr next = space->getGeodesic(pointOfReference, vector);
	//Manifold::GeodesicPtr original = next;
	//Vector3d firstVector = next->getEndPoint()->getVector();
	Manifold::GeodesicPtr current;
	do {
		current = next;
		space = current->getSpace();
		next = space->nextPiece(current);
		/*if(next) {
			std::cout << "next" << std::endl;
		}*/
	} while(next);
	//assert((current->getEndPoint()->getVector() - firstVector).squaredNorm() < EPSILON*EPSILON);	//For portals leading to themselves.
	//std::cout << "current->getEndPoint()->getVector():\n" << current->getEndPoint()->getVector() << std::endl;
	//std::cout << "firstVector:\n" << firstVector << std::endl;
	//return original;
	assert(current->getEndPoint()->isInManifold());
	return current;
}

 void quadric_simplify(Manifold& m, double keep_fraction, double singular_thresh, bool choose_optimal_positions)
 {
     gel_srand(1210);
     long int F = m.no_faces();
     VertexAttributeVector<int> mask(m.no_faces(), 0);
     long int max_work = max(static_cast<long int>(0), F- static_cast<long int>(keep_fraction * F));
     QuadricSimplifier qsim(m, mask, singular_thresh, choose_optimal_positions);
     qsim.reduce(max_work);
 }

void create_single_triangle_manifold(const Vec3f& p1, 
																		 const Vec3f& p2, 
																		 const Vec3f& p3, 
																		 Manifold& mani)
{
	// Create vector of vertices
	vector<Vec3f> vertices(3);
	vertices[0] = p1;
	vertices[1] = p2;
	vertices[2] = p3;

	// Create vector of faces. Each element corresponds to a face and tells
	// how many vertices that face contains. In the case of a triangle
	// mesh, each face has three vertices.
	vector<int> faces(1);
	faces[0] = 3;

	// Create the index vector. Each element is an index into the vertex list
	// 
	vector<int> indices(3);
	indices[0]=0;
	indices[1]=1;
	indices[2]=2;

  mani.build(3,           // Number of vertices.
	  				 reinterpret_cast<float*>(&vertices[0]),// Pointer to vertices.
						 1,           // Number of faces.
						 &faces[0],   // Pointer to faces.
						 &indices[0]);// Pointer to indices.


}

    bool obj_save(const string& filename, Manifold& m)
    {
        ofstream os(filename.data());
        if(os.bad())
            return false;

        VertexAttributeVector<int> vmap;
        int k = 0;
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
            Vec3d p = m.pos(*v);
            os << "v "<< p[0] << " " << p[1] << " " << p[2] << "\n";
            vmap[*v] = k++;
        }

        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){        
            vector<int> verts;
            for(Walker w = m.walker(*f); !w.full_circle(); w = w.circulate_face_ccw()){
                int idx = vmap[w.vertex()];			
                assert(static_cast<size_t>(idx) < m.no_vertices());
                // move subscript range from 0..size-1 to 1..size according to OBJ standards
                verts.push_back(idx + 1);
            }
			os << "f ";
            for(size_t i = 0; i < verts.size() ; ++i){
                os << verts[i] << " ";
            }
			os<<endl;
        }

        return true;
    }

void dual(HMesh::Manifold& m)
{
    FaceAttributeVector<Vec3d> face_center(m.no_faces());
    for (auto f : m.faces()) {
        // find mid point
        Vec3d mpt(0.0, 0.0, 0.0);
        int nb_p = 0;
        for (auto hw = m.walker(f); !hw.full_circle(); hw = hw.circulate_face_ccw()) {
            mpt += m.pos(hw.vertex());
            nb_p++;
        }
        mpt = mpt / nb_p;
        
        
        face_center[f] = mpt;
    }
    
    Manifold newMesh;
    for (auto v : m.vertices()) {
        vector<Vec3d> pts;
        for (auto hw = m.walker(v); !hw.full_circle(); hw = hw.circulate_vertex_ccw()) {
//            if (hw.opp().face() == InvalidFaceID)
//            {
//                pts.push_back((m.pos(hw.vertex()) + m.pos(hw.opp().vertex()))/2.0);
//                pts.push_back(face_center[hw.face()]);
//            }
//            else if(hw.face() == InvalidFaceID)
//            {
//                pts.push_back((m.pos(hw.vertex()) + m.pos(hw.opp().vertex()))/2.0);
//            }
//            else
            if (m.in_use(hw.face()))
            {
                pts.push_back(face_center[hw.face()]);
            }
        }
        
        newMesh.add_face(pts);
    }
    
    stitch_mesh(newMesh, 0.01);
    
    m = newMesh;
 
}

    bool x3d_load(const string& filename, Manifold& m) 
    {
        faces.clear();
        indices.clear();
        vertices.clear();

        Timer tim;
        tim.start();

        string baseurl;
        int idx = max(find_last_of(filename, "\\"), 
            find_last_of(filename, "/"));

        if(idx != -1)
            baseurl = string(filename.substr(0, idx+1));

        XmlDoc x3d_doc(filename.c_str());
        
        if(!x3d_doc.is_valid())
            return false;
        
        x3d_doc.add_handler("Shape", handle_Shape);    
        x3d_doc.add_handler("IndexedFaceSet", handle_IndexedFaceSet);
        x3d_doc.add_handler("Coordinate", handle_Coordinate);
        x3d_doc.process_elements();
        x3d_doc.close();
        
        cout << "vertices " << vertices.size() << endl;

        m.build(vertices.size()/3, 
            reinterpret_cast<float*>(&vertices[0]), 
            faces.size(), 
            &faces[0], 
            &indices[0]);
        
        cout << " Loading took " << tim.get_secs() << endl;
        return true;
    }

    void dual(Manifold& m)
    {
        // Create new vertices. Each face becomes a vertex whose position
        // is the centre of the face
        int i = 0;
        FaceAttributeVector<int> ftouched;
        vector<Vec3d> vertices;
        vertices.resize(m.no_faces());
        for(auto f : m.faces())
            vertices[ftouched[f] = i++] = centre(m, f);
        
        // Create new faces. Each vertex is a new face with N=valency of vertex
        // edges.
        vector<int> faces;
        vector<int> indices;
        for(auto v : m.vertices())
            if(valency(m, v) > 2 && !(boundary(m, v)))
            {
//				int N = circulate_vertex_ccw(m, v, (std::function<void(FaceID)>)[&](FaceID fid) {
//                    indices.push_back(ftouched[fid]);
//                });

                Walker w = m.walker(v);
                for(; !w.full_circle(); w = w.circulate_vertex_ccw()){
                    indices.push_back(ftouched[w.face()]);
                }
                int N = w.no_steps();
                // Insert face valency in the face vector.
                faces.push_back(N);
            }
        
        // Clear the manifold before new geometry is inserted.
        m.clear();
        
        // And build
        m.build(    vertices.size(),
                reinterpret_cast<double*>(&vertices[0]),
                faces.size(),
                &faces[0],
                &indices[0]);
    }

int main(int argc, char** argv)
{	
	// LOAD OBJ
    Manifold m;
    if(argc>1)
	{
		ArgExtracter ae(argc, argv);
		
		do_aabb = ae.extract("-A");
		do_obb = ae.extract("-O");
		ae.extract("-x", vol_dim[0]);
		ae.extract("-y", vol_dim[1]);
		ae.extract("-z", vol_dim[2]);
		do_ray_tests = ae.extract("-R");
		flip_normals = ae.extract("-f");
		string file = ae.get_last_arg();
        
        cout << "loading " << file << "... " << flush; 
		load(file, m);
        cout << " done" << endl;
	}
    else
	{
		string fn("../../data/bunny-little.x3d");
		x3d_load(fn, m);
	}
	cout << "Volume dimensions " << vol_dim << endl;
	if(!valid(m))
	{
		cout << "Not a valid manifold" << endl;
		exit(0);
	}
	triangulate_by_edge_face_split(m);
	
	Vec3d p0, p7;
	bbox(m, p0, p7);
	
	Mat4x4d T = fit_bounding_volume(p0,p7,10);
    
    cout << "Transformation " << T << endl;
	
	for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v)
		m.pos(*v) = T.mul_3D_point(m.pos(*v));
	
	
 	RGridf grid(vol_dim,FLT_MAX);
	Util::Timer tim;
	
	
	float T_build_obb=0, T_build_aabb=0, T_dist_obb=0, 
		T_dist_aabb=0, T_ray_obb=0, T_ray_aabb=0;
	
	if(do_obb)
	{
		cout << "Building OBB Tree" << endl;
		tim.start();
		OBBTree obb_tree;
		build_OBBTree(m, obb_tree);
		T_build_obb = tim.get_secs();
		
		cout << "Computing distances from OBB Tree" << endl;
		tim.start();
		DistCompCache<OBBTree> dist(&obb_tree);
		for_each_voxel(grid, dist);
		T_dist_obb = tim.get_secs();
		
		cout << "Saving distance field" << endl;
		save_raw_float("obb_dist.raw", grid);
		
		if(do_ray_tests)
		{
			cout << "Ray tests on OBB Tree" << endl;
			tim.start();
			RayCast<OBBTree> ray(&obb_tree);
			for_each_voxel(grid, ray);
			T_ray_obb = tim.get_secs();
			
			cout << "Saving ray volume" << endl;
			save_raw_float("obb_ray.raw", grid);
		}
	}
	
	if(do_aabb)
	{
		cout << "Building AABB Tree" << endl;
		tim.start();
		AABBTree aabb_tree;
		build_AABBTree(m, aabb_tree);
		T_build_aabb = tim.get_secs();
		
		cout << "Computing distances from AABB Tree" << endl;
		tim.start();
		DistCompCache<AABBTree> dist(&aabb_tree);
		for_each_voxel(grid, dist);
		T_dist_aabb = tim.get_secs();
		
		cout << "Saving distance field" << endl;
		save_raw_float("aabb_dist.raw", grid);
		
		if(do_ray_tests)
//.........这里部分代码省略.........

void mean_curvature_smooth(Manifold& m, bool implicit, double lambda)
{
    using EigMat = SparseMatrix<double>;
    using EigVec = VectorXd;
    
    int N = (int)m.no_vertices();
    VertexAttributeVector<int> indices(m.allocated_vertices());
    VertexAttributeVector<double> areas(m.allocated_vertices());
    int i=0;
    for(auto v: m.vertices()) {
        indices[v] = i++;
        areas[v] = mixed_area(m, v);
    }
    

    EigMat K(N,N);      // Sparse matrix initialized with 0
    EigVec X(N),Y(N),Z(N);
    EigVec Xp(N), Yp(N), Zp(N);
    
    //-----------------------------------------------------------
    // Student implementation
    //-----------------------------------------------------------
    double epsilon = 1e-5;
    for (auto vkey : m.vertices())
    {
        int i = indices[vkey];
        for (auto w = m.walker(vkey); !w.full_circle(); w = w.circulate_vertex_ccw())
        {
            int j = indices[w.vertex()];
            assert(i != j);
            
            if (i > j
                or w.face() == HMesh::InvalidFaceID
                or w.opp().face() == HMesh::InvalidFaceID)
            {
                continue; // Avoid recomputation
            }
            
            auto pi = m.pos(w.opp().vertex());
            auto pj = m.pos(w.vertex());
            auto pl = m.pos(w.opp().next().vertex());
            auto pk = m.pos(w.next().vertex());
            
            double cot_alpha_ij = dot(pj - pk, pi - pk) /
                                ( cross(pi - pk, pj - pk).length() + epsilon);
            double cot_beta_ij = dot(pj - pl, pi - pl) /
                                ( cross(pi - pl, pj - pl).length() + epsilon);
            
            double Ai = areas[w.opp().vertex()];
            double Aj = areas[w.vertex()];
            
            double Lij = (cot_alpha_ij + cot_beta_ij)
                        / sqrt(Ai*Aj + epsilon);
            
            K.coeffRef(i, j) = Lij;
            K.coeffRef(j, i) = Lij;
            K.coeffRef(i, i) -= Lij;
            K.coeffRef(j, j) -= Lij;
            
        }
    }
    
    EigMat I(N,N);
    for (int i = 0; i < N; i++)
    {
        I.coeffRef(i, i) = 1;
    }
    
    K = I - K*lambda;
    
    
    for (auto vkey : m.vertices())
    {
        auto p = m.pos(vkey);
        int i = indices[vkey];
        X.coeffRef(i) = p[0];
        Y.coeffRef(i) = p[1];
        Z.coeffRef(i) = p[2];
    }
    
    // Solve
    SimplicialLLT<EigMat> solver(K);
    Xp = solver.solve(X);
    Yp = solver.solve(Y);
    Zp = solver.solve(Z);
    
    // End student implementation
    //-----------------------------------------------------------
        
    for(auto v: m.vertices())
    {
        int i = indices[v];
        m.pos(v) = Vec3d(Xp[i], Yp[i], Zp[i]);
    }
}

bool Manifold::isSameType(const Manifold& other) const
{
  return this->getTypeId() == other.getTypeId();
}

#include <HMesh/triangulate.h>

#include <GLGraphics/gel_glut.h>
#include <CGLA/Vec2d.h>
#include <CGLA/Mat3x3d.h>
#include <CGLA/Mat3x3f.h>
#include <CGLA/Mat4x4f.h>

using namespace std;
using namespace CGLA;
using namespace HMesh;

// The range of the input data.
Vec2d dmin(99e99), dmax(-99e99);

Manifold m; // The triangle mesh data structure.
HalfEdgeIDIterator flipper = m.halfedges_begin(); // The halfedge we try to flip.
HalfEdgeAttributeVector<int> touched;

/*
 * Draw the triangle mesh. This function is called from GLUT (a library
 * which enables OpenGL drawing in windows.)
 */

void display()
{
	// Set up correct OpenGL projection
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluOrtho2D(dmin[0], dmax[0], dmin[1], dmax[1]);
	glMatrixMode(GL_MODELVIEW);

std::shared_ptr<Manifold> Manifold::copyManifold(const Manifold& m)
{
  std::shared_ptr<Manifold> copy = m.getNewCopy();
  return copy;
}