OpenWalnut  1.4.0
WTriangleMesh.cpp
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #include <iostream>
26 #include <list>
27 #include <map>
28 #include <set>
29 #include <sstream>
30 #include <stack>
31 #include <string>
32 #include <utility>
33 #include <vector>
34 #include <limits>
35 
36 #include <osg/io_utils>
37 
38 #include <Eigen/Eigenvalues>
39 
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"
45 
46 // init _static_ member variable and provide a linker reference to it
47 boost::shared_ptr< WPrototyped > WTriangleMesh::m_prototype = boost::shared_ptr< WPrototyped >();
48 
49 boost::shared_ptr< WPrototyped > WTriangleMesh::getPrototype()
50 {
51  if( !m_prototype )
52  {
53  m_prototype = boost::shared_ptr< WPrototyped >( new WTriangleMesh( 0, 0 ) );
54  }
55  return m_prototype;
56 }
57 
58 /**
59  * constructor that already reserves space for a given number of triangles and vertexes
60  */
61 WTriangleMesh::WTriangleMesh( size_t vertNum, size_t triangleNum )
62  : m_countVerts( 0 ),
63  m_countTriangles( 0 ),
64  m_meshDirty( true ),
65  m_autoNormal( true ),
66  m_neighborsCalculated( false ),
67  m_curvatureCalculated( false )
68 {
69  m_verts = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( vertNum ) );
70  m_textureCoordinates = 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 ) );
73 
74  m_triangles.resize( triangleNum * 3 );
75  m_triangleNormals = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( triangleNum ) );
76  m_triangleColors = osg::ref_ptr< osg::Vec4Array >( new osg::Vec4Array( triangleNum ) );
77 }
78 
79 WTriangleMesh::WTriangleMesh( osg::ref_ptr< osg::Vec3Array > vertices, const std::vector< size_t >& triangles )
80  : m_countVerts( vertices->size() ),
81  m_countTriangles( triangles.size() / 3 ),
82  m_meshDirty( true ),
83  m_autoNormal( true ),
84  m_neighborsCalculated( false ),
85  m_verts( vertices ),
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 ) )
92 {
93  WAssert( triangles.size() % 3 == 0, "Invalid triangle vector, having an invalid size (not divideable by 3)" );
94 }
95 
97 {
98 }
99 
100 void WTriangleMesh::setAutoRecalcNormals( bool autoRecalc )
101 {
102  m_autoNormal = autoRecalc;
103 }
104 
105 size_t WTriangleMesh::addVertex( float x, float y, float z )
106 {
107  return addVertex( osg::Vec3( x, y, z ) );
108 }
109 
111 {
112  return addVertex( osg::Vec3( vert[0], vert[1], vert[2] ) );
113 }
114 
115 void WTriangleMesh::addTriangle( size_t vert0, size_t vert1, size_t vert2 )
116 {
117  if( m_triangles.size() == m_countTriangles * 3 )
118  {
119  m_triangles.resize( m_countTriangles * 3 + 3 );
120  }
121  m_triangles[ m_countTriangles * 3 ] = vert0;
122  m_triangles[ m_countTriangles * 3 + 1 ] = vert1;
123  m_triangles[ m_countTriangles * 3 + 2 ] = vert2;
125 }
126 
127 void WTriangleMesh::addTriangle( osg::Vec3 vert0, osg::Vec3 vert1, osg::Vec3 vert2 )
128 {
129  addVertex( vert0 );
130  addVertex( vert1 );
131  addVertex( vert2 );
133 }
134 
135 void WTriangleMesh::setVertexNormal( size_t index, osg::Vec3 normal )
136 {
137  WAssert( index < m_countVerts, "set vertex normal: index out of range" );
138  ( *m_vertNormals )[index] = normal;
139 }
140 
141 void WTriangleMesh::setVertexNormal( size_t index, WPosition normal )
142 {
143  setVertexNormal( index, osg::Vec3( normal[0], normal[1], normal[2] ) );
144 }
145 
146 void WTriangleMesh::setVertexNormal( size_t index, float x, float y, float z )
147 {
148  setVertexNormal( index, osg::Vec3( x, y, z ) );
149 }
150 
151 void WTriangleMesh::setVertexColor( size_t index, osg::Vec4 color )
152 {
153  WAssert( index < m_countVerts, "set vertex color: index out of range" );
154  ( *m_vertColors )[index] = color;
155 }
156 
157 void WTriangleMesh::setTriangleColor( size_t index, osg::Vec4 color )
158 {
159  WAssert( index < m_countTriangles, "set triangle color: index out of range" );
160  ( *m_triangleColors )[index] = color;
161 }
162 
163 osg::ref_ptr< osg::Vec3Array >WTriangleMesh::getVertexArray()
164 {
165  return m_verts;
166 }
167 
168 osg::ref_ptr< const osg::Vec3Array >WTriangleMesh::getVertexArray() const
169 {
170  return m_verts;
171 }
172 
173 osg::ref_ptr< osg::Vec3Array > WTriangleMesh::getTextureCoordinateArray()
174 {
175  return m_textureCoordinates;
176 }
177 
178 osg::ref_ptr< const osg::Vec3Array > WTriangleMesh::getTextureCoordinateArray() const
179 {
180  return m_textureCoordinates;
181 }
182 
183 osg::ref_ptr< osg::Vec3Array >WTriangleMesh::getVertexNormalArray( bool forceRecalc )
184 {
185  if( forceRecalc || ( m_meshDirty && m_autoNormal ) )
186  {
188  }
189  return m_vertNormals;
190 }
191 
192 osg::ref_ptr< osg::Vec3Array >WTriangleMesh::getTriangleNormalArray( bool forceRecalc )
193 {
194  if( forceRecalc || ( m_meshDirty && m_autoNormal ) )
195  {
197  }
198  return m_triangleNormals;
199 }
200 
201 
202 osg::ref_ptr< osg::Vec4Array >WTriangleMesh::getVertexColorArray()
203 {
204  return m_vertColors;
205 }
206 
207 const std::vector< size_t >& WTriangleMesh::getTriangles() const
208 {
209  return m_triangles;
210 }
211 
212 osg::Vec3 WTriangleMesh::getVertex( size_t index ) const
213 {
214  WAssert( index < m_countVerts, "get vertex: index out of range" );
215  return ( *m_verts )[index];
216 }
217 
218 osg::Vec4 WTriangleMesh::getVertColor( size_t index ) const
219 {
220  WAssert( index < m_countVerts, "get vertex color: index out of range" );
221  return ( *m_vertColors )[index];
222 }
223 
225 {
226  WAssert( index < m_countVerts, "get normal as position: index out of range" );
227  if( m_meshDirty )
228  {
230  }
231  return WPosition( ( *m_vertNormals )[index][0], ( *m_vertNormals )[index][1], ( *m_vertNormals )[index][2] );
232 }
233 
234 void WTriangleMesh::removeVertex( size_t index )
235 {
236  WAssert( index < m_countVerts, "remove vertex: index out of range" );
237  if( m_vertexIsInTriangle[ index ].size() > 0 )
238  {
239  return;
240  }
241  ( *m_verts ).erase( ( *m_verts ).begin() + index );
242 
243  for( size_t i = 0; i < m_countTriangles * 3; ++i )
244  {
245  if( m_triangles[ i ] > index )
246  {
247  --m_triangles[ i ];
248  }
249  }
250  m_meshDirty = true;
251 }
252 
253 void WTriangleMesh::removeTriangle( size_t index )
254 {
255  WAssert( index < m_countTriangles, "remove triangle: index out of range" );
256  m_triangles.erase( m_triangles.begin() + index * 3, m_triangles.begin() + index * 3 + 3 );
257  m_meshDirty = true;
258 }
259 
261 {
263 
264  ( *m_vertNormals ).resize( m_countVerts );
265  ( *m_triangleNormals ).resize( m_countTriangles );
266 
267  for( size_t i = 0; i < m_countTriangles; ++i )
268  {
269  ( *m_triangleNormals )[i] = calcTriangleNormal( i );
270  }
271 
272  for( size_t vertId = 0; vertId < m_countVerts; ++vertId )
273  {
274  osg::Vec3 tempNormal( 0.0, 0.0, 0.0 );
275  for( size_t neighbour = 0; neighbour < m_vertexIsInTriangle[vertId].size(); ++neighbour )
276  {
277  tempNormal += ( *m_triangleNormals )[ m_vertexIsInTriangle[vertId][neighbour] ];
278  }
279  tempNormal *= 1./m_vertexIsInTriangle[vertId].size();
280 
281  tempNormal.normalize();
282  ( *m_vertNormals )[vertId] = tempNormal;
283  }
284 
285  m_meshDirty = false;
286 }
287 
289 {
290  m_vertexIsInTriangle.clear();
291  std::vector< size_t >v;
292  m_vertexIsInTriangle.resize( ( *m_verts ).size(), v );
293 
294  for( size_t i = 0; i < m_countTriangles; ++i )
295  {
296  m_vertexIsInTriangle[ getTriVertId0( i ) ].push_back( i );
297  m_vertexIsInTriangle[ getTriVertId1( i ) ].push_back( i );
298  m_vertexIsInTriangle[ getTriVertId2( i ) ].push_back( i );
299  }
300 }
301 
302 osg::Vec3 WTriangleMesh::calcTriangleNormal( size_t triangle )
303 {
304  osg::Vec3 v1( getTriVert( triangle, 1 ) - getTriVert( triangle, 0 ) );
305  osg::Vec3 v2( getTriVert( triangle, 2 ) - getTriVert( triangle, 0 ) );
306 
307  osg::Vec3 tempNormal( 0, 0, 0 );
308 
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];
312 
313  tempNormal.normalize();
314 
315  return tempNormal;
316 }
317 
318 osg::Vec3 WTriangleMesh::calcNormal( osg::Vec3 vert0, osg::Vec3 vert1, osg::Vec3 vert2 )
319 {
320  osg::Vec3 v1( vert1 - vert0 );
321  osg::Vec3 v2( vert2 - vert0 );
322 
323  osg::Vec3 tempNormal( 0, 0, 0 );
324 
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];
328 
329  tempNormal.normalize();
330 
331  return tempNormal;
332 }
333 
335 {
336  return m_countVerts;
337 }
338 
340 {
341  return m_countTriangles;
342 }
343 
345 {
346  std::vector<size_t> v( 3, -1 );
347  m_triangleNeighbors.resize( ( *m_triangleNormals ).size(), v );
348 
349  for( size_t triId = 0; triId < m_countTriangles; ++triId )
350  {
351  size_t coVert0 = getTriVertId0( triId );
352  size_t coVert1 = getTriVertId1( triId );
353  size_t coVert2 = getTriVertId2( triId );
354 
355  m_triangleNeighbors[triId][0] = getNeighbor( coVert0, coVert1, triId );
356  m_triangleNeighbors[triId][1] = getNeighbor( coVert1, coVert2, triId );
357  m_triangleNeighbors[triId][2] = getNeighbor( coVert2, coVert0, triId );
358  }
359  m_neighborsCalculated = true;
360 }
361 
362 size_t WTriangleMesh::getNeighbor( const size_t coVert1, const size_t coVert2, const size_t triangleNum )
363 {
364  std::vector< size_t > candidates = m_vertexIsInTriangle[coVert1];
365  std::vector< size_t > compares = m_vertexIsInTriangle[coVert2];
366 
367  for( size_t i = 0; i < candidates.size(); ++i )
368  {
369  for( size_t k = 0; k < compares.size(); ++k )
370  {
371  if( ( candidates[i] != triangleNum ) && ( candidates[i] == compares[k] ) )
372  {
373  return candidates[i];
374  }
375  }
376  }
377  return triangleNum;
378 }
379 
381 {
384 
385  ( *m_verts ).resize( m_numTriVerts * 4 );
386  m_triangles.resize( m_numTriFaces * 4 * 3 );
387 
389 
390  osg::Vec3* newVertexPositions = new osg::Vec3[m_numTriVerts];
391 
392  //std::cout << "Loop subdivision pass 1" << std::endl;
393  for( size_t i = 0; i < m_numTriVerts; ++i )
394  {
395  newVertexPositions[i] = loopCalcNewPosition( i );
396  }
397 
398  //std::cout << "Loop subdivision pass 2" << std::endl;
399  for( size_t i = 0; i < m_numTriFaces; ++i )
400  {
402  }
403  ( *m_verts ).resize( m_countVerts );
404  std::vector< size_t >v;
405  m_vertexIsInTriangle.resize( ( *m_verts ).size(), v );
406 
407  //std::cout << "Loop subdivision pass 3" << std::endl;
408  for( size_t i = 0; i < m_numTriFaces; ++i )
409  {
411  }
412 
413  //std::cout << "Loop subdivision pass 4" << std::endl;
414  for( size_t i = 0; i < m_numTriVerts; ++i )
415  {
416  ( *m_verts )[i] = newVertexPositions[i];
417  }
418 
419  delete[] newVertexPositions;
420 
421  m_vertNormals->resize( m_verts->size() );
422  m_vertColors->resize( m_verts->size() );
423  m_triangleColors->resize( m_triangles.size() / 3 );
424 
425  m_meshDirty = true;
426 }
427 
428 
429 osg::Vec3 WTriangleMesh::loopCalcNewPosition( size_t vertId )
430 {
431  std::vector< size_t > starP = m_vertexIsInTriangle[vertId];
432  int starSize = starP.size();
433 
434  osg::Vec3 oldPos = getVertex( vertId );
435  double alpha = loopGetAlpha( starSize );
436 
437  double scale = 1.0 - ( static_cast<double>( starSize ) * alpha );
438  oldPos *= scale;
439 
440  osg::Vec3 newPos;
441  int edgeV = 0;
442  for( int i = 0; i < starSize; i++ )
443  {
444  edgeV = loopGetNextVertex( starP[i], vertId );
445  osg::Vec3 translate = getVertex( edgeV );
446  newPos += translate;
447  }
448  newPos *= alpha;
449 
450  return oldPos + newPos;
451 }
452 
454 {
455  size_t edgeVerts[3];
456 
457  edgeVerts[0] = loopCalcEdgeVert( triId, getTriVertId0( triId ), getTriVertId1( triId ), getTriVertId2( triId ) );
458  edgeVerts[1] = loopCalcEdgeVert( triId, getTriVertId1( triId ), getTriVertId2( triId ), getTriVertId0( triId ) );
459  edgeVerts[2] = loopCalcEdgeVert( triId, getTriVertId2( triId ), getTriVertId0( triId ), getTriVertId1( triId ) );
460 
461  addTriangle( edgeVerts[0], edgeVerts[1], edgeVerts[2] );
462 }
463 
464 
465 size_t WTriangleMesh::loopCalcEdgeVert( size_t triId, size_t edgeV1, size_t edgeV2, size_t V3 )
466 {
467  size_t neighborVert = -1;
468  size_t neighborFaceNum = -1;
469  osg::Vec3 edgeVert;
470 
471  neighborFaceNum = getNeighbor( edgeV1, edgeV2, triId );
472 
473  if( neighborFaceNum == triId )
474  {
475  osg::Vec3 edgeVert = ( ( *m_verts )[edgeV1] + ( *m_verts )[edgeV2] ) / 2.0;
476  size_t vertId = m_countVerts;
477  addVertex( edgeVert );
478  return vertId;
479  }
480 
481  else if( neighborFaceNum > triId )
482  {
483  neighborVert = loopGetThirdVert( edgeV1, edgeV2, neighborFaceNum );
484 
485  osg::Vec3 edgePart = ( *m_verts )[edgeV1] + ( *m_verts )[edgeV2];
486  osg::Vec3 neighborPart = ( *m_verts )[neighborVert] + ( *m_verts )[V3];
487 
488  edgeVert = ( ( edgePart * ( 3.0 / 8.0 ) ) + ( neighborPart * ( 1.0 / 8.0 ) ) );
489  size_t vertId = m_countVerts;
490  addVertex( edgeVert );
491  return vertId;
492  }
493  else
494  {
495  size_t neighborCenterP = neighborFaceNum + m_numTriFaces;
496  size_t neighborP = neighborFaceNum;
497 
498  if( getTriVertId0( neighborP ) == edgeV2 )
499  {
500  return getTriVertId0( neighborCenterP );
501  }
502  else if( getTriVertId1( neighborP ) == edgeV2 )
503  {
504  return getTriVertId1( neighborCenterP );
505  }
506  else
507  {
508  return getTriVertId2( neighborCenterP );
509  }
510  }
511  return -1;
512 }
513 
515 {
516  // comment: center are twisted from the orignal vertices.
517  // original: 0, 1, 2
518  // center: a, b, c
519  // reAsgnOrig: 0, a, c
520  // addTris: 1, b, a
521  // addTris: 2, c, b
522  //
523  size_t originalTri0 = getTriVertId0( triId );
524  size_t originalTri1 = getTriVertId1( triId );
525  size_t originalTri2 = getTriVertId2( triId );
526 
527  size_t centerTri0 = getTriVertId0( triId + m_numTriFaces );
528  size_t centerTri1 = getTriVertId1( triId + m_numTriFaces );
529  size_t centerTri2 = getTriVertId2( triId + m_numTriFaces );
530 
531  addTriangle( originalTri1, centerTri1, centerTri0 );
532  addTriangle( originalTri2, centerTri2, centerTri1 );
533  loopSetTriangle( triId, originalTri0, centerTri0, centerTri2 );
534 }
535 
536 void WTriangleMesh::loopSetTriangle( size_t triId, size_t vertId1, size_t vertId2, size_t vertId3 )
537 {
538  loopEraseTriangleFromVertex( triId, getTriVertId1( triId ) );
539  loopEraseTriangleFromVertex( triId, getTriVertId2( triId ) );
540 
541  setTriVert0( triId, vertId1 );
542  setTriVert1( triId, vertId2 );
543  setTriVert2( triId, vertId3 );
544 
545  m_vertexIsInTriangle[vertId2].push_back( triId );
546  m_vertexIsInTriangle[vertId3].push_back( triId );
547 }
548 
549 void WTriangleMesh::loopEraseTriangleFromVertex( size_t triId, size_t vertId )
550 {
551  std::vector< size_t > temp;
552  for( size_t i = 0; i < m_vertexIsInTriangle[vertId].size(); ++i )
553  {
554  if( triId != m_vertexIsInTriangle[vertId][i] )
555  temp.push_back( m_vertexIsInTriangle[vertId][i] );
556  }
557  m_vertexIsInTriangle[vertId] = temp;
558 }
559 
561 {
562  double answer;
563  if( n > 3 )
564  {
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 );
567  }
568  else
569  {
570  answer = 3.0 / 16.0;
571  }
572  return answer;
573 }
574 
575 size_t WTriangleMesh::loopGetNextVertex( size_t triNum, size_t vertNum )
576 {
577  if( getTriVertId0( triNum ) == vertNum )
578  {
579  return getTriVertId1( triNum );
580  }
581  else if( getTriVertId1( triNum ) == vertNum )
582  {
583  return getTriVertId2( triNum );
584  }
585  return getTriVertId0( triNum );
586 }
587 
588 size_t WTriangleMesh::loopGetThirdVert( size_t coVert1, size_t coVert2, size_t triangleNum )
589 {
590  if( !( getTriVertId0( triangleNum ) == coVert1 ) && !( getTriVertId0( triangleNum ) == coVert2 ) )
591  {
592  return getTriVertId0( triangleNum );
593  }
594  else if( !( getTriVertId1( triangleNum ) == coVert1 ) && !( getTriVertId1( triangleNum ) == coVert2 ) )
595  {
596  return getTriVertId1( triangleNum );
597  }
598  return getTriVertId2( triangleNum );
599 }
600 
601 void WTriangleMesh::addMesh( boost::shared_ptr<WTriangleMesh> mesh, float xOff, float yOff, float zOff )
602 {
603  size_t oldVertSize = m_countVerts;
604 
605  ( *m_vertColors ).resize( oldVertSize + mesh->vertSize() );
606  for( size_t i = 0; i < mesh->vertSize(); ++i )
607  {
608  osg::Vec3 v( mesh->getVertex( i ) );
609  v[0] += xOff;
610  v[1] += yOff;
611  v[2] += zOff;
612  addVertex( v );
613  setVertexColor( oldVertSize + i, mesh->getVertColor( i ) );
614  }
615  for( size_t i = 0; i < mesh->triangleSize(); ++i )
616  {
617  addTriangle( mesh->getTriVertId0( i ) + oldVertSize, mesh->getTriVertId1( i ) + oldVertSize, mesh->getTriVertId2( i ) + oldVertSize );
618  }
619  m_meshDirty = true;
620 }
621 
622 void WTriangleMesh::translateMesh( float xOff, float yOff, float zOff )
623 {
624  osg::Vec3 t( xOff, yOff, zOff );
625  for( size_t i = 0; i < ( *m_verts ).size(); ++i )
626  {
627  ( *m_verts )[i] += t;
628  }
629 }
630 
631 void WTriangleMesh::zoomMesh( float zoom )
632 {
633  for( size_t i = 0; i < ( *m_verts ).size(); ++i )
634  {
635  ( *m_verts )[i] *= zoom;
636  }
637 }
638 
640 {
641  float maxR = 0;
642  float maxG = 0;
643  float maxB = 0;
644  for( size_t vertId = 0; vertId < m_vertColors->size(); ++vertId )
645  {
646  if( ( *m_vertColors )[vertId][0] > maxR )
647  {
648  maxR = ( *m_vertColors )[vertId][0];
649  }
650  if( ( *m_vertColors )[vertId][1] > maxG )
651  {
652  maxG = ( *m_vertColors )[vertId][1];
653  }
654  if( ( *m_vertColors )[vertId][2] > maxB )
655  {
656  maxB = ( *m_vertColors )[vertId][2];
657  }
658  }
659  for( size_t vertId = 0; vertId < m_vertColors->size(); ++vertId )
660  {
661  ( *m_vertColors )[vertId][0] /= maxR;
662  ( *m_vertColors )[vertId][1] /= maxG;
663  ( *m_vertColors )[vertId][2] /= maxB;
664  }
665 }
666 
667 std::ostream& tm_utils::operator<<( std::ostream& os, const WTriangleMesh& rhs )
668 {
669  std::stringstream ss;
670  ss << "WTriangleMesh( #vertices=" << rhs.vertSize() << " #triangles=" << rhs.triangleSize() << " )" << std::endl;
671  using string_utils::operator<<;
672  size_t count = 0;
673  ss << std::endl;
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 )
677  {
678  std::stringstream prefix;
679  prefix << "triangle: " << count << "[ ";
680  std::string indent( prefix.str().size(), ' ' );
681  using osg::operator<<; // using operator<< as defined in osg/io_utils
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;
686  }
687  return os << ss.str();
688 }
689 
690 boost::shared_ptr< std::list< boost::shared_ptr< WTriangleMesh > > > tm_utils::componentDecomposition( const WTriangleMesh& mesh )
691 {
692  boost::shared_ptr< std::list< boost::shared_ptr< WTriangleMesh > > > result( new std::list< boost::shared_ptr< WTriangleMesh > >() );
693  if( mesh.vertSize() <= 0 ) // no component possible
694  {
695  return result;
696  }
697  if( mesh.triangleSize() < 3 )
698  {
699  if( mesh.vertSize() > 0 )
700  {
701  // there are vertices but no triangles
702  WAssert( false, "Not implemented the decomposition of a TriangleMesh without any triangles" );
703  }
704  else // no component possible
705  {
706  return result;
707  }
708  }
709 
710  WUnionFind uf( mesh.vertSize() ); // idea: every vertex in own component, then successivley join in accordance with the triangles
711 
712  const std::vector< size_t >& triangles = mesh.getTriangles();
713  for( size_t vID = 0; vID <= triangles.size() - 3; vID += 3)
714  {
715  uf.merge( triangles[ vID ], triangles[ vID + 1 ] );
716  uf.merge( triangles[ vID ], triangles[ vID + 2 ] ); // uf.merge( triangles[ vID + 2 ], triangles[ vID + 1 ] ); they are already in same
717  }
718 
719  // ATTENTION: The reason for using the complex BucketType instead of pasting vertices directly into a new WTriangleMesh
720  // is performance! For example: If there are many vertices reused inside the former WTriangleMesh mesh, then we want
721  // to reuse them in the new components too. Hence we must determine if a certain vertex is already inside the new component.
722  // Since the vertices are organized in a vector, we can use std::find( v.begin, v.end(), vertexToLookUp ) which results
723  // in O(N^2) or we could use faster lookUp via key and value leading to the map and the somehow complicated BucketType.
724  typedef std::map< osg::Vec3, size_t > VertexType; // look up fast if a vertex is already inside the new mesh!
725  typedef std::vector< size_t > TriangleType;
726  typedef std::pair< VertexType, TriangleType > BucketType; // Later on the Bucket will be transformed into the new WTriangleMesh component
727  std::map< size_t, BucketType > buckets; // Key identify with the cannonical element from UnionFind the new connected component
728 
729  osg::ref_ptr< const osg::Vec3Array > vertices = mesh.getVertexArray();
730  for( size_t vID = 0; vID <= triangles.size() - 3; vID += 3 )
731  {
732  size_t component = uf.find( triangles[ vID ] );
733  if( buckets.find( component ) == buckets.end() )
734  {
735  buckets[ component ] = BucketType( VertexType(), TriangleType() ); // create new bucket
736  }
737 
738  // Note: We discard the order of the points and indices, but semantically the structure remains the same
739  VertexType& mapRef = buckets[ component ].first; // short hand alias
740  for( int i = 0; i < 3; ++i )
741  {
742  size_t id = 0;
743  const osg::Vec3& vertex = ( *vertices )[ triangles[ vID + i ] ];
744  if( mapRef.find( vertex ) == mapRef.end() )
745  {
746  id = mapRef.size(); // since size might change in next line
747  mapRef[ vertex ] = id;
748  }
749  else
750  {
751  id = mapRef[ vertex ];
752  }
753  buckets[ component ].second.push_back( id );
754  }
755  }
756 
757  for( std::map< size_t, BucketType >::const_iterator cit = buckets.begin(); cit != buckets.end(); ++cit )
758  {
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 )
762  {
763  newVertices->at( vit->second ) = vit->first; // if you are sure that vit->second is always valid replace at() call with operator[]
764  }
765  boost::shared_ptr< WTriangleMesh > newMesh( new WTriangleMesh( newVertices, cit->second.second ) );
766  result->push_back( newMesh );
767  }
768 
769  return result;
770 }
771 
772 osg::ref_ptr< osg::Vec4Array > WTriangleMesh::getTriangleColors() const
773 {
774  return m_triangleColors;
775 }
776 
777 void WTriangleMesh::performFeaturePreservingSmoothing( float sigmaDistance, float sigmaInfluence )
778 {
780 
781  calcNeighbors();
782 
783  // we perform a first smoothing pass and write the resulting vertex coords into a buffer
784  // this will only update the normals
785  performFeaturePreservingSmoothingMollificationPass( sigmaDistance / 2.0f, sigmaInfluence );
786 
787  // using the smoothed normals, we now perform a second smoothing pass, this time writing the new vertex coords
788  performFeaturePreservingSmoothingVertexPass( sigmaDistance, sigmaInfluence );
789 }
790 
791 void WTriangleMesh::performFeaturePreservingSmoothingMollificationPass( float sigmaDistance, float sigmaInfluence )
792 {
793  // calc Eq. 3 for every triangle
794  osg::ref_ptr< osg::Vec3Array > vtxArray = new osg::Vec3Array( m_verts->size() );
795 
796  for( std::size_t k = 0; k < m_verts->size(); ++k )
797  {
798  vtxArray->operator[] ( k ) = estimateSmoothedVertexPosition( k, sigmaDistance, sigmaInfluence, true );
799  }
800 
801  // calc the new normal directions - update triangle normals
802  for( std::size_t k = 0; k < m_triangles.size() / 3; ++k )
803  {
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 ] );
807 
808  m_triangleNormals->operator[] ( k ) = ( p1 - p0 ) ^ ( p2 - p1 );
809  m_triangleNormals->operator[] ( k ).normalize();
810  }
811 }
812 
813 void WTriangleMesh::performFeaturePreservingSmoothingVertexPass( float sigmaDistance, float sigmaInfluence )
814 {
815  for( std::size_t k = 0; k < m_verts->size(); ++k )
816  {
817  m_verts->operator[] ( k ) = estimateSmoothedVertexPosition( k, sigmaDistance, sigmaInfluence, false );
818  }
819 
821 }
822 
823 osg::Vec3 WTriangleMesh::estimateSmoothedVertexPosition( std::size_t vtx, float sigmaDistance, float sigmaInfluence, bool mollify )
824 {
825  std::stack< std::size_t > triStack;
826  std::set< std::size_t > triSet;
827 
828  for( std::size_t k = 0; k < m_vertexIsInTriangle[ vtx ].size(); ++k )
829  {
830  triStack.push( m_vertexIsInTriangle[ vtx ][ k ] );
831  triSet.insert( m_vertexIsInTriangle[ vtx ][ k ] );
832  }
833 
834  while( !triStack.empty() )
835  {
836  std::size_t currentTriangle = triStack.top();
837  triStack.pop();
838 
839  for( std::size_t k = 0; k < m_triangleNeighbors[ currentTriangle ].size(); ++k )
840  {
841  osg::Vec3 center = calcTriangleCenter( m_triangleNeighbors[ currentTriangle ][ k ] );
842 
843  if( ( m_verts->operator[] ( vtx ) - center ).length() > 4.0 * sigmaDistance )
844  {
845  continue;
846  }
847 
848  if( triSet.find( m_triangleNeighbors[ currentTriangle ][ k ] ) == triSet.end() )
849  {
850  triStack.push( m_triangleNeighbors[ currentTriangle ][ k ] );
851  triSet.insert( m_triangleNeighbors[ currentTriangle ][ k ] );
852  }
853  }
854  }
855 
856  double sum = 0.0;
857  osg::Vec3 res( 0.0, 0.0, 0.0 );
858 
859  for( std::set< std::size_t >::const_iterator it = triSet.begin(); it != triSet.end(); ++it )
860  {
861  osg::Vec3 center = calcTriangleCenter( *it );
862  double area = calcTriangleArea( *it );
863 
864  // calc f
865  double dist = ( m_verts->operator[] ( vtx ) - center ).length();
866  double f = 1.0 / ( sigmaDistance * sqrt( 2.0 * piDouble ) ) * exp( -0.5 * dist * dist );
867 
868  double g;
869  if( !mollify )
870  {
871  // calc prediction
872  osg::Vec3f const& p = m_verts->operator[] ( vtx );
873  osg::Vec3f const& n = m_triangleNormals->operator[] ( *it );
874  osg::Vec3f pred = p + n * ( n * ( center - p ) );
875 
876  dist = ( p - pred ).length();
877  }
878  g = 1.0 / ( sigmaInfluence * sqrt( 2.0 * piDouble ) ) * exp( -0.5 * dist * dist );
879 
880  sum += area * f * g;
881  res += center * area * f * g;
882  }
883 
884  res /= sum;
885  return res;
886 }
887 
888 osg::Vec3 WTriangleMesh::calcTriangleCenter( std::size_t triIdx ) const
889 {
890  osg::Vec3 res = m_verts->operator[] ( m_triangles[ 3 * triIdx + 0 ] );
891  res += m_verts->operator[] ( m_triangles[ 3 * triIdx + 1 ] );
892  res += m_verts->operator[] ( m_triangles[ 3 * triIdx + 2 ] );
893 
894  res /= 3.0f;
895  return res;
896 }
897 
898 float WTriangleMesh::calcTriangleArea( std::size_t triIdx ) const
899 {
900  osg::Vec3 const& p0 = m_verts->operator[] ( m_triangles[ 3 * triIdx + 0 ] );
901  osg::Vec3 const& p1 = m_verts->operator[] ( m_triangles[ 3 * triIdx + 1 ] );
902  osg::Vec3 const& p2 = m_verts->operator[] ( m_triangles[ 3 * triIdx + 2 ] );
903 
904  return ( ( p1 - p0 ) ^ ( p2 - p0 ) ).length() * 0.5;
905 }
906 
908 {
910  calcNeighbors();
911 
912  std::vector< osg::Vec3 > normals( m_verts->size() );
913 
914  // init vectors
915  m_mainNormalCurvature = boost::shared_ptr< std::vector< float > >( new std::vector< float >( m_verts->size() ) );
916  m_secondaryNormalCurvature = boost::shared_ptr< std::vector< float > >( new std::vector< float >( m_verts->size() ) );
917  m_mainCurvaturePrincipalDirection = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( m_verts->size() ) );
918  m_secondaryCurvaturePrincipalDirection = osg::ref_ptr< osg::Vec3Array >( new osg::Vec3Array( m_verts->size() ) );
919 
920  // calculate vertex normals using distance-weighted summing of neighbor-triangle normals
921  for( std::size_t vtxId = 0; vtxId < m_verts->size(); ++vtxId )
922  {
923  osg::Vec3 const& p = m_verts->operator[] ( vtxId );
924  osg::Vec3 n( 0.0, 0.0, 0.0 );
925 
926  for( std::size_t k = 0; k < m_vertexIsInTriangle[ vtxId ].size(); ++k )
927  {
928  std::size_t triId = m_vertexIsInTriangle[ vtxId ][ k ];
929 
930  osg::Vec3 center = calcTriangleCenter( triId );
931  double w = 1.0 / ( center - p ).length();
932 
933  n += m_triangleNormals->operator[] ( triId ) * w;
934  }
935 
936  WAssert( n.length() > 0.0001, "Invalid normal!" );
937 
938  n.normalize();
939  normals[ vtxId ] = n;
940  }
941 
942  // calculate curvatures for every vertex
943  for( std::size_t vtxId = 0; vtxId < m_verts->size(); ++vtxId )
944  {
945  osg::Vec3 const& p = m_verts->operator[] ( vtxId );
946 
947  osg::Vec3 const& normal = normals[ vtxId ];
948 
949  // get the set of neighbor vertices
950  std::set< std::size_t > neighbors;
951 
952  for( std::size_t k = 0; k < m_vertexIsInTriangle[ vtxId ].size(); ++k )
953  {
954  std::size_t triId = m_vertexIsInTriangle[ vtxId ][ k ];
955 
956  for( std::size_t j = 0; j < 3; ++j )
957  {
958  std::size_t e = m_triangles[ 3 * triId + j ];
959 
960  if( neighbors.find( e ) == neighbors.end() && e != vtxId )
961  {
962  neighbors.insert( e );
963  }
964  }
965  }
966 
967  WAssert( neighbors.size() > 2, "Vertex has too few neighbors! Does this mesh have holes?" );
968 
969  double maxCurvature = -std::numeric_limits< double >::infinity();
970  std::vector< double > curvatures;
971 
972  osg::Vec3 maxCurvatureTangent( 0.0, 0.0, 0.0 );
973  std::vector< osg::Vec3 > tangents;
974 
975  // part 1: get curvatures at tangents and their maximum curvature
976  for( std::set< std::size_t >::const_iterator it = neighbors.begin(); it != neighbors.end(); ++it )
977  {
978  osg::Vec3 const& neighbPos = m_verts->operator[] ( *it );
979  osg::Vec3 const& neighbNormal = normals[ *it ];
980 
981  // project ( neighbPos - p ) onto the tangent plane
982  osg::Vec3 tangent = ( neighbPos - p ) - normal * ( ( neighbPos - p ) * normal );
983  tangent.normalize();
984 
985  // approximate normal curvature in tangent direction
986  double ncurv = -1.0 * ( ( neighbPos - p ) * ( neighbNormal - normal ) ) / ( ( neighbPos - p ) * ( neighbPos - p ) );
987 
988  if( ncurv > maxCurvature )
989  {
990  maxCurvature = ncurv;
991  maxCurvatureTangent = tangent;
992  }
993 
994  tangents.push_back( tangent );
995  curvatures.push_back( ncurv );
996  }
997 
998  WAssert( maxCurvatureTangent.length() > 0.0001, "Invalid tangent length!" );
999 
1000  // part 2: choose a coordinate system in the tangent plane
1001  osg::Vec3 const& e1 = maxCurvatureTangent;
1002  osg::Vec3 e2 = e1 ^ normal;
1003 
1004  e2.normalize();
1005 
1006  bool significantCurvature = false;
1007  for( std::vector< double >::const_iterator it = curvatures.begin(); it != curvatures.end(); ++it )
1008  {
1009  if( fabs( *it ) > 0.00001 )
1010  {
1011  significantCurvature = true;
1012  break;
1013  }
1014  }
1015 
1016  if( !significantCurvature )
1017  {
1018  // curvatures were all almost zero
1019  // the mesh is flat at this point, write values accordingly
1020  m_mainNormalCurvature->operator[] ( vtxId ) = 0.0;
1021  m_mainCurvaturePrincipalDirection->operator[] ( vtxId ) = e1;
1022  m_secondaryNormalCurvature->operator[] ( vtxId ) = 0.0;
1023  m_secondaryCurvaturePrincipalDirection->operator[] ( vtxId ) = e2;
1024 
1025  continue;
1026  }
1027 
1028  // calculate coefficients of the ellipse a * cos²(theta) + b * cos(theta) * sin(theta) * c * sin²(theta)
1029  // this is done by estimating a as the largest curvature amoung the estimated curvatures in all tangent
1030  // direction belonging to points that share a triangle with the current one
1031  double const& a = maxCurvature;
1032 
1033  Eigen::Matrix< double, -1, -1 > X( tangents.size(), 2 );
1034  Eigen::Matrix< double, -1, -1 > y( tangents.size(), 1 );
1035 
1036  for( std::size_t k = 0; k < tangents.size(); ++k )
1037  {
1038  double theta = calcAngleBetweenNormalizedVectors( tangents[ k ], e1 );
1039 
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 );
1043  }
1044 
1045  // use LU decomposition to calculate rank
1046  Eigen::FullPivLU< Eigen::Matrix< double, -1, -1 > > lu( X );
1047 
1048  // we need a rank of at least 2 to calculate the coeffs
1049  if( lu.rank() < 2 )
1050  {
1051  wlog::warn( "WTriangleMesh::estimateCurvature" ) << "Rank too low, cannot estimate curvature for this vertex!";
1052 
1053  m_mainNormalCurvature->operator[] ( vtxId ) = a;
1054  m_mainCurvaturePrincipalDirection->operator[] ( vtxId ) = e1;
1055  m_secondaryNormalCurvature->operator[] ( vtxId ) = 0.0;
1056  m_secondaryCurvaturePrincipalDirection->operator[] ( vtxId ) = e2;
1057 
1058  continue;
1059  }
1060 
1061  // do least squares
1062  Eigen::Matrix< double, -1, -1 > bb = ( X.transpose() * X ).inverse() * X.transpose() * y;
1063 
1064  // the missing coeffs b and c are now:
1065  double b = bb( 0, 0 );
1066  double c = bb( 1, 0 );
1067 
1068  // this calculates the maximum and minimum normal curvatures
1069  // (as eigenvalues)
1070  double Kg = a * c - 0.25 * b * b;
1071  double H = 0.5 * ( a + c );
1072  double s = H * H - Kg;
1073 
1074  if( s < 0.0 )
1075  {
1076  s = 0.0;
1077  }
1078 
1079  double k1 = H + sqrt( s );
1080  double k2 = H - sqrt( s );
1081 
1082  if( fabs( k1 - k2 ) < 0.000001 )
1083  {
1084  // if the curvatures are equal, there is no single principal direction
1085  m_mainNormalCurvature->operator[] ( vtxId ) = a;
1086  m_mainCurvaturePrincipalDirection->operator[] ( vtxId ) = e1;
1087  }
1088  else
1089  {
1090  // if the curvatures differ, we can now find the direction of maximum curvature
1091  double temp = b / ( k2 - k1 );
1092 
1093  if( temp > 1.0 )
1094  {
1095  temp = 1.0;
1096  }
1097  if( temp < -1.0 )
1098  {
1099  temp = -1.0;
1100  }
1101 
1102  double theta = 0.5 * asin( temp );
1103 
1104  osg::Vec3 ne1 = e1 * cos( theta ) + e2 * sin( theta );
1105  osg::Vec3 ne2 = e2 * cos( theta ) - e1 * sin( theta );
1106 
1107  ne1.normalize();
1108  ne2.normalize();
1109 
1110  e2 = ne2;
1111 
1112  theta = calcAngleBetweenNormalizedVectors( ne1, e1 );
1113 
1114  m_mainNormalCurvature->operator[] ( vtxId ) = a * cos( theta ) * cos( theta )
1115  + b * cos( theta ) * sin( theta )
1116  + c * sin( theta ) * sin( theta );
1117  m_mainCurvaturePrincipalDirection->operator[] ( vtxId ) = ne1;
1118  }
1119 
1120  double theta = calcAngleBetweenNormalizedVectors( e2, e1 );
1121 
1122  m_secondaryNormalCurvature->operator[] ( vtxId ) = a * cos( theta ) * cos( theta )
1123  + b * cos( theta ) * sin( theta )
1124  + c * sin( theta ) * sin( theta );
1125  m_secondaryCurvaturePrincipalDirection->operator[] ( vtxId ) = e2;
1126  }
1127 
1128  m_curvatureCalculated = true;
1129 }
1130 
1131 double WTriangleMesh::calcAngleBetweenNormalizedVectors( osg::Vec3 const& v1, osg::Vec3 const& v2 )
1132 {
1133  // assumes vectors are normalized
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!" );
1138 
1139  double temp = v1 * v2;
1140 
1141  // avoid NaNs due to numerical errors
1142  if( temp < -1.0 )
1143  {
1144  temp = -1.0;
1145  }
1146  if( temp > 1.0 )
1147  {
1148  temp = 1.0;
1149  }
1150 
1151  return acos( temp );
1152 }
1153 
1154 double WTriangleMesh::getMainCurvature( std::size_t vtxId ) const
1155 {
1156  return m_mainNormalCurvature->operator[] ( vtxId );
1157 }
1158 
1159 double WTriangleMesh::getSecondaryCurvature( std::size_t vtxId ) const
1160 {
1161  return m_secondaryNormalCurvature->operator[] ( vtxId );
1162 }
1163 
1164 boost::shared_ptr< std::vector< float > > const& WTriangleMesh::getMainCurvatures() const
1165 {
1166  return m_mainNormalCurvature;
1167 }
1168 
1169 boost::shared_ptr< std::vector< float > > const& WTriangleMesh::getSecondaryCurvatures() const
1170 {
1172 }
1173 
1174 osg::Vec3 WTriangleMesh::getCurvatureMainPrincipalDirection( std::size_t vtxId ) const
1175 {
1176  return m_mainCurvaturePrincipalDirection->operator[] ( vtxId );
1177 }
1178 
1180 {
1181  return m_secondaryCurvaturePrincipalDirection->operator[] ( vtxId );
1182 }
1183 
1184 void WTriangleMesh::setTextureCoord( std::size_t index, osg::Vec3 texCoord )
1185 {
1186  m_textureCoordinates->operator[] ( index ) = texCoord;
1187 }
1188 
1190 {
1192 }
1193 
1195 {
1197 }
1198