MeshKit  1.0
Tri2Quad.cpp
Go to the documentation of this file.
00001 #include "meshkit/Tri2Quad.hpp"
00002 #include "meshkit/StopWatch.hpp"
00003 
00004 using namespace Jaal;
00005 
00007 
00008 int Tri2Quads::verify(Mesh *mesh, const vector<FacePair> &matching)
00009 {
00010      int relexist2 = mesh->build_relations(0, 2);
00011 
00012      // This is the Graph matching verification. It is useful to verify the
00013      // matching, if done by other algorithms.
00014 
00015      size_t numfaces = mesh->getSize(2);
00016      for (size_t i = 0; i < numfaces; i++) {
00017           Face * f = mesh->getFaceAt(i);
00018           f->setVisitMark(0);
00019      }
00020 
00021      FaceSequence faceneighs;
00022      size_t nsize = matching.size();
00023      for (size_t i = 0; i < nsize; i++) {
00024           size_t f1 = matching[i].first;
00025           size_t f2 = matching[i].second;
00026           assert(f1 != f2);
00027           if (f1 < f2) {
00028                Face *f0 = mesh->getFaceAt(f1);
00029                Face *f1 = mesh->getFaceAt(f2);
00030                f0->getRelations12( faceneighs );
00031                if (find(faceneighs.begin(), faceneighs.end(), f1) == faceneighs.end()) {
00032                     cout << "Error: Face matching error: faces not neighbors  " << endl;
00033                     exit(0);
00034                     return 1;
00035                }
00036                assert(!f0->isVisited());
00037                assert(!f1->isVisited());
00038                f0->setVisitMark(1);
00039                f1->setVisitMark(1);
00040           }
00041      }
00042 
00043      for (size_t i = 0; i < numfaces; i++) {
00044           Face * f = mesh->getFaceAt(i);
00045           assert(f->isVisited());
00046      }
00047 
00048      if( !relexist2 ) mesh->clear_relations(0,2 );
00049 
00050      return 0;
00051 
00052 }
00054 
00055 Mesh* Tri2Quads::collapse_matched_triangles(Mesh *trimesh, const vector<FacePair> &matching,
00056           int replace)
00057 {
00058 #ifdef DEBUG
00059      if(Tri2Quads::verify(trimesh, matching) != 0) return NULL;
00060      cout << " Verification Done " << endl;
00061 #endif
00062 
00063      Mesh *quadmesh = new Mesh;
00064 
00065      size_t nsize = matching.size();
00066 
00067      assert( nsize );
00068 
00069      Face *tri1, *tri2, *quad;
00070 
00071      NodeSequence  tnodes;
00072      trimesh->getNodes( tnodes );
00073      quadmesh->addNodes(tnodes );
00074      quadmesh->reserve( nsize, 2 );
00075 
00076      for (size_t i = 0; i < nsize; i++) {
00077           size_t f1 = matching[i].first;
00078           size_t f2 = matching[i].second;
00079           tri1 = trimesh->getFaceAt(f1);
00080           tri2 = trimesh->getFaceAt(f2);
00081           quad = Face::create_quad(tri1, tri2, replace);
00082           quadmesh->addFace(quad);
00083           if (replace) delete tri2;
00084      }
00085 
00086      if (replace) trimesh->emptyAll();
00087 
00088      return quadmesh;
00089 }
00091 
00092 int Tri2Quads::refine_boundary_triangle(Face *tri0)
00093 {
00094      if (tri0->getSize(0) != 3)
00095           return 1;
00096 
00097      Vertex *bv0 = NULL;
00098      Vertex *bv1 = NULL;
00099      Vertex *bv2 = NULL;
00100 
00101      for (int i = 0; i < 3; i++) {
00102           Vertex *ev1 = tri0->getNodeAt(i + 1);
00103           Vertex *ev2 = tri0->getNodeAt(i + 2);
00104           if (ev1->isBoundary() && ev2->isBoundary()) {
00105                bv0 = ev1;
00106                bv1 = ev2;
00107                bv2 = tri0->getNodeAt(i);
00108                break;
00109           }
00110      }
00111 
00112      if (bv0 == NULL || bv1 == NULL)
00113           return 2;
00114 
00115      Point3D p3d;
00116      Vertex::mid_point(bv0, bv1, p3d);
00117 
00118      Vertex *bound = Vertex::newObject();
00119      bound->setXYZCoords(p3d);
00120 
00121      trimesh->addNode(bound);
00122 
00123      NodeSequence tconnect(3);
00124      tconnect[0] = bv0;
00125      tconnect[1] = bound;
00126      tconnect[2] = bv2;
00127      tri0->setNodes(tconnect);
00128 
00129      tconnect[0] = bound;
00130      tconnect[1] = bv1;
00131      tconnect[2] = bv2;
00132      maxfaceID++;
00133      Face *tri1 = Face::newObject();
00134      tri1->setID(maxfaceID);
00135      tri1->setNodes(tconnect);
00136      trimesh->addFace(tri1);
00137 
00138      steinerNodes.push_back(bound);
00139 
00140      FacePair facepair;
00141      facepair.first = tri0->getID();
00142      facepair.second = tri1->getID();
00143      facematching.push_back(facepair);
00144 
00145      return 0;
00146 }
00147 
00149 
00150 
00151 void Tri2Quads::splitParent(BinaryNode *parent, BinaryNode *child1,
00152                             BinaryNode *child2)
00153 {
00154      Vertex* dnode = NULL;
00155      dnode = parent->getDualNode();
00156 
00157      Face *parentface = NULL;
00158      dnode->getAttribute("PrimalFace", parentface);
00159 
00160      dnode = child1->getDualNode();
00161      Face *face1 = NULL;
00162      dnode->getAttribute("PrimalFace", face1);
00163 
00164      dnode = child2->getDualNode();
00165      Face *face2 = NULL;
00166      dnode->getAttribute("PrimalFace", face2);
00167 
00168      NodeSequence connect(3);
00169 
00170      // Remove all existing vertex-face relations;
00171      Vertex *vertex;
00172      for (int i = 0; i < 3; i++) {
00173           vertex = parentface->getNodeAt(i);
00174           vertex->clearRelations(2);
00175 
00176           vertex = face1->getNodeAt(i);
00177           vertex->clearRelations(2);
00178 
00179           vertex = face2->getNodeAt(i);
00180           vertex->clearRelations(2);
00181      }
00182 
00183      // Rebuild vertex-face relations...
00184      for (int i = 0; i < 3; i++) {
00185           vertex = parentface->getNodeAt(i);
00186           vertex->addRelation(parentface);
00187 
00188           vertex = face1->getNodeAt(i);
00189           vertex->addRelation(face1);
00190 
00191           vertex = face2->getNodeAt(i);
00192           vertex->addRelation(face2);
00193      }
00194 
00195      dnode = NULL;
00196      parentface->getAttribute("DualNode", dnode);
00197      Vertex *steiner = dnode->getClone();
00198      steiner->setID(parentface->getID());
00199      trimesh->addNode(steiner);
00200      steinerNodes.push_back(steiner);
00201 
00202      int edge1, edge2, edge3;
00203 
00204      edge1 = edge2 = edge3 = -1;
00205      FaceSequence neighs;
00206      for (int i = 0; i < 3; i++) {
00207           Vertex *v0 = parentface->getNodeAt(i + 1);
00208           Vertex *v1 = parentface->getNodeAt(i + 2);
00209           Mesh::getRelations112(v0, v1, neighs);
00210 
00211           if (neighs.size() == 1)
00212                edge3 = i;
00213 
00214           if (neighs.size() == 2) {
00215                if (find(neighs.begin(), neighs.end(), face1) != neighs.end())
00216                     edge1 = i;
00217                if (find(neighs.begin(), neighs.end(), face2) != neighs.end())
00218                     edge2 = i;
00219           }
00220      }
00221 
00222      Face *qface;
00223      Vertex *ev0, *ev1;
00224 
00225      // Match Child1 and One of the Split Triangle ...
00226      maxfaceID++;
00227      ev0 = parentface->getNodeAt(edge1 + 1);
00228      ev1 = parentface->getNodeAt(edge1 + 2);
00229      connect[0] = steiner;
00230      connect[1] = ev0;
00231      connect[2] = ev1;
00232      qface = Face::newObject();
00233      qface->setNodes(connect);
00234      qface->setID(maxfaceID);
00235      Vertex *dc1 = DualGraph::getNewDualNode( qface );
00236      dc1->setID(maxfaceID);
00237      trimesh->addFace(qface);
00238      steinerFaces.push_back(qface);
00239      dnode = NULL;
00240      face1->getAttribute("DualNode", dnode);
00241      matchnodes( dnode, dc1);
00242      BinaryNode *bnode1 = new BinaryNode(dc1);
00243      bnode1->setMatchMark(1);
00244      bnode1->setParent(parent);
00245      bnode1->addChild(child1);
00246      btree->addNode(bnode1);
00247 
00248      // Match Child2 and One of the Split Triangle ...
00249      maxfaceID++;
00250      ev0 = parentface->getNodeAt(edge2 + 1);
00251      ev1 = parentface->getNodeAt(edge2 + 2);
00252      connect[0] = steiner;
00253      connect[1] = ev0;
00254      connect[2] = ev1;
00255      qface = Face::newObject();
00256      qface->setID(maxfaceID);
00257      qface->setNodes(connect);
00258      Vertex *dc2 = DualGraph::getNewDualNode( qface );
00259      dc2->setID(maxfaceID);
00260      trimesh->addFace(qface);
00261      steinerFaces.push_back(qface);
00262      dnode = NULL;
00263      face2->getAttribute( "DualNode", dnode);
00264      matchnodes( dnode, dc2);
00265 
00266      BinaryNode *bnode2 = new BinaryNode(dc2);
00267      bnode2->setMatchMark(1);
00268      bnode2->setParent(parent);
00269      bnode2->addChild(child2);
00270      btree->addNode(bnode2);
00271 
00272      // Now Parent have different connectivity ...
00273      ev0 = parentface->getNodeAt(edge3 + 1);
00274      ev1 = parentface->getNodeAt(edge3 + 2);
00275      connect[0] = steiner;
00276      connect[1] = ev0;
00277      connect[2] = ev1;
00278      parentface->setNodes(connect);
00279      Point3D p3d;
00280      parentface->getAvgPos( p3d );
00281 
00282      Vertex *dc3 = NULL;
00283      parentface->getAttribute("DualNode", dc3);
00284      dc3->setXYZCoords(p3d);
00285      parent->addChild(bnode1);
00286      parent->addChild(bnode2);
00287      modifiedFaces.push_back(parentface);
00288 
00289      for (int i = 0; i < 3; i++) {
00290           vertex = parentface->getNodeAt(i);
00291           vertex->clearRelations(2);
00292 
00293           vertex = face1->getNodeAt(i);
00294           vertex->clearRelations(2);
00295 
00296           vertex = face2->getNodeAt(i);
00297           vertex->clearRelations(2);
00298      }
00299 
00300      child1->setMatchMark(1);
00301      child2->setMatchMark(1);
00302 
00303      btree->removeNode(child1);
00304      btree->removeNode(child2);
00305 }
00306 
00308 
00309 void Tri2Quads::matchnode(BinaryNode* v)
00310 {
00311      BinaryNode *parv = v->getParent();
00312 
00313      if (parv == NULL)
00314           return;
00315 
00316      int degree = parv->getDegree();
00317 
00318      if (parv->isRoot() && degree == 2) {
00319           BinaryNode *vsib = v->getSibling();
00320           splitParent(parv, v, vsib);
00321           return;
00322      }
00323 
00324      if (degree == 1 || degree == 2) {
00325           matchnodes(v, parv);
00326           return;
00327      }
00328 
00329      if (degree == 3) {
00330 
00331           if (required_topology == ALL_QUADS) {
00332                BinaryNode *vsib = v->getSibling();
00333                splitParent(parv, v, vsib);
00334                return;
00335           }
00336 
00337           BinaryNode *vsib = v->getSibling();
00338           Vertex *d0 = v->getDualNode();
00339           Vertex *d1 = vsib->getDualNode();
00340           if (d0->getNumRelations(0) < d1->getNumRelations(0)) {
00341                matchnodes(v, parv);
00342                btree->unlinkNode(vsib);
00343           } else {
00344                matchnodes(vsib, parv);
00345                btree->unlinkNode(v);
00346           }
00347      }
00348 }
00349 
00351 
00352 BinaryNode* Tri2Quads::getChildofDegreeNParent(BNodeList &levelnodes,
00353           int nd)
00354 {
00355      BinaryNode *currnode, *parent, *child;
00356 
00357      int ncount;
00358      BNodeList::const_iterator it;
00359 
00360      for (it = levelnodes.begin(); it != levelnodes.end(); ++it) {
00361           currnode = *it;
00362           parent = currnode->getParent();
00363           if (parent) {
00364                if (!parent->isMatched()) {
00365                     ncount = 0;
00366                     if (parent->getParent())
00367                          ncount = 1;
00368                     for (int i = 0; i < parent->getNumChildren(); i++) {
00369                          child = parent->getChild(i);
00370                          if (!child->isMatched())
00371                               ncount++;
00372                     }
00373                     if (ncount == nd)
00374                          return currnode;
00375                }
00376           }
00377      }
00378 
00379      return NULL;
00380 }
00381 
00383 
00384 BinaryNode *Tri2Quads::getNextNode(BNodeList &levelnodes)
00385 {
00386      BinaryNode *currnode = NULL;
00387 
00388      if (levelnodes.empty())
00389           return currnode;
00390 
00391      BNodeList::iterator it;
00392      for (it = levelnodes.begin(); it != levelnodes.end(); ++it) {
00393           currnode = *it;
00394           if (currnode->isMatched())
00395                btree->unlinkNode(currnode);
00396      }
00397 
00398      it = remove_if(levelnodes.begin(), levelnodes.end(), already_matched);
00399      levelnodes.erase(it, levelnodes.end());
00400 
00401      BinaryNode *child = NULL;
00402 
00403      // High Priority: parent having degree = 1;
00404      child = getChildofDegreeNParent(levelnodes, 1);
00405 
00406      if (!child)
00407           child = getChildofDegreeNParent(levelnodes, 2);
00408 
00409      // Low Priority: parent having degree = 3;
00410      if (!child)
00411           child = getChildofDegreeNParent(levelnodes, 3);
00412 
00413      return child;
00414 }
00415 
00417 
00418 void Tri2Quads::prunelevel(BNodeList &levelnodes)
00419 {
00420      while (1) {
00421           BinaryNode *currnode = getNextNode(levelnodes);
00422           if (currnode == NULL) break;
00423           matchnode(currnode);
00424      }
00425 }
00426 
00428 
00429 void Tri2Quads::percolateup()
00430 {
00431      steinerNodes.clear();
00432      steinerFaces.clear();
00433 
00434      int height = btree->getHeight();
00435      BNodeList levelnodes;
00436      BNodeList::const_iterator it;
00437 
00438      //Reset all the Matching marks to 0;
00439      for (int i = 0; i < height; i++) {
00440           levelnodes = btree->getLevelNodes(height - i - 1);
00441           BinaryNode *currnode;
00442           for (it = levelnodes.begin(); it != levelnodes.end(); ++it) {
00443                currnode = *it;
00444                currnode->setMatchMark(0);
00445           }
00446      }
00447 
00448      // Start Prunning the level. At most the root will be unmatched.
00449      for (int i = 0; i < height; i++) {
00450           levelnodes = btree->getLevelNodes(height - i - 1);
00451           prunelevel(levelnodes);
00452      }
00453 
00454      size_t numfaces = trimesh->getSize(2);
00455      facematching.reserve(numfaces);
00456 
00457      for (size_t i = 0; i < numfaces; i++) {
00458           Face *face = trimesh->getFaceAt(i);
00459           Vertex *u = NULL;
00460           face->getAttribute("DualNode", u);
00461            
00462           assert(u);
00463           Vertex *v = NULL;
00464           u->getAttribute("DualMate", v);
00465           if (v) {
00466                if (v->getID() > u->getID()) {
00467                     FacePair facepair;
00468                     facepair.first = u->getID();
00469                     facepair.second = v->getID();
00470                     facematching.push_back(facepair);
00471                }
00472           }
00473      }
00474 
00475      // If the root is unmatched, bring it down to a leaf and then split the
00476      // leaf. Do this step after the triangles have been matched.
00477      BinaryNode *root = btree->getRoot();
00478      if (!root->isMatched()) {
00479 #ifdef VERBOSE
00480           cout << "Warning: Boundary Triangle modified " << endl;
00481 #endif
00482           Vertex *dnode = root->getDualNode();
00483           Face *rootface = NULL;
00484           dnode->getAttribute("PrimalFace", rootface);
00485           refine_boundary_triangle(rootface);
00486      }
00487 
00488 }
00489 
00491 
00492 void Tri2Quads::maximum_tree_matching()
00493 {
00494      // In order to insert any steiner point on the boundary triangle (at the root)
00495      // We should know which triangles and nodes are on the boundary. Therefore,
00496      // call this function to set the boundary flags. Building the relationship
00497      // at this stage is good as even the DualGraph construction require it.
00498 
00499      trimesh->build_relations(0, 2);
00500      trimesh->search_boundary();
00501 
00502 #ifdef VERBOSE
00503      cout << " Creating Dual Graph ... " << endl;
00504 #endif
00505 
00506      DualGraph *dgraph = new DualGraph;
00507      dgraph->build(trimesh);
00508 
00509      if (verbose)
00510           cout << " Building Binary Tree of Dual Graph ... " << endl;
00511 
00512      btree = new BinaryTree(dgraph);
00513      btree->build();
00514      btree->saveAs("btree");
00515 
00516 #ifdef VERBOSE
00517      cout << " Tree Matching ... " << endl;
00518 #endif
00519 
00520      percolateup();
00521 
00522      // Percolateup will remove relations, so it is necessary to clear all the
00523      // relations.
00524      trimesh->clear_relations(0, 2);
00525 
00526      btree->deleteAll();
00527      dgraph->deleteAll();
00528 
00529      delete btree;
00530      delete dgraph;
00531 }
00532 
00534 
00535 Mesh* Tri2Quads::getQuadMesh(Mesh *inmesh, int replace, int topo)
00536 {
00537 
00538 #ifdef DEBUG
00539      if (inmesh->isHomogeneous() != 3) {
00540           cout << "Warning: Input mesh is not triangular " << endl;
00541           return NULL;
00542      }
00543 
00544      if (!inmesh->isSimple()) {
00545           cout << "Warning: Input mesh is not simple, use edmonds' algorithm " << endl;
00546           return NULL;
00547      }
00548 
00549      if (inmesh->getNumOfComponents() > 1) {
00550           cout << "Warning: There are multiple components in the mesh" << endl;
00551           cout << "         Algorithm works for single component " << endl;
00552           return NULL;
00553      }
00554 
00555      int euler0 = trimesh->getEulerCharacteristic();
00556      cout << " Input Euler # : " << euler0 << endl;
00557      cout << inmesh->saveAs( "Check.dat");
00558 #endif
00559 
00560      trimesh = inmesh;
00561 
00562      required_topology = topo;
00563 
00564      trimesh->enumerate(2);
00565      maxfaceID = trimesh->getSize(2) - 1;
00566 
00568      // Generate Maximum Matching on a binary tree using Suneeta's Algorithm.
00569      // If the required topology is set to ALL_QUADS, steiner points( and new
00570      // faces) will be inserted in the input triangle mesh.  Please note that
00571      // this implementation doesn't produces "The" optimal soluation as
00572      // described in the original papers, and doesn't even guarantee that
00573      // the resulting quadrilaterals will be convex. This along with other
00574      // topological and geometric optimization are anyhow essential and
00575      // are carried out during the "Post Processing" step. Therefore, we
00576      // have sacrifised performance over quality in this implementation.
00577      // Roughly we can expect that about 4-5% steiner points are inserted in
00578      // most of the general cases.
00580      StopWatch swatch;
00581      swatch.start();
00582 
00583      maximum_tree_matching();
00584      Mesh *quadmesh = collapse_matched_triangles(trimesh, facematching, replace);
00585      swatch.stop();
00586      cout << "Info: Tri->Quad Elapsed Time " << swatch.getSeconds() << endl;
00587 
00588      if( quadmesh ) {
00589           if (!quadmesh->isSimple())
00590                cout << "Warning: Quadrilateral Mesh is not simple " << endl;
00591 
00592           quadmesh->enumerate(2);
00593           if (!quadmesh->is_consistently_oriented()) {
00594                cout << "Warning:Trying to make Quadrilateal Mesh consistently oriented " << endl;
00595                quadmesh->make_consistently_oriented();
00596                if (quadmesh->is_consistently_oriented())
00597                     cout << "Info: Quadrilateral Mesh is now consistently oriented: Very good " << endl;
00598                else
00599                     cout << "Alas ! Quadrilateral Mesh is still inconsistently oriented: Check manually " << endl;
00600           }
00601 
00602           quadmesh->enumerate(0);
00603           quadmesh->enumerate(2);
00604      }
00605 
00607      // Since Steiner points may be inserted in the mesh ( and new triangles).
00608      // Renumber all the nodes and faces for future processing.
00610 
00611      trimesh->enumerate(0);
00612      trimesh->enumerate(2);
00613 
00614      return quadmesh;
00615 }
00616 
00618 
00619 void Tri2Quads::match_tree_walk(BinaryTree *btree, BinaryNode *parent)
00620 {
00621      //
00622      // Brings all the internal unmatched nodes at the leaf.
00623      //
00624      if (parent == NULL)
00625           return;
00626      if (parent->isLeaf())
00627           return;
00628 
00629      int numChildren = parent->getNumChildren();
00630 
00631      for (int i = 0; i < numChildren; i++) {
00632           BinaryNode *child1 = parent->getChild(i);
00633           if (!btree->isMatched(parent, child1)) {
00634                int numGrandChildren = child1->getNumChildren();
00635                for (int j = 0; j < numGrandChildren; j++) {
00636                     BinaryNode *child2 = child1->getChild(j);
00637                     if (btree->isMatched(child1, child2)) {
00638                          Vertex *np = parent->getDualNode();
00639                          assert(np);
00640                          Vertex *c1 = child1->getDualNode();
00641                          assert(c1);
00642                          Vertex *c2 = child2->getDualNode();
00643                          assert(c2);
00644                          matchnodes(np, c1);
00645                          c2->setAttribute("DualMate", 0);
00646                          c2->setStatus(MeshEntity::ACTIVE);
00647                          match_tree_walk(btree, child2);
00648                          return;
00649                     }
00650                }
00651           }
00652      }
00653 
00654 }
00655 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines