The Parametric Pseudo-Manifold (PPS) Library 1.0
|
00001 00026 #ifndef PPS_H 00027 #define PPS_H 00028 00029 #include <cassert> // assert 00030 #include <cmath> // fabs, sqrt, acos, asin, cos, sin 00031 #include <vector> // std::vector 00032 00033 #include "bezier.h" // Bezier 00034 00049 namespace pps { 00050 00073 template < typename Mesh > 00074 class PPS { 00075 public: 00076 // --------------------------------------------------------------- 00077 // 00078 // Type name definitions for the generic mesh components. 00079 // 00080 // --------------------------------------------------------------- 00081 00087 typedef typename Mesh::Vertex Vertex ; 00088 00089 00095 typedef typename Mesh::Halfedge Halfedge ; 00096 00097 00103 typedef typename Mesh::Edge Edge ; 00104 00105 00111 typedef typename Mesh::Face Face ; 00112 00113 00119 typedef typename Mesh::VertexIterator VertexIterator ; 00120 00121 00127 typedef typename Mesh::EdgeIterator EdgeIterator ; 00128 00129 00135 typedef typename Mesh::FaceIterator FaceIterator ; 00136 00137 00138 // --------------------------------------------------------------- 00139 // 00140 // Public methods 00141 // 00142 // --------------------------------------------------------------- 00143 00152 PPS( Mesh* mesh ) : _mesh( mesh ) , _MYPI( acos( -1 ) ) 00153 {} 00154 00155 00161 virtual ~PPS() {} 00162 00163 00172 void build() ; 00173 00174 00195 void eval_pps( 00196 Face* face , 00197 double u , 00198 double v , 00199 double w , 00200 double& x , 00201 double& y , 00202 double& z 00203 ) 00204 const ; 00205 00206 00225 virtual void eval_surface( 00226 Face* face , 00227 double u , 00228 double v , 00229 double w , 00230 double& x , 00231 double& y , 00232 double& z 00233 ) 00234 const = 0 ; 00235 00236 00244 inline Mesh* get_mesh() const { 00245 return _mesh; 00246 } 00247 00248 00249 protected: 00250 00251 // --------------------------------------------------------------- 00252 // 00253 // Protected methods 00254 // 00255 // --------------------------------------------------------------- 00256 00257 00265 void build_shape_functions() ; 00266 00267 00287 void build_one_shape_function( Vertex* vertex ) ; 00288 00289 00305 inline unsigned get_shape_function_degree( unsigned D ) const 00306 { 00307 return ( D > 6 ) ? 7 : ( D + 1 ) ; 00308 } 00309 00310 00324 inline unsigned get_number_of_parameter_points( unsigned D ) const 00325 { 00326 return ( D << 1 ) ; 00327 } 00328 00329 00341 inline unsigned get_degree( Halfedge* h ) const 00342 { 00343 return get_degree( get_org( h ) ) ; 00344 } 00345 00346 00374 bool find_triangle( 00375 Halfedge*& h , 00376 double x , 00377 double y , 00378 double& u , 00379 double& v , 00380 double& w 00381 ) 00382 const ; 00383 00384 00411 void get_barycentric_coordinates( 00412 double x0 , 00413 double y0 , 00414 double x1 , 00415 double y1 , 00416 double x2 , 00417 double y2 , 00418 double xp , 00419 double yp , 00420 double& u , 00421 double& v , 00422 double& w 00423 ) 00424 const ; 00425 00426 00453 void eval_pps( 00454 Halfedge* h , 00455 double u , 00456 double v , 00457 double w , 00458 double& x , 00459 double& y , 00460 double& z 00461 ) 00462 const ; 00463 00464 00478 void select_pdomain( 00479 Face* face , 00480 double& u , 00481 double& v , 00482 double& w , 00483 Halfedge*& h 00484 ) 00485 const ; 00486 00487 00513 void from_barycentric_to_Cartesian( 00514 double u , 00515 double v , 00516 double w , 00517 double x0 , 00518 double y0 , 00519 double x1 , 00520 double y1 , 00521 double x2 , 00522 double y2 , 00523 double& x , 00524 double& y 00525 ) 00526 const ; 00527 00528 00553 void compute_pdomain_contribution( 00554 Halfedge* h , 00555 double u , 00556 double v , 00557 double& weight, 00558 double& x , 00559 double& y , 00560 double& z 00561 ) 00562 const ; 00563 00564 00585 void gfunction( 00586 Halfedge* h , 00587 double x , 00588 double y , 00589 double& u , 00590 double& v 00591 ) 00592 const ; 00593 00594 00611 void to_canonical_domain( 00612 Halfedge* h , 00613 double x , 00614 double y , 00615 double& u , 00616 double& v 00617 ) 00618 const ; 00619 00620 00637 void from_canonical_domain( 00638 Halfedge* h , 00639 double x , 00640 double y , 00641 double& u , 00642 double& v 00643 ) 00644 const ; 00645 00646 00662 void rot_2d( 00663 double x , 00664 double y , 00665 double ang , 00666 double& u , 00667 double& v 00668 ) 00669 const ; 00670 00671 00690 virtual unsigned int get_id( Halfedge* h ) const ; 00691 00692 00704 double weight_function( double x , double y , double R ) const ; 00705 00706 00719 double eta_function( double s , double d1 , double d2 ) const ; 00720 00721 00722 // --------------------------------------------------------------- 00723 // 00724 // Abstract methods - must be implemented by the PPS library user. 00725 // 00726 // --------------------------------------------------------------- 00727 00738 virtual bool mesh_has_boundary() const = 0 ; 00739 00740 00751 virtual bool mesh_is_simplicial() const = 0 ; 00752 00753 00765 virtual Vertex* get_org( Halfedge* h ) const = 0 ; 00766 00767 00780 virtual Vertex* get_dst( Halfedge* h ) const = 0 ; 00781 00782 00794 virtual Edge* get_edge( Halfedge* h ) const = 0 ; 00795 00796 00808 virtual Face* get_face( Halfedge* h ) const = 0 ; 00809 00810 00824 virtual Halfedge* get_prev( Halfedge* h ) const = 0 ; 00825 00826 00840 virtual Halfedge* get_next( Halfedge* h ) const = 0 ; 00841 00842 00854 virtual Halfedge* get_mate( Halfedge* h ) const = 0 ; 00855 00856 00868 virtual Halfedge* get_halfedge( Face* face ) const = 0 ; 00869 00870 00885 virtual Halfedge* get_halfedge( Vertex* vertex ) const = 0 ; 00886 00887 00900 virtual unsigned int get_degree( Vertex* vertex ) const = 0 ; 00901 00902 00915 virtual Bezier* get_shape_function( Vertex* vertex ) const = 0 ; 00916 00917 00928 virtual void set_shape_function( Vertex* vertex, Bezier* patch ) = 0 ; 00929 00930 00931 public: 00932 00933 // --------------------------------------------------------------- 00934 // 00935 // Public abstract methods - must be implemented by the user. 00936 // 00937 // --------------------------------------------------------------- 00938 00949 virtual VertexIterator vertices_begin() const = 0 ; 00950 00951 00966 virtual bool is_done( const VertexIterator& iterator ) const = 0 ; 00967 00968 00979 virtual void move_forward( VertexIterator& iterator ) const = 0 ; 00980 00981 00994 virtual Vertex* get_vertex( const VertexIterator& iterator ) const = 0 ; 00995 00996 01007 virtual EdgeIterator edges_begin() const = 0 ; 01008 01009 01024 virtual bool is_done( const EdgeIterator& iterator ) const = 0 ; 01025 01026 01036 virtual void move_forward( EdgeIterator& iterator ) const = 0 ; 01037 01038 01051 virtual Edge* get_edge( const EdgeIterator& iterator ) const = 0 ; 01052 01053 01064 virtual FaceIterator faces_begin() const = 0 ; 01065 01066 01081 virtual bool is_done( const FaceIterator& iterator ) const = 0 ; 01082 01083 01093 virtual void move_forward( FaceIterator& iterator ) const = 0 ; 01094 01095 01108 virtual Face* get_face( const FaceIterator& iterator ) const = 0 ; 01109 01110 01111 protected: 01112 01113 // --------------------------------------------------------------- 01114 // 01115 // Protected data members 01116 // 01117 // --------------------------------------------------------------- 01118 01119 Mesh* _mesh ; 01120 01121 const double _MYPI ; 01122 01123 }; 01124 01125 01126 // ----------------------------------------------------------------- 01127 // 01128 // Implementation of non-virtual public methods 01129 // 01130 // ----------------------------------------------------------------- 01131 01149 template< typename Mesh > 01150 void 01151 PPS< Mesh >::eval_pps( 01152 Face* face , 01153 double u , 01154 double v , 01155 double w , 01156 double& x , 01157 double& y , 01158 double& z 01159 ) 01160 const 01161 { 01162 /* 01163 * The barycentric coordinates must define a point in the face. 01164 */ 01165 assert( ( u >= 0 ) && ( u <= 1 ) ); 01166 assert( ( v >= 0 ) && ( v <= 1 ) ); 01167 assert( ( w >= 0 ) && ( w <= 1 ) ); 01168 assert( fabs( 1 - ( u + v + w ) ) <= 1e-15 ) ; 01169 01170 /* 01171 * We choose the p-domain whose distance from the given point is 01172 * the smallest. 01173 */ 01174 01175 Halfedge* h ; 01176 double uaux = u ; 01177 double vaux = v ; 01178 double waux = w ; 01179 01180 select_pdomain( face , uaux , vaux , waux , h ) ; 01181 01182 /* 01183 * Evaluate the PPS using the chosen p-domain. 01184 */ 01185 eval_pps( h , uaux , vaux , waux , x , y , z ) ; 01186 01187 return ; 01188 } 01189 01190 01191 // ----------------------------------------------------------------- 01192 // 01193 // Implementation of non-virtual protected methods. 01194 // 01195 // ----------------------------------------------------------------- 01196 01205 template < typename Mesh > 01206 void 01207 PPS< Mesh >::build() 01208 { 01209 /* 01210 * The given underlying mesh must be a triangle mesh. 01211 */ 01212 assert( mesh_is_simplicial() ) ; 01213 01214 /* 01215 * The given mesh must have an empty boundary. 01216 */ 01217 assert( !mesh_has_boundary() ) ; 01218 01219 /* 01220 * Compute the shape functions of this PPS. 01221 */ 01222 build_shape_functions() ; 01223 } 01224 01225 01236 template< typename Mesh > 01237 void 01238 PPS< Mesh >::build_shape_functions() 01239 { 01240 /* 01241 * For each vertex "v" of the underlying mesh, sample the 01242 * P-polygon associated with it, and then compute the 01243 * corresponding points in the generic parametric patch. These 01244 * points are then used to compute the control points of the shape 01245 * functions. 01246 */ 01247 for ( VertexIterator vi = vertices_begin() ; !is_done( vi ) ; move_forward( vi ) ) { 01248 /* 01249 * Get the current vertex. 01250 */ 01251 Vertex* vertex = get_vertex( vi ); 01252 01253 /* 01254 * Compute the shape function corresponding to the current 01255 * vertex. 01256 */ 01257 build_one_shape_function( vertex ) ; 01258 } 01259 01260 return ; 01261 } 01262 01263 01281 template < typename Mesh > 01282 void 01283 PPS< Mesh >::build_one_shape_function( Vertex* vertex ) { 01284 #ifdef DEBUGMODE 01285 01289 assert ( vertex != 0 ) ; 01290 #endif 01291 01292 /* 01293 * Get one halfedge with origin at vertex \var vertex. 01294 */ 01295 Halfedge* h = get_halfedge( vertex ) ; 01296 01297 #ifdef DEBUGMODE 01298 01302 assert ( h != 0 ) ; 01303 #endif 01304 01305 unsigned int nu = get_degree( h ) ; 01306 01307 #ifdef DEBUGMODE 01308 01312 assert ( nu >= 3 ) ; 01313 #endif 01314 01315 /* 01316 * Generate a rectangular grid with (N + 1) * (N + 1) parameter 01317 * points. 01318 */ 01319 01320 const unsigned int D = get_shape_function_degree( nu ) ; 01321 const unsigned int N = get_number_of_parameter_points( D ) ; 01322 01323 const double R = cos( _MYPI / nu ) ; 01324 01325 /* 01326 * Coordinates of the lower leftmost point of the grid. 01327 */ 01328 const double x0 = -R ; 01329 const double y0 = -R ; 01330 01331 /* 01332 * Spacing between two consecutives points in both X and Y 01333 * directions. 01334 */ 01335 const double dd = ( 2 * R ) / N ; 01336 01337 /* 01338 * Initialize the Y direction spacing counter. 01339 */ 01340 double dy = 0 ; 01341 01342 std::vector< double > param_pts ; // Array of parameter points. 01343 std::vector< double > patch_pts ; // Array of surface points. 01344 01345 for (unsigned int j = 0 ; j <= N ; j++ ) { 01346 /* Y coordinate of the point. */ 01347 double y = y0 + dy ; 01348 01349 /* 01350 * Initialize the X direction spacing counter. 01351 */ 01352 double dx = 0 ; 01353 01354 for ( unsigned int i = 0 ; i <= N ; i++ ) { 01355 /* X coordinate of the point. */ 01356 double x = x0 + dx ; 01357 01358 /* 01359 * If this parameter point is inside the P-polygon the 01360 * P-domain is inscribed in, then we can compute a point on 01361 * the given generic patch. 01362 */ 01363 01364 double u ; 01365 double v ; 01366 double w ; 01367 Halfedge* haux = h ; 01368 if ( find_triangle( haux , x , y , u , v , w ) ) { 01369 /* 01370 * Store the coordinates of the parameter point. 01371 */ 01372 param_pts.push_back( x ) ; 01373 param_pts.push_back( y ) ; 01374 01375 /* 01376 * Compute the corresponding point in the generic surface 01377 * patch. 01378 */ 01379 #ifdef DEBUGMODE 01380 assert( haux != 0 ) ; 01381 #endif 01382 01383 /* Get the face the half-edge haux belongs to. */ 01384 Face* face = get_face( haux ) ; 01385 01386 #ifdef DEBUGMODE 01387 assert( face != 0 ) ; 01388 #endif 01389 01390 double pt[ 3 ] ; 01391 01392 if ( haux == get_halfedge( face ) ) { 01393 eval_surface( face , u , v , w , pt[ 0 ] , pt[ 1 ] , pt[ 2 ] ) ; 01394 } 01395 else if ( get_next( haux ) == get_halfedge( face ) ) { 01396 eval_surface( face , v , w , u , pt[ 0 ] , pt[ 1 ] , pt[ 2 ] ) ; 01397 } 01398 else { 01399 eval_surface( face , w , u , v , pt[ 0 ] , pt[ 1 ] , pt[ 2 ] ) ; 01400 } 01401 01402 /* 01403 * Store the coordinates of the point on the generic surface 01404 * patch. 01405 */ 01406 patch_pts.push_back( pt[0] ) ; 01407 patch_pts.push_back( pt[1] ) ; 01408 patch_pts.push_back( pt[2] ) ; 01409 } 01410 01411 /* Increment X direction spacing counter. */ 01412 dx += dd ; 01413 } 01414 01415 /* Increment Y direction spacing counter. */ 01416 dy += dd ; 01417 01418 } 01419 01420 /* Creates the shape function. */ 01421 set_shape_function( 01422 vertex , 01423 new Bezier( ¶m_pts[0] , 01424 &patch_pts[0] , 01425 param_pts.size() >> 1 , 01426 D , 01427 D , 01428 -R , 01429 -R , 01430 R , 01431 R 01432 ) 01433 ) ; 01434 01435 01436 return ; 01437 01438 } 01439 01440 01468 template < typename Mesh > 01469 bool 01470 PPS< Mesh >::find_triangle( 01471 Halfedge*& h , 01472 double x , 01473 double y , 01474 double& u , 01475 double& v , 01476 double& w 01477 ) 01478 const 01479 { 01480 /* 01481 * Get the degree of the vertex associated with the P-polygon. 01482 */ 01483 unsigned int nu = get_degree( h ) ; 01484 01485 /* 01486 * Loop over all triangles of the canonical triangulation of the 01487 * P-polygon associated with the origin vertex of halfedge h. For 01488 * each triangle, checks if the given point belongs to it. If so, 01489 * compute the barycentric coordinates of the point with respect 01490 * to the triangle. 01491 */ 01492 Halfedge* haux = h ; 01493 do { 01494 01495 #ifdef DEBUGMODE 01496 assert( haux != 0 ) ; 01497 #endif 01498 01499 /* 01500 * Get the identifier of the current half-edge. 01501 */ 01502 unsigned int id = get_id( haux ) ; 01503 01504 /* 01505 * Compute the vertices of the P-polygon canonical triangulation 01506 * triangle associated with the face that contains the half-edge 01507 * haux. 01508 */ 01509 const double ang = 2 * ( _MYPI / nu ) ; 01510 01511 double x1 = cos( id * ang ) ; 01512 double y1 = sin( id * ang ) ; 01513 double x2 = cos( ( id + 1 ) * ang ) ; 01514 double y2 = sin( ( id + 1 ) * ang ) ; 01515 01516 /* 01517 * Compute the barycentric coordinates of the given point with 01518 * respect to the triangle given by the vertices ( 0 , 0 ) , ( 01519 * x1 , y1 ) , and ( x2 , y2 ). 01520 */ 01521 get_barycentric_coordinates( 0. , 0. , x1 , y1 , x2 , y2 , 01522 x , y , u , v , w ) ; 01523 01524 /* 01525 * If all barycentric coordinates are equal to or greater than 01526 * zero, then the triangle contains the given point and the 01527 * search ends. Otherwise, keep looking for a triangle that 01528 * contains the point. 01529 */ 01530 if ( 01531 ( u >= 0 ) && ( u <= 1 ) && ( v >= 0 ) && 01532 ( v <= 1 ) && ( w >= 0 ) && ( w <= 1 ) 01533 ) 01534 { 01535 h = haux; 01536 return true; 01537 } 01538 else { 01539 haux = get_mate( get_prev( haux ) ) ; 01540 } 01541 01542 } 01543 while ( haux != h ); 01544 01545 /* 01546 * If the code reached this point, then no triangle contains the 01547 * given point. 01548 */ 01549 return false; 01550 } 01551 01552 01579 template < typename Mesh > 01580 void 01581 PPS< Mesh >::get_barycentric_coordinates( 01582 double x0 , 01583 double y0 , 01584 double x1 , 01585 double y1 , 01586 double x2 , 01587 double y2 , 01588 double xp , 01589 double yp , 01590 double& u , 01591 double& v , 01592 double& w 01593 ) 01594 const 01595 { 01596 /* 01597 * Compute the determinant. 01598 */ 01599 double dd = ( x1 * y0 ) - ( x2 * y0 ) - ( x0 * y1 ) 01600 + ( x2 * y1 ) + ( x0 * y2 ) - ( x1 * y2 ) ; 01601 01602 /* 01603 * The determinant cannot be zero. 01604 */ 01605 assert( fabs( dd ) > 1e-16 ) ; 01606 01607 /* 01608 * Compute the barycentric coordinates. 01609 */ 01610 u = ( x2 * y1 ) - ( xp * y1 ) - ( x1 * y2 ) 01611 + ( xp * y2 ) + ( x1 * yp ) - ( x2 * yp ) ; 01612 01613 u /= dd ; 01614 01615 v = ( xp * y0 ) - ( x2 * y0 ) + ( x0 * y2 ) 01616 - ( xp * y2 ) - ( x0 * yp ) + ( x2 * yp ) ; 01617 01618 v /= dd ; 01619 01620 if ( fabs( u ) < 1e-14 ) { 01621 u = 0 ; 01622 } 01623 else if ( fabs( 1 - u ) < 1e-14 ) { 01624 u = 1 ; 01625 } 01626 01627 if ( fabs( v ) < 1e-14 ) { 01628 v = 0 ; 01629 } 01630 else if ( fabs( 1 - v ) < 1e-14 ) { 01631 v = 1 ; 01632 } 01633 01634 w = 1 - u - v ; 01635 01636 if ( fabs( w ) < 1e-14 ) { 01637 w = 0 ; 01638 } 01639 else if ( fabs( 1 - w ) < 1e-14 ) { 01640 w = 1 ; 01641 } 01642 01643 return ; 01644 } 01645 01646 01672 template < typename Mesh > 01673 void 01674 PPS< Mesh >::eval_pps( 01675 Halfedge* h , 01676 double u , 01677 double v , 01678 double w , 01679 double& x , 01680 double& y , 01681 double& z 01682 ) 01683 const 01684 { 01685 #ifdef DEBUGMODE 01686 01690 assert( h != 0 ) ; 01691 #endif 01692 01693 /* 01694 * Computes the Cartesian coordinates of the given point with 01695 * respect to the upper triangle of the canonical domain 01696 * (quadrilateral). 01697 */ 01698 double xc ; 01699 double yc ; 01700 from_barycentric_to_Cartesian( u , v , w , 0. , 0. , 1. , 0. , 01701 0.5 , 0.5 * sqrt( 3. ) , xc , yc ) ; 01702 01703 /* 01704 * Map the point from the canonical domain to the p-domain 01705 * associated with the origin vertex of h. 01706 */ 01707 double xr ; 01708 double yr ; 01709 from_canonical_domain( h , xc , yc , xr , yr ) ; 01710 01711 /* Get the degree of the origin vertex of the given half-edge. */ 01712 unsigned int nu = get_degree( h ) ; 01713 01714 /* 01715 * Compute the radius of the p-domain associated with the origin 01716 * vertex of the given half-edge. 01717 */ 01718 const double R = cos( _MYPI / nu ) ; 01719 01720 double ll = ( xr * xr ) + ( yr * yr ) ; 01721 01722 /* 01723 * The given point must belong to the p-domain. 01724 */ 01725 assert( ll < R * R ) ; 01726 01727 /* 01728 * Initialize the accumulator of the sum of weight function values 01729 * and the accumulator of the contribution of p-domain shape 01730 * functions. 01731 */ 01732 double sw = 0 ; 01733 double sf[ 3 ] = { 0. , 0. , 0. } ; 01734 01735 /* 01736 * Computes the contribution of the p-domain associated with the 01737 * origin vertex of h. 01738 */ 01739 01740 /* Compute the value of the weight function. */ 01741 sw = weight_function( xr , yr , R ) ; 01742 01743 /* 01744 * The weight of the given point can never be zero. 01745 */ 01746 assert( sw >= 1e-16 ) ; 01747 01749 Vertex* vertex = get_org( h ) ; 01750 01751 #ifdef DEBUGMODE 01752 assert( vertex != 0 ) ; 01753 #endif 01754 01755 /* 01756 * Get the shape function associated with the origin vertex of the 01757 * given half-edge. 01758 */ 01759 Bezier* patch = get_shape_function( vertex ) ; 01760 01761 #ifdef DEBUGMODE 01762 assert( patch != 0 ) ; 01763 #endif 01764 01765 /* 01766 * Compute a point on the patch corresponding to the shape 01767 * function. 01768 */ 01769 patch->point( xr , yr , sf[ 0 ] , sf[ 1 ] , sf[ 2 ] ) ; 01770 01771 sf[ 0 ] = sw * sf[ 0 ] ; 01772 sf[ 1 ] = sw * sf[ 1 ] ; 01773 sf[ 2 ] = sw * sf[ 2 ] ; 01774 01775 /* 01776 * Compute the contribution of the other two p-domains (if any). 01777 */ 01778 if ( ( u != 1 ) && ( v != 1 ) && ( w != 1 ) ) { 01779 /* 01780 * If the point lies on an edge, then we can consider only one 01781 * more p-domain. 01782 */ 01783 if ( w == 0 ) { 01784 double ur ; 01785 double vr ; 01786 gfunction( h , xr , yr , ur , vr ) ; 01787 01788 compute_pdomain_contribution( get_mate( h ) , 01789 ur , 01790 vr , 01791 sw , 01792 sf[ 0 ] , 01793 sf[ 1 ] , 01794 sf[ 2 ] 01795 ) ; 01796 } 01797 else if ( v == 0 ) { 01798 double ur ; 01799 double vr ; 01800 Halfedge* h2 = get_mate( get_prev( h ) ) ; 01801 01802 gfunction( h2 , xr , yr , ur , vr ) ; 01803 01804 compute_pdomain_contribution( get_mate( h2 ) , 01805 ur , 01806 vr , 01807 sw , 01808 sf[ 0 ] , 01809 sf[ 1 ] , 01810 sf[ 2 ] 01811 ) ; 01812 } 01813 else { 01814 double ur ; 01815 double vr ; 01816 gfunction( h , xr , yr , ur , vr ) ; 01817 01818 compute_pdomain_contribution( get_mate( h ) , 01819 ur , 01820 vr , 01821 sw , 01822 sf[ 0 ] , 01823 sf[ 1 ] , 01824 sf[ 2 ] 01825 ) ; 01826 01827 Halfedge* h2 = get_mate( get_prev( h ) ) ; 01828 01829 gfunction( h2 , xr , yr , ur , vr ) ; 01830 01831 compute_pdomain_contribution( get_mate( h2 ) , 01832 ur , 01833 vr , 01834 sw , 01835 sf[ 0 ] , 01836 sf[ 1 ] , 01837 sf[ 2 ] 01838 ) ; 01839 } 01840 } 01841 01842 01843 /* 01844 * Compute the point on the PPS after weighting the contribution 01845 * of the three p-domains. 01846 */ 01847 x = sf[ 0 ] / sw ; 01848 y = sf[ 1 ] / sw ; 01849 z = sf[ 2 ] / sw ; 01850 01851 return ; 01852 } 01853 01854 01868 template< typename Mesh > 01869 void 01870 PPS< Mesh >::select_pdomain( 01871 Face* face , 01872 double& u , 01873 double& v , 01874 double& w , 01875 Halfedge*& h 01876 ) 01877 const 01878 { 01879 #ifdef DEBUGMODE 01880 01884 assert( face != 0 ) ; 01885 #endif 01886 01887 /* 01888 * Find out which P-domain contains the given point after mapping 01889 * the point to the canonical quadrilateral. 01890 */ 01891 double xc ; 01892 double yc ; 01893 from_barycentric_to_Cartesian( u , v , w , 0. , 0. , 1. , 0. , 01894 0.5 , 0.5 * sqrt( 3. ) , xc , yc ) ; 01895 01896 /* 01897 * We choose the p-domain whose distance from the given point is 01898 * the smallest. 01899 */ 01900 double l1 = ( xc * xc ) + ( yc * yc ) ; 01901 01902 double x2 = xc - 1 ; 01903 double l2 = ( x2 * x2 ) + ( yc * yc ) ; 01904 01905 double x3 = xc - 0.5 ; 01906 double y3 = yc - ( 0.5 * sqrt( 3. ) ) ; 01907 01908 double l3 = ( x3 * x3 ) + ( y3 * y3 ) ; 01909 01910 if ( ( l1 <= l2 ) && ( l1 <= l3 ) ) { 01911 /* 01912 * We pick the p-domain associated with the origin vertex of the 01913 * first half-edge of the given face. 01914 */ 01915 h = get_halfedge( face ) ; 01916 } 01917 else if ( l2 <= l3 ) { 01918 /* 01919 * We pick the p-domain associated with the origin vertex of the 01920 * second half-edge of the given face. 01921 */ 01922 h = get_next( get_halfedge( face ) ) ; 01923 01924 double uaux = u ; 01925 u = v ; 01926 v = w ; 01927 w = uaux ; 01928 } 01929 else { 01930 /* 01931 * We pick the p-domain associated with the origin vertex of the 01932 * third half-edge of the given face. 01933 */ 01934 h = get_prev( get_halfedge( face ) ) ; 01935 01936 double uaux = u ; 01937 u = w ; 01938 w = v ; 01939 v = uaux ; 01940 } 01941 01942 return ; 01943 } 01944 01945 01971 template < typename Mesh > 01972 void 01973 PPS< Mesh >::from_barycentric_to_Cartesian( 01974 double u , 01975 double v , 01976 double w , 01977 double x0 , 01978 double y0 , 01979 double x1 , 01980 double y1 , 01981 double x2 , 01982 double y2 , 01983 double& x , 01984 double& y 01985 ) 01986 const 01987 { 01988 x = ( u * x0 ) + ( v * x1 ) + ( w * x2 ) ; 01989 y = ( u * y0 ) + ( v * y1 ) + ( w * y2 ) ; 01990 } 01991 01992 02017 template < typename Mesh > 02018 void 02019 PPS< Mesh >::compute_pdomain_contribution( 02020 Halfedge* h , 02021 double u , 02022 double v , 02023 double& weight, 02024 double& x , 02025 double& y , 02026 double& z 02027 ) 02028 const 02029 { 02030 /* 02031 * Get the degree of the origin vertex of the given half-edge. 02032 */ 02033 unsigned int nu = get_degree( h ) ; 02034 02035 /* 02036 * Get the radius of the boundary circumference of the p-domain 02037 * associated with the origin vertex of the given half-edge. 02038 */ 02039 const double R = cos( _MYPI / nu ) ; 02040 02041 /* 02042 * Compute the squared distance of the point (u,v) to the origin 02043 * of the local coordinate system of the p-domain associated with 02044 * the given half-edge. 02045 */ 02046 double ll = u * u + v * v ; 02047 02048 /* 02049 * If the point (u,v) is inside the p-domain, then compute the 02050 * contribution of the shape function associated with the origin 02051 * vertex of the given half-edge. 02052 */ 02053 if ( ll < R * R ) { 02054 /* 02055 * Compute the weight of the contribution. 02056 */ 02057 double w = weight_function( u , v , R ) ; 02058 weight += w ; 02059 02060 /* 02061 * Get the shape function associated with the origin vertex of 02062 * the given half-edge. 02063 */ 02064 Bezier* patch = get_shape_function( get_org( h ) ) ; 02065 02066 #ifdef DEBUGMODE 02067 assert( patch != 0 ) ; 02068 #endif 02069 02070 /* 02071 * Compute a point on the patch corresponding to the shape 02072 * function. 02073 */ 02074 double fx , fy , fz ; 02075 patch->point( u , v , fx , fy , fz ) ; 02076 02077 /* 02078 * Weight the contribution of the shape function. 02079 */ 02080 x += w * fx ; 02081 y += w * fy ; 02082 z += w * fz ; 02083 } 02084 02085 return ; 02086 } 02087 02088 02109 template < typename Mesh > 02110 void 02111 PPS< Mesh >::gfunction( 02112 Halfedge* h , 02113 double x , 02114 double y , 02115 double& u , 02116 double& v 02117 ) 02118 const 02119 { 02120 02121 #ifdef DEBUGMODE 02122 assert( h != 0 ) ; 02123 #endif 02124 02125 /* 02126 * Map the point in the source gluing domain to the canonical 02127 * lens. 02128 */ 02129 double utemp ; 02130 double vtemp ; 02131 to_canonical_domain( h , x , y , utemp , vtemp ) ; 02132 02133 /* 02134 * Apply the reflection in the canonical quadrilateral. 02135 */ 02136 utemp = 1 - utemp ; 02137 vtemp = - vtemp ; 02138 02139 /* 02140 * Map the point in the canonical lens to the target p-domain. 02141 */ 02142 from_canonical_domain( get_mate( h ) , utemp , vtemp , u , v ) ; 02143 } 02144 02145 02162 template < typename Mesh > 02163 void 02164 PPS< Mesh >::to_canonical_domain( 02165 Halfedge* h , 02166 double x , 02167 double y , 02168 double& u , 02169 double& v 02170 ) 02171 const 02172 { 02173 /* 02174 * Get the degree of the origin vertex of the given half-edge. 02175 */ 02176 const unsigned int nu = get_degree( h ) ; 02177 02178 /* 02179 * Get the rotation angle associated with the given half-edge. 02180 */ 02181 const double au = ( 2 * _MYPI ) / nu ; 02182 02183 /* 02184 * Get the identifier of the given half-edge. 02185 */ 02186 const unsigned i = get_id( h ) ; 02187 02188 /* 02189 * Perform a 2D rotation by -i * au around the origin of the 02190 * p-domain associated with the origin vertex of the given 02191 * half-edge. 02192 */ 02193 rot_2d( x , y , -( i * au ) , u , v ) ; 02194 02195 /* 02196 * Change from Cartesian to polar coordinates. 02197 */ 02198 double rr = sqrt( ( u * u ) + ( v * v ) ) ; 02199 02200 02201 if ( fabs( rr ) > 1e-16 ) { 02202 double aa ; 02203 02204 if ( u < 0 ) { 02205 if ( v >= 0 ) { 02206 aa = _MYPI - acos( fabs( u ) / rr ) ; 02207 } 02208 else { 02209 aa = -_MYPI + acos( fabs( u ) / rr ) ; 02210 } 02211 } 02212 else { 02213 if ( v >= 0 ) { 02214 aa = asin( v / rr ) ; 02215 } 02216 else { 02217 aa = -asin( fabs( v ) / rr ) ; 02218 } 02219 } 02220 02221 aa *= ( nu / 6. ) ; 02222 02223 rr *= ( cos( _MYPI / 6. ) / cos( _MYPI / nu ) ) ; 02224 02225 u = rr * cos( aa ) ; 02226 v = rr * sin( aa ) ; 02227 } 02228 else { 02229 u = v = 0. ; 02230 } 02231 02232 if ( fabs( u ) < 1e-15 ) u = 0. ; 02233 02234 if ( fabs( v ) < 1e-15 ) v = 0. ; 02235 02236 if ( fabs( 1 - u ) < 1e-15 ) u = 1. ; 02237 02238 if ( fabs( 1 - v ) < 1e-15 ) v = 1. ; 02239 02240 return ; 02241 } 02242 02243 02260 template < typename Mesh > 02261 void 02262 PPS< Mesh >::from_canonical_domain( 02263 Halfedge* h , 02264 double x , 02265 double y , 02266 double& u , 02267 double& v 02268 ) 02269 const 02270 { 02271 02272 #ifdef DEBUGMODE 02273 assert( h != 0 ) ; 02274 #endif 02275 02276 /* 02277 * Get the degree of the origin vertex of the given half-edge. 02278 */ 02279 const unsigned int nu = get_degree( h ) ; 02280 02281 /* 02282 * Get the rotation angle associated with the given half-edge. 02283 */ 02284 const double au = ( 2 * _MYPI ) / nu ; 02285 02286 /* 02287 * Get the identifier of the given half-edge. 02288 */ 02289 const unsigned i = get_id( h ) ; 02290 02291 /* 02292 * Change from Cartesian to polar coordinates. 02293 */ 02294 double rr = sqrt( ( x * x ) + ( y * y ) ) ; 02295 02296 if ( fabs( rr ) > 1e-16 ) { 02297 double aa; 02298 02299 if ( x < 0 ) { 02300 if ( y >= 0 ) { 02301 aa = _MYPI - acos( fabs( x ) / rr ) ; 02302 } 02303 else { 02304 aa = -_MYPI + acos( fabs( x ) / rr ) ; 02305 } 02306 } 02307 else { 02308 if ( y >= 0 ) { 02309 aa = asin( y / rr ) ; 02310 } 02311 else { 02312 aa = -asin( fabs( y ) / rr ) ; 02313 } 02314 } 02315 02316 aa *= ( 6. / double( nu ) ) ; 02317 02318 rr *= ( cos( _MYPI / nu ) / cos( _MYPI / 6 ) ) ; 02319 02320 u = rr * cos( aa ) ; 02321 v = rr * sin( aa ) ; 02322 } 02323 else { 02324 u = v = 0.; 02325 } 02326 02327 /* 02328 * Perform a 2D rotation by i * au around the origin of the 02329 * p-domain associated with the origin vertex of the given 02330 * half-edge. 02331 */ 02332 rot_2d( u , v , ( i * au ) , u , v ) ; 02333 02334 if ( fabs( u ) < 1e-15 ) u =0. ; 02335 02336 if ( fabs( v ) < 1e-15 ) v = 0. ; 02337 02338 if ( fabs( 1 - u ) < 1e-15 ) u = 1. ; 02339 02340 if ( fabs( 1 - v ) < 1e-15 ) v = 1. ; 02341 02342 return ; 02343 } 02344 02345 02361 template < typename Mesh > 02362 void 02363 PPS< Mesh >::rot_2d( 02364 double x , 02365 double y , 02366 double ang , 02367 double& u , 02368 double& v 02369 ) 02370 const 02371 { 02372 u = ( x * cos( ang ) ) - ( y * sin( ang ) ) ; 02373 v = ( x * sin( ang ) ) + ( y * cos( ang ) ) ; 02374 } 02375 02376 02394 template< typename Mesh > 02395 unsigned int 02396 PPS< Mesh >::get_id( Halfedge* h ) const 02397 { 02398 02399 #ifdef DEBUGMODE 02400 assert( h != 0 ) ; 02401 #endif 02402 02403 unsigned int i = 0 ; 02404 Halfedge* h2 = get_halfedge( get_org( h ) ) ; 02405 while ( h2 != h ) { 02406 ++i; 02407 h2 = get_mate( get_prev( h2 ) ) ; 02408 } 02409 02410 return i; 02411 } 02412 02413 02425 template < typename Mesh > 02426 double 02427 PPS< Mesh >::weight_function( double x , double y , double R ) const 02428 { 02429 double ll = sqrt( ( x * x ) + ( y * y ) ) ; 02430 02431 return eta_function( ll , 0.25 * R , R ) ; 02432 } 02433 02434 02447 template < typename Mesh > 02448 double 02449 PPS< Mesh >::eta_function( double s , double d1 , double d2 ) const 02450 { 02451 02452 #ifdef DEBUGMODE 02453 assert( d2 > d1 ) ; 02454 assert( d1 > 0 ) ; 02455 assert( d2 < 1 ) ; 02456 #endif 02457 02458 double res ; 02459 02460 if ( s <= d1 ) { 02461 res = 1 ; 02462 } 02463 else if ( s < d2 ) { 02464 double h1 = ( s - d1 ) / ( d2 - d1 ) ; 02465 double h2 = 1 / sqrt( 1 - h1 ) ; 02466 02467 h1 = 1 / sqrt( h1 ) ; 02468 h2 = exp( h2 - h1 ) ; 02469 02470 res = 1 / ( 1 + ( h2 * h2 ) ) ; 02471 } 02472 else { 02473 res = 0 ; 02474 } 02475 02476 return res ; 02477 } 02478 02479 } 02480 //end of group class. 02482 02483 #endif // PPS_H