2
/* osgEarth - Dynamic map generation toolkit for OpenSceneGraph
3
* Copyright 2008-2010 Pelican Mapping
6
* osgEarth is free software; you can redistribute it and/or modify
7
* it under the terms of the GNU Lesser General Public License as published by
8
* the Free Software Foundation; either version 2 of the License, or
9
* (at your option) any later version.
11
* This program is distributed in the hope that it will be useful,
12
* but WITHOUT ANY WARRANTY; without even the implied warranty of
13
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14
* GNU Lesser General Public License for more details.
16
* You should have received a copy of the GNU Lesser General Public License
17
* along with this program. If not, see <http://www.gnu.org/licenses/>
19
#include <osgEarthSymbology/MeshSubdivider>
20
#include <osgEarthSymbology/LineFunctor>
21
#include <osgEarth/GeoMath>
22
#include <osg/TriangleFunctor>
23
//#include <osgUtil/MeshOptimizers>
28
#define LC "[MeshSubdivider] "
30
using namespace osgEarth;
31
using namespace osgEarth::Symbology;
33
//------------------------------------------------------------------------
37
// convert geocenric coords to spherical geodetic coords in radians.
39
geocentricToGeodetic( const osg::Vec3d& g, osg::Vec2d& out_geod )
41
double r = g.length();
42
out_geod.set( atan2(g.y(),g.x()), acos(g.z()/r) );
45
// calculate the lat/long midpoint, taking care to use the shortest
48
geodeticMidpoint( const osg::Vec2d& g0, const osg::Vec2d& g1, osg::Vec2d& out_mid )
50
if ( fabs(g0.x()-g1.x()) < osg::PI )
51
out_mid.set( 0.5*(g0.x()+g1.x()), 0.5*(g0.y()+g1.y()) );
52
else if ( g1.x() > g0.x() )
53
out_mid.set( 0.5*((g0.x()+2*osg::PI)+g1.x()), 0.5*(g0.y()+g1.y()) );
55
out_mid.set( 0.5*(g0.x()+(g1.x()+2*osg::PI)), 0.5*(g0.y()+g1.y()) );
58
// finds the midpoint between two geocentric coordinates. We have to convert
59
// back to geographic in order to get the correct interpolation. Spherical
60
// conversion is good enough
62
geocentricMidpoint( const osg::Vec3d& v0, const osg::Vec3d& v1 )
64
// geocentric to spherical:
66
geocentricToGeodetic(v0, g0);
67
geocentricToGeodetic(v1, g1);
70
geodeticMidpoint(g0, g1, mid);
72
double size = 0.5*(v0.length() + v1.length());
74
// spherical to geocentric:
75
double sin_lat = sin(mid.y());
76
return osg::Vec3d( cos(mid.x())*sin_lat, sin(mid.x())*sin_lat, cos(mid.y()) ) * size;
79
// the approximate surface-distance between two geocentric points (spherical)
81
geocentricSurfaceDistance( const osg::Vec3d& v0, const osg::Vec3d& v1 )
84
geocentricToGeodetic(v0, g0);
85
geocentricToGeodetic(v1, g1);
86
return GeoMath::distance( v0.y(), v0.x(), v1.y(), v1.x() );
89
// returns the geocentric bisection vector
91
bisector( const osg::Vec3d& v0, const osg::Vec3d& v1 )
93
osg::Vec3d f = (v0+v1)*0.5;
95
return f * 0.5*(v0.length() + v1.length());
98
// the angle between two 3d vectors
100
angleBetween( const osg::Vec3d& v0, const osg::Vec3d& v1 )
102
osg::Vec3d v0n = v0; v0n.normalize();
103
osg::Vec3d v1n = v1; v1n.normalize();
104
return fabs( acos( v0n * v1n ) );
107
//--------------------------------------------------------------------
112
Triangle(GLuint i0, GLuint i1, GLuint i2) : _i0(i0), _i1(i1), _i2(i2) { }
113
GLuint _i0, _i1, _i2;
116
typedef std::queue<Triangle> TriangleQueue;
117
typedef std::vector<Triangle> TriangleVector;
121
typedef std::map<osg::Vec3,GLuint> VertMap;
123
osg::Vec3Array* _verts;
128
_verts = new osg::Vec3Array();
131
GLuint record( const osg::Vec3& v )
133
VertMap::iterator i = _vertMap.find(v);
134
if ( i == _vertMap.end() )
136
GLuint index = _verts->size();
137
_verts->push_back(v);
147
void operator()( const osg::Vec3& v0, const osg::Vec3& v1, const osg::Vec3& v2, bool temp )
149
_tris.push( Triangle(record(v0), record(v1), record(v2)) );
156
Edge( GLuint i0, GLuint i1 ) : _i0(i0), _i1(i1) { }
158
bool operator < (const Edge& rhs) const {
160
_i0 < rhs._i0 ? true :
161
_i0 > rhs._i0 ? false :
162
_i1 < rhs._i1 ? true :
163
_i1 > rhs._i1 ? false :
166
bool operator == (const Edge& rhs) const { return _i0==rhs._i0 && _i1==rhs._i1; }
169
typedef std::map<Edge,GLuint> EdgeMap;
172
* Populates the geometry object with a collection of index elements primitives.
174
template<typename ETYPE, typename VTYPE>
175
void populateTriangles( osg::Geometry& geom, const TriangleVector& tris, unsigned int maxElementsPerEBO )
177
unsigned int totalTris = tris.size();
178
unsigned int totalTrisWritten = 0;
179
unsigned int numElementsInCurrentEBO = maxElementsPerEBO;
183
for( TriangleVector::const_iterator i = tris.begin(); i != tris.end(); ++i )
185
if ( numElementsInCurrentEBO+2 >= maxElementsPerEBO )
189
geom.addPrimitiveSet( ebo );
192
ebo = new ETYPE( GL_TRIANGLES );
194
unsigned int trisRemaining = totalTris - totalTrisWritten;
195
ebo->reserve( osg::minimum( trisRemaining*3, maxElementsPerEBO ) );
197
numElementsInCurrentEBO = 0;
199
ebo->push_back( static_cast<VTYPE>( i->_i0 ) );
200
ebo->push_back( static_cast<VTYPE>( i->_i1 ) );
201
ebo->push_back( static_cast<VTYPE>( i->_i2 ) );
203
numElementsInCurrentEBO += 3;
207
if ( ebo && ebo->size() > 0 )
209
geom.addPrimitiveSet( ebo );
213
//----------------------------------------------------------------------
218
Line(GLuint i0, GLuint i1) : _i0(i0), _i1(i1) { }
222
typedef std::queue<Line> LineQueue;
223
typedef std::vector<Line> LineVector;
227
typedef std::map<osg::Vec3,GLuint> VertMap;
229
osg::Vec3Array* _verts;
234
_verts = new osg::Vec3Array();
237
GLuint record( const osg::Vec3& v )
239
VertMap::iterator i = _vertMap.find(v);
240
if ( i == _vertMap.end() )
242
GLuint index = _verts->size();
243
_verts->push_back(v);
253
void operator()( const osg::Vec3& v0, const osg::Vec3& v1, bool temp )
255
_lines.push( Line( record(v0), record(v1) ) );
260
* Populates the geometry object with a collection of index elements primitives.
262
template<typename ETYPE, typename VTYPE>
263
void populateLines( osg::Geometry& geom, const LineVector& lines, unsigned int maxElementsPerEBO )
265
unsigned int totalLines = lines.size();
266
unsigned int totalLinesWritten = 0;
267
unsigned int numElementsInCurrentEBO = maxElementsPerEBO;
271
for( LineVector::const_iterator i = lines.begin(); i != lines.end(); ++i )
273
if ( numElementsInCurrentEBO+2 >= maxElementsPerEBO )
277
geom.addPrimitiveSet( ebo );
280
ebo = new ETYPE( GL_LINES );
282
unsigned int linesRemaining = totalLines - totalLinesWritten;
283
ebo->reserve( osg::minimum( linesRemaining*2, maxElementsPerEBO ) );
285
numElementsInCurrentEBO = 0;
287
ebo->push_back( static_cast<VTYPE>( i->_i0 ) );
288
ebo->push_back( static_cast<VTYPE>( i->_i1 ) );
290
numElementsInCurrentEBO += 3;
294
if ( ebo && ebo->size() > 0 )
296
geom.addPrimitiveSet( ebo );
300
static const osg::Vec3d s_pole(0,0,1);
301
static const double s_maxLatAdjustment(0.75);
304
* Collects all the line segments from the geometry, coalesces them into a single
305
* line set, subdivides it according to the granularity threshold, and replaces
306
* the data in the Geometry object with the new vertex and primitive data.
311
const osg::Matrixd& W2L, // world=>local xform
312
const osg::Matrixd& L2W, // local=>world xform
313
unsigned int maxElementsPerEBO )
315
// collect all the line segments in the geometry.
316
LineFunctor<LineData> data;
319
int numLinesIn = data._lines.size();
322
done.reserve( 2 * data._lines.size() );
324
// Subdivide lines until we run out.
325
while( data._lines.size() > 0 )
327
Line line = data._lines.front();
330
osg::Vec3d v0_w = (*data._verts)[line._i0] * L2W;
331
osg::Vec3d v1_w = (*data._verts)[line._i1] * L2W;
333
double g0 = angleBetween(v0_w, v1_w);
335
if ( g0 > granularity )
337
data._verts->push_back( geocentricMidpoint(v0_w, v1_w) * W2L );
338
GLuint i = data._verts->size()-1;
340
data._lines.push( Line( line._i0, i ) );
341
data._lines.push( Line( i, line._i1 ) );
345
// line is small enough- put it on the "done" list.
346
done.push_back( line );
350
if ( done.size() > 0 )
352
while( geom.getNumPrimitiveSets() > 0 )
353
geom.removePrimitiveSet(0);
356
geom.setVertexArray( data._verts );
358
if ( data._verts->size() < 256 )
359
populateLines<osg::DrawElementsUByte,GLubyte>( geom, done, maxElementsPerEBO );
360
else if ( data._verts->size() < 65536 )
361
populateLines<osg::DrawElementsUShort,GLushort>( geom, done, maxElementsPerEBO );
363
populateLines<osg::DrawElementsUInt,GLuint>( geom, done, maxElementsPerEBO );
368
//----------------------------------------------------------------------
371
* Collects all the triangles from the geometry, coalesces them into a single
372
* triangle set, subdivides them according to the granularity threshold, and
373
* replaces the data in the Geometry object with the new vertex and primitive
376
* The subdivision algorithm is adapted from http://bit.ly/dTIagq
377
* (c) Copyright 2010 Patrick Cozzi and Deron Ohlarik, MIT License.
379
void subdivideTriangles(
382
const osg::Matrixd& W2L, // world=>local xform
383
const osg::Matrixd& L2W, // local=>world xform
384
unsigned int maxElementsPerEBO )
386
// collect all the triangled in the geometry.
387
osg::TriangleFunctor<TriangleData> data;
390
int numTrisIn = data._tris.size();
393
done.reserve(2.0 * data._tris.size());
395
// Used to make sure shared edges are not split more than once.
398
// Subdivide triangles until we run out
399
while( data._tris.size() > 0 )
401
Triangle tri = data._tris.front();
404
osg::Vec3d v0_w = (*data._verts)[tri._i0] * L2W;
405
osg::Vec3d v1_w = (*data._verts)[tri._i1] * L2W;
406
osg::Vec3d v2_w = (*data._verts)[tri._i2] * L2W;
408
double g0 = angleBetween(v0_w, v1_w);
409
double g1 = angleBetween(v1_w, v2_w);
410
double g2 = angleBetween(v2_w, v0_w);
411
double max = osg::maximum( g0, osg::maximum(g1, g2) );
413
if ( max > granularity )
417
Edge edge( osg::minimum(tri._i0, tri._i1), osg::maximum(tri._i0, tri._i1) );
419
EdgeMap::iterator ei = edges.find(edge);
421
if ( ei == edges.end() )
423
data._verts->push_back( geocentricMidpoint(v0_w, v1_w) * W2L );
424
i = data._verts->size() - 1;
432
data._tris.push( Triangle(tri._i0, i, tri._i2) );
433
data._tris.push( Triangle(i, tri._i1, tri._i2) );
435
else if ( g1 == max )
437
Edge edge( osg::minimum(tri._i1, tri._i2), osg::maximum(tri._i1,tri._i2) );
439
EdgeMap::iterator ei = edges.find(edge);
441
if ( ei == edges.end() )
443
data._verts->push_back( geocentricMidpoint(v1_w, v2_w) * W2L );
444
i = data._verts->size() - 1;
452
data._tris.push( Triangle(tri._i1, i, tri._i0) );
453
data._tris.push( Triangle(i, tri._i2, tri._i0) );
455
else if ( g2 == max )
457
Edge edge( osg::minimum(tri._i2, tri._i0), osg::maximum(tri._i2,tri._i0) );
459
EdgeMap::iterator ei = edges.find(edge);
461
if ( ei == edges.end() )
463
data._verts->push_back( geocentricMidpoint(v2_w, v0_w) * W2L );
464
i = data._verts->size() - 1;
472
data._tris.push( Triangle(tri._i2, i, tri._i1) );
473
data._tris.push( Triangle(i, tri._i0, tri._i1) );
478
// triangle is small enough- put it on the "done" list.
483
if ( done.size() > 0 )
485
// first, remove the old primitive sets.
486
while( geom.getNumPrimitiveSets() > 0 )
487
geom.removePrimitiveSet( 0 );
490
geom.setVertexArray( data._verts );
492
if ( data._verts->size() < 256 )
493
populateTriangles<osg::DrawElementsUByte,GLubyte>( geom, done, maxElementsPerEBO );
494
else if ( data._verts->size() < 65536 )
495
populateTriangles<osg::DrawElementsUShort,GLushort>( geom, done, maxElementsPerEBO );
497
populateTriangles<osg::DrawElementsUInt,GLuint>( geom, done, maxElementsPerEBO );
504
const osg::Matrixd& W2L, // world=>local xform
505
const osg::Matrixd& L2W, // local=>world xform
506
unsigned int maxElementsPerEBO )
508
GLenum mode = geom.getPrimitiveSet(0)->getMode();
510
if ( mode == GL_POINTS )
513
if ( mode == GL_LINES || mode == GL_LINE_STRIP || mode == GL_LINE_LOOP )
515
subdivideLines( granularity, geom, W2L, L2W, maxElementsPerEBO );
519
subdivideTriangles( granularity, geom, W2L, L2W, maxElementsPerEBO );
521
//osgUtil::VertexCacheVisitor cacheOptimizer;
522
//cacheOptimizer.optimizeVertices( geom );
525
//osgUtil::VertexAccessOrderVisitor orderOptimizer;
526
//orderOptimizer.optimizeOrder( geom );
530
//------------------------------------------------------------------------
532
MeshSubdivider::MeshSubdivider(const osg::Matrixd& world2local,
533
const osg::Matrixd& local2world ) :
534
_local2world(local2world),
535
_world2local(world2local),
536
_maxElementsPerEBO( INT_MAX )
538
if ( !_world2local.isIdentity() && _local2world.isIdentity() )
539
_local2world = osg::Matrixd::inverse(_world2local);
540
else if ( _world2local.isIdentity() && !_local2world.isIdentity() )
541
_world2local = osg::Matrixd::inverse(_local2world);
545
MeshSubdivider::run(double granularity, osg::Geometry& geom)
547
if ( geom.getNumPrimitiveSets() < 1 )
550
subdivide( granularity, geom, _world2local, _local2world, _maxElementsPerEBO );
added
removed
src/osgEarthSymbology/MeshSubdivider.cpp
