Skip to content

Commit 5b74bdf

Browse files
committedOct 12, 2017
Fix numerical problems in the tessellator
Due to the recent changes to allow non-horizontal surfaces, tessellator was loosing precision in float math with large numbers, so there could be artefacts on extruded buildings between roof and walls Also clean up the code a bit - move stuff out of the main method + remove some duplicate code
1 parent d97a51b commit 5b74bdf

File tree

1 file changed

+85
-80
lines changed

1 file changed

+85
-80
lines changed
 

‎src/3d/qgstessellator.cpp

Lines changed: 85 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -168,7 +168,61 @@ static QVector3D _calculateNormal( const QgsCurve *curve, double originX, double
168168
return normal;
169169
}
170170

171-
void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHeight )
171+
172+
static void _normalVectorToXYVectors( const QVector3D &pNormal, QVector3D &pXVector, QVector3D &pYVector )
173+
{
174+
// Here we define the two perpendicular vectors that define the local
175+
// 2D space on the plane. They will act as axis for which we will
176+
// calculate the projection coordinates of a 3D point to the plane.
177+
if ( pNormal.z() > 0.001 || pNormal.z() < -0.001 )
178+
{
179+
pXVector = QVector3D( 1, 0, -pNormal.x() / pNormal.z() );
180+
}
181+
else if ( pNormal.y() > 0.001 || pNormal.y() < -0.001 )
182+
{
183+
pXVector = QVector3D( 1, -pNormal.x() / pNormal.y(), 0 );
184+
}
185+
else
186+
{
187+
pXVector = QVector3D( -pNormal.y() / pNormal.x(), 1, 0 );
188+
}
189+
pXVector.normalize();
190+
pYVector = QVector3D::normal( pNormal, pXVector );
191+
}
192+
193+
194+
static void _ringToPoly2tri( const QgsCurve *ring, const QgsPoint &ptFirst, const QVector3D &pXVector, const QVector3D &pYVector, std::vector<p2t::Point *> &polyline, QHash<p2t::Point *, float> &zHash )
195+
{
196+
QgsVertexId::VertexType vt;
197+
QgsPoint pt;
198+
199+
const int pCount = ring->numPoints();
200+
double x0 = ptFirst.x(), y0 = ptFirst.y(), z0 = ( std::isnan( ptFirst.z() ) ? 0 : ptFirst.z() );
201+
202+
polyline.reserve( pCount );
203+
204+
for ( int i = 0; i < pCount - 1; ++i )
205+
{
206+
ring->pointAt( i, pt, vt );
207+
const float z = std::isnan( pt.z() ) ? 0 : pt.z();
208+
QVector3D tempPt( pt.x() - x0, pt.y() - y0, z - z0 );
209+
const float x = QVector3D::dotProduct( tempPt, pXVector );
210+
const float y = QVector3D::dotProduct( tempPt, pYVector );
211+
212+
const bool found = std::find_if( polyline.begin(), polyline.end(), [x, y]( p2t::Point *&p ) { return *p == p2t::Point( x, y ); } ) != polyline.end();
213+
214+
if ( found )
215+
{
216+
continue;
217+
}
218+
219+
p2t::Point *pt2 = new p2t::Point( x, y );
220+
polyline.push_back( pt2 );
221+
zHash[pt2] = z;
222+
}
223+
}
224+
225+
static bool _check_intersecting_rings( const QgsPolygonV2 &polygon )
172226
{
173227
// At this point we assume that input polygons are valid according to the OGC definition.
174228
// This means e.g. no duplicate points, polygons are simple (no butterfly shaped polygon with self-intersection),
@@ -198,24 +252,28 @@ void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHei
198252
for ( int j = i + 1; j < geomRings.count(); ++j )
199253
{
200254
if ( geomRings[i].intersects( geomRings[j] ) )
201-
{
202-
// skip the polygon - it would cause a crash inside poly2tri library
203-
qDebug() << "polygon rings intersect each other - skipping";
204-
return;
205-
}
255+
return false;
206256
}
207257
}
258+
return true;
259+
}
260+
261+
262+
void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHeight )
263+
{
264+
if ( !_check_intersecting_rings( polygon ) )
265+
{
266+
// skip the polygon - it would cause a crash inside poly2tri library
267+
qDebug() << "polygon rings intersect each other - skipping";
268+
return;
269+
}
208270

209271
const QgsCurve *exterior = polygon.exteriorRing();
210272

211273
QList< std::vector<p2t::Point *> > polylinesToDelete;
212274
QHash<p2t::Point *, float> z;
213275

214276
std::vector<p2t::Point *> polyline;
215-
polyline.reserve( exterior->numPoints() );
216-
217-
QgsVertexId::VertexType vt;
218-
QgsPoint pt;
219277

220278
const QVector3D pNormal = _calculateNormal( exterior, mOriginX, mOriginY );
221279
const int pCount = exterior->numPoints();
@@ -224,6 +282,7 @@ void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHei
224282
if ( pCount == 4 )
225283
{
226284
QgsPoint pt;
285+
QgsVertexId::VertexType vt;
227286
for ( int i = 0; i < 3; i++ )
228287
{
229288
exterior->pointAt( i, pt, vt );
@@ -234,65 +293,28 @@ void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHei
234293
}
235294
else
236295
{
237-
238296
if ( !qgsDoubleNear( pNormal.length(), 1, 0.001 ) )
239-
{
240-
return;
241-
}
242-
243-
const QgsPoint ptFirst( exterior->startPoint() );
244-
QVector3D pOrigin( ptFirst.x(), ptFirst.y(), std::isnan( ptFirst.z() ) ? 0 : ptFirst.z() );
245-
QVector3D pXVector;
246-
// Here we define the two perpendicular vectors that define the local
247-
// 2D space on the plane. They will act as axis for which we will
248-
// calculate the projection coordinates of a 3D point to the plane.
249-
if ( pNormal.z() > 0.001 || pNormal.z() < -0.001 )
250-
{
251-
pXVector = QVector3D( 1, 0, -pNormal.x() / pNormal.z() );
252-
}
253-
else if ( pNormal.y() > 0.001 || pNormal.y() < -0.001 )
254-
{
255-
pXVector = QVector3D( 1, -pNormal.x() / pNormal.y(), 0 );
256-
}
257-
else
258-
{
259-
pXVector = QVector3D( -pNormal.y() / pNormal.x(), 1, 0 );
260-
}
261-
pXVector.normalize();
262-
const QVector3D pYVector = QVector3D::normal( pNormal, pXVector );
297+
return; // this should not happen - pNormal should be normalized to unit length
263298

264-
for ( int i = 0; i < pCount - 1; ++i )
265-
{
266-
exterior->pointAt( i, pt, vt );
267-
QVector3D tempPt( pt.x(), pt.y(), ( std::isnan( pt.z() ) ? 0 : pt.z() ) );
268-
const float x = QVector3D::dotProduct( tempPt - pOrigin, pXVector );
269-
const float y = QVector3D::dotProduct( tempPt - pOrigin, pYVector );
270-
271-
const bool found = std::find_if( polyline.begin(), polyline.end(), [x, y]( p2t::Point *&p ) { return *p == p2t::Point( x, y ); } ) != polyline.end();
299+
QVector3D pXVector, pYVector;
300+
_normalVectorToXYVectors( pNormal, pXVector, pYVector );
272301

273-
if ( found )
274-
{
275-
continue;
276-
}
277-
278-
p2t::Point *pt2 = new p2t::Point( x, y );
279-
polyline.push_back( pt2 );
280-
281-
z[pt2] = std::isnan( pt.z() ) ? 0 : pt.z();
282-
}
302+
const QgsPoint ptFirst( exterior->startPoint() );
303+
_ringToPoly2tri( exterior, ptFirst, pXVector, pYVector, polyline, z );
283304
polylinesToDelete << polyline;
284305

285306
// TODO: robustness (no nearly duplicate points, invalid geometries ...)
286307

308+
double x0 = ptFirst.x(), y0 = ptFirst.y();
287309
if ( polyline.size() == 3 && polygon.numInteriorRings() == 0 )
288310
{
289311
for ( std::vector<p2t::Point *>::iterator it = polyline.begin(); it != polyline.end(); it++ )
290312
{
291313
p2t::Point *p = *it;
292314
const double zPt = z[p];
293-
QVector3D nPoint = pOrigin + pXVector * p->x + pYVector * p->y;
294-
const double fx = nPoint.x() - mOriginX;
295-
const double fy = nPoint.y() - mOriginY;
315+
QVector3D nPoint = pXVector * p->x + pYVector * p->y;
316+
const double fx = nPoint.x() - mOriginX + x0;
317+
const double fy = nPoint.y() - mOriginY + y0;
296318
const double fz = extrusionHeight + ( std::isnan( zPt ) ? 0 : zPt );
297319
mData << fx << fz << -fy;
298320
if ( mAddNormals )
@@ -307,31 +329,14 @@ void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHei
307329
for ( int i = 0; i < polygon.numInteriorRings(); ++i )
308330
{
309331
std::vector<p2t::Point *> holePolyline;
310-
holePolyline.reserve( exterior->numPoints() );
311332
const QgsCurve *hole = polygon.interiorRing( i );
312-
for ( int j = 0; j < hole->numPoints() - 1; ++j )
313-
{
314-
hole->pointAt( j, pt, vt );
315-
QVector3D tempPt( pt.x(), pt.y(), ( std::isnan( pt.z() ) ? 0 : pt.z() ) );
316-
317-
const float x = QVector3D::dotProduct( tempPt - pOrigin, pXVector );
318-
const float y = QVector3D::dotProduct( tempPt - pOrigin, pYVector );
319-
320-
const bool found = std::find_if( polyline.begin(), polyline.end(), [x, y]( p2t::Point *&p ) { return *p == p2t::Point( x, y ); } ) != polyline.end();
321-
322-
if ( found )
323-
{
324-
continue;
325-
}
326333

327-
p2t::Point *pt2 = new p2t::Point( x, y );
328-
holePolyline.push_back( pt2 );
334+
_ringToPoly2tri( hole, ptFirst, pXVector, pYVector, holePolyline, z );
329335

330-
z[pt2] = std::isnan( pt.z() ) ? 0 : pt.z();
331-
}
332336
cdt->AddHole( holePolyline );
333337
polylinesToDelete << holePolyline;
334338
}
339+
335340
try
336341
{
337342
cdt->Triangulate();
@@ -344,11 +349,11 @@ void QgsTessellator::addPolygon( const QgsPolygonV2 &polygon, float extrusionHei
344349
for ( int j = 0; j < 3; ++j )
345350
{
346351
p2t::Point *p = t->GetPoint( j );
347-
float zPt = z[p];
348-
QVector3D nPoint = pOrigin + pXVector * p->x + pYVector * p->y;
349-
float fx = nPoint.x() - mOriginX;
350-
float fy = nPoint.y() - mOriginY;
351-
float fz = extrusionHeight + ( std::isnan( zPt ) ? 0 : zPt );
352+
const double zPt = z[p];
353+
QVector3D nPoint = pXVector * p->x + pYVector * p->y;
354+
const double fx = nPoint.x() - mOriginX + x0;
355+
const double fy = nPoint.y() - mOriginY + y0;
356+
const double fz = extrusionHeight + ( std::isnan( zPt ) ? 0 : zPt );
352357
mData << fx << fz << -fy;
353358
if ( mAddNormals )
354359
mData << pNormal.x() << pNormal.z() << - pNormal.y();

0 commit comments

Comments
 (0)
Please sign in to comment.