The Parametric Pseudo-Manifold (PPS) Library 1.0
pps.h
Go to the documentation of this file.
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( &param_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