fix clipping for geocentric projections

This commit is contained in:
bdm-oslandia 2021-08-31 14:10:55 +02:00 committed by Nyall Dawson
parent 349eb7cb1e
commit 8e6a8fc8bb
12 changed files with 600 additions and 54 deletions

View File

@ -235,6 +235,20 @@ Returns the smallest distance between the box and the point ``point``
bool operator==( const QgsBox3d &other ) const;
void scale( double scaleFactor, const QgsPoint *c = 0 );
%Docstring
Scale the rectangle around a center :py:class:`QgsPoint`.
.. versionadded:: 3.21
%End
void scale( double scaleFactor, double centerX, double centerY, double centerZ );
%Docstring
Scale the rectangle around a center coordinates.
.. versionadded:: 3.21
%End
};
/************************************************************************

View File

@ -822,6 +822,13 @@ corresponds to the last point in the line.
%End
QgsBox3d calculateBoundingBox3d() const;
%Docstring
Calculator for the minimal 3d bounding box for the geometry.
.. seealso:: :py:func:`calculateBoundingBox`
%End
protected:
int compareToSameClass( const QgsAbstractGeometry *other ) const final;

View File

@ -11,6 +11,7 @@
%Feature ARM // Some parts are not available in sip bindings on ARM because of qreal double vs. float issues
@ -45,7 +46,9 @@ magnitude of the screen coordinates (+/- 32768, i.e. 16 bit integer).
XMax,
XMin,
YMax,
YMin
YMin,
ZMax,
ZMin
};
%If (!ARM) // Not available on ARM sip bindings because of qreal issues
@ -68,19 +71,28 @@ an open shape (linestring).
static void trimPolygon( QPolygonF &pts, const QgsRectangle &clipRect );
static QPolygonF clippedLine( const QgsCurve &curve, const QgsRectangle &clipExtent );
static void trimPolygon( QgsLineString &pts, const QgsBox3d &clipRect );
static QgsLineString clipped3dLine( const QgsCurve &curve, const QgsBox3d &clipExtent );
%Docstring
Takes a linestring and clips it to clipExtent
Takes a ``curve`` with 3D coordinates and clips it to clipExtent
:param curve: the linestring
:param clipExtent: clipping bounds
:return: clipped line coordinates
.. versionadded:: 3.21
%End
static QPolygonF clippedLine( const QPolygonF &curve, const QgsRectangle &clipExtent );
%Docstring
Takes a ``curve`` and clips it to clipExtent.
Takes a 2D ``curve`` and clips it to clipExtent.
:param curve: the linestring
:param clipExtent: clipping bounds
:return: clipped line coordinates
.. versionadded:: 3.16
%End
@ -109,6 +121,7 @@ An implementation of the 'Fast clipping' algorithm (Sobkow et al. 1987, Computer
/************************************************************************
* This file has been generated automatically from *
* *

View File

@ -133,3 +133,38 @@ bool QgsBox3d::operator==( const QgsBox3d &other ) const
qgsDoubleNear( mZmin, other.mZmin ) &&
qgsDoubleNear( mZmax, other.mZmax );
}
void QgsBox3d::scale( double scaleFactor, const QgsPoint *c )
{
// scale from the center
double centerX, centerY, centerZ;
if ( c )
{
centerX = c->x();
centerY = c->y();
centerZ = c->z();
}
else
{
centerX = ( xMinimum() + xMaximum() ) / 2;
centerY = ( yMinimum() + yMaximum() ) / 2;
centerZ = ( zMinimum() + zMaximum() ) / 2;
}
scale( scaleFactor, centerX, centerY, centerZ );
}
void QgsBox3d::scale( double scaleFactor, double centerX, double centerY, double centerZ )
{
double newX = ( xMaximum() - xMinimum() ) * scaleFactor * 0.5;
double newY = ( yMaximum() - yMinimum() ) * scaleFactor * 0.5;
double newZ = ( zMaximum() - zMinimum() ) * scaleFactor * 0.5;
setXMinimum( centerX - newX );
setXMaximum( centerX + newX );
setYMinimum( centerY - newY );
setYMaximum( centerY + newY );
setZMinimum( centerZ - newZ );
setZMaximum( centerZ + newZ );
}

View File

@ -215,6 +215,20 @@ class CORE_EXPORT QgsBox3d
bool operator==( const QgsBox3d &other ) const;
/**
* Scale the rectangle around a center QgsPoint.
*
* \since QGIS 3.21
*/
void scale( double scaleFactor, const QgsPoint *c = nullptr );
/**
* Scale the rectangle around a center coordinates.
*
* \since QGIS 3.21
*/
void scale( double scaleFactor, double centerX, double centerY, double centerZ );
private:
QgsRectangle mBounds2d;

View File

@ -34,6 +34,7 @@
#include <QDomDocument>
#include <QJsonObject>
#include "qgsbox3d.h"
/***************************************************************************
* This class is considered CRITICAL and any change MUST be accompanied with
@ -574,9 +575,15 @@ bool QgsLineString::fromWkb( QgsConstWkbPtr &wkbPtr )
}
QgsRectangle QgsLineString::calculateBoundingBox() const
{
QgsBox3d b = calculateBoundingBox3d();
return QgsRectangle( b.xMinimum(), b.yMinimum(), b.xMaximum(), b.yMaximum(), false );
}
QgsBox3d QgsLineString::calculateBoundingBox3d() const
{
if ( mX.empty() )
return QgsRectangle();
return QgsBox3d();
auto result = std::minmax_element( mX.begin(), mX.end() );
const double xmin = *result.first;
@ -584,7 +591,19 @@ QgsRectangle QgsLineString::calculateBoundingBox() const
result = std::minmax_element( mY.begin(), mY.end() );
const double ymin = *result.first;
const double ymax = *result.second;
return QgsRectangle( xmin, ymin, xmax, ymax, false );
if ( is3D() )
{
result = std::minmax_element( mZ.begin(), mZ.end() );
const double zmin = *result.first;
const double zmax = *result.second;
return QgsBox3d( xmin, ymin, zmin, xmax, ymax, zmax );
}
else
{
// TODO Nan or Inf?
return QgsBox3d( xmin, ymin, -HUGE_VAL, xmax, ymax, HUGE_VAL );
}
}
void QgsLineString::scroll( int index )

View File

@ -27,6 +27,7 @@
#include "qgscompoundcurve.h"
class QgsLineSegment2D;
class QgsBox3d;
/***************************************************************************
* This class is considered CRITICAL and any change MUST be accompanied with
@ -984,6 +985,12 @@ class CORE_EXPORT QgsLineString: public QgsCurve
#endif
/**
* Calculator for the minimal 3d bounding box for the geometry.
* \see calculateBoundingBox()
*/
QgsBox3d calculateBoundingBox3d() const;
protected:
int compareToSameClass( const QgsAbstractGeometry *other ) const final;

View File

@ -18,7 +18,6 @@
#include "qgsclipper.h"
#include "qgsgeometry.h"
#include "qgscurve.h"
#include "qgslogger.h"
// Where has all the code gone?
@ -38,9 +37,68 @@ const double QgsClipper::MIN_Y = -16000;
const double QgsClipper::SMALL_NUM = 1e-12;
QPolygonF QgsClipper::clippedLine( const QgsCurve &curve, const QgsRectangle &clipExtent )
QgsLineString QgsClipper::clipped3dLine( const QgsCurve &curve, const QgsBox3d &clipExtent )
{
return clippedLine( curve.asQPolygonF(), clipExtent );
const int nPoints = curve.numPoints();
double p0x, p0y, p0z, p1x = 0.0, p1y = 0.0, p1z = 0.0; //original coordinates
double p1x_c, p1y_c, p1z_c; //clipped end coordinates
double lastClipX = 0.0, lastClipY = 0., lastClipZ = 0.0; //last successfully clipped coords
QgsLineString line;
QgsPoint curPoint;
QgsVertexId::VertexType type;
for ( int i = 0; i < nPoints; ++i )
{
curve.pointAt( i, curPoint, type );
if ( i == 0 )
{
p1x = curPoint.x();
p1y = curPoint.y();
p1z = curPoint.z();
}
else
{
p0x = p1x;
p0y = p1y;
p0z = p1z;
p1x = curPoint.x();
p1y = curPoint.y();
p1z = curPoint.z();
p1x_c = p1x;
p1y_c = p1y;
p1z_c = p1z;
// TODO: should be in 3D
if ( clipLineSegment( clipExtent, p0x, p0y, p0z, p1x_c, p1y_c, p1z_c ) )
{
bool newLine = !line.isEmpty() && ( !qgsDoubleNear( p0x, lastClipX )
|| !qgsDoubleNear( p0y, lastClipY )
|| !qgsDoubleNear( p0z, lastClipZ ) );
if ( newLine )
{
//add edge points to connect old and new line
// TODO: should be (really) in 3D
connectSeparatedLines( lastClipX, lastClipY, lastClipZ, p0x, p0y, p0z, clipExtent, line );
}
if ( line.numPoints() == 0 || newLine )
{
//add first point
line.addVertex( QgsPoint( p0x, p0y, p0z ) );
}
//add second point
lastClipX = p1x_c;
lastClipY = p1y_c;
lastClipZ = p1z_c;
line.addVertex( QgsPoint( p1x_c, p1y_c, p1z_c ) );
}
}
}
return line;
}
QPolygonF QgsClipper::clippedLine( const QPolygonF &curve, const QgsRectangle &clipExtent )
@ -74,7 +132,7 @@ QPolygonF QgsClipper::clippedLine( const QPolygonF &curve, const QgsRectangle &c
p1x_c = p1x;
p1y_c = p1y;
if ( clipLineSegment( clipExtent.xMinimum(), clipExtent.xMaximum(), clipExtent.yMinimum(), clipExtent.yMaximum(),
p0x, p0y, p1x_c, p1y_c ) )
p0x, p0y, p1x_c, p1y_c ) )
{
const bool newLine = !line.isEmpty() && ( !qgsDoubleNear( p0x, lastClipX ) || !qgsDoubleNear( p0y, lastClipY ) );
if ( newLine )
@ -196,3 +254,105 @@ void QgsClipper::connectSeparatedLines( double x0, double y0, double x1, double
}
}
}
void QgsClipper::connectSeparatedLines( double x0, double y0, double z0, double x1, double y1, double z1,
const QgsBox3d &clipRect, QgsLineString &pts )
{
// TODO: really relevant and sufficient?
double meanZ = ( z0 + z1 ) / 2.0;
//test the different edge combinations...
if ( qgsDoubleNear( x0, clipRect.xMinimum() ) )
{
if ( qgsDoubleNear( x1, clipRect.xMinimum() ) )
{
return;
}
else if ( qgsDoubleNear( y1, clipRect.yMaximum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMinimum(), clipRect.yMaximum(), meanZ ) );
return;
}
else if ( qgsDoubleNear( x1, clipRect.xMaximum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMinimum(), clipRect.yMaximum(), meanZ ) );
pts.addVertex( QgsPoint( clipRect.xMaximum(), clipRect.yMaximum(), meanZ ) );
return;
}
else if ( qgsDoubleNear( y1, clipRect.yMinimum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMinimum(), clipRect.yMinimum(), meanZ ) );
return;
}
}
else if ( qgsDoubleNear( y0, clipRect.yMaximum() ) )
{
if ( qgsDoubleNear( y1, clipRect.yMaximum() ) )
{
return;
}
else if ( qgsDoubleNear( x1, clipRect.xMaximum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMaximum(), clipRect.yMaximum(), meanZ ) );
return;
}
else if ( qgsDoubleNear( y1, clipRect.yMinimum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMaximum(), clipRect.yMaximum(), meanZ ) );
pts.addVertex( QgsPoint( clipRect.xMaximum(), clipRect.yMinimum(), meanZ ) );
return;
}
else if ( qgsDoubleNear( x1, clipRect.xMinimum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMinimum(), clipRect.yMaximum(), meanZ ) );
return;
}
}
else if ( qgsDoubleNear( x0, clipRect.xMaximum() ) )
{
if ( qgsDoubleNear( x1, clipRect.xMaximum() ) )
{
return;
}
else if ( qgsDoubleNear( y1, clipRect.yMinimum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMaximum(), clipRect.yMinimum(), meanZ ) );
return;
}
else if ( qgsDoubleNear( x1, clipRect.xMinimum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMaximum(), clipRect.yMinimum(), meanZ ) );
pts.addVertex( QgsPoint( clipRect.xMinimum(), clipRect.yMinimum(), meanZ ) );
return;
}
else if ( qgsDoubleNear( y1, clipRect.yMaximum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMaximum(), clipRect.yMaximum(), meanZ ) );
return;
}
}
else if ( qgsDoubleNear( y0, clipRect.yMinimum() ) )
{
if ( qgsDoubleNear( y1, clipRect.yMinimum() ) )
{
return;
}
else if ( qgsDoubleNear( x1, clipRect.xMinimum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMinimum(), clipRect.yMinimum(), meanZ ) );
return;
}
else if ( qgsDoubleNear( y1, clipRect.yMaximum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMinimum(), clipRect.yMinimum(), meanZ ) );
pts.addVertex( QgsPoint( clipRect.xMinimum(), clipRect.yMaximum(), meanZ ) );
return;
}
else if ( qgsDoubleNear( x1, clipRect.xMaximum() ) )
{
pts.addVertex( QgsPoint( clipRect.xMaximum(), clipRect.yMinimum(), meanZ ) );
return;
}
}
}

View File

@ -27,7 +27,10 @@
#include <QVector>
#include <QPolygonF>
#include "qgsbox3d.h"
class QgsCurve;
class QgsLineString;
SIP_FEATURE( ARM ) // Some parts are not available in sip bindings on ARM because of qreal double vs. float issues
@ -80,7 +83,9 @@ class CORE_EXPORT QgsClipper
XMax,
XMin,
YMax,
YMin
YMin,
ZMax,
ZMin
};
SIP_IF_FEATURE( !ARM ) // Not available on ARM sip bindings because of qreal issues
@ -101,17 +106,22 @@ class CORE_EXPORT QgsClipper
static void trimPolygon( QPolygonF &pts, const QgsRectangle &clipRect );
static void trimPolygon( QgsLineString &pts, const QgsBox3d &clipRect );
/**
* Takes a linestring and clips it to clipExtent
* Takes a \a curve with 3D coordinates and clips it to clipExtent
* \param curve the linestring
* \param clipExtent clipping bounds
* \returns clipped line coordinates
* \since QGIS 3.21
*/
static QPolygonF clippedLine( const QgsCurve &curve, const QgsRectangle &clipExtent );
static QgsLineString clipped3dLine( const QgsCurve &curve, const QgsBox3d &clipExtent );
/**
* Takes a \a curve and clips it to clipExtent.
*
* Takes a 2D \a curve and clips it to clipExtent.
* \param curve the linestring
* \param clipExtent clipping bounds
* \returns clipped line coordinates
* \since QGIS 3.16
*/
static QPolygonF clippedLine( const QPolygonF &curve, const QgsRectangle &clipExtent );
@ -151,11 +161,15 @@ class CORE_EXPORT QgsClipper
static inline void trimPolygonToBoundary( const QPolygonF &inPts, QPolygonF &outPts, const QgsRectangle &rect, Boundary b, double boundaryValue );
static inline void trimPolygonToBoundary( const QgsLineString &inPts, QgsLineString &outPts, const QgsBox3d &rect, Boundary b, double boundaryValue );
// Determines if a point is inside or outside the given boundary
static inline bool inside( double x, double y, Boundary b );
static inline bool inside( QPointF pt, Boundary b, double val );
static inline bool inside( QgsPoint pt, Boundary b, double val );
// Calculates the intersection point between a line defined by a
// (x1, y1), and (x2, y2) and the given boundary
static inline QgsPointXY intersect( double x1, double y1,
@ -166,6 +180,15 @@ class CORE_EXPORT QgsClipper
QPointF pt2,
Boundary b, const QgsRectangle &rect );
static inline QgsPoint intersectRect( QgsPoint pt1,
QgsPoint pt2,
Boundary b, const QgsBox3d &rect );
//Implementation of 'Fast clipping' algorithm (Sobkow et al. 1987, Computers & Graphics Vol.11, 4, p.459-467)
static bool clipLineSegment( double xLeft, double xRight, double yBottom, double yTop, double &x0, double &y0, double &x1, double &y1 );
static bool clipLineSegment( const QgsBox3d &extent, double &x0, double &y0, double &z0, double &x1, double &y1, double &z1 );
/**
* Connects two lines split by the clip (by inserting points on the clip border)
* \param x0 x-coordinate of the first line end
@ -178,6 +201,20 @@ class CORE_EXPORT QgsClipper
static void connectSeparatedLines( double x0, double y0, double x1, double y1,
const QgsRectangle &clipRect, QPolygonF &pts );
/**
* Connects two 3D lines split by the clip (by inserting points on the clip border)
* \param x0 x-coordinate of the first line end
* \param y0 y-coordinate of the first line end
* \param z0 z-coordinate of the first line end
* \param x1 x-coordinate of the second line start
* \param y1 y-coordinate of the second line start
* \param z1 z-coordinate of the second line start
* \param clipRect clip box
* \param pts: in/out array of clipped points
*/
static void connectSeparatedLines( double x0, double y0, double z0, double x1, double y1, double z1,
const QgsBox3d &clipRect, QgsLineString &pts );
//low level clip methods for fast clip algorithm
static void clipStartTop( double &x0, double &y0, double x1, double y1, double yMax );
static void clipStartBottom( double &x0, double &y0, double x1, double y1, double yMin );
@ -189,6 +226,9 @@ class CORE_EXPORT QgsClipper
static void clipEndLeft( double x0, double y0, double &x1, double &y1, double xMin );
};
#include "qgslinestring.h"
#include "qgspoint.h"
#ifndef SIP_RUN
// The inline functions
@ -237,6 +277,26 @@ inline void QgsClipper::trimPolygon( QPolygonF &pts, const QgsRectangle &clipRec
trimPolygonToBoundary( tmpPts, pts, clipRect, YMin, clipRect.yMinimum() );
}
inline void QgsClipper::trimPolygon( QgsLineString &pts, const QgsBox3d &clipRect )
{
QgsLineString tmpPts;
trimPolygonToBoundary( pts, tmpPts, clipRect, XMax, clipRect.xMaximum() );
pts.clear();
trimPolygonToBoundary( tmpPts, pts, clipRect, YMax, clipRect.yMaximum() );
tmpPts.clear();
trimPolygonToBoundary( pts, tmpPts, clipRect, XMin, clipRect.xMinimum() );
pts.clear();
trimPolygonToBoundary( tmpPts, pts, clipRect, YMin, clipRect.yMinimum() );
if ( !clipRect.is2d() )
{
tmpPts.clear();
trimPolygonToBoundary( pts, tmpPts, clipRect, ZMax, clipRect.zMaximum() );
pts.clear();
trimPolygonToBoundary( tmpPts, pts, clipRect, ZMin, clipRect.zMinimum() );
}
}
// An auxiliary function that is part of the polygon trimming
// code. Will trim the given polygon to the given boundary and return
// the trimmed polygon in the out pointer. Uses Sutherland and
@ -342,6 +402,44 @@ inline void QgsClipper::trimPolygonToBoundary( const QPolygonF &inPts, QPolygonF
}
}
inline void QgsClipper::trimPolygonToBoundary( const QgsLineString &inPts, QgsLineString &outPts, const QgsBox3d &rect, Boundary b, double boundaryValue )
{
QgsVertexId::VertexType type;
QgsPoint inI1, inI2;
inPts.pointAt( inPts.numPoints() - 1, inI1, type ); // start with last point
// and compare to the first point initially.
for ( int i2 = 0; i2 < inPts.numPoints() ; ++i2 )
{
inPts.pointAt( i2, inI2, type );
// look at each edge of the polygon in turn
if ( inside( inI2, b, boundaryValue ) ) // end point of edge is inside boundary
{
if ( inside( inI1, b, boundaryValue ) )
{
outPts.addVertex( inI2 );
}
else
{
// edge crosses into the boundary, so trim back to the boundary, and
// store both ends of the new edge
outPts.addVertex( intersectRect( inI1, inI2, b, rect ) );
outPts.addVertex( inI2 );
}
}
else // end point of edge is outside boundary
{
// start point is in boundary, so need to trim back
if ( inside( inI1, b, boundaryValue ) )
{
outPts.addVertex( intersectRect( inI1, inI2, b, rect ) );
}
}
inPts.pointAt( i2, inI1, type );
}
}
// An auxiliary function to trimPolygonToBoundarY() that returns
// whether a point is inside or outside the given boundary.
@ -365,6 +463,8 @@ inline bool QgsClipper::inside( const double x, const double y, Boundary b )
if ( y > MIN_Y )
return true;
break;
default:
return false;
}
return false;
}
@ -381,6 +481,28 @@ inline bool QgsClipper::inside( QPointF pt, Boundary b, double val )
return ( pt.y() < val );
case YMin: // y > MIN_Y is inside
return ( pt.y() > val );
default:
return false;
}
return false;
}
inline bool QgsClipper::inside( QgsPoint pt, Boundary b, double val )
{
switch ( b )
{
case XMax: // x < MAX_X is inside
return ( pt.x() < val );
case XMin: // x > MIN_X is inside
return ( pt.x() > val );
case YMax: // y < MAX_Y is inside
return ( pt.y() < val );
case YMin: // y > MIN_Y is inside
return ( pt.y() > val );
case ZMax: // z < MAX_Z is inside
return ( pt.z() < val );
case ZMin: // z > MIN_Z is inside
return ( pt.z() > val );
}
return false;
}
@ -418,6 +540,8 @@ inline QgsPointXY QgsClipper::intersect( const double x1, const double y1,
r_n = ( y1 - MIN_Y ) * ( MAX_X - MIN_X );
r_d = -( y2 - y1 ) * ( MAX_X - MIN_X );
break;
default:
break;
}
QgsPointXY p;
@ -468,6 +592,8 @@ inline QPointF QgsClipper::intersectRect( QPointF pt1,
r_n = ( y1 - rect.yMinimum() ) * ( rect.xMaximum() - rect.xMinimum() );
r_d = -( y2 - y1 ) * ( rect.xMaximum() - rect.xMinimum() );
break;
default:
break;
}
double r = 0;
@ -478,6 +604,55 @@ inline QPointF QgsClipper::intersectRect( QPointF pt1,
return QPointF( x1 + r * ( x2 - x1 ), y1 + r * ( y2 - y1 ) );
}
inline QgsPoint QgsClipper::intersectRect( QgsPoint pt1,
QgsPoint pt2,
Boundary b, const QgsBox3d &rect )
{
// This function assumes that the two given points (x1, y1), and
// (x2, y2) cross the given boundary. Making this assumption allows
// some optimisations.
double r_n = SMALL_NUM, r_d = SMALL_NUM;
const double x1 = pt1.x(), x2 = pt2.x();
const double y1 = pt1.y(), y2 = pt2.y();
const double z1 = pt1.z(), z2 = pt2.z();
switch ( b )
{
case XMax: // x = MAX_X boundary
r_n = -( x1 - rect.xMaximum() ) * ( rect.yMaximum() - rect.yMinimum() );
r_d = ( x2 - x1 ) * ( rect.yMaximum() - rect.yMinimum() );
break;
case XMin: // x = MIN_X boundary
r_n = -( x1 - rect.xMinimum() ) * ( rect.yMaximum() - rect.yMinimum() );
r_d = ( x2 - x1 ) * ( rect.yMaximum() - rect.yMinimum() );
break;
case YMax: // y = MAX_Y boundary
r_n = ( y1 - rect.yMaximum() ) * ( rect.xMaximum() - rect.xMinimum() );
r_d = -( y2 - y1 ) * ( rect.xMaximum() - rect.xMinimum() );
break;
case YMin: // y = MIN_Y boundary
r_n = ( y1 - rect.yMinimum() ) * ( rect.xMaximum() - rect.xMinimum() );
r_d = -( y2 - y1 ) * ( rect.xMaximum() - rect.xMinimum() );
break;
case ZMax: // z = MAX_Z boundary
r_n = ( z1 - rect.zMaximum() );
r_d = -( z2 - z1 );
break;
case ZMin: // z = MIN_Z boundary
r_n = ( z1 - rect.zMinimum() );
r_d = -( z2 - z1 );
break;
}
double r = 0;
if ( !qgsDoubleNear( r_d, 0.0 ) )
{
r = r_n / r_d;
}
return QgsPoint( x1 + r * ( x2 - x1 ), y1 + r * ( y2 - y1 ), z1 + r * ( z2 - z1 ) );
}
inline void QgsClipper::clipStartTop( double &x0, double &y0, double x1, double y1, double yMax )
{
x0 += ( x1 - x0 ) * ( yMax - y0 ) / ( y1 - y0 );
@ -526,6 +701,13 @@ inline void QgsClipper::clipEndLeft( double x0, double y0, double &x1, double &y
x1 = xMin;
}
inline bool QgsClipper::clipLineSegment( const QgsBox3d &clipExtent, double &x0, double &y0, double &, double &x1, double &y1, double & )
{
// TODO: while waiting for a valid 3d solution, using the 2d one:
return clipLineSegment( clipExtent.xMinimum(), clipExtent.xMaximum(), clipExtent.yMinimum(), clipExtent.yMaximum(),
x0, y0, x1, y1 );
}
//'Fast clipping' algorithm (Sobkow et al. 1987, Computers & Graphics Vol.11, 4, p.459-467)
inline bool QgsClipper::clipLineSegment( double xLeft, double xRight, double yBottom, double yTop, double &x0, double &y0, double &x1, double &y1 )
{

View File

@ -83,11 +83,19 @@ Q_NOWARN_DEPRECATED_POP
QPolygonF QgsSymbol::_getLineString( QgsRenderContext &context, const QgsCurve &curve, bool clipToExtent )
{
printf( "_getLineString(in): %d\n", curve.numPoints() );
for ( int i = 0 ; i < 10 && i < curve.numPoints(); i++ )
{
QgsPoint p;
QgsVertexId::VertexType t;
curve.pointAt( i, p, t );
printf( "_getLineString(in): %f, %f\n", p.x(), p.y() );
}
const unsigned int nPoints = curve.numPoints();
QgsCoordinateTransform ct = context.coordinateTransform();
const QgsMapToPixel &mtp = context.mapToPixel();
QPolygonF pts;
QgsLineString pts;
//apply clipping for large lines to achieve a better rendering performance
if ( clipToExtent && nPoints > 1 && !( context.flags() & Qgis::RenderContextFlag::ApplyClipAfterReprojection ) )
@ -95,12 +103,15 @@ QPolygonF QgsSymbol::_getLineString( QgsRenderContext &context, const QgsCurve &
const QgsRectangle e = context.extent();
const double cw = e.width() / 10;
const double ch = e.height() / 10;
const QgsRectangle clipRect( e.xMinimum() - cw, e.yMinimum() - ch, e.xMaximum() + cw, e.yMaximum() + ch );
pts = QgsClipper::clippedLine( curve, clipRect );
const QgsBox3d clipRect( QgsPoint( e.xMinimum() - cw, e.yMinimum() - ch ), QgsPoint( e.xMaximum() + cw, e.yMaximum() + ch ) ); // TODO also need to be clipped according to z axis
pts = QgsClipper::clipped3dLine( curve, clipRect );
}
else
{
pts = curve.asQPolygonF();
// clone...
QgsPointSequence seq;
curve.points( seq );
pts.setPoints( seq );
}
//transform the QPolygonF to screen coordinates
@ -108,7 +119,7 @@ QPolygonF QgsSymbol::_getLineString( QgsRenderContext &context, const QgsCurve &
{
try
{
ct.transformPolygon( pts );
pts.transform( ct, QgsCoordinateTransform::TransformDirection::ForwardTransform, true );
}
catch ( QgsCsException & )
{
@ -117,11 +128,17 @@ QPolygonF QgsSymbol::_getLineString( QgsRenderContext &context, const QgsCurve &
}
// remove non-finite points, e.g. infinite or NaN points caused by reprojecting errors
pts.erase( std::remove_if( pts.begin(), pts.end(),
[]( const QPointF point )
pts.filterVertices( []( const QgsPoint point )
{
return !std::isfinite( point.x() ) || !std::isfinite( point.y() );
} ), pts.end() );
if ( point.is3D() )
{
return std::isfinite( point.x() ) && std::isfinite( point.y() ) && std::isfinite( point.z() );
}
else
{
return std::isfinite( point.x() ) && std::isfinite( point.y() );
}
} );
if ( clipToExtent && nPoints > 1 && context.flags() & Qgis::RenderContextFlag::ApplyClipAfterReprojection )
{
@ -129,25 +146,44 @@ QPolygonF QgsSymbol::_getLineString( QgsRenderContext &context, const QgsCurve &
const QgsRectangle e = context.mapExtent();
const double cw = e.width() / 10;
const double ch = e.height() / 10;
const QgsRectangle clipRect( e.xMinimum() - cw, e.yMinimum() - ch, e.xMaximum() + cw, e.yMaximum() + ch );
pts = QgsClipper::clippedLine( pts, clipRect );
const QgsBox3d clipRect( QgsPoint( e.xMinimum() - cw, e.yMinimum() - ch ), QgsPoint( e.xMaximum() + cw, e.yMaximum() + ch ) ); // TODO also need to be clipped according to z axis
pts = QgsClipper::clipped3dLine( pts, clipRect );
}
QPointF *ptr = pts.data();
for ( int i = 0; i < pts.size(); ++i, ++ptr )
QPolygonF out = pts.asQPolygonF();
QPointF *ptr = out.data();
for ( int i = 0; i < out.size(); ++i, ++ptr )
{
mtp.transformInPlace( ptr->rx(), ptr->ry() );
}
return pts;
printf( "_getLineString(out): %d\n", out.size() );
for ( int i = 0 ; i < 10 && i < out.size(); i++ )
{
printf( "_getLineString(out): %f, %f\n", out.at( i ).x(), out.at( i ).y() );
}
return out;
}
QPolygonF QgsSymbol::_getPolygonRing( QgsRenderContext &context, const QgsCurve &curve, const bool clipToExtent, const bool isExteriorRing, const bool correctRingOrientation )
{
printf( "_getPolygonRing(in): %d\n", curve.numPoints() );
for ( int i = 0 ; i < 10 && i < curve.numPoints(); i++ )
{
QgsPoint p;
QgsVertexId::VertexType t;
curve.pointAt( i, p, t );
printf( "_getPolygonRing(in): %f, %f\n", p.x(), p.y() );
}
const QgsCoordinateTransform ct = context.coordinateTransform();
const QgsMapToPixel &mtp = context.mapToPixel();
QPolygonF poly = curve.asQPolygonF();
// clone...
QgsLineString poly;
QgsPointSequence seq;
curve.points( seq );
poly.setPoints( seq );
if ( curve.numPoints() < 1 )
return QPolygonF();
@ -155,19 +191,19 @@ QPolygonF QgsSymbol::_getPolygonRing( QgsRenderContext &context, const QgsCurve
if ( correctRingOrientation )
{
// ensure consistent polygon ring orientation
if ( isExteriorRing && curve.orientation() != Qgis::AngularDirection::Clockwise )
std::reverse( poly.begin(), poly.end() );
else if ( !isExteriorRing && curve.orientation() != Qgis::AngularDirection::CounterClockwise )
std::reverse( poly.begin(), poly.end() );
if ( (isExteriorRing && curve.orientation() != Qgis::AngularDirection::Clockwise) || ( !isExteriorRing && curve.orientation() != Qgis::AngularDirection::CounterClockwise )) {
poly.reversed()->points( seq );
poly.setPoints(seq);
}
}
//clip close to view extent, if needed
if ( clipToExtent && !( context.flags() & Qgis::RenderContextFlag::ApplyClipAfterReprojection ) && !context.extent().contains( poly.boundingRect() ) )
if ( clipToExtent && !( context.flags() & Qgis::RenderContextFlag::ApplyClipAfterReprojection ) && !context.extent().contains( poly.boundingBox() ) )
{
const QgsRectangle e = context.extent();
const double cw = e.width() / 10;
const double ch = e.height() / 10;
const QgsRectangle clipRect( e.xMinimum() - cw, e.yMinimum() - ch, e.xMaximum() + cw, e.yMaximum() + ch );
const QgsBox3d clipRect( QgsPoint( e.xMinimum() - cw, e.yMinimum() - ch ), QgsPoint( e.xMaximum() + cw, e.yMaximum() + ch ) ); // TODO also need to be clipped according to z axis
QgsClipper::trimPolygon( poly, clipRect );
}
@ -176,7 +212,7 @@ QPolygonF QgsSymbol::_getPolygonRing( QgsRenderContext &context, const QgsCurve
{
try
{
ct.transformPolygon( poly );
poly.transform( ct, QgsCoordinateTransform::TransformDirection::ForwardTransform, true );
}
catch ( QgsCsException & )
{
@ -185,32 +221,45 @@ QPolygonF QgsSymbol::_getPolygonRing( QgsRenderContext &context, const QgsCurve
}
// remove non-finite points, e.g. infinite or NaN points caused by reprojecting errors
poly.erase( std::remove_if( poly.begin(), poly.end(),
[]( const QPointF point )
poly.filterVertices( []( const QgsPoint point )
{
return !std::isfinite( point.x() ) || !std::isfinite( point.y() );
} ), poly.end() );
if ( point.is3D() )
{
return std::isfinite( point.x() ) && std::isfinite( point.y() ) && std::isfinite( point.z() );
}
else
{
return std::isfinite( point.x() ) && std::isfinite( point.y() );
}
} );
if ( clipToExtent && context.flags() & Qgis::RenderContextFlag::ApplyClipAfterReprojection && !context.mapExtent().contains( poly.boundingRect() ) )
if ( clipToExtent && context.flags() & Qgis::RenderContextFlag::ApplyClipAfterReprojection && !context.mapExtent().contains( poly.boundingBox() ) )
{
// early clipping was not possible, so we have to apply it here after transformation
const QgsRectangle e = context.mapExtent();
const double cw = e.width() / 10;
const double ch = e.height() / 10;
const QgsRectangle clipRect( e.xMinimum() - cw, e.yMinimum() - ch, e.xMaximum() + cw, e.yMaximum() + ch );
const QgsBox3d clipRect( QgsPoint( e.xMinimum() - cw, e.yMinimum() - ch ), QgsPoint( e.xMaximum() + cw, e.yMaximum() + ch ) ); // TODO also need to be clipped according to z axis
QgsClipper::trimPolygon( poly, clipRect );
}
QPointF *ptr = poly.data();
for ( int i = 0; i < poly.size(); ++i, ++ptr )
QPolygonF out = poly.asQPolygonF();
QPointF *ptr = out.data();
for ( int i = 0; i < out.size(); ++i, ++ptr )
{
mtp.transformInPlace( ptr->rx(), ptr->ry() );
}
if ( !poly.empty() && !poly.isClosed() )
poly << poly.at( 0 );
if ( !out.empty() && !poly.isClosed() )
out << out.at( 0 );
return poly;
printf( "_getPolygonRing(out): %d\n", out.size() );
for ( int i = 0 ; i < 10 && i < out.size(); i++ )
{
printf( "_getPolygonRing(out): %f, %f\n", out.at( i ).x(), out.at( i ).y() );
}
return out;
}
void QgsSymbol::_getPolygon( QPolygonF &pts, QVector<QPolygonF> &holes, QgsRenderContext &context, const QgsPolygon &polygon, const bool clipToExtent, const bool correctRingOrientation )

View File

@ -562,14 +562,14 @@ bool QgsMapToolIdentify::identifyVectorLayer( QList<QgsMapToolIdentify::Identify
QgsGeometry selectionGeom = geometry;
bool isPointOrRectangle;
QgsPointXY point;
QgsPoint point;
bool isSingleClick = selectionGeom.type() == QgsWkbTypes::PointGeometry;
if ( isSingleClick )
{
isPointOrRectangle = true;
point = selectionGeom.asPoint();
point = selectionGeom.vertexAt( 0 );
commonDerivedAttributes = derivedAttributesForPoint( QgsPoint( point ) );
commonDerivedAttributes = derivedAttributesForPoint( point );
}
else
{
@ -695,18 +695,18 @@ void QgsMapToolIdentify::closestVertexAttributes( const QgsAbstractGeometry &geo
QgsPoint closestPoint = geometry.vertexAt( vId );
QgsPointXY closestPointMapCoords = mCanvas->mapSettings().layerToMapCoordinates( layer, QgsPointXY( closestPoint.x(), closestPoint.y() ) );
QgsPoint closestPointMapCoords = mCanvas->mapSettings().layerToMapCoordinates( layer, closestPoint );
derivedAttributes.insert( tr( "Closest vertex X" ), formatXCoordinate( closestPointMapCoords ) );
derivedAttributes.insert( tr( "Closest vertex Y" ), formatYCoordinate( closestPointMapCoords ) );
if ( closestPoint.is3D() )
{
str = QLocale().toString( closestPoint.z(), 'g', 10 );
str = QLocale().toString( closestPointMapCoords.z(), 'g', 10 );
derivedAttributes.insert( tr( "Closest vertex Z" ), str );
}
if ( closestPoint.isMeasure() )
{
str = QLocale().toString( closestPoint.m(), 'g', 10 );
str = QLocale().toString( closestPointMapCoords.m(), 'g', 10 );
derivedAttributes.insert( tr( "Closest vertex M" ), str );
}

View File

@ -34,8 +34,11 @@ class TestQgsClipper: public QObject
void init() {} // will be called before each testfunction is executed.
void cleanup() {} // will be called after every testfunction.
void basic();
void basicWithZ();
private:
bool checkBoundingBox( const QPolygonF &polygon, const QgsRectangle &clipRect );
bool checkBoundingBox( const QgsLineString &polygon, const QgsBox3d &clipRect );
};
void TestQgsClipper::initTestCase()
@ -43,6 +46,44 @@ void TestQgsClipper::initTestCase()
}
void TestQgsClipper::basicWithZ()
{
// QgsClipper is static only
QgsLineString polygon;
polygon.addVertex( QgsPoint( 10.4, 20.5, 10.3 ) );
polygon.addVertex( QgsPoint( 20.2, 30.2, 20.3 ) );
QgsBox3d clipRect( 10, 10, 10, 25, 30, 20 );
QgsClipper::trimPolygon( polygon, clipRect );
// Check nothing sticks out.
QVERIFY( checkBoundingBox( polygon, clipRect ) );
// Check that it didn't clip too much
QgsBox3d clipRectInner( clipRect );
clipRectInner.scale( 0.999 );
QVERIFY( ! checkBoundingBox( polygon, clipRectInner ) );
// A more complex example
polygon.clear();
polygon.addVertex( QgsPoint( 1.0, 9.0, 1.0 ) );
polygon.addVertex( QgsPoint( 11.0, 11.0, 11.0 ) );
polygon.addVertex( QgsPoint( 9.0, 1.0, 9.0 ) );
clipRect = QgsBox3d( 0.0, 0.0, 0.0, 10.0, 10.0, 10.0 );
QgsClipper::trimPolygon( polygon, clipRect );
// We should have 5 vertices now?
QCOMPARE( polygon.numPoints(), 5 );
// Check nothing sticks out.
QVERIFY( checkBoundingBox( polygon, clipRect ) );
// Check that it didn't clip too much
clipRectInner = clipRect;
clipRectInner.scale( 0.999 );
QVERIFY( ! checkBoundingBox( polygon, clipRectInner ) );
}
void TestQgsClipper::basic()
{
// QgsClipper is static only
@ -78,6 +119,11 @@ void TestQgsClipper::basic()
QVERIFY( ! checkBoundingBox( polygon, clipRectInner ) );
}
bool TestQgsClipper::checkBoundingBox( const QgsLineString &polygon, const QgsBox3d &clipRect )
{
return clipRect.contains( polygon.calculateBoundingBox3d() );
}
bool TestQgsClipper::checkBoundingBox( const QPolygonF &polygon, const QgsRectangle &clipRect )
{
const QgsRectangle bBox( polygon.boundingRect() );