36 #include <osg/io_utils>
38 #include <Eigen/Eigenvalues>
40 #include "../common/WAssert.h"
41 #include "../common/WLogger.h"
42 #include "../common/math/WMath.h"
43 #include "../common/datastructures/WUnionFind.h"
44 #include "WTriangleMesh.h"
63 m_countTriangles( 0 ),
66 m_neighborsCalculated( false ),
67 m_curvatureCalculated( false )
69 m_verts = osg::ref_ptr< osg::Vec3Array >(
new osg::Vec3Array( vertNum ) );
71 m_vertNormals = osg::ref_ptr< osg::Vec3Array >(
new osg::Vec3Array( vertNum ) );
72 m_vertColors = osg::ref_ptr< osg::Vec4Array >(
new osg::Vec4Array( vertNum ) );
75 m_triangleNormals = osg::ref_ptr< osg::Vec3Array >(
new osg::Vec3Array( triangleNum ) );
76 m_triangleColors = osg::ref_ptr< osg::Vec4Array >(
new osg::Vec4Array( triangleNum ) );
80 : m_countVerts( vertices->size() ),
81 m_countTriangles( triangles.size() / 3 ),
84 m_neighborsCalculated( false ),
86 m_textureCoordinates( new osg::Vec3Array( vertices->size() ) ),
87 m_vertNormals( new osg::Vec3Array( vertices->size() ) ),
88 m_vertColors( new osg::Vec4Array( vertices->size() ) ),
89 m_triangles( triangles ),
90 m_triangleNormals( new osg::Vec3Array( triangles.size() / 3 ) ),
91 m_triangleColors( new osg::Vec4Array( triangles.size() / 3 ) )
93 WAssert( triangles.size() % 3 == 0,
"Invalid triangle vector, having an invalid size (not divideable by 3)" );
107 return addVertex( osg::Vec3( x, y, z ) );
112 return addVertex( osg::Vec3( vert[0], vert[1], vert[2] ) );
137 WAssert( index <
m_countVerts,
"set vertex normal: index out of range" );
138 ( *m_vertNormals )[index] = normal;
143 setVertexNormal( index, osg::Vec3( normal[0], normal[1], normal[2] ) );
153 WAssert( index <
m_countVerts,
"set vertex color: index out of range" );
154 ( *m_vertColors )[index] = color;
159 WAssert( index <
m_countTriangles,
"set triangle color: index out of range" );
160 ( *m_triangleColors )[index] = color;
214 WAssert( index <
m_countVerts,
"get vertex: index out of range" );
220 WAssert( index <
m_countVerts,
"get vertex color: index out of range" );
226 WAssert( index <
m_countVerts,
"get normal as position: index out of range" );
236 WAssert( index <
m_countVerts,
"remove vertex: index out of range" );
241 ( *m_verts ).erase( ( *m_verts ).begin() + index );
272 for(
size_t vertId = 0; vertId <
m_countVerts; ++vertId )
274 osg::Vec3 tempNormal( 0.0, 0.0, 0.0 );
281 tempNormal.normalize();
282 ( *m_vertNormals )[vertId] = tempNormal;
291 std::vector< size_t >v;
307 osg::Vec3 tempNormal( 0, 0, 0 );
309 tempNormal[0] = v1[1] * v2[2] - v1[2] * v2[1];
310 tempNormal[1] = v1[2] * v2[0] - v1[0] * v2[2];
311 tempNormal[2] = v1[0] * v2[1] - v1[1] * v2[0];
313 tempNormal.normalize();
320 osg::Vec3 v1( vert1 - vert0 );
321 osg::Vec3 v2( vert2 - vert0 );
323 osg::Vec3 tempNormal( 0, 0, 0 );
325 tempNormal[0] = v1[1] * v2[2] - v1[2] * v2[1];
326 tempNormal[1] = v1[2] * v2[0] - v1[0] * v2[2];
327 tempNormal[2] = v1[0] * v2[1] - v1[1] * v2[0];
329 tempNormal.normalize();
346 std::vector<size_t> v( 3, -1 );
367 for(
size_t i = 0; i < candidates.size(); ++i )
369 for(
size_t k = 0; k < compares.size(); ++k )
371 if( ( candidates[i] != triangleNum ) && ( candidates[i] == compares[k] ) )
373 return candidates[i];
390 osg::Vec3* newVertexPositions =
new osg::Vec3[
m_numTriVerts];
404 std::vector< size_t >v;
416 ( *m_verts )[i] = newVertexPositions[i];
419 delete[] newVertexPositions;
432 int starSize = starP.size();
437 double scale = 1.0 - (
static_cast<double>( starSize ) * alpha );
442 for(
int i = 0; i < starSize; i++ )
445 osg::Vec3 translate =
getVertex( edgeV );
450 return oldPos + newPos;
461 addTriangle( edgeVerts[0], edgeVerts[1], edgeVerts[2] );
467 size_t neighborVert = -1;
468 size_t neighborFaceNum = -1;
471 neighborFaceNum =
getNeighbor( edgeV1, edgeV2, triId );
473 if( neighborFaceNum == triId )
475 osg::Vec3 edgeVert = ( ( *m_verts )[edgeV1] + ( *m_verts )[edgeV2] ) / 2.0;
481 else if( neighborFaceNum > triId )
485 osg::Vec3 edgePart = ( *m_verts )[edgeV1] + ( *m_verts )[edgeV2];
486 osg::Vec3 neighborPart = ( *m_verts )[neighborVert] + ( *m_verts )[V3];
488 edgeVert = ( ( edgePart * ( 3.0 / 8.0 ) ) + ( neighborPart * ( 1.0 / 8.0 ) ) );
496 size_t neighborP = neighborFaceNum;
531 addTriangle( originalTri1, centerTri1, centerTri0 );
532 addTriangle( originalTri2, centerTri2, centerTri1 );
551 std::vector< size_t > temp;
565 double center = ( 0.375 + ( 0.25 * cos( ( 2.0 * 3.14159265358979 ) / static_cast<double>( n ) ) ) );
566 answer = ( 0.625 - ( center * center ) ) /
static_cast<double>( n );
605 ( *m_vertColors ).resize( oldVertSize + mesh->vertSize() );
606 for(
size_t i = 0; i < mesh->vertSize(); ++i )
608 osg::Vec3 v( mesh->getVertex( i ) );
615 for(
size_t i = 0; i < mesh->triangleSize(); ++i )
617 addTriangle( mesh->getTriVertId0( i ) + oldVertSize, mesh->getTriVertId1( i ) + oldVertSize, mesh->getTriVertId2( i ) + oldVertSize );
624 osg::Vec3 t( xOff, yOff, zOff );
625 for(
size_t i = 0; i < ( *m_verts ).size(); ++i )
627 ( *m_verts )[i] += t;
633 for(
size_t i = 0; i < ( *m_verts ).size(); ++i )
635 ( *m_verts )[i] *= zoom;
644 for(
size_t vertId = 0; vertId <
m_vertColors->size(); ++vertId )
648 maxR = ( *m_vertColors )[vertId][0];
652 maxG = ( *m_vertColors )[vertId][1];
656 maxB = ( *m_vertColors )[vertId][2];
659 for(
size_t vertId = 0; vertId <
m_vertColors->size(); ++vertId )
661 ( *m_vertColors )[vertId][0] /= maxR;
662 ( *m_vertColors )[vertId][1] /= maxG;
663 ( *m_vertColors )[vertId][2] /= maxB;
669 std::stringstream ss;
670 ss <<
"WTriangleMesh( #vertices=" << rhs.
vertSize() <<
" #triangles=" << rhs.
triangleSize() <<
" )" << std::endl;
671 using string_utils::operator<<;
674 const std::vector< size_t >& triangles = rhs.
getTriangles();
675 osg::ref_ptr< const osg::Vec3Array > vertices = rhs.
getVertexArray();
676 for(
size_t vID = 0 ; vID <= triangles.size() - 3; ++count )
678 std::stringstream prefix;
679 prefix <<
"triangle: " << count <<
"[ ";
680 std::string indent( prefix.str().size(),
' ' );
681 using osg::operator<<;
682 ss << prefix.str() << vertices->at( triangles[ vID++ ] ) << std::endl;
683 ss << indent << vertices->at( triangles[ vID++ ] ) << std::endl;
684 ss << indent << vertices->at( triangles[ vID++ ] ) << std::endl;
685 ss << std::string( indent.size() - 2,
' ' ) <<
"]" << std::endl;
687 return os << ss.str();
692 boost::shared_ptr< std::list< boost::shared_ptr< WTriangleMesh > > > result(
new std::list< boost::shared_ptr< WTriangleMesh > >() );
702 WAssert(
false,
"Not implemented the decomposition of a TriangleMesh without any triangles" );
712 const std::vector< size_t >& triangles = mesh.
getTriangles();
713 for(
size_t vID = 0; vID <= triangles.size() - 3; vID += 3)
715 uf.merge( triangles[ vID ], triangles[ vID + 1 ] );
716 uf.merge( triangles[ vID ], triangles[ vID + 2 ] );
724 typedef std::map< osg::Vec3, size_t > VertexType;
725 typedef std::vector< size_t > TriangleType;
726 typedef std::pair< VertexType, TriangleType > BucketType;
727 std::map< size_t, BucketType > buckets;
729 osg::ref_ptr< const osg::Vec3Array > vertices = mesh.
getVertexArray();
730 for(
size_t vID = 0; vID <= triangles.size() - 3; vID += 3 )
732 size_t component = uf.find( triangles[ vID ] );
733 if( buckets.find( component ) == buckets.end() )
735 buckets[ component ] = BucketType( VertexType(), TriangleType() );
739 VertexType& mapRef = buckets[ component ].first;
740 for(
int i = 0; i < 3; ++i )
743 const osg::Vec3& vertex = ( *vertices )[ triangles[ vID + i ] ];
744 if( mapRef.find( vertex ) == mapRef.end() )
747 mapRef[ vertex ] = id;
751 id = mapRef[ vertex ];
753 buckets[ component ].second.push_back(
id );
757 for( std::map< size_t, BucketType >::const_iterator cit = buckets.begin(); cit != buckets.end(); ++cit )
759 osg::ref_ptr< osg::Vec3Array > newVertices(
new osg::Vec3Array );
760 newVertices->resize( cit->second.first.size() );
761 for( VertexType::const_iterator vit = cit->second.first.begin(); vit != cit->second.first.end(); ++vit )
763 newVertices->at( vit->second ) = vit->first;
765 boost::shared_ptr< WTriangleMesh > newMesh(
new WTriangleMesh( newVertices, cit->second.second ) );
766 result->push_back( newMesh );
794 osg::ref_ptr< osg::Vec3Array > vtxArray =
new osg::Vec3Array(
m_verts->size() );
796 for( std::size_t k = 0; k <
m_verts->size(); ++k )
802 for( std::size_t k = 0; k <
m_triangles.size() / 3; ++k )
804 osg::Vec3
const& p0 = vtxArray->operator[](
m_triangles[ 3 * k + 0 ] );
805 osg::Vec3
const& p1 = vtxArray->operator[](
m_triangles[ 3 * k + 1 ] );
806 osg::Vec3
const& p2 = vtxArray->operator[](
m_triangles[ 3 * k + 2 ] );
815 for( std::size_t k = 0; k <
m_verts->size(); ++k )
825 std::stack< std::size_t > triStack;
826 std::set< std::size_t > triSet;
834 while( !triStack.empty() )
836 std::size_t currentTriangle = triStack.top();
843 if( (
m_verts->operator[] ( vtx ) - center ).length() > 4.0 * sigmaDistance )
857 osg::Vec3 res( 0.0, 0.0, 0.0 );
859 for( std::set< std::size_t >::const_iterator it = triSet.begin(); it != triSet.end(); ++it )
865 double dist = (
m_verts->operator[] ( vtx ) - center ).length();
866 double f = 1.0 / ( sigmaDistance * sqrt( 2.0 * piDouble ) ) * exp( -0.5 * dist * dist );
872 osg::Vec3f
const& p =
m_verts->operator[] ( vtx );
874 osg::Vec3f pred = p + n * ( n * ( center - p ) );
876 dist = ( p - pred ).length();
878 g = 1.0 / ( sigmaInfluence * sqrt( 2.0 * piDouble ) ) * exp( -0.5 * dist * dist );
881 res += center * area * f * g;
904 return ( ( p1 - p0 ) ^ ( p2 - p0 ) ).length() * 0.5;
912 std::vector< osg::Vec3 > normals(
m_verts->size() );
921 for( std::size_t vtxId = 0; vtxId <
m_verts->size(); ++vtxId )
923 osg::Vec3
const& p =
m_verts->operator[] ( vtxId );
924 osg::Vec3 n( 0.0, 0.0, 0.0 );
931 double w = 1.0 / ( center - p ).length();
936 WAssert( n.length() > 0.0001,
"Invalid normal!" );
939 normals[ vtxId ] = n;
943 for( std::size_t vtxId = 0; vtxId <
m_verts->size(); ++vtxId )
945 osg::Vec3
const& p =
m_verts->operator[] ( vtxId );
947 osg::Vec3
const& normal = normals[ vtxId ];
950 std::set< std::size_t > neighbors;
956 for( std::size_t j = 0; j < 3; ++j )
960 if( neighbors.find( e ) == neighbors.end() && e != vtxId )
962 neighbors.insert( e );
967 WAssert( neighbors.size() > 2,
"Vertex has too few neighbors! Does this mesh have holes?" );
969 double maxCurvature = -std::numeric_limits< double >::infinity();
970 std::vector< double > curvatures;
972 osg::Vec3 maxCurvatureTangent( 0.0, 0.0, 0.0 );
973 std::vector< osg::Vec3 > tangents;
976 for( std::set< std::size_t >::const_iterator it = neighbors.begin(); it != neighbors.end(); ++it )
978 osg::Vec3
const& neighbPos =
m_verts->operator[] ( *it );
979 osg::Vec3
const& neighbNormal = normals[ *it ];
982 osg::Vec3 tangent = ( neighbPos - p ) - normal * ( ( neighbPos - p ) * normal );
986 double ncurv = -1.0 * ( ( neighbPos - p ) * ( neighbNormal - normal ) ) / ( ( neighbPos - p ) * ( neighbPos - p ) );
988 if( ncurv > maxCurvature )
990 maxCurvature = ncurv;
991 maxCurvatureTangent = tangent;
994 tangents.push_back( tangent );
995 curvatures.push_back( ncurv );
998 WAssert( maxCurvatureTangent.length() > 0.0001,
"Invalid tangent length!" );
1001 osg::Vec3
const& e1 = maxCurvatureTangent;
1002 osg::Vec3 e2 = e1 ^ normal;
1006 bool significantCurvature =
false;
1007 for( std::vector< double >::const_iterator it = curvatures.begin(); it != curvatures.end(); ++it )
1009 if( fabs( *it ) > 0.00001 )
1011 significantCurvature =
true;
1016 if( !significantCurvature )
1020 m_mainNormalCurvature->operator[] ( vtxId ) = 0.0;
1031 double const& a = maxCurvature;
1033 Eigen::Matrix< double, -1, -1 > X( tangents.size(), 2 );
1034 Eigen::Matrix< double, -1, -1 > y( tangents.size(), 1 );
1036 for( std::size_t k = 0; k < tangents.size(); ++k )
1040 X( k, 0 ) = cos( theta ) * sin( theta );
1041 X( k, 1 ) = sin( theta ) * sin( theta );
1042 y( k, 0 ) = curvatures[ k ] - a * cos( theta ) * cos( theta );
1046 Eigen::FullPivLU< Eigen::Matrix< double, -1, -1 > > lu( X );
1051 wlog::warn(
"WTriangleMesh::estimateCurvature" ) <<
"Rank too low, cannot estimate curvature for this vertex!";
1053 m_mainNormalCurvature->operator[] ( vtxId ) = a;
1062 Eigen::Matrix< double, -1, -1 > bb = ( X.transpose() * X ).inverse() * X.transpose() * y;
1065 double b = bb( 0, 0 );
1066 double c = bb( 1, 0 );
1070 double Kg = a * c - 0.25 * b * b;
1071 double H = 0.5 * ( a + c );
1072 double s = H * H - Kg;
1079 double k1 = H + sqrt( s );
1080 double k2 = H - sqrt( s );
1082 if( fabs( k1 - k2 ) < 0.000001 )
1085 m_mainNormalCurvature->operator[] ( vtxId ) = a;
1091 double temp = b / ( k2 - k1 );
1102 double theta = 0.5 * asin( temp );
1104 osg::Vec3 ne1 = e1 * cos( theta ) + e2 * sin( theta );
1105 osg::Vec3 ne2 = e2 * cos( theta ) - e1 * sin( theta );
1114 m_mainNormalCurvature->operator[] ( vtxId ) = a * cos( theta ) * cos( theta )
1115 + b * cos( theta ) * sin( theta )
1116 + c * sin( theta ) * sin( theta );
1123 + b * cos( theta ) * sin( theta )
1124 + c * sin( theta ) * sin( theta );
1134 WAssert( v1.length() < 1.0001,
"Vector is not normalized!" );
1135 WAssert( v1.length() > -1.0001,
"Vector is not normalized!" );
1136 WAssert( v2.length() < 1.0001,
"Vector is not normalized!" );
1137 WAssert( v2.length() > -1.0001,
"Vector is not normalized!" );
1139 double temp = v1 * v2;
1151 return acos( temp );