QuOc

 

triMesh.cpp

Go to the documentation of this file.
00001 /* Copyright (c) 2001-2010 AG Rumpf, INS, Universitaet Bonn                      */
00002 /*                                                                               */
00003 /* The contents of this file are subject to the terms of the Common Development  */
00004 /* and Distribution License Version 1.0 (the "License"); you may not use         */
00005 /* this file except in compliance with the License. You may obtain a copy of     */
00006 /* the License at http://www.sun.com/cddl/                                       */
00007 /*                                                                               */
00008 /* Software distributed under the License is distributed on an "AS IS" basis,    */
00009 /* WITHOUT WARRANTY OF ANY KIND, either express or implied.                      */
00010 /*                                                                               */
00011 #ifndef USE_EXTERNAL_OPENMESH
00012 #error OpenMesh external is needed to compile module openmesh
00013 #endif
00014 
00015 #include <triMesh.h>
00016 #include <openMeshOps.h>
00017 #include <preconditioner.h>
00018 
00019 namespace om {
00020 
00021 template <typename RealType>
00022 RealType TriMesh<RealType>::enclosedVolume( ) const {
00023   RealType v = 0.;
00024   int num = 0;
00025 
00026   for ( ElementIteratorType it = *this; it.notAtEnd(); ++it ) {
00027     aol::Vec3<RealType> n, c;
00028     it->barycenter ( c );
00029     it->normal ( n );
00030     v += ( n * c ) * it->vol();
00031     num++;
00032   }
00033   return v / 3.;
00034 }
00035 
00036 template <typename RealType>
00037 RealType TriMesh<RealType>::area( ) const {
00038   RealType a = 0.;
00039   int num = 0;
00040   for ( ElementIteratorType it = *this; it.notAtEnd(); ++it ) {
00041     a += it->vol();
00042     num++;
00043   }
00044 #ifdef VERBOSE
00045   std::cerr << "anum = " << num << std::endl;
00046 #endif
00047   return a;
00048 }
00049 
00050 template <typename RealType>
00051 void TriMesh<RealType>::normalizeTo01( bool translate ) {
00052   const int size = getNumVertices();
00053   aol::Vec3<RealType> min, max;
00054 
00055   RealType factor = getMeshBoundingBox ( min, max );
00056   Point offset;
00057   for ( int i = 0; i < 3; ++i)
00058     offset[i] = min[i];
00059 #ifdef VERBOSE
00060   std::cerr << "maximum diam = " << factor << endl;
00061 #endif
00062   if ( (factor > 1.0) || translate )
00063     for ( int i = 0; i < size; i++ ) {
00064       Point &p = getOpenMeshObject().point ( getOpenMeshObject().vertex_handle ( i ) );
00065       if (translate)
00066         p -= offset;
00067       if (factor > 1.0)
00068         for ( int j = 0; j < 3; j++ )
00069           p[j] /= factor;
00070     }
00071 }
00072 
00073 
00074 template <typename RealType>
00075 void TriMesh<RealType>::getFaceNormal ( const FaceHandle &fh, aol::Vec3<RealType> &wn, bool unit) const  {
00076 
00077   aol::Mat<3,2,RealType> mat;
00078   aol::Vec3<RealType> dx1, dx2;
00079   getDx ( fh,  mat );
00080   mat.getCol(0, dx1);
00081   mat.getCol(1, dx2);
00082   wn.makeCrossProduct ( dx1, dx2 );
00083   if(unit)
00084     wn.normalize();
00085 }
00086 
00087 
00088 template <typename RealType>
00089 void TriMesh<RealType>::getDx(const FaceHandle &fh,  aol::Mat<3,2,RealType> &Dx ) const {
00090 
00091   // Dx = [x1-x0 | x2-x0] = [ e_1 | -e_0]
00092   Dx.setZero();
00093   aol::Vec3<RealType> x0; 
00094   int i = 0;
00095   for(ConstFaceVertexIter fv_it = getOpenMeshObject().cfv_iter(fh); fv_it; ++fv_it ){
00096     const Point &p = getOpenMeshObject().point ( fv_it.handle() );
00097     if( i == 0 )
00098       for ( int j = 0; j < 3; j++ ) 
00099         x0[j] = p[j];
00100     else
00101       for ( int j = 0; j < 3; j++ ) 
00102         Dx.set( j, i-1,  p[j] - x0[j] );
00103     ++i;
00104   }
00105 }
00106 
00107 
00108 template <typename RealType>
00109 void TriMesh<RealType>::getFirstFundForm(const FaceHandle &fh, aol::Matrix22<RealType> &g ) const{
00110   aol::Mat<3,2,RealType> Dx;
00111   getDx(fh,Dx);
00112   g.makeProductAtransposedB(Dx,Dx);    
00113 }
00114 
00115 
00116 template <typename RealType>
00117 void TriMesh<RealType>::toVector ( aol::MultiVector<RealType> &dst ) const {
00118   const int size = getNumVertices();
00119   dst.reallocate( 3, size );
00120   for ( int i = 0; i < size; i++ ) {
00121     const Point &p = this->point ( getOpenMeshObject().vertex_handle ( i ) );
00122     for ( int j = 0; j < 3; j++ ) {
00123       dst[j][i] = p[j];
00124     }
00125   }
00126 }
00127 
00128 template <typename RealType>
00129 void TriMesh<RealType>::fromVector ( const aol::MultiVector<RealType> &src ) {
00130 
00131   if (  src.numComponents() != 3 )
00132     throw aol::Exception ( "source must have 3 components", __FILE__, __LINE__ );
00133 
00134   const int size = getNumVertices();
00135   if( src[0].size() != size )
00136     throw aol::Exception ( "sizes do not match", __FILE__, __LINE__ );
00137 
00138   for ( int i = 0; i < size; i++ ) {
00139     Point &p = this->point ( getOpenMeshObject().vertex_handle ( i ) );
00140     for ( int j = 0; j < 3; j++ )
00141       p[j] = src[j][i];
00142   }
00143 }
00144 
00145 
00146 template <typename RealType>
00147 void TriMesh<RealType>::regularizeMesh ( ) {
00148   for ( NodeIteratorType vIter = *this; vIter.notAtEnd(); ++vIter ) {
00149     int faces = 0;
00150     aol::Vec3<RealType> wbary, wn;
00151     RealType totalvol = 0;
00152 
00153     wbary.setZero();
00154     wn.setZero();
00155     for ( VertexFaceIter vf_it = getOpenMeshObject().vf_iter ( vIter.vertexHandle() ); vf_it; ++vf_it ) {
00156 
00157       // compute the barycenter
00158       aol::Vec3<RealType> bary, coords[3], v[2], n;
00159       bary.setZero();
00160       int i = 0;
00161       for ( FaceVertexIter iv_it = getOpenMeshObject().fv_iter ( vf_it.handle() ); iv_it; ++iv_it ) {
00162         const Point &p = this->point ( iv_it.handle() );
00163         for ( int j = 0; j < 3; j++ ) { coords[i][j] = p[j]; }
00164         bary += coords[i];
00165         ++i;
00166       }
00167       bary /= 3.;
00168       v[0] = coords[1]; v[0] -= coords[0];
00169       v[1] = coords[2]; v[1] -= coords[0];
00170       n = v[0].crossProduct ( v[1] );
00171       n *= 0.5;
00172       totalvol += n.norm();
00173       wn += n;
00174 
00175       bary *= n.norm();
00176       wbary += bary;
00177       ++faces;
00178     }
00179     wn /= totalvol;
00180     //wbary /= (RealType)faces;
00181     wbary /= totalvol;
00182 
00183     // second cycle: compute t
00184     RealType nom = 0, denom = 0;
00185     for ( VertexFaceIter vf_it = getOpenMeshObject().vf_iter ( vIter.vertexHandle() ); vf_it; ++vf_it ) {
00186 
00187       // compute the barycenter
00188       aol::Vec3<RealType> coords[3];
00189       int center = -1;
00190       int i = 0;
00191       for ( FaceVertexIter iv_it = getOpenMeshObject().fv_iter ( vf_it.handle() ); iv_it; ++iv_it ) {
00192         if ( vIter.vertexHandle() == iv_it.handle() ) {
00193           center = i;
00194         }
00195         ++i;
00196       }
00197 #ifdef VERBOSE
00198       std::cerr << "center = " << center << endl;
00199 #endif
00200       i = 0;
00201       for ( FaceVertexIter iv_it = getOpenMeshObject().fv_iter ( vf_it.handle() ); iv_it; ++iv_it ) {
00202         const Point &p = this->point ( iv_it.handle() );
00203         for ( int j = 0; j < 3; j++ ) { coords[ ( i+3-center ) %3][j] = p[j]; }
00204         ++i;
00205       }
00206       for ( int j = 0; j < 3; j++ ) {
00207         coords[j] -= wbary;
00208       }
00209       aol::Vec3<RealType> cross;
00210       cross = coords[0].crossProduct ( coords[1] );
00211       nom += ( cross * coords[2] );
00212       cross = wn.crossProduct ( coords[1] );
00213       denom += cross * coords[2];
00214     }
00215 
00216     RealType t = nom / denom;
00217     Point &p = this->point ( vIter.vertexHandle() );
00218     // wn *=  t / (1 + aol::Abs(t*100));
00219     wn *= 0.8 * t;
00220     //std::cerr << " t = " << t ;
00221     wbary += wn;
00222 
00223     for ( int j = 0; j < 3; j++ ) {
00224       p[j] = wbary[j];
00225     }
00226   }
00227 }
00228 
00229 
00230 template <typename RealType>
00231 RealType TriMesh<RealType>::getMeshBoundingBox ( aol::Vec3<RealType> &min,
00232                                                     aol::Vec3<RealType> &max ) const {
00233   const int size = getNumVertices();
00234   const Point &p = this->point ( getOpenMeshObject().vertex_handle ( 0 ) );
00235   min[0] = p[0]; min[1] = p[1]; min[2] = p[2];
00236   max = min;
00237 
00238   for ( int i = 1; i < size; i++ ) {
00239     const Point &p = getOpenMeshObject().point ( getOpenMeshObject().vertex_handle ( i ) );
00240     for ( int j = 0; j < 3; j++ ) {
00241       max[j] = aol::Max ( ( RealType ) p[j], max[j] );
00242       min[j] = aol::Min ( ( RealType ) p[j], min[j] );
00243     }
00244   }
00245   return aol::Max ( aol::Max ( max[0] - min[0], max[1] - min[1] ),  max[2] - min[2] );
00246 }
00247 
00248 
00249 template <typename RealType>
00250 void TriMesh<RealType>::loadFromUDPLY ( const std::string filename ) {
00251 
00252   aol::TriangMesh<RealType> triangMesh;
00253   triangMesh.loadFromUDPLY ( filename.c_str() );
00254   importFromTriangMesh(triangMesh);
00255 }
00256 
00257 template <typename RealType>
00258 void TriMesh<RealType>::loadFromPLY ( const std::string filename ) {
00259 
00260   aol::TriangMesh<RealType> triangMesh;
00261   triangMesh.loadFromPLY ( filename.c_str() );
00262   importFromTriangMesh(triangMesh);
00263 }
00264 
00265 // convert this to TriangMesh
00266 template <typename RealType>
00267 void TriMesh<RealType>::convertToTriangMesh ( aol::TriangMesh<RealType> & triangMesh ) const {
00268 
00269   const int numVertex = getNumVertices();
00270   const int numTriang = getNumFaces();
00271   triangMesh.resize ( numVertex, numTriang );
00272 
00273   int nIndex = 0;
00274   for ( NodeIteratorType nIter = *this; nIter.notAtEnd(); ++nIter, ++nIndex ) {
00275     const Point &p = getOpenMeshObject().point ( nIter.vertexHandle() );
00276 
00277     aol::Vec3<RealType> vertex ( p[0], p[1], p[2] );
00278     triangMesh.setVertex (nIndex, vertex);
00279   }
00280 
00281   int nFace = 0;
00282   for ( ElementIteratorType fIter = *this; fIter.notAtEnd(); ++fIter, ++nFace ) {
00283 
00284     int n = 0;
00285     aol::Vec3<int> triangle;
00286     for ( ConstFaceVertexIter v_it = getOpenMeshObject().cfv_iter ( fIter.faceHandle() ); v_it; ++v_it ) {
00287       triangle[n] = v_it.handle().idx();
00288       n++;
00289     }
00290     triangMesh.setTriang ( nFace, triangle );
00291   }
00292 }
00293 
00294 // import from aol::TriangMesh
00295 template <typename RealType>
00296 TriMesh<RealType>& TriMesh<RealType>::importFromTriangMesh(const aol::TriangMesh<RealType>& triangMesh){
00297   this->clear();
00298 
00299   const int numVertex = triangMesh.getNumVertices();
00300   const int numTriang = triangMesh.getNumTriangs();
00301   
00302   for ( int i = 0; i < numVertex; i++ )
00303     addVertex ( triangMesh.getVertex ( i ) );
00304 
00305   for ( int i = 0; i < numTriang; i++ )
00306     addFace ( triangMesh.getTriang ( i ) );
00307   
00308   return (*this);
00309 
00310 }
00311 
00312 //=======================================================================================================
00313 //  TriMeshWithEdgeNormals
00314 //=======================================================================================================
00315 
00316 template <typename RealType>
00317 void TriMeshWithEdgeNormals<RealType>::getEdgeNormal ( const FaceHandle &fh, int localCoord, aol::Vec3<RealType> &wn, bool normalized ) const{
00318   if( localCoord > 2 )
00319     throw aol::Exception ( "unvalid local coordinate (must be 0,1,2)!", __FILE__, __LINE__ );
00320 
00321   getFaceNormal( fh, wn );
00322 
00323   typename TriMesh<RealType>::ConstFaceFaceIter ffIter = this->getOpenMeshObject().cff_iter( fh );
00324   for(int i = 0; i < localCoord; i++)
00325     ++ffIter;
00326 
00327   // edge on boundary?
00328   if(ffIter.handle().idx() != -1 ){
00329     aol::Vec3<RealType> temp;
00330     getFaceNormal( ffIter.handle(), temp );
00331     wn += temp;
00332   }
00333 
00334   if( normalized )
00335     wn.normalize();
00336 }
00337 
00338 
00339 template <typename RealType>
00340 void TriMeshWithEdgeNormals<RealType>::getEdgeNormal ( const EdgeHandle &eh, aol::Vec3<RealType> &wn, bool normalized ) const{
00341   if(!_edgeNormalProp)
00342     throw aol::Exception ( "edge normal property has not been set yet!", __FILE__, __LINE__ );
00343   wn = this->getOpenMeshObject().property( _edgeNormals, eh );
00344   if( normalized )
00345     wn.normalize();
00346 }
00347 
00348 template <typename RealType>
00349 void TriMeshWithEdgeNormals<RealType>::releaseEdgeNormals (){
00350   if( _edgeNormalCount > 0 ){
00351     _edgeNormalCount--;
00352     if ( _edgeNormalCount == 0 ) {
00353       this->getOpenMeshObject().remove_property ( _edgeNormals );
00354       _edgeNormalProp  = false;
00355     }
00356   }
00357 }
00358 
00359 template <typename RealType>
00360 bool TriMeshWithEdgeNormals<RealType>::hasEdgeNormals() const {
00361   return _edgeNormalProp;
00362 }
00363 
00364 template <typename RealType>
00365 int  TriMeshWithEdgeNormals<RealType>::updateEdgeNormals (){
00366   requestEdgeNormals ( );
00367   return _edgeNormalCount;
00368 }
00369 
00370 
00371 template <typename RealType>
00372 void TriMeshWithEdgeNormals<RealType>::requestEdgeNormals(  ) {
00373 
00374   if(_edgeNormalProp==false)
00375     this->add_property(_edgeNormals);
00376 
00377   // iterate over all faces
00378   for ( typename TriMesh<RealType>::ElementIteratorType fIter = *this; fIter.notAtEnd(); ++fIter ) {
00379 
00380     // get face normal
00381     aol::Vec3<RealType> faceNormal;
00382     this->getFaceNormal ( fIter.faceHandle(), faceNormal );
00383 
00384     // iterate over edges
00385     typename TriMesh<RealType>::ConstFaceEdgeIter feIter = this->getOpenMeshObject().cfe_iter( fIter.faceHandle() );
00386     for(; feIter; ++feIter)
00387       this->getOpenMeshObject().property( _edgeNormals, feIter.handle() ) += faceNormal;
00388 
00389   }
00390 
00391   _edgeNormalProp = true;
00392   _edgeNormalCount++;
00393 }
00394 
00395 template <typename RealType>
00396 void TriMeshWithEdgeNormals<RealType>::getSecondFundForm( const FaceHandle &fh, aol::Matrix22<RealType> &h ) const{
00397 
00398   // c = 1/(2|\omega|^2) where \omega represents the unit reference element in \R^2 having vertices at (0,0), (1,0), (0,1)
00399   RealType c = 2.;
00400 
00401   // get edges e_i
00402   aol::Vec3<RealType> edges[3];
00403   typename TriMesh<RealType>::Point points[3];
00404   int i = 0;
00405   for( typename TriMesh<RealType>::ConstFaceVertexIter fv_it = this->getOpenMeshObject().cfv_iter(fh); fv_it; ++fv_it )
00406     points[i++] = this->getOpenMeshObject().point ( fv_it.handle() );
00407 
00408   // in Openmesh vertex i has e_i as incoming edge and e_{(i+1)%3} as outgoing edge
00409   for( int j = 0; j<3; j++ )
00410    for( int k = 0; k<3; k++ )
00411      edges[j][k] = points[j][k] - points[ (j+2)%3 ][k];
00412 
00413   // edge normals
00414   aol::Vec3<RealType> normals[3];
00415   i = 0;
00416   for( typename TriMesh<RealType>::ConstFaceEdgeIter fe_it = this->getOpenMeshObject().cfe_iter(fh); fe_it; ++fe_it ){
00417     getEdgeNormal ( fh, i, normals[i] );
00418     i++;
00419   }
00420 
00421   // h = c \sum_{i=0}^2 ( N_i, e_{\pi(i)} ) B_i, where N_i is the normal on edge e_i, \pi(i):= (i+2)%3, i is the local coordinate
00422   // rotation of e_i by 90 degrees leads to t_i; set B_i := t_i t_i^T
00423   // t_0 = (-1, 0) => B_0 = (1 0; 0 0), t_1 = (0 , -1) =>  B_1 = (0 0; 0 1), t_2 = (1, 1) => B_2 = (1 1; 1 1)
00424   RealType aux = normals[2]*edges[1];
00425   h.set (0, 0, aux + normals[0]*edges[2]);
00426   h.set (0, 1, aux);
00427   h.set (1, 0, aux);
00428   h.set (1, 1, aux + normals[1]*edges[0]);
00429   h *= c;
00430 
00431 }
00432 
00433 template <typename RealType>
00434 void TriMeshWithEdgeNormals<RealType>::getSecondFundFormFast( const FaceHandle &fh, aol::Matrix22<RealType> &h ) const{
00435   if(!_edgeNormalProp)
00436     throw aol::Exception ( "edge normal property has not been set yet!", __FILE__, __LINE__ );
00437 
00438   // c = 1/(2|\omega|^2) where \omega represents the unit reference element in \R^2 having vertices at (0,0), (1,0), (0,1)
00439   RealType c = 2.;
00440 
00441   // get edges e_i
00442   aol::Vec3<RealType> edges[3];
00443   typename TriMesh<RealType>::Point points[3];
00444   int i = 0;
00445   for( typename TriMesh<RealType>::ConstFaceVertexIter fv_it = this->getOpenMeshObject().cfv_iter(fh); fv_it; ++fv_it )
00446     points[i++] = this->getOpenMeshObject().point ( fv_it.handle() );
00447 
00448   // in Openmesh vertex i has e_i as incoming edge and e_{(i+1)%3} as outgoing edge
00449   for( int j = 0; j<3; j++ )
00450    for( int k = 0; k<3; k++ )
00451      edges[j][k] = points[j][k] - points[ (j+2)%3 ][k];
00452 
00453   // edge normals
00454   aol::Vec3<RealType> normals[3];
00455   i = 0;
00456   for( typename TriMesh<RealType>::ConstFaceEdgeIter fe_it = this->getOpenMeshObject().cfe_iter(fh); fe_it; ++fe_it )
00457     getEdgeNormal ( fe_it.handle(), normals[i++] );
00458 
00459   // h = c \sum_{i=0}^2 ( N_i, e_{\pi(i)} ) B_i, where N_i is the normal on edge e_i, \pi(i):= (i+2)%3, i is the local coordinate
00460   // rotation of e_i by 90 degrees leads to t_i; set B_i := t_i t_i^T
00461   // t_0 = (0, -1) => B_0 = (1 0; 0 0), t_1 = (1 , 1) =>  B_1 = (1 1; 1 1), t_2 = (-1, 0) => B_2 = (0 0; 0 1)
00462   RealType aux = normals[1]*edges[0];
00463   h.set (0, 0, aux + normals[0]*edges[2]);
00464   h.set (0, 1, aux);
00465   h.set (1, 0, aux);
00466   h.set (1, 1, aux + normals[2]*edges[1]);
00467   h *= c;
00468 }
00469 
00470 template <typename RealType>
00471 void TriMeshWithEdgeNormals<RealType>::getShapeOperatorFast( const FaceHandle &fh, aol::Matrix22<RealType> &S ) const {
00472     aol::Matrix22<RealType> temp;
00473     getFirstFundForm(fh,temp);
00474     S.makeInverse(temp);
00475     getSecondFundFormFast(fh,temp);
00476     S*=temp;
00477 }
00478 
00479 //=======================================================================================================
00480 //=======================================================================================================
00481 
00482 template <typename RealType>
00483 void TriMeshWithVertexNormals<RealType>::getVertexNormal ( const VertexHandle &vh, aol::Vec3<RealType> &wn, bool unit) const {
00484 
00485   wn.setZero();
00486   int faces = 0;
00487   for ( typename TriMesh<RealType>::ConstVertexFaceIter vf_it = this->getOpenMeshObject().cvf_iter ( vh ); vf_it; ++vf_it ) {
00488     // compute the barycenter
00489     aol::Vec3<RealType> coords[3], v[2], n;
00490     int i = 0;
00491     for ( typename TriMesh<RealType>::ConstFaceVertexIter iv_it = this->getOpenMeshObject().cfv_iter ( vf_it.handle() ); iv_it; ++iv_it ) {
00492       const typename TriMesh<RealType>::Point &p = this->getOpenMeshObject().point ( iv_it.handle() );
00493       for ( int j = 0; j < 3; j++ ) { coords[i][j] = p[j]; }
00494       ++i;
00495     }
00496     v[0] = coords[1]; v[0] -= coords[0];
00497     v[1] = coords[2]; v[1] -= coords[0];
00498     n = v[0].crossProduct ( v[1] );
00499     n.normalize();
00500     wn += n;
00501 
00502     ++faces;
00503   }
00504   
00505   if(unit)
00506     wn.normalize();
00507   else
00508     wn /= static_cast<RealType>(faces);
00509 }
00510 
00511 template <typename RealType>
00512 void TriMeshWithVertexNormals<RealType>::getDn(const FaceHandle &fh, aol::Mat<3,2,RealType> &Dn ) const {
00513 
00514   Dn.setZero();
00515   aol::Vec3<RealType> wn[3]; 
00516   int i =0;
00517   for( typename TriMesh<RealType>::ConstFaceVertexIter fv_it = this->getOpenMeshObject().cfv_iter(fh); fv_it; ++fv_it )
00518     getVertexNormal(fv_it.handle(),wn[i++]);
00519 
00520   // Dn = [n1-n0 | n2-n0]
00521   aol::Vec<2,RealType> row;
00522   for(int i =0; i<3; ++i){
00523     row[0] = wn[1][i] - wn[0][i];
00524     row[1] = wn[2][i] - wn[0][i];
00525     Dn[i]  = row;
00526   }
00527 }
00528 
00529 template <typename RealType>
00530 void TriMeshWithVertexNormals<RealType>::getSecondFundForm( const FaceHandle &fh, aol::Matrix22<RealType> &h ) const{
00531   aol::Mat<3,2,RealType> Dx, Dn;
00532   getDx(fh,Dx);
00533   getDn(fh,Dn);
00534   h.makeSymmetricProduct(Dn,Dx); 
00535 }
00536 
00537 template <typename RealType>
00538 void TriMeshWithVertexNormals<RealType>::noise ( RealType factor ) {
00539   for ( typename TriMesh<RealType>::NodeIteratorType vIter = *this; vIter.notAtEnd(); ++vIter ) {
00540     aol::Vec3<RealType> wn;
00541     getVertexNormal ( vIter.vertexHandle(), wn );
00542 
00543     typename TriMesh<RealType>::Point &p = this->point ( vIter.vertexHandle() );
00544     RealType r = ( static_cast<RealType> ( rand() ) / static_cast<RealType> ( RAND_MAX ) * 2. - 1. ) * factor;
00545     for ( int j = 0;j < 3;j++ ) {
00546       p[j] += r * wn[j];
00547     }
00548   }
00549 }
00550 //=======================================================================================================
00551 //=======================================================================================================
00552 
00553 template <typename RealType, typename Imp>
00554 void TriMeshWithGeomProps<RealType, Imp>::getCurvature( const FaceHandle &fh, aol::Vec2<RealType> &curv ) const{
00555   aol::Matrix22<RealType> shapeOp;
00556   getShapeOperator(fh, shapeOp );
00557   shapeOp.eigenValues(curv); 
00558 }
00559 
00560 template <typename RealType, typename Imp>
00561 void TriMeshWithGeomProps<RealType, Imp>::writePrincipalCurvaturesTo ( aol::MultiVector<RealType> & pc ) {
00562 
00563   if ( pc.numComponents() != 2 || ! pc.allDimsEqual ( this->getNumFaces() ) )
00564     throw aol::DimensionMismatchException ( "om::TrimMesh::writePrincipalCurvaturesTo(): "
00565                                             "passed MultiVector has not the right size.",
00566                                             __FILE__, __LINE__ );
00567 
00568   aol::Vec2<RealType> temp;
00569   for ( typename TriMesh<RealType>::ElementIteratorType fIter = *this; fIter.notAtEnd(); ++fIter ) {
00570 
00571     getCurvature( fIter.faceHandle(), temp );
00572 
00573     pc[0][fIter.getIndex()] = temp[0];
00574     pc[1][fIter.getIndex()] = temp[1];
00575   }
00576 }
00577 
00578 //! exports mesh with curvatures as vtk-file
00579 template <typename RealType, typename Imp>
00580 void TriMeshWithGeomProps<RealType, Imp>::exportCurvature( const std::string constfilename ){
00581 
00582   int n = this->getNumFaces();
00583   aol::MultiVector<RealType> curv ( 2, n );
00584   writePrincipalCurvaturesTo ( curv );
00585 
00586   aol::Vector<RealType> meanCurv  ( n ),
00587                         gaussCurv ( n ),
00588                         big20     ( n ),
00589                         big50     ( n ),
00590                         big100    ( n );
00591 
00592   for ( int i = 0; i < n; ++i ) {
00593     meanCurv[i]  = curv[0][i] + curv[1][i];
00594     gaussCurv[i] = curv[0][i] * curv[1][i];
00595     big20[i]     = std::max(curv[0][i], curv[1][i])> 20  ? 1 : 0;
00596     big50[i]     = std::max(curv[0][i], curv[1][i])> 50  ? 1 : 0;
00597     big100[i]    = std::max(curv[0][i], curv[1][i])> 100 ? 1 : 0;
00598   }
00599 
00600   // Does file end with ".vtk"?
00601   string filename = constfilename.substr( 0, constfilename.rfind(".") ) + ".vtk";
00602 
00603   aol::MeshWithData<TriMesh<RealType> > ( *this )
00604     .addData ( curv[0],   "k1", aol::FACE_DATA )
00605     .addData ( curv[1],   "k2", aol::FACE_DATA )
00606     .addData ( gaussCurv, "K",  aol::FACE_DATA )
00607     .addData ( meanCurv,  "H",  aol::FACE_DATA )
00608     .addData ( big20,     "big20",   aol::FACE_DATA )
00609     .addData ( big50,     "big50",   aol::FACE_DATA )
00610     .addData ( big100,    "big100",  aol::FACE_DATA )
00611     .saveAsLegacyVTK ( filename );
00612 }
00613 
00614 //=======================================================================================================
00615 //=======================================================================================================
00616 
00617 template class TriMesh<float>;
00618 template class TriMesh<double>;
00619 template class TriMesh<long double>;
00620 
00621 template class TriMeshWithVertexNormals<float>;
00622 template class TriMeshWithVertexNormals<double>;
00623 template class TriMeshWithVertexNormals<long double>;
00624 
00625 template class TriMeshWithEdgeNormals<float>;
00626 template class TriMeshWithEdgeNormals<double>;
00627 template class TriMeshWithEdgeNormals<long double>;
00628 
00629 template class TriMeshWithGeomProps< float,  TriMeshWithEdgeNormals<float> >;
00630 template class TriMeshWithGeomProps< float,  TriMeshWithVertexNormals<float> >;
00631 template class TriMeshWithGeomProps< double, TriMeshWithEdgeNormals<double> >;
00632 template class TriMeshWithGeomProps< double, TriMeshWithVertexNormals<double> >;
00633 
00634 } // end of namespace om.

Generated on Fri Sep 9 2011 21:08:05 for QuocMesh by doxygen 1.7.1