The Parametric Pseudo-Manifold (PPS) Library 1.0
pntriangle.cpp
Go to the documentation of this file.
00001 
00025 #include "pntriangle.h"
00026 
00027 #include <cassert>  // assert
00028 #include <cmath>    // fabs, sqrt
00029 
00030 
00045 namespace ppsfrompnt {
00046 
00047   
00063   PNTriangle::PNTriangle(
00064                          double p0[ 3 ] ,
00065                          double p1[ 3 ] ,
00066                          double p2[ 3 ] ,
00067                          double n0[ 3 ] ,
00068                          double n1[ 3 ] ,
00069                          double n2[ 3 ]
00070                         )
00071   {
00072     // ---------------------------------------------------------------
00073     //
00074     // Compute control points b300, b030, and b003.
00075     //
00076     // ---------------------------------------------------------------
00077 
00081     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00082       _b[ 0 ][ i ] = p0[ i ] ; 
00083       _b[ 1 ][ i ] = p1[ i ] ;
00084       _b[ 2 ][ i ] = p2[ i ] ;
00085     }
00086 
00087     // ---------------------------------------------------------------
00088     //
00089     // Compute control points b210, b120, b201, b021, b102, and b012.
00090     //
00091     // ---------------------------------------------------------------
00092 
00100     double w210 = 0 ;
00101     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00102       w210 += ( p1[ i ] - p0[ i ] ) * n0[ i ] ;
00103     }
00104 
00108     double w120 = 0 ;
00109     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00110       w120 += ( p0[ i ] - p1[ i ] ) * n1[ i ] ;
00111     }
00112 
00116     double w201 = 0 ;
00117     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00118       w201 += ( p2[ i ] - p0[ i ] ) * n0[ i ] ;
00119     }
00120 
00124     double w021 = 0 ;
00125     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00126       w021 += ( p2[ i ] - p1[ i ] ) * n1[ i ] ;
00127     }
00128 
00132     double w102 = 0 ;
00133     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00134       w102 += ( p0[ i ] - p2[ i ] ) * n2[ i ] ;
00135     }    
00136 
00140     double w012 = 0 ;
00141     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00142       w012 += ( p1[ i ] - p2[ i ] ) * n2[ i ] ;
00143     }
00144 
00148     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00149       _b[ 3 ][ i ]  = ( 2 * p0[ i ] ) + p1[ i ] - ( w210 * n0[ i ] ) ;
00150       _b[ 3 ][ i ] /= 3. ;
00151     }
00152 
00156     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00157       _b[ 4 ][ i ]  = ( 2 * p1[ i ] ) + p0[ i ] - ( w120 * n1[ i ] ) ;
00158       _b[ 4 ][ i ] /=  3. ;
00159     }
00160 
00164     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00165       _b[ 5 ][ i ]  = ( 2 * p0[ i ] ) + p2[ i ] - ( w201 * n0[ i ] ) ;
00166       _b[ 5 ][ i ] /= 3. ;
00167     }
00168 
00172     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00173       _b[ 6 ][ i ]  = ( 2 * p1[ i ] ) + p2[ i ] - ( w021 * n1[ i ] ) ;
00174       _b[ 6 ][ i ] /= 3. ;
00175     }
00176 
00180     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00181       _b[ 7 ][ i ]  = ( 2 * p2[ i ] ) + p0[ i ] - ( w102 * n2[ i ] ) ;
00182       _b[ 7 ][ i ] /= 3. ;
00183     }
00184 
00188     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00189       _b[ 8 ][ i ]  = ( 2 * p2[ i ] ) + p1[ i ] - ( w012 * n2[ i ] ) ;
00190       _b[ 8 ][ i ] /= 3. ;
00191     }
00192 
00193     // ---------------------------------------------------------------
00194     //
00195     // Compute control point b111.
00196     //
00197     // ---------------------------------------------------------------
00198 
00211     double e[ 3 ] ;
00212     double v[ 3 ] ;
00213     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00214       e[ i ] = 0.0 ;
00215       for ( unsigned j = 3 ; j < 9 ; j++ ) {
00216         e[ i ] += _b[ j ][ i ] ;
00217       }
00218 
00219       e[ i ] /= 6. ;
00220 
00221       v[ i ] = ( p0[ i ] + p1[ i ] + p2[ i ] ) / 3. ;
00222 
00223       _b[ 9 ][ i ] = e[ i ] + ( 0.5 * ( e[ i ] - v[ i ] ) ) ;
00224     }
00225 
00226     // ---------------------------------------------------------------
00227     // 
00228     // Compute unit normals n200, n020, and n002.
00229     //
00230     // ---------------------------------------------------------------
00231 
00235     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00236       _n[ 0 ][ i ] = n0[ i ] ; 
00237       _n[ 1 ][ i ] = n1[ i ] ;
00238       _n[ 2 ][ i ] = n2[ i ] ;
00239     }
00240 
00241     // ---------------------------------------------------------------
00242     // 
00243     // Compute unit normals n110, n101, and n011.
00244     //
00245     // ---------------------------------------------------------------
00246 
00259     double mv0[ 3 ] ;
00260     double mv1[ 3 ] ;
00261     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00262       mv0[ i ] = p1[ i ] - p0[ i ] ;
00263       mv1[ i ] = n0[ i ] + n1[ i ] ;
00264     }
00265 
00270     double d110 = 0 ;
00271     double w110 = 0 ;
00272     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00273       d110 += mv0[ i ] * mv0[ i ] ;
00274       w110 += mv0[ i ] * mv1[ i ] ;
00275     }
00276 
00277     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00278       _n[ 3 ][ i ] = mv1[ i ] - ( ( 2 * w110 ) / d110 ) * ( mv0[ i ] ) ;
00279     }
00280 
00289     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00290       mv0[ i ] = p2[ i ] - p0[ i ] ;
00291       mv1[ i ] = n0[ i ] + n2[ i ] ;
00292     }
00293 
00298     double d101 = 0 ;
00299     double w101 = 0 ;
00300     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00301       d101 += mv0[ i ] * mv0[ i ] ;
00302       w101 += mv0[ i ] * mv1[ i ] ;
00303     }
00304 
00305     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00306       _n[ 4 ][ i ] = mv1[ i ] - ( ( 2 * w101 ) / d101 ) * ( mv0[ i ] ) ;
00307     }
00308 
00317     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00318       mv0[ i ] = p2[ i ] - p1[ i ] ;
00319       mv1[ i ] = n1[ i ] + n2[ i ] ;
00320     }
00321 
00326     double d011 = 0 ;
00327     double w011 = 0 ;
00328     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00329       d011 += mv0[ i ] * mv0[ i ] ;
00330       w011 += mv0[ i ] * mv1[ i ] ;
00331     }
00332 
00333     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00334       _n[ 5 ][ i ] = mv1[ i ] - ( ( 2 * w011 ) / d011 ) * ( mv0[ i ] ) ;
00335     }
00336 
00340     double l[ 3 ] ;
00341     for ( unsigned i = 3 ; i < 6 ; i++ ) {
00342       l[ i - 3 ] = 0.0 ;
00343 
00344       for ( unsigned j = 0 ; j < 3 ; j++ ) {
00345         l[ i - 3 ] += _n[ i ][ j ] * _n[ i ][ j ] ; 
00346       }
00347 
00348       double sqrtl = sqrt( l[ i - 3 ] ) ;
00349 
00350       for ( unsigned j = 0 ; j < 3 ; j++ ) {
00351         _n[ i ][ j ] /= sqrtl ;
00352       }      
00353     }
00354 
00355     return ;
00356   }
00357 
00358 
00366   PNTriangle::PNTriangle( const PNTriangle& pnt ) 
00367   {
00368     if ( &pnt != this ) {
00369       for ( unsigned i = 0 ; i < 10 ; i++ ) {
00370         _b[ i ][ 0 ] = pnt._b[ i ][ 0 ] ;
00371         _b[ i ][ 1 ] = pnt._b[ i ][ 1 ] ;
00372         _b[ i ][ 2 ] = pnt._b[ i ][ 2 ] ;
00373         
00374         if ( i < 6 ) {
00375           _n[ i ][ 0 ] = pnt._n[ i ][ 0 ] ;
00376           _n[ i ][ 1 ] = pnt._n[ i ][ 1 ] ;
00377           _n[ i ][ 2 ] = pnt._n[ i ][ 2 ] ;
00378         }         
00379       }
00380     }
00381     
00382     return ;
00383   }
00384 
00385 
00398   void 
00399   PNTriangle::point(
00400                     double u ,
00401                     double v ,
00402                     double& x ,
00403                     double& y ,
00404                     double& z
00405                     )
00406     const
00407   {
00408     if ( fabs( u ) <= 1e-15 ) {
00409       u = 0 ;
00410     }
00411     else if ( fabs( 1 - u ) <= 1e-15 ) {
00412       u = 1 ;
00413     }
00414 
00415     if ( fabs( v ) <= 1e-15 ) {
00416       v = 0 ;
00417     }
00418     else if ( fabs( 1 - v ) <= 1e-15 ) {
00419       v = 1 ;
00420     }
00421 
00422     double w = 1 - ( u + v ) ;
00423 
00424     if ( fabs( w ) <= 1e-15 ) {
00425       w = 0 ;
00426     }
00427     else if ( fabs( 1 - w ) <= 1e-15 ) {
00428       w = 1 ;
00429     }
00430 
00431     double u2 = u * u ;
00432     double v2 = v * v ;
00433     double w2 = w * w ;
00434 
00435     double u3 = u * u2 ;
00436     double v3 = v * v2 ;
00437     double w3 = w * w2 ;
00438 
00439     double pt[ 3 ] ;
00440     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00441       pt[ i ] =            u3       * _b[ 0 ][ i ]
00442               +            v3       * _b[ 1 ][ i ]
00443               +            w3       * _b[ 2 ][ i ]
00444               + ( 3 * u2 * v  )     * _b[ 3 ][ i ] 
00445               + ( 3 * u  * v2 )     * _b[ 4 ][ i ]
00446               + ( 3 * u2 * w  )     * _b[ 5 ][ i ] 
00447               + ( 3 * v2 * w  )     * _b[ 6 ][ i ]
00448               + ( 3 * u  * w2 )     * _b[ 7 ][ i ]
00449               + ( 3 * v  * w2 )     * _b[ 8 ][ i ]
00450               + ( 6 * u  * v  * w ) * _b[ 9 ][ i ] ; 
00451     }
00452 
00453     x = pt[ 0 ] ;
00454     y = pt[ 1 ] ;
00455     z = pt[ 2 ] ;
00456 
00457     return ;
00458   }
00459 
00460 
00473   void
00474   PNTriangle::normal(
00475                       double u ,
00476                       double v ,
00477                       double& x ,
00478                       double& y ,
00479                       double& z
00480                      )
00481     const
00482   {
00483 
00484     if ( fabs( u ) <= 1e-15 ) {
00485       u = 0 ;
00486     }
00487     else if ( fabs( 1 - u ) <= 1e-15 ) {
00488       u = 1 ;
00489     }
00490 
00491     if ( fabs( v ) <= 1e-15 ) {
00492       v = 0 ;
00493     }
00494     else if ( fabs( 1 - v ) <= 1e-15 ) {
00495       v = 1 ;
00496     }
00497 
00498     double w = 1 - ( u + v ) ;
00499 
00500     if ( fabs( w ) <= 1e-15 ) {
00501       w = 0 ;
00502     }
00503     else if ( fabs( 1 - w ) <= 1e-15 ) {
00504       w = 1 ;
00505     }
00506 
00507     assert( 
00508            ( u >= 0 ) && ( u <= 1 ) && ( v >= 0 ) && 
00509            ( v <= 1 ) && ( w >= 0 ) && ( w <= 1 )
00510           ) ;
00511 
00512     double u2 = u * u ;
00513     double v2 = v * v ;
00514     double w2 = w * w ;
00515 
00516     double pt[ 3 ] ;
00517     for ( unsigned i = 0 ; i < 3 ; i++ ) {
00518       pt[ i ] =   u2       * _n[ 0 ][ i ]
00519               +   v2       * _n[ 1 ][ i ]
00520               +   w2       * _n[ 2 ][ i ]
00521               + ( u  * v ) * _n[ 3 ][ i ] 
00522               + ( u  * w ) * _n[ 4 ][ i ]
00523               + ( v  * w ) * _n[ 5 ][ i ] ;
00524      }
00525 
00526     x = pt[ 0 ] ;
00527     y = pt[ 1 ] ;
00528     z = pt[ 2 ] ;
00529 
00530     double l = sqrt( ( x * x ) + ( y * y ) + ( z * z ) ) ;
00531 
00532     assert( l > 0 ) ;
00533 
00534     x /= l ;
00535     y /= l ;
00536     z /= l ;
00537 
00538     return ; 
00539   }
00540 
00541 }
00542  //end of group class.