MeshKit  1.0
QuadCleanUp.hpp
Go to the documentation of this file.
00001 #ifndef QUADCLEAN_H
00002 #define QUADCLEAN_H
00003 
00005 //                      Quad-Cleanup
00006 //
00007 // Developed by:  Chaman Singh Verma
00008 //                Department of Computer Sciences.
00009 //                The University of Wisconsin, Madison
00010 //
00011 // Work Supported by:
00012 //                 Dr. Tim Tautges
00013 //                 Argonne National Lab, Chicago
00014 //
00015 //
00016 // Objective:  Given a quadrilateral mesh, this class implements various strategies
00017 // to improve the quadrilateral mesh both geometrically and topologically. The
00018 // Laplacian ( local and global ) is used for geometric quality improvement, and for
00019 // topological improvements various operations are used.
00020 // The two basis operations for topological improvements are
00021 //  1)   Face close
00022 //  2)   doublet insertion and removal.
00023 //
00024 
00025 // Reference Papers:
00026 //  1) Topological Improvement Procedures for Quadrilateral Finite Element Meshes
00027 //     S.A. Canann,  S.N. Muthikrishnan and R.K. Phillips
00028 
00029 //  2) Automated All Quadrilateral Mesh Adaptation Through Refinment and Coarsening
00030 //     Bret Dallas Anderson
00031 //     Master Thesis, Brigham Young University.
00032 //
00033 //  3) Non-Local Topological Clean-Up ( The idea of yring  is from this paper)
00034 //     Guy Bunin.
00035 //
00036 // For suggestios, bugs and criticisms, please send e-mail to
00037 //                      csverma@cs.wisc.edu
00038 //
00039 // Last Date update:  16th Feb 2010.
00040 //
00042 
00043 #include "Mesh.hpp"
00044 #include "Tri2Quad.hpp"
00045 #include "basic_math.hpp"
00046 #include "StopWatch.hpp"
00047 #include "tfiblend.hpp"
00048 
00049 #include "DijkstraShortestPath.hpp"
00050 #include "QuadPatchRemesh.hpp"
00051 
00052 extern double area_of_poly3d(int n, double *x, double *y, double *z);
00053 
00054 namespace Jaal {
00055 
00056 struct FirstIrregularNode : public MeshFilter {
00057      bool pass( const Vertex *vertex ) const {
00058           if( vertex->isBoundary() ) return 1;
00059           if( vertex->getNumRelations(2) != 4 ) return 0;
00060           return 1;
00061      }
00062 };
00063 
00065 // Diamond:  An element whose at least one of the opposite vertex is surrounded by
00066 //           three faces. In many cases, diamonds are essential in the quadrilateral
00067 // mesh and they can not be removed, Finding the minimum number of diamonds is hard,
00068 // and we are working towards it.
00070 
00071 class FaceClose {
00072 public:
00073 
00074      FaceClose(Mesh *m, Face *f, Vertex *v0, Vertex *v1) {
00075           mesh = m;
00076           face = f;
00077           vertex0 = v0;
00078           vertex2 = v1;
00079           replacedNode = NULL;
00080      }
00081 
00082      ~FaceClose() {
00083           if( replacedNode ) delete replacedNode;
00084      }
00085 
00086      int remove();
00087 
00088      int build();
00089      int commit();
00090 
00091      Mesh *mesh;
00092      Face *face;
00093      Vertex *vertex0, *vertex2;
00094      Vertex *replacedNode;
00095 private:
00096      bool isSafe() const;
00097      int  backup();
00098      int  rollback();
00099 };
00100 
00102 
00103 struct Diamond {
00104 
00105      Diamond(Mesh *m, Face *f, int p) {
00106           mesh = m;
00107           face = f;
00108           position = p;
00109           faceclose = NULL;
00110           vertex0   = NULL;
00111           vertex2   = NULL;
00112 
00113           if (position == 0 || position == 2) {
00114                vertex0 = face->getNodeAt(0);
00115                vertex2 = face->getNodeAt(2);
00116           }
00117 
00118           if (position == 1 || position == 3) {
00119                vertex0 = face->getNodeAt(1);
00120                vertex2 = face->getNodeAt(3);
00121           }
00122      }
00123      ~Diamond() {
00124           if( faceclose) delete faceclose;
00125      }
00126 
00127      int remove();
00128      int commit();
00129      int isSafe();
00130      int makeShield();
00131 
00132      bool operator<(const Diamond & rhs) const {
00133           return face->getArea() < rhs.face->getArea();
00134      }
00135 
00136      Vertex * getNewNode() const {
00137           if (faceclose) return faceclose->replacedNode;
00138           return NULL;
00139      }
00140 
00141      double getDiagonalRatio() const {
00142           Vertex *v0 = face->getNodeAt((position + 0) % 4);
00143           Vertex *v1 = face->getNodeAt((position + 1) % 4);
00144           Vertex *v2 = face->getNodeAt((position + 2) % 4);
00145           Vertex *v3 = face->getNodeAt((position + 3) % 4);
00146           double len0 = Vertex::length(v0, v2);
00147           double len1 = Vertex::length(v1, v3);
00148           return len0 / len1;
00149      }
00150      int build();
00151 
00152      Face *face;
00153 private:
00154      Vertex *vertex0, *vertex2;
00155      Mesh *mesh;
00156      int position;
00157      FaceClose *faceclose;
00158 };
00159 
00161 
00162 struct Doublet {
00163 
00164      Doublet(Mesh *m, Vertex * v) {
00165           mesh = m;
00166           vertex = v;
00167           replacedFace = NULL;
00168      }
00169 
00170      bool isSafe() const;
00171      void makeShield();
00172      int remove();
00173 
00174      Mesh *mesh;
00175      Vertex *vertex;
00176      Face *replacedFace;
00177      Face * shield[2];
00178 };
00179 
00181 
00182 struct Singlet {
00183      Singlet(Mesh *m, Vertex * v) {
00184           mesh = m;
00185           vertex = v;
00186           active = 1;
00187      }
00188 
00189      int remove();
00190 
00191 private:
00192      Mesh *mesh;
00193      Vertex *vertex;
00194      bool active;
00195      NodeSequence oldNodes, newNodes;
00196      FaceSequence oldFaces, newFaces;
00197 
00198      int remove_by_refinement();
00199      int remove_by_swapping();
00200 
00201      int commit();
00202      void clear();
00203 };
00204 
00206 
00207 class OneDefectPatch {
00208 public:
00209      static size_t MAX_FACES_ALLOWED;
00210 
00211      static size_t num_boundaries;
00212      static double exec_time;
00213      static size_t num_3_patches;
00214      static size_t num_4_patches;
00215      static size_t num_5_patches;
00216      static size_t disk_remeshable;
00217 
00218      OneDefectPatch( Mesh *m ) {
00219           mesh = m;
00220           apex = NULL;
00221           quad_splitting_node = NULL;
00222           new_defective_node  = NULL;
00223      }
00224 
00225      OneDefectPatch( Mesh *m, Vertex *v) {
00226           mesh = m;
00227           apex = v;
00228           quad_splitting_node = NULL;
00229           new_defective_node  = NULL;
00230      }
00231 
00232      void set_initial_path( const NodeSequence &sq) {
00233           nodepath = sq;
00234      }
00235 
00236      size_t getSize(int e) const {
00237           if( e == 0) return inner_nodes.size() + bound_nodes.size();
00238           if( e == 2) return faces.size();
00239           return 0;
00240      }
00241 
00242      const NodeSequence &get_irregular_nodes_removed() {
00243           return irregular_nodes_removed;
00244      }
00245 
00246      bool isBoundaryEven() const {
00247           if( bound_nodes.size()%2 == 0 ) return 1;
00248           return 1;
00249      }
00250 
00251      size_t count_irregular_nodes(int where);
00252 
00253      int  build_remeshable_boundary();
00254 
00255      double get_isoperimetic_quotient() const {
00256           // This definiation is taken from Wikipedia..
00257           double A = getArea();
00258           double L = getPerimeter();
00259           double q = 4*M_PI*A/(L*L);
00260           return q;
00261      }
00262 
00263      void getFaces( FaceSequence &result) const {
00264           size_t nSize = faces.size();
00265 
00266           result.clear();
00267           if( nSize == 0 ) return;
00268 
00269           result.resize( nSize );
00270 
00271           int index = 0;
00272           FaceSet::const_iterator it;
00273           for( it = faces.begin(); it != faces.end(); ++it)
00274                result[index++] = *it;
00275      }
00276 
00277      const NodeSequence &getBoundaryNodes() const {
00278           return bound_nodes;
00279      }
00280 
00281      bool operator < ( const OneDefectPatch &rhs) const {
00282           return getSize(2) < rhs.getSize(2);
00283      }
00284 
00285      Vertex *get_new_defective_node() {
00286           return new_defective_node;
00287      }
00288 
00289      int  remesh();
00290 
00291      void setAttributes();
00292 
00293      void clear() {
00294           nodepath.clear();
00295           new_defective_node  = NULL;
00296           quad_splitting_node = NULL;
00297           corners.clear();
00298           inner_nodes.clear();
00299           bound_nodes.clear();
00300           faces.clear();
00301           inner_faces.clear();
00302           irregular_nodes_removed.clear();
00303           boundary.clear();
00304           cornerPos.clear();
00305           segSize.clear();
00306           newnodes.clear();
00307           newfaces.clear();
00308      }
00309 
00310 private:
00311      bool isSafe();
00312 
00313      // Input data.
00314      Mesh   *mesh;
00315      Vertex *apex;               // Seed: Irregular vertex to start from.
00316      NodeSequence  nodepath;     // Initial joining two irregular nodes..
00317 
00318      FaceSet faces;
00319      FaceSet inner_faces;
00320 
00321 #ifdef USE_HASHMAP
00322      std::tr1::unordered_map<Vertex*, FaceSet> relations02;
00323      std::tr1::unordered_map<Vertex*, FaceSet>::iterator miter;
00324 //   std::tr1::unordered_set<Face*> inner_faces;
00325 #else
00326      std::map<Vertex*, FaceSet> relations02;
00327      std::map<Vertex*, FaceSet>::iterator miter;
00328 #endif
00329 
00330      Vertex *new_defective_node;
00331 
00332      // Local data ...
00333      Vertex *quad_splitting_node;     // One special node that splits a quad loop
00334      int quad_splitting_node_degree;  // Valence of the splitting node.
00335 
00336      NodeSet corners;                 // Corners of the blob
00337      NodeSet nodes;                   // All the nodes (inner + boundary)
00338      NodeSequence inner_nodes;        // Inner nodes (not on the boundary ) of the blob
00339      NodeSequence bound_nodes;        // Boundary nodes
00340      NodeSequence irregular_nodes_removed;
00341 
00342      vector<Edge> boundary;          // boundary of the blob.
00343      vector<int>  cornerPos;         // Positions of the corners in the bound_nodes.
00344      vector<int>  segSize;
00345      int   partSegments[10];
00346 
00347      TriRemeshTemplate    template3;
00348      QuadRemeshTemplate   template4;
00349      PentaRemeshTemplate  template5;
00350 
00351      // Variable used in 3-5 sided patch...
00352      NodeSequence anodes, bnodes, cnodes, dnodes, enodes;  // Nodes on each segment.
00353 
00354      // Variables used in 4 sided patch..
00355      NodeSequence a1nodes, a2nodes, b1nodes, c0nodes, c1nodes, c2nodes,
00356      abnodes, canodes, bcnodes, d1nodes;
00357 
00358      // New nodes and faces in the patch...
00359      NodeSequence  newnodes, nnodes;
00360      FaceSequence  newfaces, nfaces;
00361 
00362      // Get the position on the boundary ...
00363      int getPosOf( const Vertex *v) const {
00364           size_t nSize = bound_nodes.size();
00365           for (size_t i = 0; i <  nSize; i++)
00366                if (bound_nodes[i] == v) return i;
00367           return -1;
00368      }
00369      // Return nodes within the range (src, dst)
00370      void get_bound_nodes( const Vertex *src, const Vertex *dst, NodeSequence &s);
00371 
00372      // randomly select one irregular node
00373      bool  has_irregular_node_on_first_segment() const;
00374 
00375      // re-orient boundary nodes so that it starts from a given vertex.
00376      void start_boundary_loop_from (Vertex *v);
00377 
00378      // re-orient loops ...
00379      int reorient_4_sided_loop();
00380 
00381      // Patch creation functions...
00382      int   init_blob();
00383      int   update_boundary();
00384      int   finalize_boundary();
00385      int   expand_blob(Vertex *v);
00386      int   expand_blob();
00387      int   get_topological_outer_angle( Vertex *v);
00388      bool  is_simply_connected();
00389 
00390      bool  is_quad_breakable_at( const Vertex *v);
00391      // Query for the validity of 3-4-5 sided patches.
00392      bool  is_4_sided_convex_loop_quad_meshable();
00393 
00394      // Set the boundary pattern string.
00395      void set_boundary_segments();
00396 
00397      // If the resulting mesh is invalid for some reasons, revert back to
00398      // original and restore all information.
00399      void rollback();
00400 
00401      void pre_remesh();  // Before we start remeshing, do some clean-up
00402      int  remesh_3_sided_patch();
00403      int  remesh_4_sided_patch();
00404      int  remesh_5_sided_patch();
00405      void local_smoothing();
00406      void post_remesh(); // After successful remeshing, do some clean-up
00407 
00408      double getArea() const {
00409           double a = 0.0;
00410           FaceSet::const_iterator it;
00411           assert( faces.size() );
00412           for( it = faces.begin(); it != faces.end(); ++it) {
00413                a += fabs( (*it)->getArea() );
00414           }
00415           return a;
00416      }
00417 
00418      double getPerimeter() const {
00419           double l = 0.0;
00420           int nSize = bound_nodes.size();
00421           for( int i = 0; i < nSize; i++)
00422                l += Vertex::length( bound_nodes[i], bound_nodes[(i+1)%nSize] );
00423           return l;
00424      }
00425 
00426 };
00427 
00429 
00430 class QuadCleanUp {
00431 public:
00432      static bool isDoublet(const Vertex *v);
00433      static bool isSinglet(const Vertex *v);
00434      static bool isRegular( const Vertex *v);
00435      static bool hasSinglet(const Face *f);
00436      static bool isTunnel(const Edge *e);
00437      static bool isEdge33(const Edge *e);
00438      static bool isEdge35(const Edge *e);
00439      static bool isDiamond(Face *f, int &pos, int type = 33);
00440 
00441      QuadCleanUp() {
00442           djkpath = NULL;
00443           defective_patch = NULL;
00444      }
00445 
00446      QuadCleanUp(Mesh *m) {
00447           setMesh(m);
00448           djkpath = NULL;
00449           defective_patch = NULL;
00450      }
00451 
00452      ~QuadCleanUp() {
00453           /*
00454                     if( lapweight ) delete lapweight;
00455                     if( lapsmooth ) delete lapsmooth;
00456           */
00457           if( djkpath   ) delete djkpath;
00458           if( defective_patch ) delete defective_patch;
00459      }
00460 
00461      void setMesh( Mesh *m ) {
00462           mesh = m;
00463           /*
00464                     lapsmooth = new LaplaceSmoothing(mesh);
00465                     lapweight = new LaplaceLengthWeight;
00466                     lapsmooth->setWeight( lapweight );
00467           */
00468      }
00469 
00470      // Query methods ...
00471      NodeSequence  search_restricted_nodes();
00472      FaceSequence  search_restricted_faces();
00473      FaceSequence  search_flat_quads();
00474 
00475      vector<Diamond> search_diamonds(int type = 33 );
00476      vector<Singlet> search_boundary_singlets();
00477      vector<Doublet> search_interior_doublets();
00478      vector<Edge>    search_tunnels();
00479      vector<OneDefectPatch> search_one_defect_patches();
00480 
00481      OneDefectPatch* build_one_defect_patch(Vertex *vertex = NULL );
00482 
00483      int  degree_5_dominated();
00484 
00485      // Global Cleanup methods ..
00486      int remesh_defective_patches();
00487 
00488      // Local Cleanup methods ..
00489      int reduce_degree( Vertex *v );
00490      int vertex_degree_reduction();
00491 
00492      int swap_concave_faces();
00493 
00494      // Removal Methods ...
00495      int  remove_diamonds();
00496      int  remove_tunnels();
00497      int  remove_interior_doublets();
00498      int  remove_boundary_singlets();
00499      int  remove_bridges();
00500      int  shift_irregular_nodes();
00501 
00502 //  int irregular_nodes_clustering();
00503 
00504      //  void remove_ynodes();
00505      int  clean_layer(int id);
00506      void cleanup_boundary(double cutOffAngle = 100.0);
00507      void advancing_front_cleanup();
00508      void advancing_front_edges_swap();
00509 
00510      int  automatic();
00511      void report();
00512 
00513      int atomic_op_swap_edge( Vertex *v0, Vertex *v1);
00514      int atomic_op_face_close( Face *f);
00515 
00516      // Some Feature that may be obsolete in the next version...
00517      Vertex* insert_doublet(Face *face);
00518      Vertex* insert_boundary_doublet(Face *face);
00519      Vertex* insert_doublet(Face *face, Vertex *v0, Vertex *v2);
00520 
00521 //   vector<Edge33>  search_edges33();
00522 //   int remove_edges35();
00523 
00524      int refine_restricted_node(Vertex *resnode, Vertex *bndnode);
00525      int refine_degree3_faces();
00526      int refine_bridges_face();
00527      // Utility functions ...
00528      void get_strips(Face *face, FaceSequence &strip1, FaceSequence strip2);
00529 
00530      int  reduce_internal_vertex_degree(Vertex *v);
00531      int  reduce_boundary_vertex_degree(Vertex *v);
00532 
00533 private:
00534      Mesh *mesh;
00535 
00536      MeshOptimization mopt;
00537      /*
00538           LaplaceSmoothing *lapsmooth;
00539           LaplaceWeight *lapweight;
00540      */
00541 
00542      DijkstraShortestPath *djkpath; // Used in one defect remeshing ....
00543      OneDefectPatch* defective_patch;
00544 
00545      int  has_interior_nodes_degree_345();
00546 
00547      NodeSequence  irregular_nodes;
00548 
00549      vector<OneDefectPatch>  vDefectPatches;
00550      vector<Doublet> vDoublets;
00551      vector<Singlet> vSinglets;
00552      vector<Diamond> vDiamonds; // Diamonds in the mesh;
00553      vector<Diamond> search_diamonds_in_layer(int l);
00554 
00555      int clean_layer_once(int id);
00556      int face_close(Face *face, Vertex *v0, Vertex *v2);
00557      int diamond_collapse(FaceClose &d);
00558      int remove_interior_doublet(Doublet &d);
00559      int remove_boundary_singlet_type1(const Singlet &s);
00560      int remove_boundary_singlet_type2(const Singlet &s);
00561      int remove_boundary_singlets_once();
00562      int remove_bridges_in_layer( int l);
00563      int remove_bridges_once();
00564      int remove_diamonds_once();
00565      int remove_diamonds_in_layer( int l);
00566      int advance_front_edges_swap_once(int layerid);
00567 
00568      int apply_advance_front_bridge_rule( Vertex *v0, Vertex *v1);
00569      int apply_advance_front_excess_rule( Vertex *v);
00570      int apply_advance_front_triplet_rule( Vertex *v);
00571      int apply_advance_front_singlet_rule( Vertex *v);
00572 
00573      int remove_doublets_once();
00574      int remove_interior_doublets_once();
00575 
00576      int boundary_vertex_degree_reduction_once();
00577      int internal_vertex_degree_reduction_once();
00578 
00579      // High level utility function composed of basic functions...
00580      void cleanup_internal_boundary_face();
00581 
00582      // May become obsolere
00583      int refine_3434_pattern( Face *face, int pos);
00584      int refine_3454_pattern( Face *face, int pos);
00585      int refine_3444_pattern( Face *face, int pos);
00586      int apply_shift_node3_rule( Vertex *vertex);
00587 };
00588 
00590 
00591 inline bool
00592 QuadCleanUp::isRegular (const Vertex *v)
00593 {
00594      assert(v);
00595      // Any interior vertex having four nodes( or faces ) is a regular node.
00596      if (!v->isBoundary () && (v->getNumRelations(2) == 4)) return 1;
00597      return 0;
00598 }
00599 
00601 
00602 inline bool
00603 QuadCleanUp::isDoublet (const Vertex *v)
00604 {
00605      assert(v);
00606      // Any interior node having two neighboring face is a  doublet node.
00607      if (!v->isBoundary () && (v->getNumRelations(2) == 2)) return 1;
00608      return 0;
00609 }
00610 
00612 
00613 inline bool
00614 QuadCleanUp::isSinglet (const Vertex *v)
00615 {
00616      assert( v );
00617      // Any boundary node having only one neigbour cell is a singlet node ...
00618      int numfaces = v->getNumRelations(2);
00619      assert (numfaces >= 0);
00620      if (v->isBoundary () && (numfaces == 1)) return 1;
00621      return 0;
00622 }
00623 
00625 
00626 inline bool
00627 QuadCleanUp::hasSinglet (const Face *face)
00628 {
00629      assert( face );
00630 
00631      if( face->isRemoved() ) return 0;
00632 
00633      for (int i = 0; i < face->getSize (0); i++) {
00634           if (isSinglet (face->getNodeAt (i))) return 1;
00635      }
00636      return 0;
00637 }
00638 
00640 
00641 void set_singlet_tag(Mesh *m, const string &s = "Singlet" );
00642 void set_doublet_tag(Mesh *m, const string &s = "Doublet" );
00643 void set_diamond_tag(Mesh *mesh, const string &s = "Diamond" );
00644 
00645 } // namespace Jaal
00646 
00647 #endif
00648 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines