OpenWalnut  1.4.0
WTriangleMesh.cpp
00001 //---------------------------------------------------------------------------
00002 //
00003 // Project: OpenWalnut ( http://www.openwalnut.org )
00004 //
00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
00006 // For more information see http://www.openwalnut.org/copying
00007 //
00008 // This file is part of OpenWalnut.
00009 //
00010 // OpenWalnut is free software: you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as published by
00012 // the Free Software Foundation, either version 3 of the License, or
00013 // (at your option) any later version.
00014 //
00015 // OpenWalnut is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 // GNU Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public License
00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
00022 //
00023 //---------------------------------------------------------------------------
00024 
00025 #include <iostream>
00026 #include <list>
00027 #include <map>
00028 #include <set>
00029 #include <sstream>
00030 #include <stack>
00031 #include <string>
00032 #include <utility>
00033 #include <vector>
00034 #include <limits>
00035 
00036 #include <osg/io_utils>
00037 
00038 #include <Eigen/Eigenvalues>
00039 
00040 #include "../common/WAssert.h"
00041 #include "../common/WLogger.h"
00042 #include "../common/math/WMath.h"
00043 #include "../common/datastructures/WUnionFind.h"
00044 #include "WTriangleMesh.h"
00045 
00046 // init _static_ member variable and provide a linker reference to it
00047 boost::shared_ptr< WPrototyped > WTriangleMesh::m_prototype = boost::shared_ptr< WPrototyped >();
00048 
00049 boost::shared_ptr< WPrototyped > WTriangleMesh::getPrototype()
00050 {
00051     if( !m_prototype )
00052     {
00053          m_prototype = boost::shared_ptr< WPrototyped >( new WTriangleMesh( 0, 0 ) );
00054     }
00055     return m_prototype;
00056 }
00057 
00058 /**
00059  * constructor that already reserves space for a given number of triangles and vertexes
00060  */
00061 WTriangleMesh::WTriangleMesh( size_t vertNum, size_t triangleNum )
00062     : m_countVerts( 0 ),
00063       m_countTriangles( 0 ),
00064       m_meshDirty( true ),
00065       m_autoNormal( true ),
00066       m_neighborsCalculated( false ),
00067       m_curvatureCalculated( false )
00068 {
00069     m_verts = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( vertNum ) );
00070     m_textureCoordinates = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( vertNum ) );
00071     m_vertNormals = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( vertNum ) );
00072     m_vertColors = osg::ref_ptr< osg::Vec4Array >( new osg::Vec4Array( vertNum ) );
00073 
00074     m_triangles.resize( triangleNum * 3 );
00075     m_triangleNormals = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( triangleNum ) );
00076     m_triangleColors = osg::ref_ptr< osg::Vec4Array >( new osg::Vec4Array( triangleNum ) );
00077 }
00078 
00079 WTriangleMesh::WTriangleMesh( osg::ref_ptr< osg::Vec3Array > vertices, const std::vector< size_t >& triangles )
00080     : m_countVerts( vertices->size() ),
00081       m_countTriangles( triangles.size() / 3 ),
00082       m_meshDirty( true ),
00083       m_autoNormal( true ),
00084       m_neighborsCalculated( false ),
00085       m_verts( vertices ),
00086       m_textureCoordinates( new osg::Vec3Array( vertices->size() ) ),
00087       m_vertNormals( new osg::Vec3Array( vertices->size() ) ),
00088       m_vertColors( new osg::Vec4Array( vertices->size() ) ),
00089       m_triangles( triangles ),
00090       m_triangleNormals( new osg::Vec3Array( triangles.size() / 3 ) ),
00091       m_triangleColors( new osg::Vec4Array( triangles.size() / 3 ) )
00092 {
00093     WAssert( triangles.size() % 3 == 0, "Invalid triangle vector, having an invalid size (not divideable by 3)" );
00094 }
00095 
00096 WTriangleMesh::~WTriangleMesh()
00097 {
00098 }
00099 
00100 void WTriangleMesh::setAutoRecalcNormals( bool autoRecalc )
00101 {
00102     m_autoNormal = autoRecalc;
00103 }
00104 
00105 size_t WTriangleMesh::addVertex( float x, float y, float z )
00106 {
00107     return addVertex( osg::Vec3( x, y, z ) );
00108 }
00109 
00110 size_t WTriangleMesh::addVertex( WPosition vert )
00111 {
00112     return addVertex( osg::Vec3( vert[0], vert[1], vert[2] ) );
00113 }
00114 
00115 void WTriangleMesh::addTriangle( size_t vert0, size_t vert1, size_t vert2 )
00116 {
00117     if( m_triangles.size() == m_countTriangles * 3 )
00118     {
00119         m_triangles.resize( m_countTriangles * 3 + 3 );
00120     }
00121     m_triangles[ m_countTriangles * 3 ] = vert0;
00122     m_triangles[ m_countTriangles * 3 + 1 ] = vert1;
00123     m_triangles[ m_countTriangles * 3 + 2 ] = vert2;
00124     ++m_countTriangles;
00125 }
00126 
00127 void WTriangleMesh::addTriangle( osg::Vec3 vert0, osg::Vec3 vert1, osg::Vec3 vert2 )
00128 {
00129     addVertex( vert0 );
00130     addVertex( vert1 );
00131     addVertex( vert2 );
00132     addTriangle( m_countVerts - 3, m_countVerts - 2, m_countVerts - 1 );
00133 }
00134 
00135 void WTriangleMesh::setVertexNormal( size_t index, osg::Vec3 normal )
00136 {
00137     WAssert( index < m_countVerts, "set vertex normal: index out of range" );
00138     ( *m_vertNormals )[index] = normal;
00139 }
00140 
00141 void WTriangleMesh::setVertexNormal( size_t index, WPosition normal )
00142 {
00143     setVertexNormal( index, osg::Vec3( normal[0], normal[1], normal[2] ) );
00144 }
00145 
00146 void WTriangleMesh::setVertexNormal( size_t index, float x, float y, float z )
00147 {
00148     setVertexNormal( index, osg::Vec3( x, y, z ) );
00149 }
00150 
00151 void WTriangleMesh::setVertexColor( size_t index, osg::Vec4 color )
00152 {
00153     WAssert( index < m_countVerts, "set vertex color: index out of range" );
00154     ( *m_vertColors )[index] = color;
00155 }
00156 
00157 void WTriangleMesh::setTriangleColor( size_t index, osg::Vec4 color )
00158 {
00159     WAssert( index < m_countTriangles, "set triangle color: index out of range" );
00160     ( *m_triangleColors )[index] = color;
00161 }
00162 
00163 osg::ref_ptr< osg::Vec3Array >WTriangleMesh::getVertexArray()
00164 {
00165     return m_verts;
00166 }
00167 
00168 osg::ref_ptr< const osg::Vec3Array >WTriangleMesh::getVertexArray() const
00169 {
00170     return m_verts;
00171 }
00172 
00173 osg::ref_ptr< osg::Vec3Array > WTriangleMesh::getTextureCoordinateArray()
00174 {
00175     return m_textureCoordinates;
00176 }
00177 
00178 osg::ref_ptr< const osg::Vec3Array > WTriangleMesh::getTextureCoordinateArray() const
00179 {
00180     return m_textureCoordinates;
00181 }
00182 
00183 osg::ref_ptr< osg::Vec3Array >WTriangleMesh::getVertexNormalArray( bool forceRecalc )
00184 {
00185     if( forceRecalc || ( m_meshDirty && m_autoNormal ) )
00186     {
00187         recalcVertNormals();
00188     }
00189     return m_vertNormals;
00190 }
00191 
00192 osg::ref_ptr< osg::Vec3Array >WTriangleMesh::getTriangleNormalArray( bool forceRecalc )
00193 {
00194     if( forceRecalc || ( m_meshDirty && m_autoNormal ) )
00195     {
00196         recalcVertNormals();
00197     }
00198     return m_triangleNormals;
00199 }
00200 
00201 
00202 osg::ref_ptr< osg::Vec4Array >WTriangleMesh::getVertexColorArray()
00203 {
00204     return m_vertColors;
00205 }
00206 
00207 const std::vector< size_t >& WTriangleMesh::getTriangles() const
00208 {
00209     return m_triangles;
00210 }
00211 
00212 osg::Vec3 WTriangleMesh::getVertex( size_t index ) const
00213 {
00214     WAssert( index < m_countVerts, "get vertex: index out of range" );
00215     return ( *m_verts )[index];
00216 }
00217 
00218 osg::Vec4 WTriangleMesh::getVertColor( size_t index ) const
00219 {
00220     WAssert( index < m_countVerts, "get vertex color: index out of range" );
00221     return ( *m_vertColors )[index];
00222 }
00223 
00224 WVector3d WTriangleMesh::getNormal( size_t index )
00225 {
00226     WAssert( index < m_countVerts, "get normal as position: index out of range" );
00227     if( m_meshDirty )
00228     {
00229         recalcVertNormals();
00230     }
00231     return WPosition( ( *m_vertNormals )[index][0], ( *m_vertNormals )[index][1], ( *m_vertNormals )[index][2] );
00232 }
00233 
00234 void WTriangleMesh::removeVertex( size_t index )
00235 {
00236     WAssert( index < m_countVerts, "remove vertex: index out of range" );
00237     if( m_vertexIsInTriangle[ index ].size() > 0 )
00238     {
00239         return;
00240     }
00241     ( *m_verts ).erase( ( *m_verts ).begin() + index );
00242 
00243     for( size_t i = 0; i < m_countTriangles * 3; ++i )
00244     {
00245         if( m_triangles[ i ] > index )
00246         {
00247             --m_triangles[ i ];
00248         }
00249     }
00250     m_meshDirty = true;
00251 }
00252 
00253 void WTriangleMesh::removeTriangle( size_t index )
00254 {
00255     WAssert( index < m_countTriangles, "remove triangle: index out of range" );
00256     m_triangles.erase( m_triangles.begin() + index * 3, m_triangles.begin() + index * 3 + 3 );
00257     m_meshDirty = true;
00258 }
00259 
00260 void WTriangleMesh::recalcVertNormals()
00261 {
00262     updateVertsInTriangles();
00263 
00264     ( *m_vertNormals ).resize( m_countVerts );
00265     ( *m_triangleNormals ).resize( m_countTriangles );
00266 
00267     for( size_t i = 0; i < m_countTriangles; ++i )
00268     {
00269         ( *m_triangleNormals )[i] = calcTriangleNormal( i );
00270     }
00271 
00272     for( size_t vertId = 0; vertId < m_countVerts; ++vertId )
00273     {
00274         osg::Vec3 tempNormal( 0.0, 0.0, 0.0 );
00275         for( size_t neighbour = 0; neighbour < m_vertexIsInTriangle[vertId].size(); ++neighbour )
00276         {
00277             tempNormal += ( *m_triangleNormals )[ m_vertexIsInTriangle[vertId][neighbour] ];
00278         }
00279         tempNormal *= 1./m_vertexIsInTriangle[vertId].size();
00280 
00281         tempNormal.normalize();
00282         ( *m_vertNormals )[vertId] = tempNormal;
00283     }
00284 
00285     m_meshDirty = false;
00286 }
00287 
00288 void WTriangleMesh::updateVertsInTriangles()
00289 {
00290     m_vertexIsInTriangle.clear();
00291     std::vector< size_t >v;
00292     m_vertexIsInTriangle.resize( ( *m_verts ).size(), v );
00293 
00294     for( size_t i = 0; i < m_countTriangles; ++i )
00295     {
00296         m_vertexIsInTriangle[ getTriVertId0( i ) ].push_back( i );
00297         m_vertexIsInTriangle[ getTriVertId1( i ) ].push_back( i );
00298         m_vertexIsInTriangle[ getTriVertId2( i ) ].push_back( i );
00299     }
00300 }
00301 
00302 osg::Vec3 WTriangleMesh::calcTriangleNormal( size_t triangle )
00303 {
00304     osg::Vec3 v1( getTriVert( triangle, 1 ) - getTriVert( triangle, 0 ) );
00305     osg::Vec3 v2( getTriVert( triangle, 2 ) - getTriVert( triangle, 0 ) );
00306 
00307     osg::Vec3 tempNormal( 0, 0, 0 );
00308 
00309     tempNormal[0] = v1[1] * v2[2] - v1[2] * v2[1];
00310     tempNormal[1] = v1[2] * v2[0] - v1[0] * v2[2];
00311     tempNormal[2] = v1[0] * v2[1] - v1[1] * v2[0];
00312 
00313     tempNormal.normalize();
00314 
00315     return tempNormal;
00316 }
00317 
00318 osg::Vec3 WTriangleMesh::calcNormal( osg::Vec3 vert0, osg::Vec3 vert1, osg::Vec3 vert2 )
00319 {
00320     osg::Vec3 v1( vert1 - vert0 );
00321     osg::Vec3 v2( vert2 - vert0 );
00322 
00323     osg::Vec3 tempNormal( 0, 0, 0 );
00324 
00325     tempNormal[0] = v1[1] * v2[2] - v1[2] * v2[1];
00326     tempNormal[1] = v1[2] * v2[0] - v1[0] * v2[2];
00327     tempNormal[2] = v1[0] * v2[1] - v1[1] * v2[0];
00328 
00329     tempNormal.normalize();
00330 
00331     return tempNormal;
00332 }
00333 
00334 size_t WTriangleMesh::vertSize() const
00335 {
00336     return m_countVerts;
00337 }
00338 
00339 size_t WTriangleMesh::triangleSize() const
00340 {
00341     return m_countTriangles;
00342 }
00343 
00344 void WTriangleMesh::calcNeighbors()
00345 {
00346     std::vector<size_t> v( 3, -1 );
00347     m_triangleNeighbors.resize( ( *m_triangleNormals ).size(), v );
00348 
00349     for( size_t triId = 0; triId < m_countTriangles; ++triId )
00350     {
00351         size_t coVert0 = getTriVertId0( triId );
00352         size_t coVert1 = getTriVertId1( triId );
00353         size_t coVert2 = getTriVertId2( triId );
00354 
00355         m_triangleNeighbors[triId][0] = getNeighbor( coVert0, coVert1, triId );
00356         m_triangleNeighbors[triId][1] = getNeighbor( coVert1, coVert2, triId );
00357         m_triangleNeighbors[triId][2] = getNeighbor( coVert2, coVert0, triId );
00358     }
00359     m_neighborsCalculated = true;
00360 }
00361 
00362 size_t WTriangleMesh::getNeighbor( const size_t coVert1, const size_t coVert2, const size_t triangleNum )
00363 {
00364     std::vector< size_t > candidates = m_vertexIsInTriangle[coVert1];
00365     std::vector< size_t > compares = m_vertexIsInTriangle[coVert2];
00366 
00367     for( size_t i = 0; i < candidates.size(); ++i )
00368     {
00369         for( size_t k = 0; k < compares.size(); ++k )
00370         {
00371             if( ( candidates[i] != triangleNum ) && ( candidates[i] == compares[k] ) )
00372             {
00373                 return candidates[i];
00374             }
00375         }
00376     }
00377     return triangleNum;
00378 }
00379 
00380 void WTriangleMesh::doLoopSubD()
00381 {
00382     m_numTriVerts = m_countVerts;
00383     m_numTriFaces = m_countTriangles;
00384 
00385     ( *m_verts ).resize( m_numTriVerts * 4 );
00386     m_triangles.resize( m_numTriFaces * 4 * 3 );
00387 
00388     updateVertsInTriangles();
00389 
00390     osg::Vec3* newVertexPositions = new osg::Vec3[m_numTriVerts];
00391 
00392     //std::cout << "Loop subdivision pass 1" << std::endl;
00393     for( size_t i = 0; i < m_numTriVerts; ++i )
00394     {
00395         newVertexPositions[i] = loopCalcNewPosition( i );
00396     }
00397 
00398     //std::cout << "Loop subdivision pass 2" << std::endl;
00399     for( size_t i = 0; i < m_numTriFaces; ++i )
00400     {
00401         loopInsertCenterTriangle( i );
00402     }
00403     ( *m_verts ).resize( m_countVerts );
00404     std::vector< size_t >v;
00405     m_vertexIsInTriangle.resize( ( *m_verts ).size(), v );
00406 
00407     //std::cout << "Loop subdivision pass 3" << std::endl;
00408     for( size_t i = 0; i < m_numTriFaces; ++i )
00409     {
00410         loopInsertCornerTriangles( i );
00411     }
00412 
00413     //std::cout << "Loop subdivision pass 4" << std::endl;
00414     for( size_t i = 0; i < m_numTriVerts; ++i )
00415     {
00416         ( *m_verts )[i] = newVertexPositions[i];
00417     }
00418 
00419     delete[] newVertexPositions;
00420 
00421     m_vertNormals->resize( m_verts->size() );
00422     m_vertColors->resize( m_verts->size() );
00423     m_triangleColors->resize( m_triangles.size() / 3 );
00424 
00425     m_meshDirty = true;
00426 }
00427 
00428 
00429 osg::Vec3 WTriangleMesh::loopCalcNewPosition( size_t vertId )
00430 {
00431     std::vector< size_t > starP = m_vertexIsInTriangle[vertId];
00432     int starSize = starP.size();
00433 
00434     osg::Vec3 oldPos = getVertex( vertId );
00435     double alpha = loopGetAlpha( starSize );
00436 
00437     double scale = 1.0 - ( static_cast<double>( starSize ) * alpha );
00438     oldPos *= scale;
00439 
00440     osg::Vec3 newPos;
00441     int edgeV = 0;
00442     for( int i = 0; i < starSize; i++ )
00443     {
00444         edgeV = loopGetNextVertex( starP[i], vertId );
00445         osg::Vec3 translate = getVertex( edgeV );
00446         newPos += translate;
00447     }
00448     newPos *= alpha;
00449 
00450     return oldPos + newPos;
00451 }
00452 
00453 void WTriangleMesh::loopInsertCenterTriangle( size_t triId )
00454 {
00455     size_t edgeVerts[3];
00456 
00457     edgeVerts[0] = loopCalcEdgeVert( triId, getTriVertId0( triId ), getTriVertId1( triId ), getTriVertId2( triId ) );
00458     edgeVerts[1] = loopCalcEdgeVert( triId, getTriVertId1( triId ), getTriVertId2( triId ), getTriVertId0( triId ) );
00459     edgeVerts[2] = loopCalcEdgeVert( triId, getTriVertId2( triId ), getTriVertId0( triId ), getTriVertId1( triId ) );
00460 
00461     addTriangle( edgeVerts[0], edgeVerts[1], edgeVerts[2] );
00462 }
00463 
00464 
00465 size_t WTriangleMesh::loopCalcEdgeVert( size_t triId, size_t edgeV1, size_t edgeV2, size_t V3 )
00466 {
00467     size_t neighborVert = -1;
00468     size_t neighborFaceNum = -1;
00469     osg::Vec3 edgeVert;
00470 
00471     neighborFaceNum = getNeighbor( edgeV1, edgeV2, triId );
00472 
00473     if( neighborFaceNum == triId )
00474     {
00475         osg::Vec3 edgeVert = ( ( *m_verts )[edgeV1] + ( *m_verts )[edgeV2] ) / 2.0;
00476         size_t vertId = m_countVerts;
00477         addVertex( edgeVert );
00478         return vertId;
00479     }
00480 
00481     else if( neighborFaceNum > triId )
00482     {
00483         neighborVert = loopGetThirdVert( edgeV1, edgeV2, neighborFaceNum );
00484 
00485         osg::Vec3 edgePart = ( *m_verts )[edgeV1] + ( *m_verts )[edgeV2];
00486         osg::Vec3 neighborPart = ( *m_verts )[neighborVert] + ( *m_verts )[V3];
00487 
00488         edgeVert = ( ( edgePart * ( 3.0 / 8.0 ) ) + ( neighborPart * ( 1.0 / 8.0 ) ) );
00489         size_t vertId = m_countVerts;
00490         addVertex( edgeVert );
00491         return vertId;
00492     }
00493     else
00494     {
00495         size_t neighborCenterP = neighborFaceNum + m_numTriFaces;
00496         size_t neighborP = neighborFaceNum;
00497 
00498         if( getTriVertId0( neighborP ) == edgeV2 )
00499         {
00500             return getTriVertId0( neighborCenterP );
00501         }
00502         else if( getTriVertId1( neighborP ) == edgeV2 )
00503         {
00504             return getTriVertId1( neighborCenterP );
00505         }
00506         else
00507         {
00508             return getTriVertId2( neighborCenterP );
00509         }
00510     }
00511     return -1;
00512 }
00513 
00514 void WTriangleMesh::loopInsertCornerTriangles( size_t triId )
00515 {
00516     // comment:             center are twisted from the orignal vertices.
00517     // original:    0, 1, 2
00518     // center:              a, b, c
00519     // reAsgnOrig:  0, a, c
00520     // addTris:             1, b, a
00521     // addTris:             2, c, b
00522     //
00523     size_t originalTri0 = getTriVertId0( triId );
00524     size_t originalTri1 = getTriVertId1( triId );
00525     size_t originalTri2 = getTriVertId2( triId );
00526 
00527     size_t centerTri0 = getTriVertId0( triId + m_numTriFaces );
00528     size_t centerTri1 = getTriVertId1( triId + m_numTriFaces );
00529     size_t centerTri2 = getTriVertId2( triId + m_numTriFaces );
00530 
00531     addTriangle( originalTri1, centerTri1, centerTri0 );
00532     addTriangle( originalTri2, centerTri2, centerTri1 );
00533     loopSetTriangle( triId, originalTri0, centerTri0, centerTri2 );
00534 }
00535 
00536 void WTriangleMesh::loopSetTriangle( size_t triId, size_t vertId1, size_t vertId2, size_t vertId3 )
00537 {
00538     loopEraseTriangleFromVertex( triId, getTriVertId1( triId ) );
00539     loopEraseTriangleFromVertex( triId, getTriVertId2( triId ) );
00540 
00541     setTriVert0( triId, vertId1 );
00542     setTriVert1( triId, vertId2 );
00543     setTriVert2( triId, vertId3 );
00544 
00545     m_vertexIsInTriangle[vertId2].push_back( triId );
00546     m_vertexIsInTriangle[vertId3].push_back( triId );
00547 }
00548 
00549 void WTriangleMesh::loopEraseTriangleFromVertex( size_t triId, size_t vertId )
00550 {
00551     std::vector< size_t > temp;
00552     for( size_t i = 0; i < m_vertexIsInTriangle[vertId].size(); ++i )
00553     {
00554         if( triId != m_vertexIsInTriangle[vertId][i] )
00555             temp.push_back( m_vertexIsInTriangle[vertId][i] );
00556     }
00557     m_vertexIsInTriangle[vertId] = temp;
00558 }
00559 
00560 double WTriangleMesh::loopGetAlpha( int n )
00561 {
00562     double answer;
00563     if( n > 3 )
00564     {
00565         double center = ( 0.375 + ( 0.25 * cos( ( 2.0 * 3.14159265358979 ) / static_cast<double>( n ) ) ) );
00566         answer = ( 0.625 - ( center * center ) ) / static_cast<double>( n );
00567     }
00568     else
00569     {
00570         answer = 3.0 / 16.0;
00571     }
00572     return answer;
00573 }
00574 
00575 size_t WTriangleMesh::loopGetNextVertex( size_t triNum, size_t vertNum )
00576 {
00577     if( getTriVertId0( triNum ) == vertNum )
00578     {
00579         return getTriVertId1( triNum );
00580     }
00581     else if( getTriVertId1( triNum ) == vertNum )
00582     {
00583         return getTriVertId2( triNum );
00584     }
00585     return getTriVertId0( triNum );
00586 }
00587 
00588 size_t WTriangleMesh::loopGetThirdVert( size_t coVert1, size_t coVert2, size_t triangleNum )
00589 {
00590     if( !( getTriVertId0( triangleNum ) == coVert1 ) && !( getTriVertId0( triangleNum ) == coVert2 ) )
00591     {
00592         return getTriVertId0( triangleNum );
00593     }
00594     else if( !( getTriVertId1( triangleNum ) == coVert1 ) && !( getTriVertId1( triangleNum ) == coVert2 ) )
00595     {
00596         return getTriVertId1( triangleNum );
00597     }
00598     return getTriVertId2( triangleNum );
00599 }
00600 
00601 void WTriangleMesh::addMesh( boost::shared_ptr<WTriangleMesh> mesh, float xOff, float yOff, float zOff )
00602 {
00603     size_t oldVertSize = m_countVerts;
00604 
00605     ( *m_vertColors ).resize( oldVertSize + mesh->vertSize() );
00606     for( size_t i = 0; i < mesh->vertSize(); ++i )
00607     {
00608         osg::Vec3 v( mesh->getVertex( i ) );
00609         v[0] += xOff;
00610         v[1] += yOff;
00611         v[2] += zOff;
00612         addVertex( v );
00613         setVertexColor( oldVertSize + i, mesh->getVertColor( i ) );
00614     }
00615     for( size_t i = 0; i < mesh->triangleSize(); ++i )
00616     {
00617         addTriangle( mesh->getTriVertId0( i ) + oldVertSize, mesh->getTriVertId1( i ) + oldVertSize, mesh->getTriVertId2( i ) + oldVertSize );
00618     }
00619     m_meshDirty = true;
00620 }
00621 
00622 void WTriangleMesh::translateMesh( float xOff, float yOff, float zOff )
00623 {
00624     osg::Vec3 t( xOff, yOff, zOff );
00625     for( size_t i = 0; i < ( *m_verts ).size(); ++i )
00626     {
00627         ( *m_verts )[i] += t;
00628     }
00629 }
00630 
00631 void WTriangleMesh::zoomMesh( float zoom )
00632 {
00633     for( size_t i = 0; i < ( *m_verts ).size(); ++i )
00634     {
00635         ( *m_verts )[i] *= zoom;
00636     }
00637 }
00638 
00639 void WTriangleMesh::rescaleVertexColors()
00640 {
00641     float maxR = 0;
00642     float maxG = 0;
00643     float maxB = 0;
00644     for( size_t vertId = 0; vertId < m_vertColors->size(); ++vertId )
00645     {
00646         if( ( *m_vertColors )[vertId][0] > maxR )
00647         {
00648             maxR = ( *m_vertColors )[vertId][0];
00649         }
00650         if( ( *m_vertColors )[vertId][1] > maxG )
00651         {
00652             maxG = ( *m_vertColors )[vertId][1];
00653         }
00654         if( ( *m_vertColors )[vertId][2] > maxB )
00655         {
00656             maxB = ( *m_vertColors )[vertId][2];
00657         }
00658     }
00659     for( size_t vertId = 0; vertId < m_vertColors->size(); ++vertId )
00660     {
00661         ( *m_vertColors )[vertId][0] /= maxR;
00662         ( *m_vertColors )[vertId][1] /= maxG;
00663         ( *m_vertColors )[vertId][2] /= maxB;
00664     }
00665 }
00666 
00667 std::ostream& tm_utils::operator<<( std::ostream& os, const WTriangleMesh& rhs )
00668 {
00669     std::stringstream ss;
00670     ss << "WTriangleMesh( #vertices=" << rhs.vertSize() << " #triangles=" << rhs.triangleSize() << " )" << std::endl;
00671     using string_utils::operator<<;
00672     size_t count = 0;
00673     ss << std::endl;
00674     const std::vector< size_t >& triangles = rhs.getTriangles();
00675     osg::ref_ptr< const osg::Vec3Array > vertices = rhs.getVertexArray();
00676     for( size_t vID = 0 ; vID <= triangles.size() - 3; ++count )
00677     {
00678         std::stringstream prefix;
00679         prefix << "triangle: " << count << "[ ";
00680         std::string indent( prefix.str().size(), ' ' );
00681         using osg::operator<<; // using operator<< as defined in osg/io_utils
00682         ss << prefix.str() << vertices->at( triangles[ vID++ ] ) << std::endl;
00683         ss << indent << vertices->at( triangles[ vID++ ] ) << std::endl;
00684         ss << indent << vertices->at( triangles[ vID++ ] ) << std::endl;
00685         ss << std::string( indent.size() - 2, ' ' ) << "]" << std::endl;
00686     }
00687     return os << ss.str();
00688 }
00689 
00690 boost::shared_ptr< std::list< boost::shared_ptr< WTriangleMesh > > > tm_utils::componentDecomposition( const WTriangleMesh& mesh )
00691 {
00692     boost::shared_ptr< std::list< boost::shared_ptr< WTriangleMesh > > > result( new std::list< boost::shared_ptr< WTriangleMesh > >() );
00693     if( mesh.vertSize() <= 0 ) // no component possible
00694     {
00695         return result;
00696     }
00697     if( mesh.triangleSize() < 3 )
00698     {
00699         if( mesh.vertSize() > 0 )
00700         {
00701             // there are vertices but no triangles
00702             WAssert( false, "Not implemented the decomposition of a TriangleMesh without any triangles" );
00703         }
00704         else // no component possible
00705         {
00706             return result;
00707         }
00708     }
00709 
00710     WUnionFind uf( mesh.vertSize() ); // idea: every vertex in own component, then successivley join in accordance with the triangles
00711 
00712     const std::vector< size_t >& triangles = mesh.getTriangles();
00713     for( size_t vID = 0; vID <= triangles.size() - 3; vID += 3)
00714     {
00715         uf.merge( triangles[ vID ], triangles[ vID + 1 ] );
00716         uf.merge( triangles[ vID ], triangles[ vID + 2 ] ); // uf.merge( triangles[ vID + 2 ], triangles[ vID + 1 ] ); they are already in same
00717     }
00718 
00719     // ATTENTION: The reason for using the complex BucketType instead of pasting vertices directly into a new WTriangleMesh
00720     // is performance! For example: If there are many vertices reused inside the former WTriangleMesh mesh, then we want
00721     // to reuse them in the new components too. Hence we must determine if a certain vertex is already inside the new component.
00722     // Since the vertices are organized in a vector, we can use std::find( v.begin, v.end(), vertexToLookUp ) which results
00723     // in O(N^2) or we could use faster lookUp via key and value leading to the map and the somehow complicated BucketType.
00724     typedef std::map< osg::Vec3, size_t > VertexType; // look up fast if a vertex is already inside the new mesh!
00725     typedef std::vector< size_t > TriangleType;
00726     typedef std::pair< VertexType, TriangleType > BucketType; // Later on the Bucket will be transformed into the new WTriangleMesh component
00727     std::map< size_t, BucketType > buckets; // Key identify with the cannonical element from UnionFind the new connected component
00728 
00729     osg::ref_ptr< const osg::Vec3Array > vertices = mesh.getVertexArray();
00730     for( size_t vID = 0; vID <= triangles.size() - 3; vID += 3 )
00731     {
00732         size_t component = uf.find( triangles[ vID ] );
00733         if( buckets.find( component ) == buckets.end() )
00734         {
00735             buckets[ component ] = BucketType( VertexType(), TriangleType() ); // create new bucket
00736         }
00737 
00738         // Note: We discard the order of the points and indices, but semantically the structure remains the same
00739         VertexType& mapRef = buckets[ component ].first; // short hand alias
00740         for( int i = 0; i < 3; ++i )
00741         {
00742             size_t id = 0;
00743             const osg::Vec3& vertex = ( *vertices )[ triangles[ vID + i ] ];
00744             if( mapRef.find( vertex ) == mapRef.end() )
00745             {
00746                 id = mapRef.size(); // since size might change in next line
00747                 mapRef[ vertex ] = id;
00748             }
00749             else
00750             {
00751                 id = mapRef[ vertex ];
00752             }
00753             buckets[ component ].second.push_back( id );
00754         }
00755     }
00756 
00757     for( std::map< size_t, BucketType >::const_iterator cit = buckets.begin(); cit != buckets.end(); ++cit )
00758     {
00759         osg::ref_ptr< osg::Vec3Array > newVertices( new osg::Vec3Array );
00760         newVertices->resize( cit->second.first.size() );
00761         for( VertexType::const_iterator vit = cit->second.first.begin(); vit != cit->second.first.end(); ++vit )
00762         {
00763             newVertices->at( vit->second ) = vit->first; // if you are sure that vit->second is always valid replace at() call with operator[]
00764         }
00765         boost::shared_ptr< WTriangleMesh > newMesh( new WTriangleMesh( newVertices, cit->second.second ) );
00766         result->push_back( newMesh );
00767     }
00768 
00769     return result;
00770 }
00771 
00772 osg::ref_ptr< osg::Vec4Array > WTriangleMesh::getTriangleColors() const
00773 {
00774     return m_triangleColors;
00775 }
00776 
00777 void WTriangleMesh::performFeaturePreservingSmoothing( float sigmaDistance, float sigmaInfluence )
00778 {
00779     updateVertsInTriangles();
00780 
00781     calcNeighbors();
00782 
00783     // we perform a first smoothing pass and write the resulting vertex coords into a buffer
00784     // this will only update the normals
00785     performFeaturePreservingSmoothingMollificationPass( sigmaDistance / 2.0f, sigmaInfluence );
00786 
00787     // using the smoothed normals, we now perform a second smoothing pass, this time writing the new vertex coords
00788     performFeaturePreservingSmoothingVertexPass( sigmaDistance, sigmaInfluence );
00789 }
00790 
00791 void WTriangleMesh::performFeaturePreservingSmoothingMollificationPass( float sigmaDistance, float sigmaInfluence )
00792 {
00793     // calc Eq. 3 for every triangle
00794     osg::ref_ptr< osg::Vec3Array > vtxArray = new osg::Vec3Array( m_verts->size() );
00795 
00796     for( std::size_t k = 0; k < m_verts->size(); ++k )
00797     {
00798         vtxArray->operator[] ( k ) = estimateSmoothedVertexPosition( k, sigmaDistance, sigmaInfluence, true );
00799     }
00800 
00801     // calc the new normal directions - update triangle normals
00802     for( std::size_t k = 0; k < m_triangles.size() / 3; ++k )
00803     {
00804         osg::Vec3 const& p0 = vtxArray->operator[]( m_triangles[ 3 * k + 0 ] );
00805         osg::Vec3 const& p1 = vtxArray->operator[]( m_triangles[ 3 * k + 1 ] );
00806         osg::Vec3 const& p2 = vtxArray->operator[]( m_triangles[ 3 * k + 2 ] );
00807 
00808         m_triangleNormals->operator[] ( k ) = ( p1 - p0 ) ^ ( p2 - p1 );
00809         m_triangleNormals->operator[] ( k ).normalize();
00810     }
00811 }
00812 
00813 void WTriangleMesh::performFeaturePreservingSmoothingVertexPass( float sigmaDistance, float sigmaInfluence )
00814 {
00815     for( std::size_t k = 0; k < m_verts->size(); ++k )
00816     {
00817         m_verts->operator[] ( k ) = estimateSmoothedVertexPosition( k, sigmaDistance, sigmaInfluence, false );
00818     }
00819 
00820     recalcVertNormals();
00821 }
00822 
00823 osg::Vec3 WTriangleMesh::estimateSmoothedVertexPosition( std::size_t vtx, float sigmaDistance, float sigmaInfluence, bool mollify )
00824 {
00825     std::stack< std::size_t > triStack;
00826     std::set< std::size_t > triSet;
00827 
00828     for( std::size_t k = 0; k < m_vertexIsInTriangle[ vtx ].size(); ++k )
00829     {
00830         triStack.push( m_vertexIsInTriangle[ vtx ][ k ] );
00831         triSet.insert( m_vertexIsInTriangle[ vtx ][ k ] );
00832     }
00833 
00834     while( !triStack.empty() )
00835     {
00836         std::size_t currentTriangle = triStack.top();
00837         triStack.pop();
00838 
00839         for( std::size_t k = 0; k < m_triangleNeighbors[ currentTriangle ].size(); ++k )
00840         {
00841             osg::Vec3 center = calcTriangleCenter( m_triangleNeighbors[ currentTriangle ][ k ] );
00842 
00843             if( ( m_verts->operator[] ( vtx ) - center ).length() > 4.0 * sigmaDistance )
00844             {
00845                 continue;
00846             }
00847 
00848             if( triSet.find( m_triangleNeighbors[ currentTriangle ][ k ] ) == triSet.end() )
00849             {
00850                 triStack.push( m_triangleNeighbors[ currentTriangle ][ k ] );
00851                 triSet.insert( m_triangleNeighbors[ currentTriangle ][ k ] );
00852             }
00853         }
00854     }
00855 
00856     double sum = 0.0;
00857     osg::Vec3 res( 0.0, 0.0, 0.0 );
00858 
00859     for( std::set< std::size_t >::const_iterator it = triSet.begin(); it != triSet.end(); ++it )
00860     {
00861         osg::Vec3 center = calcTriangleCenter( *it );
00862         double area = calcTriangleArea( *it );
00863 
00864         // calc f
00865         double dist = ( m_verts->operator[] ( vtx ) - center ).length();
00866         double f = 1.0 / ( sigmaDistance * sqrt( 2.0 * piDouble ) ) * exp( -0.5 * dist * dist );
00867 
00868         double g;
00869         if( !mollify )
00870         {
00871             // calc prediction
00872             osg::Vec3f const& p = m_verts->operator[] ( vtx );
00873             osg::Vec3f const& n = m_triangleNormals->operator[] ( *it );
00874             osg::Vec3f pred = p + n * ( n * ( center - p ) );
00875 
00876             dist = ( p - pred ).length();
00877         }
00878         g = 1.0 / ( sigmaInfluence * sqrt( 2.0 * piDouble ) ) * exp( -0.5 * dist * dist );
00879 
00880         sum += area * f * g;
00881         res += center * area * f * g;
00882     }
00883 
00884     res /= sum;
00885     return res;
00886 }
00887 
00888 osg::Vec3 WTriangleMesh::calcTriangleCenter( std::size_t triIdx ) const
00889 {
00890     osg::Vec3 res = m_verts->operator[] ( m_triangles[ 3 * triIdx + 0 ] );
00891     res += m_verts->operator[] ( m_triangles[ 3 * triIdx + 1 ] );
00892     res += m_verts->operator[] ( m_triangles[ 3 * triIdx + 2 ] );
00893 
00894     res /= 3.0f;
00895     return res;
00896 }
00897 
00898 float WTriangleMesh::calcTriangleArea( std::size_t triIdx ) const
00899 {
00900     osg::Vec3 const& p0 = m_verts->operator[] ( m_triangles[ 3 * triIdx + 0 ] );
00901     osg::Vec3 const& p1 = m_verts->operator[] ( m_triangles[ 3 * triIdx + 1 ] );
00902     osg::Vec3 const& p2 = m_verts->operator[] ( m_triangles[ 3 * triIdx + 2 ] );
00903 
00904     return ( ( p1 - p0 ) ^ ( p2 - p0 ) ).length() * 0.5;
00905 }
00906 
00907 void WTriangleMesh::estimateCurvature()
00908 {
00909     updateVertsInTriangles();
00910     calcNeighbors();
00911 
00912     std::vector< osg::Vec3 > normals( m_verts->size() );
00913 
00914     // init vectors
00915     m_mainNormalCurvature = boost::shared_ptr< std::vector< float > >( new std::vector< float >( m_verts->size() ) );
00916     m_secondaryNormalCurvature = boost::shared_ptr< std::vector< float > >( new std::vector< float >( m_verts->size() ) );
00917     m_mainCurvaturePrincipalDirection = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( m_verts->size() ) );
00918     m_secondaryCurvaturePrincipalDirection = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( m_verts->size() ) );
00919 
00920     // calculate vertex normals using distance-weighted summing of neighbor-triangle normals
00921     for( std::size_t vtxId = 0; vtxId < m_verts->size(); ++vtxId )
00922     {
00923         osg::Vec3 const& p = m_verts->operator[] ( vtxId );
00924         osg::Vec3 n( 0.0, 0.0, 0.0 );
00925 
00926         for( std::size_t k = 0; k < m_vertexIsInTriangle[ vtxId ].size(); ++k )
00927         {
00928             std::size_t triId = m_vertexIsInTriangle[ vtxId ][ k ];
00929 
00930             osg::Vec3 center = calcTriangleCenter( triId );
00931             double w = 1.0 / ( center - p ).length();
00932 
00933             n += m_triangleNormals->operator[] ( triId ) * w;
00934         }
00935 
00936         WAssert( n.length() > 0.0001, "Invalid normal!" );
00937 
00938         n.normalize();
00939         normals[ vtxId ] = n;
00940     }
00941 
00942     // calculate curvatures for every vertex
00943     for( std::size_t vtxId = 0; vtxId < m_verts->size(); ++vtxId )
00944     {
00945         osg::Vec3 const& p = m_verts->operator[] ( vtxId );
00946 
00947         osg::Vec3 const& normal = normals[ vtxId ];
00948 
00949         // get the set of neighbor vertices
00950         std::set< std::size_t > neighbors;
00951 
00952         for( std::size_t k = 0; k < m_vertexIsInTriangle[ vtxId ].size(); ++k )
00953         {
00954             std::size_t triId = m_vertexIsInTriangle[ vtxId ][ k ];
00955 
00956             for( std::size_t j = 0; j < 3; ++j )
00957             {
00958                 std::size_t e = m_triangles[ 3 * triId + j ];
00959 
00960                 if( neighbors.find( e ) == neighbors.end() && e != vtxId )
00961                 {
00962                     neighbors.insert( e );
00963                 }
00964             }
00965         }
00966 
00967         WAssert( neighbors.size() > 2, "Vertex has too few neighbors! Does this mesh have holes?" );
00968 
00969         double maxCurvature = -std::numeric_limits< double >::infinity();
00970         std::vector< double > curvatures;
00971 
00972         osg::Vec3 maxCurvatureTangent( 0.0, 0.0, 0.0 );
00973         std::vector< osg::Vec3 > tangents;
00974 
00975         // part 1: get curvatures at tangents and their maximum curvature
00976         for( std::set< std::size_t >::const_iterator it = neighbors.begin(); it != neighbors.end(); ++it )
00977         {
00978             osg::Vec3 const& neighbPos = m_verts->operator[] ( *it );
00979             osg::Vec3 const& neighbNormal = normals[ *it ];
00980 
00981             // project ( neighbPos - p ) onto the tangent plane
00982             osg::Vec3 tangent = ( neighbPos - p ) - normal * ( ( neighbPos - p ) * normal );
00983             tangent.normalize();
00984 
00985             // approximate normal curvature in tangent direction
00986             double ncurv = -1.0 * ( ( neighbPos - p ) * ( neighbNormal - normal ) ) / ( ( neighbPos - p ) * ( neighbPos - p ) );
00987 
00988             if( ncurv > maxCurvature )
00989             {
00990                 maxCurvature = ncurv;
00991                 maxCurvatureTangent = tangent;
00992             }
00993 
00994             tangents.push_back( tangent );
00995             curvatures.push_back( ncurv );
00996         }
00997 
00998         WAssert( maxCurvatureTangent.length() > 0.0001, "Invalid tangent length!" );
00999 
01000         // part 2: choose a coordinate system in the tangent plane
01001         osg::Vec3 const& e1 = maxCurvatureTangent;
01002         osg::Vec3 e2 = e1 ^ normal;
01003 
01004         e2.normalize();
01005 
01006         bool significantCurvature = false;
01007         for( std::vector< double >::const_iterator it = curvatures.begin(); it != curvatures.end(); ++it )
01008         {
01009             if( fabs( *it ) > 0.00001 )
01010             {
01011                 significantCurvature = true;
01012                 break;
01013             }
01014         }
01015 
01016         if( !significantCurvature )
01017         {
01018             // curvatures were all almost zero
01019             // the mesh is flat at this point, write values accordingly
01020             m_mainNormalCurvature->operator[] ( vtxId ) = 0.0;
01021             m_mainCurvaturePrincipalDirection->operator[] ( vtxId ) = e1;
01022             m_secondaryNormalCurvature->operator[] ( vtxId ) = 0.0;
01023             m_secondaryCurvaturePrincipalDirection->operator[] ( vtxId ) = e2;
01024 
01025             continue;
01026         }
01027 
01028         // calculate coefficients of the ellipse a * cos²(theta) + b * cos(theta) * sin(theta) * c * sin²(theta)
01029         // this is done by estimating a as the largest curvature amoung the estimated curvatures in all tangent
01030         // direction belonging to points that share a triangle with the current one
01031         double const& a = maxCurvature;
01032 
01033         Eigen::Matrix< double, -1, -1 > X( tangents.size(), 2 );
01034         Eigen::Matrix< double, -1, -1 > y( tangents.size(), 1 );
01035 
01036         for( std::size_t k = 0; k < tangents.size(); ++k )
01037         {
01038             double theta = calcAngleBetweenNormalizedVectors( tangents[ k ], e1 );
01039 
01040             X( k, 0 ) = cos( theta ) * sin( theta );
01041             X( k, 1 ) = sin( theta ) * sin( theta );
01042             y( k, 0 ) = curvatures[ k ] - a * cos( theta ) * cos( theta );
01043         }
01044 
01045         // use LU decomposition to calculate rank
01046         Eigen::FullPivLU< Eigen::Matrix< double, -1, -1 > > lu( X );
01047 
01048         // we need a rank of at least 2 to calculate the coeffs
01049         if( lu.rank() < 2 )
01050         {
01051             wlog::warn( "WTriangleMesh::estimateCurvature" ) << "Rank too low, cannot estimate curvature for this vertex!";
01052 
01053             m_mainNormalCurvature->operator[] ( vtxId ) = a;
01054             m_mainCurvaturePrincipalDirection->operator[] ( vtxId ) = e1;
01055             m_secondaryNormalCurvature->operator[] ( vtxId ) = 0.0;
01056             m_secondaryCurvaturePrincipalDirection->operator[] ( vtxId ) = e2;
01057 
01058             continue;
01059         }
01060 
01061         // do least squares
01062         Eigen::Matrix< double, -1, -1 > bb = ( X.transpose() * X ).inverse() * X.transpose() * y;
01063 
01064         // the missing coeffs b and c are now:
01065         double b = bb( 0, 0 );
01066         double c = bb( 1, 0 );
01067 
01068         // this calculates the maximum and minimum normal curvatures
01069         // (as eigenvalues)
01070         double Kg = a * c - 0.25 * b * b;
01071         double H = 0.5 * ( a + c );
01072         double s = H * H - Kg;
01073 
01074         if( s < 0.0 )
01075         {
01076             s = 0.0;
01077         }
01078 
01079         double k1 = H + sqrt( s );
01080         double k2 = H - sqrt( s );
01081 
01082         if( fabs( k1 - k2 ) < 0.000001 )
01083         {
01084             // if the curvatures are equal, there is no single principal direction
01085             m_mainNormalCurvature->operator[] ( vtxId ) = a;
01086             m_mainCurvaturePrincipalDirection->operator[] ( vtxId ) = e1;
01087         }
01088         else
01089         {
01090             // if the curvatures differ, we can now find the direction of maximum curvature
01091             double temp = b / ( k2 - k1 );
01092 
01093             if( temp > 1.0 )
01094             {
01095                 temp = 1.0;
01096             }
01097             if( temp < -1.0 )
01098             {
01099                 temp = -1.0;
01100             }
01101 
01102             double theta = 0.5 * asin( temp );
01103 
01104             osg::Vec3 ne1 = e1 * cos( theta ) + e2 * sin( theta );
01105             osg::Vec3 ne2 = e2 * cos( theta ) - e1 * sin( theta );
01106 
01107             ne1.normalize();
01108             ne2.normalize();
01109 
01110             e2 = ne2;
01111 
01112             theta = calcAngleBetweenNormalizedVectors( ne1, e1 );
01113 
01114             m_mainNormalCurvature->operator[] ( vtxId ) = a * cos( theta ) * cos( theta )
01115                                                         + b * cos( theta ) * sin( theta )
01116                                                         + c * sin( theta ) * sin( theta );
01117             m_mainCurvaturePrincipalDirection->operator[] ( vtxId ) = ne1;
01118         }
01119 
01120         double theta = calcAngleBetweenNormalizedVectors( e2, e1 );
01121 
01122         m_secondaryNormalCurvature->operator[] ( vtxId ) = a * cos( theta ) * cos( theta )
01123                                                          + b * cos( theta ) * sin( theta )
01124                                                          + c * sin( theta ) * sin( theta );
01125         m_secondaryCurvaturePrincipalDirection->operator[] ( vtxId ) = e2;
01126     }
01127 
01128     m_curvatureCalculated = true;
01129 }
01130 
01131 double WTriangleMesh::calcAngleBetweenNormalizedVectors( osg::Vec3 const& v1, osg::Vec3 const& v2 )
01132 {
01133     // assumes vectors are normalized
01134     WAssert( v1.length() <  1.0001, "Vector is not normalized!" );
01135     WAssert( v1.length() > -1.0001, "Vector is not normalized!" );
01136     WAssert( v2.length() <  1.0001, "Vector is not normalized!" );
01137     WAssert( v2.length() > -1.0001, "Vector is not normalized!" );
01138 
01139     double temp = v1 * v2;
01140 
01141     // avoid NaNs due to numerical errors
01142     if( temp < -1.0 )
01143     {
01144         temp = -1.0;
01145     }
01146     if( temp > 1.0 )
01147     {
01148         temp = 1.0;
01149     }
01150 
01151     return acos( temp );
01152 }
01153 
01154 double WTriangleMesh::getMainCurvature( std::size_t vtxId ) const
01155 {
01156     return m_mainNormalCurvature->operator[] ( vtxId );
01157 }
01158 
01159 double WTriangleMesh::getSecondaryCurvature( std::size_t vtxId ) const
01160 {
01161     return m_secondaryNormalCurvature->operator[] ( vtxId );
01162 }
01163 
01164 boost::shared_ptr< std::vector< float > > const& WTriangleMesh::getMainCurvatures() const
01165 {
01166     return m_mainNormalCurvature;
01167 }
01168 
01169 boost::shared_ptr< std::vector< float > > const& WTriangleMesh::getSecondaryCurvatures() const
01170 {
01171     return m_secondaryNormalCurvature;
01172 }
01173 
01174 osg::Vec3 WTriangleMesh::getCurvatureMainPrincipalDirection( std::size_t vtxId ) const
01175 {
01176     return m_mainCurvaturePrincipalDirection->operator[] ( vtxId );
01177 }
01178 
01179 osg::Vec3 WTriangleMesh::getCurvatureSecondaryPrincipalDirection( std::size_t vtxId ) const
01180 {
01181     return m_secondaryCurvaturePrincipalDirection->operator[] ( vtxId );
01182 }
01183 
01184 void WTriangleMesh::setTextureCoord( std::size_t index, osg::Vec3 texCoord )
01185 {
01186     m_textureCoordinates->operator[] ( index ) = texCoord;
01187 }
01188 
01189 osg::ref_ptr< osg::Vec3Array > WTriangleMesh::getMainPrincipalCurvatureDirectionArray()
01190 {
01191     return m_mainCurvaturePrincipalDirection;
01192 }
01193 
01194 osg::ref_ptr< osg::Vec3Array > WTriangleMesh::getSecondaryPrincipalCurvatureDirectionArray()
01195 {
01196     return m_secondaryCurvaturePrincipalDirection;
01197 }
01198