Mesh.cpp

Go to the documentation of this file.
00001 
00006 #include "Mesh.h"
00007 #include <iostream>
00008 #include <fstream>
00009 #include "util.h"
00010 #include <stack>
00011 #include <iomanip>
00012 #include <algorithm>
00013 
00014 using namespace std;
00015 
00016 Mesh::Mesh(void)
00017 {
00018    numBoundaryVertices = 0;
00019 }
00020 
00021 Mesh::~Mesh(void)
00022 {
00023    clear();
00024 }
00025 
00026 void Mesh::clear() {
00027    VertexIterator vIter = vertexIterator();
00028    FaceIterator fIter = faceIterator();
00029    HalfEdgeIterator eIter = halfEdgeIterator();
00030 
00031    while (!vIter.end()) {
00032       delete vIter.vertex();
00033       vIter++;
00034    }
00035    while (!fIter.end()) {
00036       delete fIter.face();
00037       fIter++;
00038    }
00039    while (!eIter.end()) {
00040       delete eIter.half_edge();
00041       eIter++;
00042    }
00043 }
00044 
00045 Vertex * Mesh::addVertex(const Vector3 & _p) {
00046    Vertex * v = new Vertex();
00047    v->p = _p;
00048    v->ID = (int)vertices.size();
00049    vertices.push_back(v);
00050    return v;
00051 }
00052 
00053 Face * Mesh::addFace(std::vector<int> faceVerts) {
00054    Face * f = new Face();
00055    f->ID = (int)faces.size();
00056 
00057    Edge * firstEdge = NULL;
00058    Edge * last = NULL;
00059    Edge * current = NULL;
00060 
00061    unsigned int i;
00062 
00063    for (i = 0; i < faceVerts.size()-1; i++) {
00064       current = addEdge(faceVerts[i], faceVerts[i+1]);
00065    
00066       check_error (!current, "Problem while parsing mesh faces.");
00067 
00068       current->face = f;
00069 
00070       if (last)
00071           last->next = current;
00072       else
00073           firstEdge = current;
00074       last = current;
00075    }
00076 
00077    current = addEdge (faceVerts[i], faceVerts[0]);
00078    check_error (!current, "Problem while parsing mesh faces.");
00079    
00080    current->face = f;
00081 
00082    last->next = current;
00083    current->next = firstEdge;
00084 
00085    f->edge = firstEdge;
00086    faces.push_back(f);
00087       
00088    return f;
00089 }
00090 
00091 Edge * Mesh::addEdge (int i, int j) {
00092    Edge eTmp;
00093    eTmp.vertex = vertices[i];
00094 
00095    Edge eTmpPair;
00096    eTmpPair.vertex = vertices[j];
00097    eTmp.pair = &eTmpPair;
00098 
00099    Mesh::EdgeSet::iterator eIn = edges.find(&eTmp);
00100    Edge * e;
00101 
00102    if (eIn != edges.end()){
00103       e = *eIn;
00104       if (e->face != NULL)
00105           return NULL;
00106    }
00107    else {
00108       e = new Edge();
00109       Edge * pair = new Edge();
00110 
00111       e->vertex = vertices[i];
00112       pair->vertex = vertices[j];
00113 
00114       pair->pair = e;
00115       e->pair = pair;
00116 
00117       edges.insert(e);
00118       edges.insert(pair);
00119 
00120       pair->vertex->edge = pair;
00121    }   
00122    return e;
00123 }
00124 
00134 bool Mesh::cutAlongEdge(Edge * forward) {
00135    if (!forward) return false;
00136    if (forward->isBoundary()) return false;
00137 
00138    Edge * backward = forward->pair;
00139    
00140    Edge * backwardPair = new Edge();
00141    Edge * forwardPair = new Edge();
00142 
00143    backwardPair->pair = backward;
00144    forwardPair->pair = forward;
00145    backwardPair->assignData(*backward);
00146    forwardPair->assignData(*forward);
00147    
00148    Vertex * v = forward->vertex;
00149    Vertex * endV = backward->vertex;
00150 
00151    if (v->isBoundary()) {
00152       Vertex * newV = addVertex(v->p);
00153       newV->assignData(*v);
00154       Vertex::EdgeAroundIterator aroundV = v->iterator(forward);
00155       do {
00156           aroundV++;
00157           aroundV.vertex() = newV;
00158       }
00159       while (aroundV.edge_out()->pair->face);
00160 
00161       forwardPair->next = aroundV.edge_out()->pair->next;
00162       aroundV.edge_out()->pair->next = backwardPair;
00163 
00164       backwardPair->vertex = newV;
00165       newV->edge = backwardPair;
00166    }
00167    else {
00168       forwardPair->next = backwardPair;
00169       backwardPair->vertex = v;
00170       v->edge = backwardPair;
00171    }
00172 
00173    if (endV->isBoundary()) {
00174       Vertex * newV = addVertex(endV->p);
00175       newV->assignData(*endV);
00176       Vertex::EdgeAroundIterator aroundV = endV->iterator(backward);
00177       do {
00178           aroundV++;
00179           aroundV.vertex() = newV;
00180       }
00181       while (aroundV.edge_out()->pair->face);
00182 
00183       backwardPair->next = aroundV.edge_out()->pair->next;
00184       aroundV.edge_out()->pair->next = forwardPair;
00185 
00186       forwardPair->vertex = newV;
00187       newV->edge = forwardPair;
00188    }
00189    else {
00190       backwardPair->next = forwardPair;
00191       forwardPair->vertex = endV;
00192       endV->edge = forwardPair;
00193    }
00194 
00195    forward->pair = forwardPair;
00196    backward->pair = backwardPair;
00197 
00198    return (addEdge(forwardPair) && addEdge(backwardPair));
00199 }
00200 
00202 void Mesh::computeLengths() {
00203    for (Mesh::EdgeIterator eIter = edgeIterator(); !eIter.end(); eIter++) {
00204       Edge * e = eIter.edge();
00205       e->length = (e->vertex->p - e->next->vertex->p).abs();
00206       e->pair->length = e->length;
00207    }
00208 }
00209 
00211 void Mesh::linkBoundary() {
00212    HalfEdgeIterator hIter = halfEdgeIterator();
00213 
00214    for (; !hIter.end(); hIter++) {
00215       if (!hIter.half_edge()->face) 
00216           hIter.half_edge()->vertex->edge = hIter.half_edge();
00217    }
00218 
00219    VertexIterator vIter = vertexIterator();
00220    Edge *next, *beg;
00221    while (!vIter.end()) {
00222       if (vIter.vertex()->isBoundary()) {
00223           beg = vIter.vertex()->edge;
00224           next = beg;
00225           while (beg->pair->face)
00226              beg = beg->pair->next;
00227           beg = beg->pair;
00228           beg->next = next;
00229           numBoundaryVertices++;
00230       }
00231       vIter++;
00232    }
00233 
00234 }
00235 
00245 void Mesh::computeMeshInfo() {
00246    cout << "Topological information about the mesh:" << endl;
00247    // Number of boundary loops
00248    Mesh::HalfEdgeIterator hIter = halfEdgeIterator();
00249    for (; !hIter.end(); hIter++) {
00250       hIter.half_edge()->check = false;
00251    }
00252    numBoundaryLoops = 0;
00253    for (hIter.reset(); !hIter.end(); hIter++) {
00254       Edge * e = hIter.half_edge();
00255       if (e->face)
00256           e->check = true;
00257       else if (!e->check) {
00258           Edge * beg = NULL;
00259           while (e != beg) {
00260              if (!beg) beg = e;
00261              check_error(!e, "Building the mesh failed, probem with input format.");
00262              e->check = true;
00263              e = e->next;
00264           }
00265           numBoundaryLoops++;
00266       }
00267    }
00268    cout << "Mesh has " << numBoundaryLoops << " boundary loops." << endl;
00269    // Number of connected components
00270    numConnectedComponents = 0;
00271    Mesh::FaceIterator fIter = faceIterator();
00272    for (; !fIter.end(); fIter++) {
00273       fIter.face()->check = false;
00274    }
00275    stack<Edge *> toVisit;
00276    for (fIter.reset(); !fIter.end(); fIter++) {
00277       if (!fIter.face()->check) {
00278           numConnectedComponents++;
00279           toVisit.push(fIter.face()->edge);
00280           while (!toVisit.empty()) {
00281              Face * fIn = toVisit.top()->face; 
00282              toVisit.pop();
00283              fIn->check = true;     
00284              Face::EdgeAroundIterator iter = fIn->iterator();
00285              for (; !iter.end(); iter++) 
00286                 if (iter.edge()->pair->face && !iter.edge()->pair->face->check)
00287                     toVisit.push(iter.edge()->pair);
00288           }
00289       }
00290    }
00291    cout << "Mesh has " << numConnectedComponents << " connected components." << endl;
00292    // Genus number
00293    check_error(numConnectedComponents == 0, "The mesh read is empty.");
00294    numGenus = 
00295       (1 - (numVertices() - numEdges() + numFaces() + numBoundaryLoops ) / 2) / numConnectedComponents;
00296    cout << "Mesh is genus " << numGenus << "." << endl;
00297 }
00298 
00300 bool Mesh::checkVertexConection() {
00301    Mesh::FaceIterator fIter = faceIterator();
00302    Mesh::VertexIterator vIter = vertexIterator();
00303    bool conectedVert = true;
00304 
00305    for (;!vIter.end(); vIter++)
00306       vIter.vertex()->check = false;
00307 
00308    for (fIter.reset(); !fIter.end(); fIter++) {
00309       Face::EdgeAroundIterator around = fIter.face()->iterator();
00310       for (;!around.end();around++)
00311           around.vertex()->check = true;
00312    }
00313    for (vIter.reset(); !vIter.end(); vIter++) {
00314       if (!vIter.vertex()->check) {
00315           cerr << "Vertex " << vIter.vertex()->ID << " is not connected." << endl;
00316           conectedVert = false;
00317       }
00318    }
00319 
00320    return conectedVert;
00321 }
00322 
00324 bool Mesh::checkManifold() {
00325    Mesh::HalfEdgeIterator eIter = halfEdgeIterator();
00326    Mesh::VertexIterator vIter = vertexIterator();
00327    bool manifold = true;
00328 
00329    for (;!eIter.end(); eIter++)
00330       eIter.half_edge()->check = false;
00331 
00332    for (vIter.reset(); !vIter.end(); vIter++) {
00333       Vertex::EdgeAroundIterator around = vIter.vertex()->iterator();
00334       for (;!around.end();around++)
00335           around.edge_out()->check = true;
00336    }
00337 
00338    for (eIter.reset(); !eIter.end(); eIter++) {
00339       if (!eIter.half_edge()->check) {
00340           cerr << "Mesh is not manifold - more then one disk at vertex " 
00341              << eIter.half_edge()->vertex->ID << endl;
00342           manifold = false;
00343           break;
00344       }
00345    }
00346 
00347    return manifold;
00348 }
00349 
00356 void Mesh::computeInitAngles() {
00357    Mesh::HalfEdgeIterator eIter = halfEdgeIterator();
00358    for (; !eIter.end(); eIter++) {
00359       Edge * e = eIter.half_edge();
00360       if (e->face) {
00361           double l1 = e->length;
00362           double l2 = e->next->length;
00363           double l3 = e->next->next->length;
00364           e->alphaOposite = acos((l2*l2 + l3*l3 - l1*l1)/(2.0*l2*l3));
00365       }
00366       else {
00367           e->alphaOposite = 0;
00368       }
00369    }
00370    for (eIter.reset(); !eIter.end(); eIter++) {
00371       Edge * e = eIter.half_edge();
00372       e->setTheta(M_PI - e->alphaOposite - e->pair->alphaOposite);
00373    }
00374 }
00375 
00390 void Mesh::readOBJ(const char * obj_file) {
00391    string front;
00392    string v = "v", vt = "vt", f = "f";
00393    Vector3 vert;
00394    vector<int> verts;
00395    vector<Vector3> uvVec;
00396    vector<int> uvs;
00397    char etc;
00398    int id;
00399 
00400    ifstream in(obj_file);
00401 
00402    check_error(!in, "Cannot find input obj file.");
00403 
00404    bool hasUV = false;
00405 
00406    while(!in.eof() || in.peek() != EOF) {
00407       in >> ws;
00408       if (in.eof() || in.peek() == EOF)
00409           break;
00410       if (in.peek() == '#') {
00411           in.ignore(300, '\n');
00412       }
00413       else {
00414           in >> front;
00415           if (front == v) {
00416              in >> vert.x() >> vert.y() >> vert.z();
00417              addVertex(vert);
00418           }
00419           else if (front == vt){
00420              in >> vert.x() >> vert.y();
00421              vert.z() = 0;
00422              uvVec.push_back(vert);
00423              hasUV = true;
00424           }
00425           else if (front == f) {
00426              verts.clear();
00427              uvs.clear();
00428              while (in >> id) {
00429 
00430                 check_error(id > numVertices(), "Problem with input OBJ file.");
00431 
00432                 verts.push_back(id-1);
00433                 bool vtNow = true;
00434                 if (in.peek() == '/'){
00435                     in >> etc;
00436                     if (in.peek() != '/') {
00437                        in >> id;
00438                        check_warn(id > numVertices(), "Texture coordinate index is greater then number of vertices.");
00439                        if (id < numVertices() && hasUV) {
00440                            uvs.push_back(id-1);
00441                            vtNow = false;
00442                        }
00443                     }
00444                 }
00445                 if (in.peek() == '/'){
00446                     int tmp;
00447                     in >> etc;
00448                     in >> tmp;
00449                 }
00450                 if (hasUV && vtNow) {
00451                     uvs.push_back(id-1);
00452                 }
00453              }
00454              in.clear(in.rdstate() & ~ios::failbit);
00455              Face * f = addFace(verts);
00456 
00457              if (hasUV && uvs.size() != 0){
00458                 int k = 0;
00459                 for (Face::EdgeAroundIterator e = f->iterator(); !e.end(); e++, k++)
00460                     e.vertex()->uv = uvVec[uvs[k]];
00461              }
00462           }
00463           else {
00464              string line;
00465              getline(in, line);
00466              cout << "Warning, line: " << line << " ignored." << endl;  
00467           }
00468       }
00469    }
00470 
00471    in.close();
00472 
00473    // Finnish building the mesh, should be called after each parse.
00474    finishMesh();
00475 }
00476 
00477 /* Write mesh in OBJ format to obj_file parameter */
00478 void Mesh::writeOBJ(const char * obj_file) {
00479    ofstream out(obj_file);
00480    check_error (!out, "Cannot find output file.");
00481 
00482    Mesh::VertexIterator vItr = vertexIterator();
00483    for (vItr.reset(); !vItr.end(); vItr++)
00484       out << "v " << vItr.vertex()->p << endl;
00485 
00486    for (vItr.reset(); !vItr.end(); vItr++)
00487       out << "vt " << vItr.vertex()->uv.x() << " " << vItr.vertex()->uv.y() << endl;
00488 
00489    Mesh::FaceIterator fItr = faceIterator();
00490    for (fItr.reset(); !fItr.end(); fItr++) {
00491       Face::EdgeAroundIterator around = fItr.face()->iterator();
00492       out << "f";
00493       for ( ; !around.end(); around++)
00494           out << " " << (around.vertex()->ID + 1) << "/" << (around.vertex()->ID + 1);
00495       out << endl;
00496    }
00497 }
00498 
00499 /* Write mesh in VT format (only text coodrinates) to vt_file parameter */
00500 void Mesh::writeVT(const char * vt_file) {
00501    ofstream out(vt_file);
00502    check_error (!out, "Cannot find output file.");
00503 
00504    Mesh::VertexIterator vItr = vertexIterator();
00505    for (vItr.reset(); !vItr.end(); vItr++)
00506       out << "vt " << vItr.vertex()->uv.x() << " " << vItr.vertex()->uv.y() << endl;
00507 }
00508 
00528 void Mesh::readCON(const char * conn_file) {
00529    string front;
00530    string v = "v", f = "f";
00531    Vector3 vert;
00532    vector<int> verts;
00533    int id;
00534 
00535    ifstream in(conn_file);
00536 
00537    check_error(!in, "Cannot find input obj file.");
00538 
00539    int NumFaces, NumVertices;
00540 
00541    in >> NumVertices;
00542    in >> NumFaces;
00543 
00544    for (int i = 0; i < NumVertices; i++) {
00545       addVertex(vert);
00546    }
00547    while(!in.eof() || in.peek() != EOF) {
00548       in >> ws;
00549       if (in.eof() || in.peek() == EOF)
00550           break;
00551       if (in.peek() == '#') {
00552           in.ignore(300, '\n');
00553       }
00554       else {
00555           in >> front;
00556           if (front == f) {
00557              verts.clear();
00558              
00559              for (int i = 0; i < 3; i++) {
00560                 in >> id;
00561                 check_error(id > numVertices(), "Problem with input CON file.");
00562                 verts.push_back(id-1);
00563              }
00564              in.clear(in.rdstate() & ~ios::failbit);
00565              Face * f = addFace(verts);
00566 
00567              for (Face::EdgeAroundIterator e = f->iterator(); !e.end(); e++) {
00568                 in >> e.edge()->length;
00569              }
00570           }
00571           else {
00572              string line;
00573              getline(in, line);
00574              cout << "Warning, line: " << line << " ignored." << endl;  
00575           }
00576       }
00577    }
00578 
00579    in.close();
00580    // Finnish building the mesh, should be called after each parse.
00581    finishMesh();
00582 }
00583 
00588 void Mesh::writeCON(const char * conn_file) {
00589    ofstream out(conn_file);
00590    check_error (!out, "Cannot find output file.");
00591 
00592    out << numVertices() << endl;
00593    out << numFaces() << endl;
00594 
00595    Mesh::FaceIterator fItr = faceIterator();
00596    for (fItr.reset(); !fItr.end(); fItr++) {
00597       Face::EdgeAroundIterator around = fItr.face()->iterator();
00598       out << "f";
00599       for ( ; !around.end(); around++)
00600           out << " " << (around.vertex()->ID + 1);
00601       out << endl;
00602 
00603       for (around.reset() ; !around.end(); around++) 
00604           out << " " << setprecision(8) << around.edge()->length;
00605       out << endl;
00606    }
00607 
00608    out.close();
00609 }
00610 
00624 void Mesh::readCUT_EDGES(const char * edge_file) {
00625    ifstream in(edge_file);
00626    check_error(!in, "Cannot find input file for edge cuts.");
00627 
00628    vector<Edge *> edges_to_cut;
00629 
00630    int i, j, f;
00631    while (in >> i) {
00632       in >> j >> f;
00633 
00634       if ((i < 0) || (i >= numVertices()) || (j < 0) || (j >= numVertices()))
00635           cout << "ID of a vertex in the edges to cut file is out of bounds." << endl;
00636       else {
00637 
00638           i--; j--; f--;
00639    
00640           Edge eTmp;
00641          eTmp.vertex = vertices[i];          
00642           Edge eTmpPair;
00643           eTmpPair.vertex = vertices[j];
00644          eTmp.pair = &eTmpPair;
00645          eTmpPair.pair = &eTmp;
00646 
00647          Edge * e = NULL;
00648           Mesh::EdgeSet::iterator eIn = edges.find(&eTmp);
00649           if (eIn != edges.end())
00650              e = *eIn;
00651           while (eIn != edges.end()) {
00652              if ((*eIn)->face && (*eIn)->face->ID == f) {
00653                 e = *eIn;
00654                 break;
00655              }
00656              eIn++;
00657           }
00658           if (e)
00659              edges_to_cut.push_back(e);
00660       }
00661    }
00662 
00663    while (!edges_to_cut.empty()) {
00664       Edge * e = edges_to_cut.back();
00665       if (e) {
00666              cutAlongEdge(e);
00667       }
00668       edges_to_cut.pop_back();
00669    }
00670 
00671    EdgeSet newedges(edges.begin(), edges.end());
00672    edges = newedges;
00673 }
00674 
00684 void Mesh::readCONE_VERT(const char * vert_file) {
00685    ifstream in(vert_file);
00686    check_error(!in, "Cannot find input file for cone singularities.");
00687 
00688    int id = 0;
00689    double minC = 0, maxC = 0;
00690 
00691    while (in >> id) {
00692       in >> minC >> maxC;
00693       if ((id < 0) || (id >= numVertices()))
00694           cout << "ID of a vertex in the cone singularities file is out of bounds." << endl;
00695       else {
00696           id--;
00697          double valenceAlow = (double)vertices[id]->getValence();
00698           check_warn ( maxC <=0, "Invalid cone angle. Cone angle read is less or equal to zero.");
00699           check_warn ( minC >= valenceAlow, 
00700              "Invalid cone angle. Cone angle read is greater or equal to what valence of vertex allows.");
00701           vertices[id]->min_cone_angle = minC;
00702           vertices[id]->max_cone_angle = maxC;
00703           vertices[id]->constrained = true;
00704       }
00705    }
00706 }
00707 
00708 
00711 void Mesh::checkGaussianCurvature() {
00712    double minConeGC = 0; 
00713    double maxConeGC = 0;
00714   
00715    for (int i = 0; i < numVertices(); i++)   {
00716       if (vertices[i]->constrained) {
00717           if (vertices[i]->isBoundary()) {
00718              maxConeGC += (M_PI - vertices[i]->min_cone_angle*M_PI);
00719              minConeGC += (M_PI - vertices[i]->max_cone_angle*M_PI);
00720           }
00721           else {
00722              maxConeGC += (2*M_PI - vertices[i]->min_cone_angle*M_PI);
00723              minConeGC += (2*M_PI - vertices[i]->max_cone_angle*M_PI);
00724           }
00725       }
00726       else if (vertices[i]->isBoundary()) {
00727           maxConeGC +=  M_PI;
00728           minConeGC += -M_PI;;
00729       }
00730    }
00731 
00732    double eNum1 = minConeGC / (2*M_PI);
00733    double eNum2 = maxConeGC / (2*M_PI);
00734    double eNumCurr = (double)(2 - 2*numGenus - numBoundaryLoops)/(double)numConnectedComponents;
00735    cout << "Checking Gauss-Bonet Theorem:" << endl;
00736    cout << "Minimum euler number for this set up is " << eNum1 << endl;
00737    cout << "Maximum euler number for this set up is " << eNum2 << endl;
00738    cout << "Euler number of the mesh is " << eNumCurr << "\n" << endl;
00739 
00740    check_warn((eNumCurr < eNum1) || (eNumCurr > eNum2), 
00741       "The mesh cannot be developed with current set of cone singularities.\n");
00742 }
00743 

Generated on Sat Jun 3 13:33:41 2006 for CirclePatterns by  doxygen 1.4.5