sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-04-13 00:03:09 -04:00
Code Issues Packages Projects Releases Wiki Activity

Remove Point3D

Browse Source 
... in favor of QgsPoint
Reference https://github.com/qgis/qgis3.0_api/issues/71
This commit is contained in:
Matthias Kuhn 2017-06-02 20:52:57 +02:00
parent bf9021ae92
commit 0f2a49d9b5
44 changed files with 619 additions and 808 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
python/analysis/analysis.sip
View File

@ -23,7 +23,6 @@
%Include interpolation/Node.sip
%Include interpolation/NormVecDecorator.sip
%Include interpolation/ParametricLine.sip
%Include interpolation/Point3D.sip
%Include interpolation/TriangleInterpolator.sip
%Include interpolation/Triangulation.sip
%Include interpolation/TriDecorator.sip

12
python/analysis/interpolation/Bezier3D.sip
View File

@ -8,7 +8,7 @@ class Bezier3D : ParametricLine
/** Default constructor*/
Bezier3D();
/** Constructor, par is a pointer to the parent, controlpoly a controlpolygon*/
Bezier3D( ParametricLine *par, QVector<Point3D*> *controlpoly );
Bezier3D( ParametricLine *par, QVector<QgsPoint*> *controlpoly );
/** Destructor*/
virtual ~Bezier3D();
/** Do not use this method, since a Bezier curve does not consist of other curves*/
@ -17,9 +17,9 @@ class Bezier3D : ParametricLine
virtual void calcFirstDer( float t, Vector3D *v );
/** Calculates the second derivative and assigns it to v*/
virtual void calcSecDer( float t, Vector3D *v );
//virtual Point3D calcPoint(float t);
//virtual QgsPoint calcPoint(float t);
/** Calculates the point on the curve and assigns it to p*/
virtual void calcPoint( float t, Point3D *p );
virtual void calcPoint( float t, QgsPoint *p );
/** Changes the order of control points*/
virtual void changeDirection();
//virtual void draw(QPainter *p);
@ -27,9 +27,9 @@ class Bezier3D : ParametricLine
/** Do not use this method, since a Bezier curve does not consist of other curves*/
virtual void remove( int i );
/** Returns a control point*/
virtual const Point3D *getControlPoint( int number ) const;
virtual const QgsPoint *getControlPoint( int number ) const;
/** Returns a pointer to the control polygon*/
virtual const QVector<Point3D*> *getControlPoly() const;
virtual const QVector<QgsPoint*> *getControlPoly() const;
/** Returns the degree of the curve*/
virtual int getDegree() const;
/** Returns the parent*/
@ -37,6 +37,6 @@ class Bezier3D : ParametricLine
/** Sets the parent*/
virtual void setParent( ParametricLine *par );
/** Sets the control polygon*/
virtual void setControlPoly( QVector<Point3D*> *cp );
virtual void setControlPoly( QVector<QgsPoint*> *cp );
};

2
python/analysis/interpolation/CloughTocherInterpolator.sip
View File

@ -19,6 +19,6 @@ class CloughTocherInterpolator : TriangleInterpolator
/** Calculates the normal vector and assigns it to vec (not implemented at the moment)*/
virtual bool calcNormVec( double x, double y, Vector3D *result );
/** Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point*/
virtual bool calcPoint( double x, double y, Point3D *result );
virtual bool calcPoint( double x, double y, QgsPoint *result );
virtual void setTriangulation( NormVecDecorator *tin );
};

14
python/analysis/interpolation/DualEdgeTriangulation.sip
View File

@ -11,23 +11,23 @@ class DualEdgeTriangulation: Triangulation
/** Adds a line (e.g. a break-, structure- or an isoline) to the triangulation. The class takes ownership of the line object and its points*/
void addLine( Line3D *line /Transfer/, bool breakline );
/** Adds a point to the triangulation and returns the number of this point in case of success or -100 in case of failure*/
int addPoint( Point3D *p /Transfer/);
int addPoint( QgsPoint *p /Transfer/);
/** Performs a consistency check, remove this later*/
virtual void performConsistencyTest();
/** Calculates the normal at a point on the surface*/
virtual bool calcNormal( double x, double y, Vector3D *result );
/** Calculates x-, y and z-value of the point on the surface*/
virtual bool calcPoint( double x, double y, Point3D *result );
virtual bool calcPoint( double x, double y, QgsPoint *result );
/** Draws the points, edges and the forced lines*/
//virtual void draw(QPainter *p, double xlowleft, double ylowleft, double xupright, double yupright, double width, double height) const;
/** Returns a pointer to the point with number i*/
virtual Point3D *getPoint( unsigned int i ) const;
virtual QgsPoint *getPoint( unsigned int i ) const;
/** Returns the number of the point opposite to the triangle points p1, p2 (which have to be on a halfedge)*/
int getOppositePoint( int p1, int p2 );
/** Finds out, in which triangle the point with coordinates x and y is and assigns the numbers of the vertices to 'n1', 'n2' and 'n3' and the vertices to 'p1', 'p2' and 'p3'*/
bool getTriangle( double x, double y, Point3D *p1, int *n1 /Out/, Point3D *p2 /Out/, int *n2 /Out/, Point3D *p3 /Out/, int *n3 /Out/ );
bool getTriangle( double x, double y, QgsPoint *p1, int *n1 /Out/, QgsPoint *p2 /Out/, int *n2 /Out/, QgsPoint *p3 /Out/, int *n3 /Out/ );
/** Finds out, in which triangle the point with coordinates x and y is and assigns addresses to the points at the vertices to 'p1', 'p2' and 'p3*/
bool getTriangle( double x, double y, Point3D *p1 /Out/, Point3D *p2 /Out/, Point3D *p3 /Out/ );
bool getTriangle( double x, double y, QgsPoint *p1 /Out/, QgsPoint *p2 /Out/, QgsPoint *p3 /Out/ );
/** Returns a pointer to a value list with the information of the triangles surrounding (counterclockwise) a point. Four integer values describe a triangle, the first three are the number of the half edges of the triangle and the fourth is -10, if the third (and most counterclockwise) edge is a breakline, and -20 otherwise. The value list has to be deleted by the code which called the method*/
QList<int> *getSurroundingTriangles( int pointno );
/** Returns the largest x-coordinate value of the bounding box*/
@ -80,7 +80,7 @@ class DualEdgeTriangulation: Triangulation
/** Returns the number of an edge which points to the point with number 'point' or -1 if there is an error*/
int baseEdgeOfPoint( int point );
/** Returns the number of a HalfEdge from a triangle in which 'point' is in. If the number -10 is returned, this means, that 'point' is outside the convex hull. If -5 is returned, then numerical problems with the leftOfTest occurred (and the value of the possible edge is stored in the variable 'mUnstableEdge'. -20 means, that the inserted point is exactly on an edge (the number is stored in the variable 'mEdgeWithPoint'). -25 means, that the point is already in the triangulation (the number of the point is stored in the member 'mTwiceInsPoint'. If -100 is returned, this means that something else went wrong*/
int baseEdgeOfTriangle( Point3D *point );
int baseEdgeOfTriangle( QgsPoint *point );
/** Checks, if 'edge' has to be swapped because of the empty circle criterion. If so, doSwap(...) is called.*/
bool checkSwap( unsigned int edge, unsigned int recursiveDeep );
/** Swaps 'edge' and test recursively for other swaps (delaunay criterion)*/
@ -100,5 +100,5 @@ class DualEdgeTriangulation: Triangulation
/** Returns true, if a half edge is on the convex hull and false otherwise*/
bool edgeOnConvexHull( int edge );
/** Function needed for the ruppert algorithm. Tests, if point is in the circle through both endpoints of edge and the endpoint of edge->dual->next->point. If so, the function calls itself recursively for edge->next and edge->next->next. Stops, if it finds a forced edge or a convex hull edge*/
void evaluateInfluenceRegion( Point3D *point, int edge, QSet<int> &set );
void evaluateInfluenceRegion( QgsPoint *point, int edge, QSet<int> &set );
};

2
python/analysis/interpolation/LinTriangleInterpolator.sip
View File

@ -13,7 +13,7 @@ class LinTriangleInterpolator : TriangleInterpolator
/** Calculates the normal vector and assigns it to vec*/
virtual bool calcNormVec( double x, double y, Vector3D *result );
/** Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point*/
virtual bool calcPoint( double x, double y, Point3D *result );
virtual bool calcPoint( double x, double y, QgsPoint *result );
/** Returns a pointer to the current Triangulation object*/
virtual DualEdgeTriangulation *getTriangulation() const;
/** Sets a Triangulation*/

6
python/analysis/interpolation/Line3D.sip
View File

@ -6,14 +6,14 @@ class Line3D /NoDefaultCtors/
public:
Line3D();
~Line3D();
/** Returns true, if the Line contains no Point3D, otherwise false*/
/** Returns true, if the Line contains no QgsPoint, otherwise false*/
bool empty() const;
/** Inserts a node behind the current position and sets the current position to this new node*/
void insertPoint( Point3D *p );
void insertPoint( QgsPoint *p );
/** Removes the point behind the current position*/
void removePoint();
/** Gets the point at the current position*/
Point3D *getPoint() const;
QgsPoint *getPoint() const;
/** Returns the current position*/
unsigned int getCurrent() const;
/** Returns the size of the line (the numbero of inserted Nodes without 'head' and 'z'*/

8
python/analysis/interpolation/Node.sip
View File

@ -11,10 +11,10 @@ class Node
// Node &operator=( const Node &n );
/** Returns a pointer to the next element in the linked list*/
Node *getNext() const;
/** Returns a pointer to the Point3D object associated with the node*/
Point3D *getPoint() const;
/** Returns a pointer to the QgsPoint object associated with the node*/
QgsPoint *getPoint() const;
/** Sets the pointer to the next node*/
void setNext( Node *n );
/** Sets a new pointer to an associated Point3D object*/
void setPoint( Point3D *p );
/** Sets a new pointer to an associated QgsPoint object*/
void setPoint( QgsPoint *p );
};

8
python/analysis/interpolation/NormVecDecorator.sip
View File

@ -11,13 +11,13 @@ class NormVecDecorator : TriDecorator
NormVecDecorator( Triangulation *tin );
virtual ~NormVecDecorator();
/** Adds a point to the triangulation*/
int addPoint( Point3D *p );
int addPoint( QgsPoint *p );
/** Calculates the normal at a point on the surface and assigns it to 'result'. Returns true in case of success and false in case of failure*/
bool calcNormal( double x, double y, Vector3D *result );
/** Calculates the normal of a triangle-point for the point with coordinates x and y. This is needed, if a point is on a break line and there is no unique normal stored in 'mNormVec'. Returns false, it something went wrong and true otherwise*/
bool calcNormalForPoint( double x, double y, int point, Vector3D *result );
/** Calculates x-, y and z-value of the point on the surface and assigns it to 'result'. Returns true in case of success and flase in case of failure*/
bool calcPoint( double x, double y, Point3D *result );
bool calcPoint( double x, double y, QgsPoint *result );
/** Eliminates the horizontal triangles by swapping or by insertion of new points. If alreadyestimated is true, a re-estimation of the normals will be done*/
virtual void eliminateHorizontalTriangles();
/** Estimates the first derivative a point. Return true in case of succes and false otherwise*/
@ -27,11 +27,11 @@ class NormVecDecorator : TriDecorator
/** Returns a pointer to the normal vector for the point with the number n*/
Vector3D *getNormal( int n ) const;
/** Finds out, in which triangle a point with coordinates x and y is and assigns the triangle points to p1, p2, p3 and the estimated normals to v1, v2, v3. The vectors are normally taken from 'mNormVec', except if p1, p2 or p3 is a point on a breakline. In this case, the normal is calculated on-the-fly. Returns false, if something went wrong and true otherwise*/
bool getTriangle( double x, double y, Point3D *p1 /Out/, Vector3D *v1 /Out/, Point3D *p2 /Out/, Vector3D *v2 /Out/, Point3D *p3 /Out/, Vector3D *v3 /Out/ );
bool getTriangle( double x, double y, QgsPoint *p1 /Out/, Vector3D *v1 /Out/, QgsPoint *p2 /Out/, Vector3D *v2 /Out/, QgsPoint *p3 /Out/, Vector3D *v3 /Out/ );
/** This function behaves similar to the one above. Additionally, the numbers of the points are returned (ptn1, ptn2, ptn3) as well as the PointStates of the triangle points (state1, state2, state3)
* @note not available in Python bindings
*/
//bool getTriangle( double x, double y, Point3D *p1 /Out/, int *ptn1 /Out/, Vector3D *v1 /Out/, PointState *state1 /Out/, Point3D *p2 /Out/, int *ptn2 /Out/, Vector3D *v2 /Out/, PointState *state2 /Out/, Point3D *p3 /Out/, int *ptn3 /Out/, Vector3D *v3 /Out/, PointState *state3 /Out/);
//bool getTriangle( double x, double y, QgsPoint *p1 /Out/, int *ptn1 /Out/, Vector3D *v1 /Out/, PointState *state1 /Out/, QgsPoint *p2 /Out/, int *ptn2 /Out/, Vector3D *v2 /Out/, PointState *state2 /Out/, QgsPoint *p3 /Out/, int *ptn3 /Out/, Vector3D *v3 /Out/, PointState *state3 /Out/);
/** Returns the state of the point with the number 'pointno'*/
PointState getState( int pointno ) const;
/** Sets an interpolator*/

12
python/analysis/interpolation/ParametricLine.sip
View File

@ -8,22 +8,22 @@ class ParametricLine
ParametricLine();
/** Constructor, par is a pointer to the parent object, controlpoly the controlpolygon
*/
ParametricLine( ParametricLine *par /Transfer/, QVector<Point3D*> *controlpoly );
ParametricLine( ParametricLine *par /Transfer/, QVector<QgsPoint*> *controlpoly );
/** Destructor*/
virtual ~ParametricLine();
virtual void add( ParametricLine *pl ) = 0;
virtual void calcFirstDer( float t, Vector3D *v ) = 0;
virtual void calcSecDer( float t, Vector3D *v ) = 0;
//virtual Point3D calcPoint(float t);
virtual void calcPoint( float t, Point3D* ) = 0;
//virtual QgsPoint calcPoint(float t);
virtual void calcPoint( float t, QgsPoint* ) = 0;
virtual void changeDirection() = 0;
//virtual void draw(QPainter *p);
virtual const Point3D *getControlPoint( int number ) const = 0;
virtual const QVector<Point3D*> *getControlPoly() const = 0;
virtual const QgsPoint *getControlPoint( int number ) const = 0;
virtual const QVector<QgsPoint*> *getControlPoly() const = 0;
virtual int getDegree() const = 0;
virtual ParametricLine *getParent() const = 0;
//virtual bool intersects(ParametricLine *pal);
virtual void remove( int i ) = 0;
virtual void setControlPoly( QVector<Point3D*> *cp ) = 0;
virtual void setControlPoly( QVector<QgsPoint*> *cp ) = 0;
virtual void setParent( ParametricLine *paral ) = 0;
};

29
python/analysis/interpolation/Point3D.sip
View File

@ -1,29 +0,0 @@
class Point3D
{
%TypeHeaderCode
#include <Point3D.h>
%End
public:
Point3D();
/** Constructor with the x-, y- and z-coordinate as arguments*/
Point3D( double x, double y, double z );
Point3D( const Point3D &p );
~Point3D();
// Point3D &operator=( const Point3D &p );
bool operator==( const Point3D &p ) const;
bool operator!=( const Point3D &p ) const;
/** Calculates the three-dimensional distance to another point*/
double dist3D( Point3D *p ) const;
/** Returns the x-coordinate of the point*/
double getX() const;
/** Returns the y-coordinate of the point*/
double getY() const;
/** Returns the z-coordinate of the point*/
double getZ() const;
/** Sets the x-coordinate of the point*/
void setX( double x );
/** Sets the y-coordinate of the point*/
void setY( double y );
/** Sets the z-coordinate of the point*/
void setZ( double z );
};

10
python/analysis/interpolation/TriDecorator.sip
View File

@ -10,17 +10,17 @@ class TriDecorator : Triangulation
explicit TriDecorator( Triangulation *t );
virtual ~TriDecorator();
virtual void addLine( Line3D *line, bool breakline );
virtual int addPoint( Point3D *p );
virtual int addPoint( QgsPoint *p );
/** Adds an association to a triangulation*/
virtual void addTriangulation( Triangulation *t );
/** Performs a consistency check, remove this later*/
virtual void performConsistencyTest();
virtual bool calcNormal( double x, double y, Vector3D *result );
virtual bool calcPoint( double x, double y, Point3D *result );
virtual Point3D *getPoint( unsigned int i ) const;
virtual bool calcPoint( double x, double y, QgsPoint *result );
virtual QgsPoint *getPoint( unsigned int i ) const;
virtual int getNumberOfPoints() const;
bool getTriangle( double x, double y, Point3D *p1, int *n1 /Out/, Point3D *p2 /Out/, int *n2 /Out/, Point3D *p3 /Out/, int *n3 /Out/ );
bool getTriangle( double x, double y, Point3D *p1 /Out/, Point3D *p2 /Out/, Point3D *p3 /Out/ );
bool getTriangle( double x, double y, QgsPoint *p1, int *n1 /Out/, QgsPoint *p2 /Out/, int *n2 /Out/, QgsPoint *p3 /Out/, int *n3 /Out/ );
bool getTriangle( double x, double y, QgsPoint *p1 /Out/, QgsPoint *p2 /Out/, QgsPoint *p3 /Out/ );
virtual int getOppositePoint( int p1, int p2 );
virtual QList<int> *getSurroundingTriangles( int pointno );
virtual double getXMax() const;

2
python/analysis/interpolation/TriangleInterpolator.sip
View File

@ -8,5 +8,5 @@ class TriangleInterpolator
/** Calculates the normal vector and assigns it to vec*/
virtual bool calcNormVec( double x, double y, Vector3D *result ) = 0;
/** Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point*/
virtual bool calcPoint( double x, double y, Point3D *result ) = 0;
virtual bool calcPoint( double x, double y, QgsPoint *result ) = 0;
};

10
python/analysis/interpolation/Triangulation.sip
View File

@ -27,7 +27,7 @@ class Triangulation
* Adds a point to the triangulation
* Ownership is transferred to this class
*/
virtual int addPoint( Point3D *p ) = 0;
virtual int addPoint( QgsPoint *p ) = 0;
/**
* Calculates the normal at a point on the surface and assigns it to 'result'.
@ -42,18 +42,18 @@ class Triangulation
* Calculates x-, y and z-value of the point on the surface and assigns it to 'result'.
* Returns true in case of success and flase in case of failure
*/
virtual bool calcPoint( double x, double y, Point3D *result ) = 0;
virtual bool calcPoint( double x, double y, QgsPoint *result ) = 0;
/** Returns a pointer to the point with number i. Any virtual points must have the number -1*/
virtual Point3D *getPoint( unsigned int i ) const = 0;
virtual QgsPoint *getPoint( unsigned int i ) const = 0;
/** Finds out in which triangle the point with coordinates x and y is and
* assigns the numbers of the vertices to 'n1', 'n2' and 'n3' and the vertices to 'p1', 'p2' and 'p3'
*/
virtual bool getTriangle( double x, double y, Point3D *p1, int *n1 /Out/, Point3D *p2 /Out/, int *n2 /Out/, Point3D *p3 /Out/, int *n3 /Out/ ) = 0;
virtual bool getTriangle( double x, double y, QgsPoint *p1, int *n1 /Out/, QgsPoint *p2 /Out/, int *n2 /Out/, QgsPoint *p3 /Out/, int *n3 /Out/ ) = 0;
/** Finds out, in which triangle the point with coordinates x and y is and assigns the points at the vertices to 'p1', 'p2' and 'p3*/
virtual bool getTriangle( double x, double y, Point3D *p1 /Out/, Point3D *p2 /Out/, Point3D *p3 /Out/ ) = 0;
virtual bool getTriangle( double x, double y, QgsPoint *p1 /Out/, QgsPoint *p2 /Out/, QgsPoint *p3 /Out/ ) = 0;
/** Returns the number of the point opposite to the triangle points p1, p2 (which have to be on a halfedge)*/
virtual int getOppositePoint( int p1, int p2 ) = 0;

1
python/auto_sip.blacklist
View File

@ -212,7 +212,6 @@ analysis/interpolation/Line3D.sip
analysis/interpolation/Node.sip
analysis/interpolation/NormVecDecorator.sip
analysis/interpolation/ParametricLine.sip
analysis/interpolation/Point3D.sip
analysis/interpolation/TriangleInterpolator.sip
analysis/interpolation/Triangulation.sip
analysis/interpolation/TriDecorator.sip

2
src/analysis/CMakeLists.txt
View File

@ -16,7 +16,6 @@ SET(QGIS_ANALYSIS_SRCS
interpolation/NormVecDecorator.cc
interpolation/Node.cc
interpolation/ParametricLine.cc
interpolation/Point3D.cc
interpolation/TriangleInterpolator.cc
interpolation/Triangulation.cc
interpolation/TriDecorator.cc
@ -129,7 +128,6 @@ SET(QGIS_ANALYSIS_HDRS
interpolation/ParametricLine.h
interpolation/CloughTocherInterpolator.h
interpolation/TriangleInterpolator.h
interpolation/Point3D.h
interpolation/Vector3D.h
interpolation/DualEdgeTriangulation.h
interpolation/Node.h

22
src/analysis/interpolation/Bezier3D.cc
View File

@ -35,9 +35,9 @@ void Bezier3D::calcFirstDer( float t, Vector3D *v )
for ( int n = 1; n <= int( mControlPoly->count() - 1 ); n++ )
{
double bernst = MathUtils::calcBernsteinPoly( mControlPoly->count() - 2, n - 1, t );
v->setX( v->getX() + ( ( *mControlPoly )[n]->getX() - ( *mControlPoly )[n - 1]->getX() )*bernst );
v->setY( v->getY() + ( ( *mControlPoly )[n]->getY() - ( *mControlPoly )[n - 1]->getY() )*bernst );
v->setZ( v->getZ() + ( ( *mControlPoly )[n]->getZ() - ( *mControlPoly )[n - 1]->getZ() )*bernst );
v->setX( v->getX() + ( ( *mControlPoly )[n]->x() - ( *mControlPoly )[n - 1]->x() )*bernst );
v->setY( v->getY() + ( ( *mControlPoly )[n]->y() - ( *mControlPoly )[n - 1]->y() )*bernst );
v->setZ( v->getZ() + ( ( *mControlPoly )[n]->z() - ( *mControlPoly )[n - 1]->z() )*bernst );
}
v->setX( v->getX() * ( mControlPoly->count() - 1 ) );
v->setY( v->getY() * ( mControlPoly->count() - 1 ) );
@ -50,7 +50,7 @@ void Bezier3D::calcFirstDer( float t, Vector3D *v )
}
}
void Bezier3D::calcPoint( float t, Point3D *p )
void Bezier3D::calcPoint( float t, QgsPoint *p )
{
if ( p && mControlPoly )
@ -62,9 +62,9 @@ void Bezier3D::calcPoint( float t, Point3D *p )
for ( int n = 1; n <= int( mControlPoly->count() ); n++ )
{
double bernst = MathUtils::calcBernsteinPoly( mControlPoly->count() - 1, n - 1, t );
p->setX( p->getX() + ( *mControlPoly )[n - 1]->getX()*bernst );
p->setY( p->getY() + ( *mControlPoly )[n - 1]->getY()*bernst );
p->setZ( p->getZ() + ( *mControlPoly )[n - 1]->getZ()*bernst );
p->setX( p->x() + ( *mControlPoly )[n - 1]->x()*bernst );
p->setY( p->y() + ( *mControlPoly )[n - 1]->y()*bernst );
p->setZ( p->z() + ( *mControlPoly )[n - 1]->z()*bernst );
}
}
@ -91,9 +91,9 @@ void Bezier3D::calcSecDer( float t, Vector3D *v )
for ( int n = 1; n <= int( nodes - 2 ); n++ )
{
double bernst = MathUtils::calcBernsteinPoly( nodes - 3, n - 1, t );
v->setX( v->getX() + ( ( *mControlPoly )[n + 1]->getX() - 2 * ( *mControlPoly )[n]->getX() + ( *mControlPoly )[n - 1]->getX() )*bernst );
v->setY( v->getY() + ( ( *mControlPoly )[n + 1]->getY() - 2 * ( *mControlPoly )[n]->getY() + ( *mControlPoly )[n - 1]->getY() )*bernst );
v->setZ( v->getZ() + ( ( *mControlPoly )[n + 1]->getZ() - 2 * ( *mControlPoly )[n]->getZ() + ( *mControlPoly )[n - 1]->getZ() )*bernst );
v->setX( v->getX() + ( ( *mControlPoly )[n + 1]->x() - 2 * ( *mControlPoly )[n]->x() + ( *mControlPoly )[n - 1]->x() )*bernst );
v->setY( v->getY() + ( ( *mControlPoly )[n + 1]->y() - 2 * ( *mControlPoly )[n]->y() + ( *mControlPoly )[n - 1]->y() )*bernst );
v->setZ( v->getZ() + ( ( *mControlPoly )[n + 1]->z() - 2 * ( *mControlPoly )[n]->z() + ( *mControlPoly )[n - 1]->z() )*bernst );
}
v->setX( v->getX()*MathUtils::faculty( nodes - 1 ) / MathUtils::faculty( nodes - 3 ) );
v->setY( v->getY()*MathUtils::faculty( nodes - 1 ) / MathUtils::faculty( nodes - 3 ) );
@ -111,7 +111,7 @@ void Bezier3D::changeDirection()//does this work correctly? more testing is need
{
if ( mControlPoly )
{
Point3D **pointer = new Point3D*[mControlPoly->count()];//create an array to temporarily store pointer to the control points
QgsPoint **pointer = new QgsPoint*[mControlPoly->count()];//create an array to temporarily store pointer to the control points
for ( int i = 0; i < mControlPoly->count(); i++ )//store the points
{
pointer[i] = ( *mControlPoly )[i];

20
src/analysis/interpolation/Bezier3D.h
View File

@ -31,7 +31,7 @@ class ANALYSIS_EXPORT Bezier3D: public ParametricLine
//! Default constructor
Bezier3D();
//! Constructor, par is a pointer to the parent, controlpoly a controlpolygon
Bezier3D( ParametricLine *par, QVector<Point3D *> *controlpoly );
Bezier3D( ParametricLine *par, QVector<QgsPoint *> *controlpoly );
virtual ~Bezier3D();
//! Do not use this method, since a Bezier curve does not consist of other curves
@ -40,9 +40,9 @@ class ANALYSIS_EXPORT Bezier3D: public ParametricLine
virtual void calcFirstDer( float t, Vector3D *v ) override;
//! Calculates the second derivative and assigns it to v
virtual void calcSecDer( float t, Vector3D *v ) override;
//virtual Point3D calcPoint(float t);
//virtual QgsPoint calcPoint(float t);
//! Calculates the point on the curve and assigns it to p
virtual void calcPoint( float t, Point3D *p ) override;
virtual void calcPoint( float t, QgsPoint *p ) override;
//! Changes the order of control points
virtual void changeDirection() override;
//virtual void draw(QPainter* p);
@ -50,9 +50,9 @@ class ANALYSIS_EXPORT Bezier3D: public ParametricLine
//! Do not use this method, since a Bezier curve does not consist of other curves
virtual void remove( int i ) override;
//! Returns a control point
virtual const Point3D *getControlPoint( int number ) const override;
virtual const QgsPoint *getControlPoint( int number ) const override;
//! Returns a pointer to the control polygon
virtual const QVector<Point3D *> *getControlPoly() const override;
virtual const QVector<QgsPoint *> *getControlPoly() const override;
//! Returns the degree of the curve
virtual int getDegree() const override;
//! Returns the parent
@ -60,7 +60,7 @@ class ANALYSIS_EXPORT Bezier3D: public ParametricLine
//! Sets the parent
virtual void setParent( ParametricLine *par ) override;
//! Sets the control polygon
virtual void setControlPoly( QVector<Point3D *> *cp ) override;
virtual void setControlPoly( QVector<QgsPoint *> *cp ) override;
};
@ -71,7 +71,7 @@ inline Bezier3D::Bezier3D() : ParametricLine()//default constructor
}
inline Bezier3D::Bezier3D( ParametricLine *parent, QVector<Point3D *> *controlpoly ) : ParametricLine( parent, controlpoly )
inline Bezier3D::Bezier3D( ParametricLine *parent, QVector<QgsPoint *> *controlpoly ) : ParametricLine( parent, controlpoly )
{
mDegree = mControlPoly->count() - 1;
}
@ -97,12 +97,12 @@ inline void Bezier3D::remove( int i )
//-----------------------------------------------setters and getters---------------------------------------------------------------
inline const Point3D *Bezier3D::getControlPoint( int number ) const
inline const QgsPoint *Bezier3D::getControlPoint( int number ) const
{
return ( *mControlPoly )[number - 1];
}
inline const QVector<Point3D *> *Bezier3D::getControlPoly() const
inline const QVector<QgsPoint *> *Bezier3D::getControlPoly() const
{
return mControlPoly;
}
@ -122,7 +122,7 @@ inline void Bezier3D::setParent( ParametricLine *par )
mParent = par;
}
inline void Bezier3D::setControlPoly( QVector<Point3D *> *cp )
inline void Bezier3D::setControlPoly( QVector<QgsPoint *> *cp )
{
mControlPoly = cp;
mDegree = mControlPoly->count() - 1;

264
src/analysis/interpolation/CloughTocherInterpolator.cc
View File

@ -73,25 +73,25 @@ bool CloughTocherInterpolator::calcNormVec( double x, double y, Vector3D *result
if ( result )
{
init( x, y );
Point3D barycoord( 0, 0, 0 );//barycentric coordinates of (x,y) with respect to the triangle
Point3D endpointUXY( 0, 0, 0 );//endpoint of the derivative in u-direction (in xy-coordinates)
Point3D endpointV( 0, 0, 0 );//endpoint of the derivative in v-direction (in barycentric coordinates)
Point3D endpointVXY( 0, 0, 0 );//endpoint of the derivative in v-direction (in xy-coordinates)
QgsPoint barycoord( 0, 0, 0 );//barycentric coordinates of (x,y) with respect to the triangle
QgsPoint endpointUXY( 0, 0, 0 );//endpoint of the derivative in u-direction (in xy-coordinates)
QgsPoint endpointV( 0, 0, 0 );//endpoint of the derivative in v-direction (in barycentric coordinates)
QgsPoint endpointVXY( 0, 0, 0 );//endpoint of the derivative in v-direction (in xy-coordinates)
//find out, in which subtriangle the point (x,y) is
//is the point in the first subtriangle (point1,point2,cp10)?
MathUtils::calcBarycentricCoordinates( x, y, &point1, &point2, &cp10, &barycoord );
if ( barycoord.getX() >= -mEdgeTolerance && barycoord.getX() <= ( 1 + mEdgeTolerance ) && barycoord.getY() >= -mEdgeTolerance && barycoord.getY() <= ( 1 + mEdgeTolerance ) )
if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
{
double zu = point1.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp1.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp2.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp3.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp4.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp7.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double zv = cp1.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp2.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + point2.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp4.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp5.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp8.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double zw = cp3.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp4.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp5.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp7.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp8.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp10.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
MathUtils::BarycentricToXY( barycoord.getX() + 1, barycoord.getY() - 1, barycoord.getZ(), &point1, &point2, &cp10, &endpointUXY );
double zu = point1.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp1.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp2.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp4.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
double zv = cp1.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp2.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point2.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp4.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
double zw = cp3.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp4.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point1, &point2, &cp10, &endpointUXY );
endpointUXY.setZ( 3 * ( zu - zv ) );
MathUtils::BarycentricToXY( barycoord.getX(), barycoord.getY() + 1, barycoord.getZ() - 1, &point1, &point2, &cp10, &endpointVXY );
MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point1, &point2, &cp10, &endpointVXY );
endpointVXY.setZ( 3 * ( zv - zw ) );
Vector3D v1( endpointUXY.getX() - x, endpointUXY.getY() - y, endpointUXY.getZ() );
Vector3D v2( endpointVXY.getX() - x, endpointVXY.getY() - y, endpointVXY.getZ() );
Vector3D v1( endpointUXY.x() - x, endpointUXY.y() - y, endpointUXY.z() );
Vector3D v2( endpointVXY.x() - x, endpointVXY.y() - y, endpointVXY.z() );
result->setX( v1.getY()*v2.getZ() - v1.getZ()*v2.getY() );
result->setY( v1.getZ()*v2.getX() - v1.getX()*v2.getZ() );
result->setZ( v1.getX()*v2.getY() - v1.getY()*v2.getX() );
@ -100,17 +100,17 @@ bool CloughTocherInterpolator::calcNormVec( double x, double y, Vector3D *result
}
//is the point in the second subtriangle (point2,point3,cp10)?
MathUtils::calcBarycentricCoordinates( x, y, &point2, &point3, &cp10, &barycoord );
if ( barycoord.getX() >= -mEdgeTolerance && barycoord.getX() <= ( 1 + mEdgeTolerance ) && barycoord.getY() >= -mEdgeTolerance && barycoord.getY() <= ( 1 + mEdgeTolerance ) )
if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
{
double zu = point2.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp9.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp16.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp5.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp13.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp8.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double zv = cp9.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp16.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + point3.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp13.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp15.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp12.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double zw = cp5.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp13.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp15.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp8.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp12.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp10.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
MathUtils::BarycentricToXY( barycoord.getX() + 1, barycoord.getY() - 1, barycoord.getZ(), &point2, &point3, &cp10, &endpointUXY );
double zu = point2.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp9.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp16.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp13.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
double zv = cp9.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp16.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point3.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp13.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
double zw = cp5.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp13.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point2, &point3, &cp10, &endpointUXY );
endpointUXY.setZ( 3 * ( zu - zv ) );
MathUtils::BarycentricToXY( barycoord.getX(), barycoord.getY() + 1, barycoord.getZ() - 1, &point2, &point3, &cp10, &endpointVXY );
MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point2, &point3, &cp10, &endpointVXY );
endpointVXY.setZ( 3 * ( zv - zw ) );
Vector3D v1( endpointUXY.getX() - x, endpointUXY.getY() - y, endpointUXY.getZ() );
Vector3D v2( endpointVXY.getX() - x, endpointVXY.getY() - y, endpointVXY.getZ() );
Vector3D v1( endpointUXY.x() - x, endpointUXY.y() - y, endpointUXY.z() );
Vector3D v2( endpointVXY.x() - x, endpointVXY.y() - y, endpointVXY.z() );
result->setX( v1.getY()*v2.getZ() - v1.getZ()*v2.getY() );
result->setY( v1.getZ()*v2.getX() - v1.getX()*v2.getZ() );
result->setZ( v1.getX()*v2.getY() - v1.getY()*v2.getX() );
@ -120,17 +120,17 @@ bool CloughTocherInterpolator::calcNormVec( double x, double y, Vector3D *result
}
//is the point in the third subtriangle (point3,point1,cp10)?
MathUtils::calcBarycentricCoordinates( x, y, &point3, &point1, &cp10, &barycoord );
if ( barycoord.getX() >= -mEdgeTolerance && barycoord.getX() <= ( 1 + mEdgeTolerance ) && barycoord.getY() >= -mEdgeTolerance && barycoord.getY() <= ( 1 + mEdgeTolerance ) )
if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
{
double zu = point3.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp14.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp6.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp15.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp11.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp12.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double zv = cp14.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp6.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + point1.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp11.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp3.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp7.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double zw = cp15.getZ() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp11.getZ() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp3.getZ() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp12.getZ() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp7.getZ() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp10.getZ() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
MathUtils::BarycentricToXY( barycoord.getX() + 1, barycoord.getY() - 1, barycoord.getZ(), &point3, &point1, &cp10, &endpointUXY );
double zu = point3.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp14.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp6.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp11.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
double zv = cp14.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp6.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point1.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp11.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
double zw = cp15.z() * calcBernsteinPoly( 2, 2, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp11.z() * calcBernsteinPoly( 2, 1, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 2, 0, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 2, 1, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 2, 0, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 2, 0, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() );
MathUtils::BarycentricToXY( barycoord.x() + 1, barycoord.y() - 1, barycoord.z(), &point3, &point1, &cp10, &endpointUXY );
endpointUXY.setZ( 3 * ( zu - zv ) );
MathUtils::BarycentricToXY( barycoord.getX(), barycoord.getY() + 1, barycoord.getZ() - 1, &point3, &point1, &cp10, &endpointVXY );
MathUtils::BarycentricToXY( barycoord.x(), barycoord.y() + 1, barycoord.z() - 1, &point3, &point1, &cp10, &endpointVXY );
endpointVXY.setZ( 3 * ( zv - zw ) );
Vector3D v1( endpointUXY.getX() - x, endpointUXY.getY() - y, endpointUXY.getZ() );
Vector3D v2( endpointVXY.getX() - x, endpointVXY.getY() - y, endpointVXY.getZ() );
Vector3D v1( endpointUXY.x() - x, endpointUXY.y() - y, endpointUXY.z() );
Vector3D v2( endpointVXY.x() - x, endpointVXY.y() - y, endpointVXY.z() );
result->setX( v1.getY()*v2.getZ() - v1.getZ()*v2.getY() );
result->setY( v1.getZ()*v2.getX() - v1.getX()*v2.getZ() );
result->setZ( v1.getX()*v2.getY() - v1.getY()*v2.getX() );
@ -139,7 +139,7 @@ bool CloughTocherInterpolator::calcNormVec( double x, double y, Vector3D *result
}
//the point is in none of the subtriangles, test if point has the same position as one of the vertices
if ( x == point1.getX() && y == point1.getY() )
if ( x == point1.x() && y == point1.y() )
{
result->setX( -der1X );
result->setY( -der1Y );
@ -147,7 +147,7 @@ bool CloughTocherInterpolator::calcNormVec( double x, double y, Vector3D *result
result->standardise();
return true;
}
else if ( x == point2.getX() && y == point2.getY() )
else if ( x == point2.x() && y == point2.y() )
{
result->setX( -der2X );
result->setY( -der2Y );
@ -155,7 +155,7 @@ bool CloughTocherInterpolator::calcNormVec( double x, double y, Vector3D *result
result->standardise();
return true;
}
else if ( x == point3.getX() && y == point3.getY() )
else if ( x == point3.x() && y == point3.y() )
{
result->setX( -der3X );
result->setY( -der3Y );
@ -178,18 +178,18 @@ bool CloughTocherInterpolator::calcNormVec( double x, double y, Vector3D *result
}
}
bool CloughTocherInterpolator::calcPoint( double x, double y, Point3D *result )
bool CloughTocherInterpolator::calcPoint( double x, double y, QgsPoint *result )
{
if ( result )
{
init( x, y );
//find out, in which subtriangle the point (x,y) is
Point3D barycoord( 0, 0, 0 );
QgsPoint barycoord( 0, 0, 0 );
//is the point in the first subtriangle (point1,point2,cp10)?
MathUtils::calcBarycentricCoordinates( x, y, &point1, &point2, &cp10, &barycoord );
if ( barycoord.getX() >= -mEdgeTolerance && barycoord.getX() <= ( 1 + mEdgeTolerance ) && barycoord.getY() >= -mEdgeTolerance && barycoord.getY() <= ( 1 + mEdgeTolerance ) )
if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
{
double z = point1.getZ() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp1.getZ() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp2.getZ() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + point2.getZ() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp3.getZ() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp4.getZ() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp5.getZ() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp7.getZ() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp8.getZ() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp10.getZ() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double z = point1.z() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp1.z() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp2.z() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point2.z() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp4.z() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.x(), barycoord.y(), barycoord.z() );
result->setX( x );
result->setY( y );
result->setZ( z );
@ -197,9 +197,9 @@ bool CloughTocherInterpolator::calcPoint( double x, double y, Point3D *result )
}
//is the point in the second subtriangle (point2,point3,cp10)?
MathUtils::calcBarycentricCoordinates( x, y, &point2, &point3, &cp10, &barycoord );
if ( barycoord.getX() >= -mEdgeTolerance && barycoord.getX() <= ( 1 + mEdgeTolerance ) && barycoord.getY() >= -mEdgeTolerance && barycoord.getY() <= ( 1 + mEdgeTolerance ) )
if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
{
double z = cp10.getZ() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp8.getZ() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp5.getZ() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + point2.getZ() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp12.getZ() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp13.getZ() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp9.getZ() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp15.getZ() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp16.getZ() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + point3.getZ() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double z = cp10.z() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.x(), barycoord.y(), barycoord.z() ) + cp8.z() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp5.z() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + point2.z() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp13.z() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp9.z() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp16.z() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + point3.z() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.x(), barycoord.y(), barycoord.z() );
result->setX( x );
result->setY( y );
result->setZ( z );
@ -207,9 +207,9 @@ bool CloughTocherInterpolator::calcPoint( double x, double y, Point3D *result )
}
//is the point in the third subtriangle (point3,point1,cp10)?
MathUtils::calcBarycentricCoordinates( x, y, &point3, &point1, &cp10, &barycoord );
if ( barycoord.getX() >= -mEdgeTolerance && barycoord.getX() <= ( 1 + mEdgeTolerance ) && barycoord.getY() >= -mEdgeTolerance && barycoord.getY() <= ( 1 + mEdgeTolerance ) )
if ( barycoord.x() >= -mEdgeTolerance && barycoord.x() <= ( 1 + mEdgeTolerance ) && barycoord.y() >= -mEdgeTolerance && barycoord.y() <= ( 1 + mEdgeTolerance ) )
{
double z = point1.getZ() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp3.getZ() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp7.getZ() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp10.getZ() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp6.getZ() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp11.getZ() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp12.getZ() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp14.getZ() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + cp15.getZ() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.getX(), barycoord.getY(), barycoord.getZ() ) + point3.getZ() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.getX(), barycoord.getY(), barycoord.getZ() );
double z = point1.z() * calcBernsteinPoly( 3, 0, 3, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp3.z() * calcBernsteinPoly( 3, 0, 2, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp7.z() * calcBernsteinPoly( 3, 0, 1, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp10.z() * calcBernsteinPoly( 3, 0, 0, 3, barycoord.x(), barycoord.y(), barycoord.z() ) + cp6.z() * calcBernsteinPoly( 3, 1, 2, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp11.z() * calcBernsteinPoly( 3, 1, 1, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + cp12.z() * calcBernsteinPoly( 3, 1, 0, 2, barycoord.x(), barycoord.y(), barycoord.z() ) + cp14.z() * calcBernsteinPoly( 3, 2, 1, 0, barycoord.x(), barycoord.y(), barycoord.z() ) + cp15.z() * calcBernsteinPoly( 3, 2, 0, 1, barycoord.x(), barycoord.y(), barycoord.z() ) + point3.z() * calcBernsteinPoly( 3, 3, 0, 0, barycoord.x(), barycoord.y(), barycoord.z() );
result->setX( x );
result->setY( y );
result->setZ( z );
@ -217,25 +217,25 @@ bool CloughTocherInterpolator::calcPoint( double x, double y, Point3D *result )
}
//the point is in none of the subtriangles, test if point has the same position as one of the vertices
if ( x == point1.getX() && y == point1.getY() )
if ( x == point1.x() && y == point1.y() )
{
result->setX( x );
result->setY( y );
result->setZ( point1.getZ() );
result->setZ( point1.z() );
return true;
}
else if ( x == point2.getX() && y == point2.getY() )
else if ( x == point2.x() && y == point2.y() )
{
result->setX( x );
result->setY( y );
result->setZ( point2.getZ() );
result->setZ( point2.z() );
return true;
}
else if ( x == point3.getX() && y == point3.getY() )
else if ( x == point3.x() && y == point3.y() )
{
result->setX( x );
result->setY( y );
result->setZ( point3.getZ() );
result->setZ( point3.z() );
return true;
}
else
@ -280,33 +280,33 @@ void CloughTocherInterpolator::init( double x, double y )//version, which has th
der3Y = -v3.getY() / v3.getZ();
//calculate the control points
cp1.setX( point1.getX() + ( point2.getX() - point1.getX() ) / 3 );
cp1.setY( point1.getY() + ( point2.getY() - point1.getY() ) / 3 );
cp1.setZ( point1.getZ() + ( cp1.getX() - point1.getX() )*der1X + ( cp1.getY() - point1.getY() )*der1Y );
cp1.setX( point1.x() + ( point2.x() - point1.x() ) / 3 );
cp1.setY( point1.y() + ( point2.y() - point1.y() ) / 3 );
cp1.setZ( point1.z() + ( cp1.x() - point1.x() )*der1X + ( cp1.y() - point1.y() )*der1Y );
cp2.setX( point2.getX() + ( point1.getX() - point2.getX() ) / 3 );
cp2.setY( point2.getY() + ( point1.getY() - point2.getY() ) / 3 );
cp2.setZ( point2.getZ() + ( cp2.getX() - point2.getX() )*der2X + ( cp2.getY() - point2.getY() )*der2Y );
cp2.setX( point2.x() + ( point1.x() - point2.x() ) / 3 );
cp2.setY( point2.y() + ( point1.y() - point2.y() ) / 3 );
cp2.setZ( point2.z() + ( cp2.x() - point2.x() )*der2X + ( cp2.y() - point2.y() )*der2Y );
cp9.setX( point2.getX() + ( point3.getX() - point2.getX() ) / 3 );
cp9.setY( point2.getY() + ( point3.getY() - point2.getY() ) / 3 );
cp9.setZ( point2.getZ() + ( cp9.getX() - point2.getX() )*der2X + ( cp9.getY() - point2.getY() )*der2Y );
cp9.setX( point2.x() + ( point3.x() - point2.x() ) / 3 );
cp9.setY( point2.y() + ( point3.y() - point2.y() ) / 3 );
cp9.setZ( point2.z() + ( cp9.x() - point2.x() )*der2X + ( cp9.y() - point2.y() )*der2Y );
cp16.setX( point3.getX() + ( point2.getX() - point3.getX() ) / 3 );
cp16.setY( point3.getY() + ( point2.getY() - point3.getY() ) / 3 );
cp16.setZ( point3.getZ() + ( cp16.getX() - point3.getX() )*der3X + ( cp16.getY() - point3.getY() )*der3Y );
cp16.setX( point3.x() + ( point2.x() - point3.x() ) / 3 );
cp16.setY( point3.y() + ( point2.y() - point3.y() ) / 3 );
cp16.setZ( point3.z() + ( cp16.x() - point3.x() )*der3X + ( cp16.y() - point3.y() )*der3Y );
cp14.setX( point3.getX() + ( point1.getX() - point3.getX() ) / 3 );
cp14.setY( point3.getY() + ( point1.getY() - point3.getY() ) / 3 );
cp14.setZ( point3.getZ() + ( cp14.getX() - point3.getX() )*der3X + ( cp14.getY() - point3.getY() )*der3Y );
cp14.setX( point3.x() + ( point1.x() - point3.x() ) / 3 );
cp14.setY( point3.y() + ( point1.y() - point3.y() ) / 3 );
cp14.setZ( point3.z() + ( cp14.x() - point3.x() )*der3X + ( cp14.y() - point3.y() )*der3Y );
cp6.setX( point1.getX() + ( point3.getX() - point1.getX() ) / 3 );
cp6.setY( point1.getY() + ( point3.getY() - point1.getY() ) / 3 );
cp6.setZ( point1.getZ() + ( cp6.getX() - point1.getX() )*der1X + ( cp6.getY() - point1.getY() )*der1Y );
cp6.setX( point1.x() + ( point3.x() - point1.x() ) / 3 );
cp6.setY( point1.y() + ( point3.y() - point1.y() ) / 3 );
cp6.setZ( point1.z() + ( cp6.x() - point1.x() )*der1X + ( cp6.y() - point1.y() )*der1Y );
//set x- and y-coordinates of the central point
cp10.setX( ( point1.getX() + point2.getX() + point3.getX() ) / 3 );
cp10.setY( ( point1.getY() + point2.getY() + point3.getY() ) / 3 );
cp10.setX( ( point1.x() + point2.x() + point3.x() ) / 3 );
cp10.setY( ( point1.y() + point2.y() + point3.y() ) / 3 );
//set the derivatives of the points to new values if they are on a breakline
if ( state1 == NormVecDecorator::BreakLine )
@ -320,12 +320,12 @@ void CloughTocherInterpolator::init( double x, double y )//version, which has th
if ( state2 == NormVecDecorator::Normal )
{
//recalculate cp1
cp1.setZ( point1.getZ() + ( cp1.getX() - point1.getX() )*der1X + ( cp1.getY() - point1.getY() )*der1Y );
cp1.setZ( point1.z() + ( cp1.x() - point1.x() )*der1X + ( cp1.y() - point1.y() )*der1Y );
}
if ( state3 == NormVecDecorator::Normal )
{
//recalculate cp6
cp6.setZ( point1.getZ() + ( cp6.getX() - point1.getX() )*der1X + ( cp6.getY() - point1.getY() )*der1Y );
cp6.setZ( point1.z() + ( cp6.x() - point1.x() )*der1X + ( cp6.y() - point1.y() )*der1Y );
}
}
}
@ -341,12 +341,12 @@ void CloughTocherInterpolator::init( double x, double y )//version, which has th
if ( state1 == NormVecDecorator::Normal )
{
//recalculate cp2
cp2.setZ( point2.getZ() + ( cp2.getX() - point2.getX() )*der2X + ( cp2.getY() - point2.getY() )*der2Y );
cp2.setZ( point2.z() + ( cp2.x() - point2.x() )*der2X + ( cp2.y() - point2.y() )*der2Y );
}
if ( state3 == NormVecDecorator::Normal )
{
//recalculate cp9
cp9.setZ( point2.getZ() + ( cp9.getX() - point2.getX() )*der2X + ( cp9.getY() - point2.getY() )*der2Y );
cp9.setZ( point2.z() + ( cp9.x() - point2.x() )*der2X + ( cp9.y() - point2.y() )*der2Y );
}
}
}
@ -362,131 +362,131 @@ void CloughTocherInterpolator::init( double x, double y )//version, which has th
if ( state1 == NormVecDecorator::Normal )
{
//recalculate cp14
cp14.setZ( point3.getZ() + ( cp14.getX() - point3.getX() )*der3X + ( cp14.getY() - point3.getY() )*der3Y );
cp14.setZ( point3.z() + ( cp14.x() - point3.x() )*der3X + ( cp14.y() - point3.y() )*der3Y );
}
if ( state2 == NormVecDecorator::Normal )
{
//recalculate cp16
cp16.setZ( point3.getZ() + ( cp16.getX() - point3.getX() )*der3X + ( cp16.getY() - point3.getY() )*der3Y );
cp16.setZ( point3.z() + ( cp16.x() - point3.x() )*der3X + ( cp16.y() - point3.y() )*der3Y );
}
}
}
cp3.setX( point1.getX() + ( cp10.getX() - point1.getX() ) / 3 );
cp3.setY( point1.getY() + ( cp10.getY() - point1.getY() ) / 3 );
cp3.setZ( point1.getZ() + ( cp3.getX() - point1.getX() )*der1X + ( cp3.getY() - point1.getY() )*der1Y );
cp3.setX( point1.x() + ( cp10.x() - point1.x() ) / 3 );
cp3.setY( point1.y() + ( cp10.y() - point1.y() ) / 3 );
cp3.setZ( point1.z() + ( cp3.x() - point1.x() )*der1X + ( cp3.y() - point1.y() )*der1Y );
cp5.setX( point2.getX() + ( cp10.getX() - point2.getX() ) / 3 );
cp5.setY( point2.getY() + ( cp10.getY() - point2.getY() ) / 3 );
cp5.setZ( point2.getZ() + ( cp5.getX() - point2.getX() )*der2X + ( cp5.getY() - point2.getY() )*der2Y );
cp5.setX( point2.x() + ( cp10.x() - point2.x() ) / 3 );
cp5.setY( point2.y() + ( cp10.y() - point2.y() ) / 3 );
cp5.setZ( point2.z() + ( cp5.x() - point2.x() )*der2X + ( cp5.y() - point2.y() )*der2Y );
cp15.setX( point3.getX() + ( cp10.getX() - point3.getX() ) / 3 );
cp15.setY( point3.getY() + ( cp10.getY() - point3.getY() ) / 3 );
cp15.setZ( point3.getZ() + ( cp15.getX() - point3.getX() )*der3X + ( cp15.getY() - point3.getY() )*der3Y );
cp15.setX( point3.x() + ( cp10.x() - point3.x() ) / 3 );
cp15.setY( point3.y() + ( cp10.y() - point3.y() ) / 3 );
cp15.setZ( point3.z() + ( cp15.x() - point3.x() )*der3X + ( cp15.y() - point3.y() )*der3Y );
cp4.setX( ( point1.getX() + cp10.getX() + point2.getX() ) / 3 );
cp4.setY( ( point1.getY() + cp10.getY() + point2.getY() ) / 3 );
cp4.setX( ( point1.x() + cp10.x() + point2.x() ) / 3 );
cp4.setY( ( point1.y() + cp10.y() + point2.y() ) / 3 );
Point3D midpoint3( ( cp1.getX() + cp2.getX() ) / 2, ( cp1.getY() + cp2.getY() ) / 2, ( cp1.getZ() + cp2.getZ() ) / 2 );
Vector3D cp1cp2( cp2.getX() - cp1.getX(), cp2.getY() - cp1.getY(), cp2.getZ() - cp1.getZ() );
QgsPoint midpoint3( ( cp1.x() + cp2.x() ) / 2, ( cp1.y() + cp2.y() ) / 2, ( cp1.z() + cp2.z() ) / 2 );
Vector3D cp1cp2( cp2.x() - cp1.x(), cp2.y() - cp1.y(), cp2.z() - cp1.z() );
Vector3D odir3( 0, 0, 0 );//direction orthogonal to the line from point1 to point2
if ( ( point2.getY() - point1.getY() ) != 0 ) //avoid division through zero
if ( ( point2.y() - point1.y() ) != 0 ) //avoid division through zero
{
odir3.setX( 1 );
odir3.setY( -( point2.getX() - point1.getX() ) / ( point2.getY() - point1.getY() ) );
odir3.setY( -( point2.x() - point1.x() ) / ( point2.y() - point1.y() ) );
odir3.setZ( ( der1X + der2X ) / 2 + ( der1Y + der2Y ) / 2 * odir3.getY() ); //take the linear interpolated cross-derivative
}
else
{
odir3.setY( 1 );
odir3.setX( -( point2.getY() - point1.getY() ) / ( point2.getX() - point1.getX() ) );
odir3.setX( -( point2.y() - point1.y() ) / ( point2.x() - point1.x() ) );
odir3.setZ( ( der1X + der2X ) / 2 * odir3.getX() + ( der1Y + der2Y ) / 2 );
}
Vector3D midpoint3cp4( 0, 0, 0 );
MathUtils::derVec( &cp1cp2, &odir3, &midpoint3cp4, cp4.getX() - midpoint3.getX(), cp4.getY() - midpoint3.getY() );
cp4.setZ( midpoint3.getZ() + midpoint3cp4.getZ() );
MathUtils::derVec( &cp1cp2, &odir3, &midpoint3cp4, cp4.x() - midpoint3.x(), cp4.y() - midpoint3.y() );
cp4.setZ( midpoint3.z() + midpoint3cp4.getZ() );
cp13.setX( ( point2.getX() + cp10.getX() + point3.getX() ) / 3 );
cp13.setY( ( point2.getY() + cp10.getY() + point3.getY() ) / 3 );
cp13.setX( ( point2.x() + cp10.x() + point3.x() ) / 3 );
cp13.setY( ( point2.y() + cp10.y() + point3.y() ) / 3 );
//find the point in the middle of the bezier curve between point2 and point3
Point3D midpoint1( ( cp9.getX() + cp16.getX() ) / 2, ( cp9.getY() + cp16.getY() ) / 2, ( cp9.getZ() + cp16.getZ() ) / 2 );
Vector3D cp9cp16( cp16.getX() - cp9.getX(), cp16.getY() - cp9.getY(), cp16.getZ() - cp9.getZ() );
QgsPoint midpoint1( ( cp9.x() + cp16.x() ) / 2, ( cp9.y() + cp16.y() ) / 2, ( cp9.z() + cp16.z() ) / 2 );
Vector3D cp9cp16( cp16.x() - cp9.x(), cp16.y() - cp9.y(), cp16.z() - cp9.z() );
Vector3D odir1( 0, 0, 0 );//direction orthogonal to the line from point2 to point3
if ( ( point3.getY() - point2.getY() ) != 0 ) //avoid division through zero
if ( ( point3.y() - point2.y() ) != 0 ) //avoid division through zero
{
odir1.setX( 1 );
odir1.setY( -( point3.getX() - point2.getX() ) / ( point3.getY() - point2.getY() ) );
odir1.setY( -( point3.x() - point2.x() ) / ( point3.y() - point2.y() ) );
odir1.setZ( ( der2X + der3X ) / 2 + ( der2Y + der3Y ) / 2 * odir1.getY() ); //take the linear interpolated cross-derivative
}
else
{
odir1.setY( 1 );
odir1.setX( -( point3.getY() - point2.getY() ) / ( point3.getX() - point2.getX() ) );
odir1.setX( -( point3.y() - point2.y() ) / ( point3.x() - point2.x() ) );
odir1.setZ( ( der2X + der3X ) / 2 * odir1.getX() + ( der2Y + der3Y ) / 2 );
}
Vector3D midpoint1cp13( 0, 0, 0 );
MathUtils::derVec( &cp9cp16, &odir1, &midpoint1cp13, cp13.getX() - midpoint1.getX(), cp13.getY() - midpoint1.getY() );
cp13.setZ( midpoint1.getZ() + midpoint1cp13.getZ() );
MathUtils::derVec( &cp9cp16, &odir1, &midpoint1cp13, cp13.x() - midpoint1.x(), cp13.y() - midpoint1.y() );
cp13.setZ( midpoint1.z() + midpoint1cp13.getZ() );
cp11.setX( ( point3.getX() + cp10.getX() + point1.getX() ) / 3 );
cp11.setY( ( point3.getY() + cp10.getY() + point1.getY() ) / 3 );
cp11.setX( ( point3.x() + cp10.x() + point1.x() ) / 3 );
cp11.setY( ( point3.y() + cp10.y() + point1.y() ) / 3 );
//find the point in the middle of the bezier curve between point3 and point2
Point3D midpoint2( ( cp14.getX() + cp6.getX() ) / 2, ( cp14.getY() + cp6.getY() ) / 2, ( cp14.getZ() + cp6.getZ() ) / 2 );
Vector3D cp14cp6( cp6.getX() - cp14.getX(), cp6.getY() - cp14.getY(), cp6.getZ() - cp14.getZ() );
QgsPoint midpoint2( ( cp14.x() + cp6.x() ) / 2, ( cp14.y() + cp6.y() ) / 2, ( cp14.z() + cp6.z() ) / 2 );
Vector3D cp14cp6( cp6.x() - cp14.x(), cp6.y() - cp14.y(), cp6.z() - cp14.z() );
Vector3D odir2( 0, 0, 0 );//direction orthogonal to the line from point 1 to point3
if ( ( point3.getY() - point1.getY() ) != 0 ) //avoid division through zero
if ( ( point3.y() - point1.y() ) != 0 ) //avoid division through zero
{
odir2.setX( 1 );
odir2.setY( -( point3.getX() - point1.getX() ) / ( point3.getY() - point1.getY() ) );
odir2.setY( -( point3.x() - point1.x() ) / ( point3.y() - point1.y() ) );
odir2.setZ( ( der3X + der1X ) / 2 + ( der3Y + der1Y ) / 2 * odir2.getY() ); //take the linear interpolated cross-derivative
}
else
{
odir2.setY( 1 );
odir2.setX( -( point3.getY() - point1.getY() ) / ( point3.getX() - point1.getX() ) );
odir2.setX( -( point3.y() - point1.y() ) / ( point3.x() - point1.x() ) );
odir2.setZ( ( der3X + der1X ) / 2 * odir2.getX() + ( der3Y + der1Y ) / 2 );
}
Vector3D midpoint2cp11( 0, 0, 0 );
MathUtils::derVec( &cp14cp6, &odir2, &midpoint2cp11, cp11.getX() - midpoint2.getX(), cp11.getY() - midpoint2.getY() );
cp11.setZ( midpoint2.getZ() + midpoint2cp11.getZ() );
MathUtils::derVec( &cp14cp6, &odir2, &midpoint2cp11, cp11.x() - midpoint2.x(), cp11.y() - midpoint2.y() );
cp11.setZ( midpoint2.z() + midpoint2cp11.getZ() );
cp7.setX( cp10.getX() + ( point1.getX() - cp10.getX() ) / 3 );
cp7.setY( cp10.getY() + ( point1.getY() - cp10.getY() ) / 3 );
cp7.setX( cp10.x() + ( point1.x() - cp10.x() ) / 3 );
cp7.setY( cp10.y() + ( point1.y() - cp10.y() ) / 3 );
//cp7 has to be in the same plane as cp4, cp3 and cp11
Vector3D cp4cp3( cp3.getX() - cp4.getX(), cp3.getY() - cp4.getY(), cp3.getZ() - cp4.getZ() );
Vector3D cp4cp11( cp11.getX() - cp4.getX(), cp11.getY() - cp4.getY(), cp11.getZ() - cp4.getZ() );
Vector3D cp4cp3( cp3.x() - cp4.x(), cp3.y() - cp4.y(), cp3.z() - cp4.z() );
Vector3D cp4cp11( cp11.x() - cp4.x(), cp11.y() - cp4.y(), cp11.z() - cp4.z() );
Vector3D cp4cp7( 0, 0, 0 );
MathUtils::derVec( &cp4cp3, &cp4cp11, &cp4cp7, cp7.getX() - cp4.getX(), cp7.getY() - cp4.getY() );
cp7.setZ( cp4.getZ() + cp4cp7.getZ() );
MathUtils::derVec( &cp4cp3, &cp4cp11, &cp4cp7, cp7.x() - cp4.x(), cp7.y() - cp4.y() );
cp7.setZ( cp4.z() + cp4cp7.getZ() );
cp8.setX( cp10.getX() + ( point2.getX() - cp10.getX() ) / 3 );
cp8.setY( cp10.getY() + ( point2.getY() - cp10.getY() ) / 3 );
cp8.setX( cp10.x() + ( point2.x() - cp10.x() ) / 3 );
cp8.setY( cp10.y() + ( point2.y() - cp10.y() ) / 3 );
//cp8 has to be in the same plane as cp4, cp5 and cp13
Vector3D cp4cp5( cp5.getX() - cp4.getX(), cp5.getY() - cp4.getY(), cp5.getZ() - cp4.getZ() );
Vector3D cp4cp13( cp13.getX() - cp4.getX(), cp13.getY() - cp4.getY(), cp13.getZ() - cp4.getZ() );
Vector3D cp4cp5( cp5.x() - cp4.x(), cp5.y() - cp4.y(), cp5.z() - cp4.z() );
Vector3D cp4cp13( cp13.x() - cp4.x(), cp13.y() - cp4.y(), cp13.z() - cp4.z() );
Vector3D cp4cp8( 0, 0, 0 );
MathUtils::derVec( &cp4cp5, &cp4cp13, &cp4cp8, cp8.getX() - cp4.getX(), cp8.getY() - cp4.getY() );
cp8.setZ( cp4.getZ() + cp4cp8.getZ() );
MathUtils::derVec( &cp4cp5, &cp4cp13, &cp4cp8, cp8.x() - cp4.x(), cp8.y() - cp4.y() );
cp8.setZ( cp4.z() + cp4cp8.getZ() );
cp12.setX( cp10.getX() + ( point3.getX() - cp10.getX() ) / 3 );
cp12.setY( cp10.getY() + ( point3.getY() - cp10.getY() ) / 3 );
cp12.setX( cp10.x() + ( point3.x() - cp10.x() ) / 3 );
cp12.setY( cp10.y() + ( point3.y() - cp10.y() ) / 3 );
//cp12 has to be in the same plane as cp13, cp15 and cp11
Vector3D cp13cp11( cp11.getX() - cp13.getX(), cp11.getY() - cp13.getY(), cp11.getZ() - cp13.getZ() );
Vector3D cp13cp15( cp15.getX() - cp13.getX(), cp15.getY() - cp13.getY(), cp15.getZ() - cp13.getZ() );
Vector3D cp13cp11( cp11.x() - cp13.x(), cp11.y() - cp13.y(), cp11.z() - cp13.z() );
Vector3D cp13cp15( cp15.x() - cp13.x(), cp15.y() - cp13.y(), cp15.z() - cp13.z() );
Vector3D cp13cp12( 0, 0, 0 );
MathUtils::derVec( &cp13cp11, &cp13cp15, &cp13cp12, cp12.getX() - cp13.getX(), cp12.getY() - cp13.getY() );
cp12.setZ( cp13.getZ() + cp13cp12.getZ() );
MathUtils::derVec( &cp13cp11, &cp13cp15, &cp13cp12, cp12.x() - cp13.x(), cp12.y() - cp13.y() );
cp12.setZ( cp13.z() + cp13cp12.getZ() );
//cp10 has to be in the same plane as cp7, cp8 and cp12
Vector3D cp7cp8( cp8.getX() - cp7.getX(), cp8.getY() - cp7.getY(), cp8.getZ() - cp7.getZ() );
Vector3D cp7cp12( cp12.getX() - cp7.getX(), cp12.getY() - cp7.getY(), cp12.getZ() - cp7.getZ() );
Vector3D cp7cp8( cp8.x() - cp7.x(), cp8.y() - cp7.y(), cp8.z() - cp7.z() );
Vector3D cp7cp12( cp12.x() - cp7.x(), cp12.y() - cp7.y(), cp12.z() - cp7.z() );
Vector3D cp7cp10( 0, 0, 0 );
MathUtils::derVec( &cp7cp8, &cp7cp12, &cp7cp10, cp10.getX() - cp7.getX(), cp10.getY() - cp7.getY() );
cp10.setZ( cp7.getZ() + cp7cp10.getZ() );
MathUtils::derVec( &cp7cp8, &cp7cp12, &cp7cp10, cp10.x() - cp7.x(), cp10.y() - cp7.y() );
cp10.setZ( cp7.z() + cp7cp10.getZ() );
lpoint1 = point1;
lpoint2 = point2;
@ -659,7 +659,7 @@ void CloughTocherInterpolator::init( double x, double y )//version which has uni
cp4.setX( ( point1.getX() + cp10.getX() + point2.getX() ) / 3 );
cp4.setY( ( point1.getY() + cp10.getY() + point2.getY() ) / 3 );
Point3D midpoint3( ( cp1.getX() + cp2.getX() ) / 2, ( cp1.getY() + cp2.getY() ) / 2, ( cp1.getZ() + cp2.getZ() ) / 2 );
QgsPoint midpoint3( ( cp1.getX() + cp2.getX() ) / 2, ( cp1.getY() + cp2.getY() ) / 2, ( cp1.getZ() + cp2.getZ() ) / 2 );
Vector3D cp1cp2( cp2.getX() - cp1.getX(), cp2.getY() - cp1.getY(), cp2.getZ() - cp1.getZ() );
Vector3D odir3( 0, 0, 0 );//direction orthogonal to the line from point1 to point2
if ( ( point2.getY() - point1.getY() ) != 0 ) //avoid division through zero
@ -681,7 +681,7 @@ void CloughTocherInterpolator::init( double x, double y )//version which has uni
cp13.setX( ( point2.getX() + cp10.getX() + point3.getX() ) / 3 );
cp13.setY( ( point2.getY() + cp10.getY() + point3.getY() ) / 3 );
//find the point in the middle of the bezier curve between point2 and point3
Point3D midpoint1( ( cp9.getX() + cp16.getX() ) / 2, ( cp9.getY() + cp16.getY() ) / 2, ( cp9.getZ() + cp16.getZ() ) / 2 );
QgsPoint midpoint1( ( cp9.getX() + cp16.getX() ) / 2, ( cp9.getY() + cp16.getY() ) / 2, ( cp9.getZ() + cp16.getZ() ) / 2 );
Vector3D cp9cp16( cp16.getX() - cp9.getX(), cp16.getY() - cp9.getY(), cp16.getZ() - cp9.getZ() );
Vector3D odir1( 0, 0, 0 );//direction orthogonal to the line from point2 to point3
if ( ( point3.getY() - point2.getY() ) != 0 ) //avoid division through zero
@ -704,7 +704,7 @@ void CloughTocherInterpolator::init( double x, double y )//version which has uni
cp11.setX( ( point3.getX() + cp10.getX() + point1.getX() ) / 3 );
cp11.setY( ( point3.getY() + cp10.getY() + point1.getY() ) / 3 );
//find the point in the middle of the bezier curve between point3 and point2
Point3D midpoint2( ( cp14.getX() + cp6.getX() ) / 2, ( cp14.getY() + cp6.getY() ) / 2, ( cp14.getZ() + cp6.getZ() ) / 2 );
QgsPoint midpoint2( ( cp14.getX() + cp6.getX() ) / 2, ( cp14.getY() + cp6.getY() ) / 2, ( cp14.getZ() + cp6.getZ() ) / 2 );
Vector3D cp14cp6( cp6.getX() - cp14.getX(), cp6.getY() - cp14.getY(), cp6.getZ() - cp14.getZ() );
Vector3D odir2( 0, 0, 0 );//direction orthogonal to the line from point 1 to point3
if ( ( point3.getY() - point1.getY() ) != 0 ) //avoid division through zero

48
src/analysis/interpolation/CloughTocherInterpolator.h
View File

@ -18,7 +18,7 @@
#define CLOUGHTOCHERINTERPOLATOR_H
#include "TriangleInterpolator.h"
#include "Point3D.h"
#include "qgspoint.h"
#include "qgis_analysis.h"
class NormVecDecorator;
@ -33,27 +33,27 @@ class ANALYSIS_EXPORT CloughTocherInterpolator : public TriangleInterpolator
//! Tolerance of the barycentric coordinates at the borders of the triangles (to prevent errors because of very small negativ baricentric coordinates)
double mEdgeTolerance;
//! First point of the triangle in x-,y-,z-coordinates
Point3D point1;
QgsPoint point1;
//! Second point of the triangle in x-,y-,z-coordinates
Point3D point2;
QgsPoint point2;
//! Third point of the triangle in x-,y-,z-coordinates
Point3D point3;
Point3D cp1;
Point3D cp2;
Point3D cp3;
Point3D cp4;
Point3D cp5;
Point3D cp6;
Point3D cp7;
Point3D cp8;
Point3D cp9;
Point3D cp10;
Point3D cp11;
Point3D cp12;
Point3D cp13;
Point3D cp14;
Point3D cp15;
Point3D cp16;
QgsPoint point3;
QgsPoint cp1;
QgsPoint cp2;
QgsPoint cp3;
QgsPoint cp4;
QgsPoint cp5;
QgsPoint cp6;
QgsPoint cp7;
QgsPoint cp8;
QgsPoint cp9;
QgsPoint cp10;
QgsPoint cp11;
QgsPoint cp12;
QgsPoint cp13;
QgsPoint cp14;
QgsPoint cp15;
QgsPoint cp16;
//! Derivative in x-direction at point1
double der1X;
//! Derivative in y-direction at point1
@ -67,11 +67,11 @@ class ANALYSIS_EXPORT CloughTocherInterpolator : public TriangleInterpolator
//! Derivative in y-direction at point3
double der3Y;
//! Stores point1 of the last run
Point3D lpoint1;
QgsPoint lpoint1;
//! Stores point2 of the last run
Point3D lpoint2;
QgsPoint lpoint2;
//! Stores point3 of the last run
Point3D lpoint3;
QgsPoint lpoint3;
//! Finds out, in which triangle the point with the coordinates x and y is
void init( double x, double y );
//! Calculates the Bernsteinpolynomials to calculate the Beziertriangle. 'n' is three in the cubical case, 'i', 'j', 'k' are the indices of the controllpoint and 'u', 'v', 'w' are the barycentric coordinates of the point
@ -86,7 +86,7 @@ class ANALYSIS_EXPORT CloughTocherInterpolator : public TriangleInterpolator
//! Calculates the normal vector and assigns it to vec (not implemented at the moment)
virtual bool calcNormVec( double x, double y, Vector3D *result ) override;
//! Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point
virtual bool calcPoint( double x, double y, Point3D *result ) override;
virtual bool calcPoint( double x, double y, QgsPoint *result ) override;
virtual void setTriangulation( NormVecDecorator *tin );
};

346
src/analysis/interpolation/DualEdgeTriangulation.cc
View File

@ -80,8 +80,8 @@ void DualEdgeTriangulation::addLine( Line3D *line, bool breakline )
{
line->goToNext();
// Use copy ctor since line can be deleted as well as its
// associated Node and Point3D
actpoint = mDecorator->addPoint( new Point3D( *line->getPoint() ) );
// associated Node and QgsPoint
actpoint = mDecorator->addPoint( new QgsPoint( *line->getPoint() ) );
if ( actpoint != -100 )
{
i++;
@ -98,7 +98,7 @@ void DualEdgeTriangulation::addLine( Line3D *line, bool breakline )
for ( ; i < line->getSize(); i++ )
{
line->goToNext();
currentpoint = mDecorator->addPoint( new Point3D( *line->getPoint() ) );
currentpoint = mDecorator->addPoint( new QgsPoint( *line->getPoint() ) );
if ( currentpoint != -100 && actpoint != -100 && currentpoint != actpoint )//-100 is the return value if the point could not be not inserted
{
insertForcedSegment( actpoint, currentpoint, breakline );
@ -109,7 +109,7 @@ void DualEdgeTriangulation::addLine( Line3D *line, bool breakline )
delete line;
}
int DualEdgeTriangulation::addPoint( Point3D *p )
int DualEdgeTriangulation::addPoint( QgsPoint *p )
{
if ( p )
{
@ -118,30 +118,30 @@ int DualEdgeTriangulation::addPoint( Point3D *p )
//first update the bounding box
if ( mPointVector.isEmpty() )//update bounding box when the first point is inserted
{
xMin = ( *p ).getX();
yMin = ( *p ).getY();
xMax = ( *p ).getX();
yMax = ( *p ).getY();
xMin = ( *p ).x();
yMin = ( *p ).y();
xMax = ( *p ).x();
yMax = ( *p ).y();
}
else//update bounding box else
{
if ( ( *p ).getX() < xMin )
if ( ( *p ).x() < xMin )
{
xMin = ( *p ).getX();
xMin = ( *p ).x();
}
else if ( ( *p ).getX() > xMax )
else if ( ( *p ).x() > xMax )
{
xMax = ( *p ).getX();
xMax = ( *p ).x();
}
if ( ( *p ).getY() < yMin )
if ( ( *p ).y() < yMin )
{
yMin = ( *p ).getY();
yMin = ( *p ).y();
}
else if ( ( *p ).getY() > yMax )
else if ( ( *p ).y() > yMax )
{
yMax = ( *p ).getY();
yMax = ( *p ).y();
}
}
@ -161,10 +161,10 @@ int DualEdgeTriangulation::addPoint( Point3D *p )
else if ( mPointVector.count() == 2 )//insert the second point into the triangulation
{
//test, if it is the same point as the first point
if ( p->getX() == mPointVector[0]->getX() && p->getY() == mPointVector[0]->getY() )
if ( p->x() == mPointVector[0]->x() && p->y() == mPointVector[0]->y() )
{
QgsDebugMsg( "second point is the same as the first point, it thus has not been inserted" );
Point3D *p = mPointVector[1];
QgsPoint *p = mPointVector[1];
mPointVector.remove( 1 );
delete p;
return -100;
@ -377,7 +377,7 @@ int DualEdgeTriangulation::addPoint( Point3D *p )
else if ( number == -100 || number == -5 )//this means unknown problems or a numerical error occurred in 'baseEdgeOfTriangle'
{
// QgsDebugMsg("point has not been inserted because of unknown problems");
Point3D *p = mPointVector[mPointVector.count() - 1];
QgsPoint *p = mPointVector[mPointVector.count() - 1];
mPointVector.remove( mPointVector.count() - 1 );
delete p;
return -100;
@ -385,9 +385,9 @@ int DualEdgeTriangulation::addPoint( Point3D *p )
else if ( number == -25 )//this means that the point has already been inserted in the triangulation
{
//Take the higher z-Value in case of two equal points
Point3D *newPoint = mPointVector[mPointVector.count() - 1];
Point3D *existingPoint = mPointVector[mTwiceInsPoint];
existingPoint->setZ( qMax( newPoint->getZ(), existingPoint->getZ() ) );
QgsPoint *newPoint = mPointVector[mPointVector.count() - 1];
QgsPoint *existingPoint = mPointVector[mTwiceInsPoint];
existingPoint->setZ( qMax( newPoint->z(), existingPoint->z() ) );
mPointVector.remove( mPointVector.count() - 1 );
delete newPoint;
@ -476,7 +476,7 @@ int DualEdgeTriangulation::baseEdgeOfPoint( int point )
}
}
int DualEdgeTriangulation::baseEdgeOfTriangle( Point3D *point )
int DualEdgeTriangulation::baseEdgeOfTriangle( QgsPoint *point )
{
unsigned int actedge = mEdgeInside;//start with an edge which does not point to the virtual point (usually number 3)
int counter = 0;//number of consecutive successful left-of-tests
@ -596,12 +596,12 @@ int DualEdgeTriangulation::baseEdgeOfTriangle( Point3D *point )
nr1 = mHalfEdge[actedge]->getPoint();
nr2 = mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint();
nr3 = mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext()]->getPoint();
double x1 = mPointVector[nr1]->getX();
double y1 = mPointVector[nr1]->getY();
double x2 = mPointVector[nr2]->getX();
double y2 = mPointVector[nr2]->getY();
double x3 = mPointVector[nr3]->getX();
double y3 = mPointVector[nr3]->getY();
double x1 = mPointVector[nr1]->x();
double y1 = mPointVector[nr1]->y();
double x2 = mPointVector[nr2]->x();
double y2 = mPointVector[nr2]->y();
double x3 = mPointVector[nr3]->x();
double y3 = mPointVector[nr3]->y();
//now make sure that always the same edge is returned
if ( x1 < x2 && x1 < x3 )//return the edge which points to the point with the lowest x-coordinate
@ -668,7 +668,7 @@ bool DualEdgeTriangulation::calcNormal( double x, double y, Vector3D *result )
}
}
bool DualEdgeTriangulation::calcPoint( double x, double y, Point3D *result )
bool DualEdgeTriangulation::calcPoint( double x, double y, QgsPoint *result )
{
if ( result && mTriangleInterpolator )
{
@ -687,10 +687,10 @@ bool DualEdgeTriangulation::checkSwap( unsigned int edge, unsigned int recursive
{
if ( swapPossible( edge ) )
{
Point3D *pta = mPointVector[mHalfEdge[edge]->getPoint()];
Point3D *ptb = mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()];
Point3D *ptc = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint()];
Point3D *ptd = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint()];
QgsPoint *pta = mPointVector[mHalfEdge[edge]->getPoint()];
QgsPoint *ptb = mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()];
QgsPoint *ptc = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint()];
QgsPoint *ptd = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint()];
if ( MathUtils::inCircle( ptd, pta, ptb, ptc ) && recursiveDeep < 100 )//empty circle criterion violated
{
doSwap( edge, recursiveDeep );//swap the edge (recursive)
@ -963,7 +963,7 @@ QList<int> *DualEdgeTriangulation::getSurroundingTriangles( int pointno )
}
bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, int *n1, Point3D *p2, int *n2, Point3D *p3, int *n3 )
bool DualEdgeTriangulation::getTriangle( double x, double y, QgsPoint *p1, int *n1, QgsPoint *p2, int *n2, QgsPoint *p3, int *n3 )
{
if ( mPointVector.size() < 3 )
{
@ -972,7 +972,7 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, int *n
if ( p1 && p2 && p3 )
{
Point3D point( x, y, 0 );
QgsPoint point( x, y, 0 );
int edge = baseEdgeOfTriangle( &point );
if ( edge == -10 )//the point is outside the convex hull
{
@ -985,15 +985,15 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, int *n
int ptnr1 = mHalfEdge[edge]->getPoint();
int ptnr2 = mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint();
int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint();
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
p2->setX( mPointVector[ptnr2]->getX() );
p2->setY( mPointVector[ptnr2]->getY() );
p2->setZ( mPointVector[ptnr2]->getZ() );
p3->setX( mPointVector[ptnr3]->getX() );
p3->setY( mPointVector[ptnr3]->getY() );
p3->setZ( mPointVector[ptnr3]->getZ() );
p1->setX( mPointVector[ptnr1]->x() );
p1->setY( mPointVector[ptnr1]->y() );
p1->setZ( mPointVector[ptnr1]->z() );
p2->setX( mPointVector[ptnr2]->x() );
p2->setY( mPointVector[ptnr2]->y() );
p2->setZ( mPointVector[ptnr2]->z() );
p3->setX( mPointVector[ptnr3]->x() );
p3->setY( mPointVector[ptnr3]->y() );
p3->setZ( mPointVector[ptnr3]->z() );
( *n1 ) = ptnr1;
( *n2 ) = ptnr2;
( *n3 ) = ptnr3;
@ -1008,15 +1008,15 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, int *n
{
return false;
}
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
p2->setX( mPointVector[ptnr2]->getX() );
p2->setY( mPointVector[ptnr2]->getY() );
p2->setZ( mPointVector[ptnr2]->getZ() );
p3->setX( mPointVector[ptnr3]->getX() );
p3->setY( mPointVector[ptnr3]->getY() );
p3->setZ( mPointVector[ptnr3]->getZ() );
p1->setX( mPointVector[ptnr1]->x() );
p1->setY( mPointVector[ptnr1]->y() );
p1->setZ( mPointVector[ptnr1]->z() );
p2->setX( mPointVector[ptnr2]->x() );
p2->setY( mPointVector[ptnr2]->y() );
p2->setZ( mPointVector[ptnr2]->z() );
p3->setX( mPointVector[ptnr3]->x() );
p3->setY( mPointVector[ptnr3]->y() );
p3->setZ( mPointVector[ptnr3]->z() );
( *n1 ) = ptnr1;
( *n2 ) = ptnr2;
( *n3 ) = ptnr3;
@ -1030,15 +1030,15 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, int *n
int ptnr1 = mHalfEdge[edge1]->getPoint();
int ptnr2 = mHalfEdge[edge2]->getPoint();
int ptnr3 = mHalfEdge[edge3]->getPoint();
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
p2->setX( mPointVector[ptnr2]->getX() );
p2->setY( mPointVector[ptnr2]->getY() );
p2->setZ( mPointVector[ptnr2]->getZ() );
p3->setX( mPointVector[ptnr3]->getX() );
p3->setY( mPointVector[ptnr3]->getY() );
p3->setZ( mPointVector[ptnr3]->getZ() );
p1->setX( mPointVector[ptnr1]->x() );
p1->setY( mPointVector[ptnr1]->y() );
p1->setZ( mPointVector[ptnr1]->z() );
p2->setX( mPointVector[ptnr2]->x() );
p2->setY( mPointVector[ptnr2]->y() );
p2->setZ( mPointVector[ptnr2]->z() );
p3->setX( mPointVector[ptnr3]->x() );
p3->setY( mPointVector[ptnr3]->y() );
p3->setZ( mPointVector[ptnr3]->z() );
( *n1 ) = ptnr1;
( *n2 ) = ptnr2;
( *n3 ) = ptnr3;
@ -1053,15 +1053,15 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, int *n
{
return false;
}
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
p2->setX( mPointVector[ptnr2]->getX() );
p2->setY( mPointVector[ptnr2]->getY() );
p2->setZ( mPointVector[ptnr2]->getZ() );
p3->setX( mPointVector[ptnr3]->getX() );
p3->setY( mPointVector[ptnr3]->getY() );
p3->setZ( mPointVector[ptnr3]->getZ() );
p1->setX( mPointVector[ptnr1]->x() );
p1->setY( mPointVector[ptnr1]->y() );
p1->setZ( mPointVector[ptnr1]->z() );
p2->setX( mPointVector[ptnr2]->x() );
p2->setY( mPointVector[ptnr2]->y() );
p2->setZ( mPointVector[ptnr2]->z() );
p3->setX( mPointVector[ptnr3]->x() );
p3->setY( mPointVector[ptnr3]->y() );
p3->setZ( mPointVector[ptnr3]->z() );
( *n1 ) = ptnr1;
( *n2 ) = ptnr2;
( *n3 ) = ptnr3;
@ -1081,7 +1081,7 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, int *n
}
}
bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, Point3D *p2, Point3D *p3 )
bool DualEdgeTriangulation::getTriangle( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 )
{
if ( mPointVector.size() < 3 )
{
@ -1090,7 +1090,7 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, Point3
if ( p1 && p2 && p3 )
{
Point3D point( x, y, 0 );
QgsPoint point( x, y, 0 );
int edge = baseEdgeOfTriangle( &point );
if ( edge == -10 )//the point is outside the convex hull
{
@ -1102,15 +1102,15 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, Point3
int ptnr1 = mHalfEdge[edge]->getPoint();
int ptnr2 = mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint();
int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint();
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
p2->setX( mPointVector[ptnr2]->getX() );
p2->setY( mPointVector[ptnr2]->getY() );
p2->setZ( mPointVector[ptnr2]->getZ() );
p3->setX( mPointVector[ptnr3]->getX() );
p3->setY( mPointVector[ptnr3]->getY() );
p3->setZ( mPointVector[ptnr3]->getZ() );
p1->setX( mPointVector[ptnr1]->x() );
p1->setY( mPointVector[ptnr1]->y() );
p1->setZ( mPointVector[ptnr1]->z() );
p2->setX( mPointVector[ptnr2]->x() );
p2->setY( mPointVector[ptnr2]->y() );
p2->setZ( mPointVector[ptnr2]->z() );
p3->setX( mPointVector[ptnr3]->x() );
p3->setY( mPointVector[ptnr3]->y() );
p3->setZ( mPointVector[ptnr3]->z() );
return true;
}
else if ( edge == -20 )//the point is exactly on an edge
@ -1122,15 +1122,15 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, Point3
{
return false;
}
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
p2->setX( mPointVector[ptnr2]->getX() );
p2->setY( mPointVector[ptnr2]->getY() );
p2->setZ( mPointVector[ptnr2]->getZ() );
p3->setX( mPointVector[ptnr3]->getX() );
p3->setY( mPointVector[ptnr3]->getY() );
p3->setZ( mPointVector[ptnr3]->getZ() );
p1->setX( mPointVector[ptnr1]->x() );
p1->setY( mPointVector[ptnr1]->y() );
p1->setZ( mPointVector[ptnr1]->z() );
p2->setX( mPointVector[ptnr2]->x() );
p2->setY( mPointVector[ptnr2]->y() );
p2->setZ( mPointVector[ptnr2]->z() );
p3->setX( mPointVector[ptnr3]->x() );
p3->setY( mPointVector[ptnr3]->y() );
p3->setZ( mPointVector[ptnr3]->z() );
return true;
}
else if ( edge == -25 )//x and y are the coordinates of an existing point
@ -1145,15 +1145,15 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, Point3
{
return false;
}
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
p2->setX( mPointVector[ptnr2]->getX() );
p2->setY( mPointVector[ptnr2]->getY() );
p2->setZ( mPointVector[ptnr2]->getZ() );
p3->setX( mPointVector[ptnr3]->getX() );
p3->setY( mPointVector[ptnr3]->getY() );
p3->setZ( mPointVector[ptnr3]->getZ() );
p1->setX( mPointVector[ptnr1]->x() );
p1->setY( mPointVector[ptnr1]->y() );
p1->setZ( mPointVector[ptnr1]->z() );
p2->setX( mPointVector[ptnr2]->x() );
p2->setY( mPointVector[ptnr2]->y() );
p2->setZ( mPointVector[ptnr2]->z() );
p3->setX( mPointVector[ptnr3]->x() );
p3->setY( mPointVector[ptnr3]->y() );
p3->setZ( mPointVector[ptnr3]->z() );
return true;
}
else if ( edge == -5 )//numerical problems in 'baseEdgeOfTriangle'
@ -1165,15 +1165,15 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D *p1, Point3
{
return false;
}
p1->setX( mPointVector[ptnr1]->getX() );
p1->setY( mPointVector[ptnr1]->getY() );
p1->setZ( mPointVector[ptnr1]->getZ() );
p2->setX( mPointVector[ptnr2]->getX() );
p2->setY( mPointVector[ptnr2]->getY() );
p2->setZ( mPointVector[ptnr2]->getZ() );
p3->setX( mPointVector[ptnr3]->getX() );
p3->setY( mPointVector[ptnr3]->getY() );
p3->setZ( mPointVector[ptnr3]->getZ() );
p1->setX( mPointVector[ptnr1]->x() );
p1->setY( mPointVector[ptnr1]->y() );
p1->setZ( mPointVector[ptnr1]->z() );
p2->setX( mPointVector[ptnr2]->x() );
p2->setY( mPointVector[ptnr2]->y() );
p2->setZ( mPointVector[ptnr2]->z() );
p3->setX( mPointVector[ptnr3]->x() );
p3->setY( mPointVector[ptnr3]->y() );
p3->setZ( mPointVector[ptnr3]->z() );
return true;
}
else//problems
@ -1249,7 +1249,7 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
}
//test, if the forced segment is a multiple of actedge and if the direction is the same
else if ( /*lines are parallel*/( mPointVector[p2]->getY() - mPointVector[p1]->getY() ) / ( mPointVector[mHalfEdge[actedge]->getPoint()]->getY() - mPointVector[p1]->getY() ) == ( mPointVector[p2]->getX() - mPointVector[p1]->getX() ) / ( mPointVector[mHalfEdge[actedge]->getPoint()]->getX() - mPointVector[p1]->getX() ) && ( ( mPointVector[p2]->getY() - mPointVector[p1]->getY() ) >= 0 ) == ( ( mPointVector[mHalfEdge[actedge]->getPoint()]->getY() - mPointVector[p1]->getY() ) > 0 ) && ( ( mPointVector[p2]->getX() - mPointVector[p1]->getX() ) >= 0 ) == ( ( mPointVector[mHalfEdge[actedge]->getPoint()]->getX() - mPointVector[p1]->getX() ) > 0 ) )
else if ( /*lines are parallel*/( mPointVector[p2]->y() - mPointVector[p1]->y() ) / ( mPointVector[mHalfEdge[actedge]->getPoint()]->y() - mPointVector[p1]->y() ) == ( mPointVector[p2]->x() - mPointVector[p1]->x() ) / ( mPointVector[mHalfEdge[actedge]->getPoint()]->x() - mPointVector[p1]->x() ) && ( ( mPointVector[p2]->y() - mPointVector[p1]->y() ) >= 0 ) == ( ( mPointVector[mHalfEdge[actedge]->getPoint()]->y() - mPointVector[p1]->y() ) > 0 ) && ( ( mPointVector[p2]->x() - mPointVector[p1]->x() ) >= 0 ) == ( ( mPointVector[mHalfEdge[actedge]->getPoint()]->x() - mPointVector[p1]->x() ) > 0 ) )
{
//mark actedge and Dual(actedge) as forced, reset p1 and start the method from the beginning
mHalfEdge[actedge]->setForced( true );
@ -1270,13 +1270,13 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
{
if ( mHalfEdge[mHalfEdge[actedge]->getNext()]->getForced() && mForcedCrossBehavior == Triangulation::SnappingTypeVertex )//if the crossed edge is a forced edge, we have to snap the forced line to the next node
{
Point3D crosspoint;
QgsPoint crosspoint;
int p3, p4;
p3 = mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getDual()]->getPoint();
MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p4], &crosspoint );
double dista = sqrt( ( crosspoint.getX() - mPointVector[p3]->getX() ) * ( crosspoint.getX() - mPointVector[p3]->getX() ) + ( crosspoint.getY() - mPointVector[p3]->getY() ) * ( crosspoint.getY() - mPointVector[p3]->getY() ) );
double distb = sqrt( ( crosspoint.getX() - mPointVector[p4]->getX() ) * ( crosspoint.getX() - mPointVector[p4]->getX() ) + ( crosspoint.getY() - mPointVector[p4]->getY() ) * ( crosspoint.getY() - mPointVector[p4]->getY() ) );
double dista = sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
double distb = sqrt( ( crosspoint.x() - mPointVector[p4]->x() ) * ( crosspoint.x() - mPointVector[p4]->x() ) + ( crosspoint.y() - mPointVector[p4]->y() ) * ( crosspoint.y() - mPointVector[p4]->y() ) );
if ( dista <= distb )
{
insertForcedSegment( p1, p3, breakline );
@ -1292,13 +1292,13 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
}
else if ( mHalfEdge[mHalfEdge[actedge]->getNext()]->getForced() && mForcedCrossBehavior == Triangulation::InsertVertex )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
{
Point3D crosspoint;
QgsPoint crosspoint;
int p3, p4;
p3 = mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getDual()]->getPoint();
MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p4], &crosspoint );
double distpart = sqrt( ( crosspoint.getX() - mPointVector[p4]->getX() ) * ( crosspoint.getX() - mPointVector[p4]->getX() ) + ( crosspoint.getY() - mPointVector[p4]->getY() ) * ( crosspoint.getY() - mPointVector[p4]->getY() ) );
double disttot = sqrt( ( mPointVector[p3]->getX() - mPointVector[p4]->getX() ) * ( mPointVector[p3]->getX() - mPointVector[p4]->getX() ) + ( mPointVector[p3]->getY() - mPointVector[p4]->getY() ) * ( mPointVector[p3]->getY() - mPointVector[p4]->getY() ) );
double distpart = sqrt( ( crosspoint.x() - mPointVector[p4]->x() ) * ( crosspoint.x() - mPointVector[p4]->x() ) + ( crosspoint.y() - mPointVector[p4]->y() ) * ( crosspoint.y() - mPointVector[p4]->y() ) );
double disttot = sqrt( ( mPointVector[p3]->x() - mPointVector[p4]->x() ) * ( mPointVector[p3]->x() - mPointVector[p4]->x() ) + ( mPointVector[p3]->y() - mPointVector[p4]->y() ) * ( mPointVector[p3]->y() - mPointVector[p4]->y() ) );
float frac = distpart / disttot;
if ( frac == 0 || frac == 1 )//just in case MathUtils::lineIntersection does not work as expected
@ -1358,13 +1358,13 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
{
if ( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getForced() && mForcedCrossBehavior == Triangulation::SnappingTypeVertex )//if the crossed edge is a forced edge and mForcedCrossBehavior is SnappingType_VERTICE, we have to snap the forced line to the next node
{
Point3D crosspoint;
QgsPoint crosspoint;
int p3, p4;
p3 = mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p4], &crosspoint );
double dista = sqrt( ( crosspoint.getX() - mPointVector[p3]->getX() ) * ( crosspoint.getX() - mPointVector[p3]->getX() ) + ( crosspoint.getY() - mPointVector[p3]->getY() ) * ( crosspoint.getY() - mPointVector[p3]->getY() ) );
double distb = sqrt( ( crosspoint.getX() - mPointVector[p4]->getX() ) * ( crosspoint.getX() - mPointVector[p4]->getX() ) + ( crosspoint.getY() - mPointVector[p4]->getY() ) * ( crosspoint.getY() - mPointVector[p4]->getY() ) );
double dista = sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
double distb = sqrt( ( crosspoint.x() - mPointVector[p4]->x() ) * ( crosspoint.x() - mPointVector[p4]->x() ) + ( crosspoint.y() - mPointVector[p4]->y() ) * ( crosspoint.y() - mPointVector[p4]->y() ) );
if ( dista <= distb )
{
insertForcedSegment( p1, p3, breakline );
@ -1380,13 +1380,13 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
}
else if ( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getForced() && mForcedCrossBehavior == Triangulation::InsertVertex )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
{
Point3D crosspoint;
QgsPoint crosspoint;
int p3, p4;
p3 = mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p4], &crosspoint );
double distpart = sqrt( ( crosspoint.getX() - mPointVector[p3]->getX() ) * ( crosspoint.getX() - mPointVector[p3]->getX() ) + ( crosspoint.getY() - mPointVector[p3]->getY() ) * ( crosspoint.getY() - mPointVector[p3]->getY() ) );
double disttot = sqrt( ( mPointVector[p3]->getX() - mPointVector[p4]->getX() ) * ( mPointVector[p3]->getX() - mPointVector[p4]->getX() ) + ( mPointVector[p3]->getY() - mPointVector[p4]->getY() ) * ( mPointVector[p3]->getY() - mPointVector[p4]->getY() ) );
double distpart = sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
double disttot = sqrt( ( mPointVector[p3]->x() - mPointVector[p4]->x() ) * ( mPointVector[p3]->x() - mPointVector[p4]->x() ) + ( mPointVector[p3]->y() - mPointVector[p4]->y() ) * ( mPointVector[p3]->y() - mPointVector[p4]->y() ) );
float frac = distpart / disttot;
if ( frac == 0 || frac == 1 )
{
@ -1405,13 +1405,13 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
{
if ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getForced() && mForcedCrossBehavior == Triangulation::SnappingTypeVertex )//if the crossed edge is a forced edge and mForcedCrossBehavior is SnappingType_VERTICE, we have to snap the forced line to the next node
{
Point3D crosspoint;
QgsPoint crosspoint;
int p3, p4;
p3 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint();
MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p4], &crosspoint );
double dista = sqrt( ( crosspoint.getX() - mPointVector[p3]->getX() ) * ( crosspoint.getX() - mPointVector[p3]->getX() ) + ( crosspoint.getY() - mPointVector[p3]->getY() ) * ( crosspoint.getY() - mPointVector[p3]->getY() ) );
double distb = sqrt( ( crosspoint.getX() - mPointVector[p4]->getX() ) * ( crosspoint.getX() - mPointVector[p4]->getX() ) + ( crosspoint.getY() - mPointVector[p4]->getY() ) * ( crosspoint.getY() - mPointVector[p4]->getY() ) );
double dista = sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
double distb = sqrt( ( crosspoint.x() - mPointVector[p4]->x() ) * ( crosspoint.x() - mPointVector[p4]->x() ) + ( crosspoint.y() - mPointVector[p4]->y() ) * ( crosspoint.y() - mPointVector[p4]->y() ) );
if ( dista <= distb )
{
insertForcedSegment( p1, p3, breakline );
@ -1427,13 +1427,13 @@ int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
}
else if ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getForced() && mForcedCrossBehavior == Triangulation::InsertVertex )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
{
Point3D crosspoint;
QgsPoint crosspoint;
int p3, p4;
p3 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint();
MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p4], &crosspoint );
double distpart = sqrt( ( crosspoint.getX() - mPointVector[p3]->getX() ) * ( crosspoint.getX() - mPointVector[p3]->getX() ) + ( crosspoint.getY() - mPointVector[p3]->getY() ) * ( crosspoint.getY() - mPointVector[p3]->getY() ) );
double disttot = sqrt( ( mPointVector[p3]->getX() - mPointVector[p4]->getX() ) * ( mPointVector[p3]->getX() - mPointVector[p4]->getX() ) + ( mPointVector[p3]->getY() - mPointVector[p4]->getY() ) * ( mPointVector[p3]->getY() - mPointVector[p4]->getY() ) );
double distpart = sqrt( ( crosspoint.x() - mPointVector[p3]->x() ) * ( crosspoint.x() - mPointVector[p3]->x() ) + ( crosspoint.y() - mPointVector[p3]->y() ) * ( crosspoint.y() - mPointVector[p3]->y() ) );
double disttot = sqrt( ( mPointVector[p3]->x() - mPointVector[p4]->x() ) * ( mPointVector[p3]->x() - mPointVector[p4]->x() ) + ( mPointVector[p3]->y() - mPointVector[p4]->y() ) * ( mPointVector[p3]->y() - mPointVector[p4]->y() ) );
float frac = distpart / disttot;
if ( frac == 0 || frac == 1 )
{
@ -1636,24 +1636,24 @@ void DualEdgeTriangulation::eliminateHorizontalTriangles()
}
double el1, el2, el3;//elevations of the points
el1 = mPointVector[p1]->getZ();
el2 = mPointVector[p2]->getZ();
el3 = mPointVector[p3]->getZ();
el1 = mPointVector[p1]->z();
el2 = mPointVector[p2]->z();
el3 = mPointVector[p3]->z();
if ( el1 == el2 && el2 == el3 )//we found a horizonal triangle
{
//swap edges if it is possible, if it would remove the horizontal triangle and if the minimum angle generated by the swap is high enough
if ( swapPossible( ( uint )e1 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e1]->getDual()]->getNext()]->getPoint()]->getZ() != el1 && swapMinAngle( e1 ) > minangle )
if ( swapPossible( ( uint )e1 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e1]->getDual()]->getNext()]->getPoint()]->z() != el1 && swapMinAngle( e1 ) > minangle )
{
doOnlySwap( ( uint )e1 );
swapped = true;
}
else if ( swapPossible( ( uint )e2 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e2]->getDual()]->getNext()]->getPoint()]->getZ() != el2 && swapMinAngle( e2 ) > minangle )
else if ( swapPossible( ( uint )e2 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e2]->getDual()]->getNext()]->getPoint()]->z() != el2 && swapMinAngle( e2 ) > minangle )
{
doOnlySwap( ( uint )e2 );
swapped = true;
}
else if ( swapPossible( ( uint )e3 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e3]->getDual()]->getNext()]->getPoint()]->getZ() != el3 && swapMinAngle( e3 ) > minangle )
else if ( swapPossible( ( uint )e3 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e3]->getDual()]->getNext()]->getPoint()]->z() != el3 && swapMinAngle( e3 ) > minangle )
{
doOnlySwap( ( uint )e3 );
swapped = true;
@ -1769,7 +1769,7 @@ void DualEdgeTriangulation::ruppertRefinement()
int minedgenext;
int minedgenextnext;
Point3D circumcenter;
QgsPoint circumcenter;
while ( !edge_angle.empty() )
{
@ -1795,11 +1795,11 @@ void DualEdgeTriangulation::ruppertRefinement()
continue;
}
if ( !pointInside( circumcenter.getX(), circumcenter.getY() ) )
if ( !pointInside( circumcenter.x(), circumcenter.y() ) )
{
//put all three edges to dontexamine and remove them from the other maps
QgsDebugMsg( QString( "put circumcenter %1//%2 on dontexamine list because it is outside the convex hull" )
.arg( circumcenter.getX() ).arg( circumcenter.getY() ) );
.arg( circumcenter.x() ).arg( circumcenter.y() ) );
dontexamine.insert( minedge );
edge_angle.erase( minedge );
std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
@ -2205,8 +2205,8 @@ void DualEdgeTriangulation::ruppertRefinement()
/*******otherwise, try to add the circumcenter to the triangulation************************************************************************************************/
Point3D *p = new Point3D();
mDecorator->calcPoint( circumcenter.getX(), circumcenter.getY(), p );
QgsPoint *p = new QgsPoint();
mDecorator->calcPoint( circumcenter.x(), circumcenter.y(), p );
int pointno = mDecorator->addPoint( p );
if ( pointno == -100 || pointno == mTwiceInsPoint )
@ -2214,17 +2214,17 @@ void DualEdgeTriangulation::ruppertRefinement()
if ( pointno == -100 )
{
QgsDebugMsg( QString( "put circumcenter %1//%2 on dontexamine list because of numerical instabilities" )
.arg( circumcenter.getX() ).arg( circumcenter.getY() ) );
.arg( circumcenter.x() ).arg( circumcenter.y() ) );
}
else if ( pointno == mTwiceInsPoint )
{
QgsDebugMsg( QString( "put circumcenter %1//%2 on dontexamine list because it is already inserted" )
.arg( circumcenter.getX() ).arg( circumcenter.getY() ) );
.arg( circumcenter.x() ).arg( circumcenter.y() ) );
//test, if the point is present in the triangulation
bool flag = false;
for ( int i = 0; i < mPointVector.count(); i++ )
{
if ( mPointVector[i]->getX() == circumcenter.getX() && mPointVector[i]->getY() == circumcenter.getY() )
if ( mPointVector[i]->x() == circumcenter.x() && mPointVector[i]->y() == circumcenter.y() )
{
flag = true;
}
@ -2435,10 +2435,10 @@ bool DualEdgeTriangulation::swapPossible( unsigned int edge )
return false;
}
//then, test, if the edge is in the middle of a not convex quad
Point3D *pta = mPointVector[mHalfEdge[edge]->getPoint()];
Point3D *ptb = mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()];
Point3D *ptc = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint()];
Point3D *ptd = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint()];
QgsPoint *pta = mPointVector[mHalfEdge[edge]->getPoint()];
QgsPoint *ptb = mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()];
QgsPoint *ptc = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getNext()]->getNext()]->getPoint()];
QgsPoint *ptd = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint()];
if ( MathUtils::leftOf( ptc, pta, ptb ) > leftOfTresh )
{
return false;
@ -2597,7 +2597,7 @@ void DualEdgeTriangulation::triangulatePolygon( QList<int> *poly, QList<int> *fr
bool DualEdgeTriangulation::pointInside( double x, double y )
{
Point3D point( x, y, 0 );
QgsPoint point( x, y, 0 );
unsigned int actedge = mEdgeInside;//start with an edge which does not point to the virtual point
int counter = 0;//number of consecutive successful left-of-tests
int nulls = 0;//number of left-of-tests, which returned 0. 1 means, that the point is on a line, 2 means that it is on an existing point
@ -2860,7 +2860,7 @@ bool DualEdgeTriangulation::readFromTAFF( QString filename )
textstream >> y;
textstream >> z;
Point3D *p = new Point3D( x, y, z );
QgsPoint *p = new QgsPoint( x, y, z );
// QgsDebugMsg( QString( "inserting point %1" ).arg( i ) );
mPointVector.insert( i, p );
@ -2951,7 +2951,7 @@ bool DualEdgeTriangulation::saveToTAFF( QString filename ) const
for ( int i = 0; i < getNumberOfPoints(); i++ )
{
Point3D *p = mPointVector[i];
QgsPoint *p = mPointVector[i];
outstream << p->getX() << " " << p->getY() << " " << p->getZ() << " ";
}
outstream << std::endl << std::flush;
@ -2963,14 +2963,14 @@ bool DualEdgeTriangulation::saveToTAFF( QString filename ) const
bool DualEdgeTriangulation::swapEdge( double x, double y )
{
Point3D p( x, y, 0 );
QgsPoint p( x, y, 0 );
int edge1 = baseEdgeOfTriangle( &p );
if ( edge1 >= 0 )
{
int edge2, edge3;
Point3D *point1 = nullptr;
Point3D *point2 = nullptr;
Point3D *point3 = nullptr;
QgsPoint *point1 = nullptr;
QgsPoint *point2 = nullptr;
QgsPoint *point3 = nullptr;
edge2 = mHalfEdge[edge1]->getNext();
edge3 = mHalfEdge[edge2]->getNext();
point1 = getPoint( mHalfEdge[edge1]->getPoint() );
@ -3024,15 +3024,15 @@ bool DualEdgeTriangulation::swapEdge( double x, double y )
QList<int> *DualEdgeTriangulation::getPointsAroundEdge( double x, double y )
{
Point3D p( x, y, 0 );
QgsPoint p( x, y, 0 );
int p1, p2, p3, p4;
int edge1 = baseEdgeOfTriangle( &p );
if ( edge1 >= 0 )
{
int edge2, edge3;
Point3D *point1 = nullptr;
Point3D *point2 = nullptr;
Point3D *point3 = nullptr;
QgsPoint *point1 = nullptr;
QgsPoint *point2 = nullptr;
QgsPoint *point3 = nullptr;
edge2 = mHalfEdge[edge1]->getNext();
edge3 = mHalfEdge[edge2]->getNext();
point1 = getPoint( mHalfEdge[edge1]->getPoint() );
@ -3133,11 +3133,11 @@ bool DualEdgeTriangulation::saveAsShapefile( const QString &fileName ) const
QgsFeature edgeLineFeature;
//geometry
Point3D *p1 = mPointVector[currentEdge->getPoint()];
Point3D *p2 = mPointVector[mHalfEdge[currentEdge->getDual()]->getPoint()];
QgsPoint *p1 = mPointVector[currentEdge->getPoint()];
QgsPoint *p2 = mPointVector[mHalfEdge[currentEdge->getDual()]->getPoint()];
QgsPolyline lineGeom;
lineGeom.push_back( QgsPointXY( p1->getX(), p1->getY() ) );
lineGeom.push_back( QgsPointXY( p2->getX(), p2->getY() ) );
lineGeom.push_back( QgsPointXY( p1->x(), p1->y() ) );
lineGeom.push_back( QgsPointXY( p2->x(), p2->y() ) );
edgeLineFeature.setGeometry( QgsGeometry::fromPolyline( lineGeom ) );
edgeLineFeature.initAttributes( 1 );
@ -3168,10 +3168,10 @@ bool DualEdgeTriangulation::saveAsShapefile( const QString &fileName ) const
double DualEdgeTriangulation::swapMinAngle( int edge ) const
{
Point3D *p1 = getPoint( mHalfEdge[edge]->getPoint() );
Point3D *p2 = getPoint( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() );
Point3D *p3 = getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() );
Point3D *p4 = getPoint( mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() );
QgsPoint *p1 = getPoint( mHalfEdge[edge]->getPoint() );
QgsPoint *p2 = getPoint( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() );
QgsPoint *p3 = getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() );
QgsPoint *p4 = getPoint( mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() );
//search for the minimum angle (it is important, which directions the lines have!)
double minangle;
@ -3215,19 +3215,19 @@ int DualEdgeTriangulation::splitHalfEdge( int edge, float position )
}
//create the new point on the heap
Point3D *p = new Point3D( mPointVector[mHalfEdge[edge]->getPoint()]->getX()*position + mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()]->getX() * ( 1 - position ), mPointVector[mHalfEdge[edge]->getPoint()]->getY()*position + mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()]->getY() * ( 1 - position ), 0 );
QgsPoint *p = new QgsPoint( mPointVector[mHalfEdge[edge]->getPoint()]->x()*position + mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()]->x() * ( 1 - position ), mPointVector[mHalfEdge[edge]->getPoint()]->y()*position + mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()]->y() * ( 1 - position ), 0 );
//calculate the z-value of the point to insert
Point3D zvaluepoint;
mDecorator->calcPoint( p->getX(), p->getY(), &zvaluepoint );
p->setZ( zvaluepoint.getZ() );
QgsPoint zvaluepoint;
mDecorator->calcPoint( p->x(), p->y(), &zvaluepoint );
p->setZ( zvaluepoint.z() );
//insert p into mPointVector
if ( mPointVector.count() >= mPointVector.size() )
{
mPointVector.resize( mPointVector.count() + 1 );
}
QgsDebugMsg( QString( "inserting point nr. %1, %2//%3//%4" ).arg( mPointVector.count() ).arg( p->getX() ).arg( p->getY() ).arg( p->getZ() ) );
QgsDebugMsg( QString( "inserting point nr. %1, %2//%3//%4" ).arg( mPointVector.count() ).arg( p->x() ).arg( p->y() ).arg( p->z() ) );
mPointVector.insert( mPointVector.count(), p );
//insert the six new halfedges
@ -3256,7 +3256,7 @@ int DualEdgeTriangulation::splitHalfEdge( int edge, float position )
checkSwap( mHalfEdge[dualedge]->getNext(), 0 );
checkSwap( mHalfEdge[edge3]->getNext(), 0 );
mDecorator->addPoint( new Point3D( p->getX(), p->getY(), 0 ) );//dirty hack to enforce update of decorators
mDecorator->addPoint( new QgsPoint( p->x(), p->y(), 0 ) );//dirty hack to enforce update of decorators
return mPointVector.count() - 1;
}
@ -3266,7 +3266,7 @@ bool DualEdgeTriangulation::edgeOnConvexHull( int edge )
return ( mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() == -1 );
}
void DualEdgeTriangulation::evaluateInfluenceRegion( Point3D *point, int edge, QSet<int> &set )
void DualEdgeTriangulation::evaluateInfluenceRegion( QgsPoint *point, int edge, QSet<int> &set )
{
if ( set.find( edge ) == set.end() )
{

34
src/analysis/interpolation/DualEdgeTriangulation.h
View File

@ -46,24 +46,24 @@ class ANALYSIS_EXPORT DualEdgeTriangulation: public Triangulation
//! Adds a line (e.g. a break-, structure- or an isoline) to the triangulation. The class takes ownership of the line object and its points
void addLine( Line3D *line SIP_TRANSFER, bool breakline ) override;
//! Adds a point to the triangulation and returns the number of this point in case of success or -100 in case of failure
int addPoint( Point3D *p ) override;
int addPoint( QgsPoint *p ) override;
//! Performs a consistency check, remove this later
virtual void performConsistencyTest() override;
//! Calculates the normal at a point on the surface
virtual bool calcNormal( double x, double y, Vector3D *result ) override;
//! Calculates x-, y and z-value of the point on the surface
virtual bool calcPoint( double x, double y, Point3D *result ) override;
virtual bool calcPoint( double x, double y, QgsPoint *result ) override;
//! Draws the points, edges and the forced lines
//virtual void draw(QPainter* p, double xlowleft, double ylowleft, double xupright, double yupright, double width, double height) const;
//! Returns a pointer to the point with number i
virtual Point3D *getPoint( unsigned int i ) const override;
virtual QgsPoint *getPoint( unsigned int i ) const override;
//! Returns the number of the point opposite to the triangle points p1, p2 (which have to be on a halfedge)
int getOppositePoint( int p1, int p2 ) override;
//! Finds out, in which triangle the point with coordinates x and y is and assigns the numbers of the vertices to 'n1', 'n2' and 'n3' and the vertices to 'p1', 'p2' and 'p3'
//! \note not available in Python bindings
virtual bool getTriangle( double x, double y, Point3D *p1, int *n1, Point3D *p2, int *n2, Point3D *p3, int *n3 ) override SIP_SKIP;
virtual bool getTriangle( double x, double y, QgsPoint *p1, int *n1, QgsPoint *p2, int *n2, QgsPoint *p3, int *n3 ) override SIP_SKIP;
//! Finds out, in which triangle the point with coordinates x and y is and assigns addresses to the points at the vertices to 'p1', 'p2' and 'p3
virtual bool getTriangle( double x, double y, Point3D *p1 SIP_OUT, Point3D *p2 SIP_OUT, Point3D *p3 SIP_OUT ) override;
virtual bool getTriangle( double x, double y, QgsPoint *p1 SIP_OUT, QgsPoint *p2 SIP_OUT, QgsPoint *p3 SIP_OUT ) override;
//! Returns a pointer to a value list with the information of the triangles surrounding (counterclockwise) a point. Four integer values describe a triangle, the first three are the number of the half edges of the triangle and the fourth is -10, if the third (and most counterclockwise) edge is a breakline, and -20 otherwise. The value list has to be deleted by the code which called the method
QList<int> *getSurroundingTriangles( int pointno ) override;
//! Returns the largest x-coordinate value of the bounding box
@ -117,7 +117,7 @@ class ANALYSIS_EXPORT DualEdgeTriangulation: public Triangulation
//! Default value for the number of storable points at the beginning
static const unsigned int DEFAULT_STORAGE_FOR_POINTS = 100000;
//! Stores pointers to all points in the triangulations (including the points contained in the lines)
QVector<Point3D *> mPointVector;
QVector<QgsPoint *> mPointVector;
//! Default value for the number of storable HalfEdges at the beginning
static const unsigned int DEFAULT_STORAGE_FOR_HALF_EDGES = 300006;
//! Stores pointers to the HalfEdges
@ -145,7 +145,7 @@ class ANALYSIS_EXPORT DualEdgeTriangulation: public Triangulation
//! Returns the number of an edge which points to the point with number 'point' or -1 if there is an error
int baseEdgeOfPoint( int point );
//! Returns the number of a HalfEdge from a triangle in which 'point' is in. If the number -10 is returned, this means, that 'point' is outside the convex hull. If -5 is returned, then numerical problems with the leftOfTest occurred (and the value of the possible edge is stored in the variable 'mUnstableEdge'. -20 means, that the inserted point is exactly on an edge (the number is stored in the variable 'mEdgeWithPoint'). -25 means, that the point is already in the triangulation (the number of the point is stored in the member 'mTwiceInsPoint'. If -100 is returned, this means that something else went wrong
int baseEdgeOfTriangle( Point3D *point );
int baseEdgeOfTriangle( QgsPoint *point );
//! Checks, if 'edge' has to be swapped because of the empty circle criterion. If so, doSwap(...) is called.
bool checkSwap( unsigned int edge, unsigned int recursiveDeep );
//! Swaps 'edge' and test recursively for other swaps (delaunay criterion)
@ -175,7 +175,7 @@ class ANALYSIS_EXPORT DualEdgeTriangulation: public Triangulation
//! Returns true, if a half edge is on the convex hull and false otherwise
bool edgeOnConvexHull( int edge );
//! Function needed for the ruppert algorithm. Tests, if point is in the circle through both endpoints of edge and the endpoint of edge->dual->next->point. If so, the function calls itself recursively for edge->next and edge->next->next. Stops, if it finds a forced edge or a convex hull edge
void evaluateInfluenceRegion( Point3D *point, int edge, QSet<int> &set );
void evaluateInfluenceRegion( QgsPoint *point, int edge, QSet<int> &set );
};
inline DualEdgeTriangulation::DualEdgeTriangulation()
@ -225,7 +225,7 @@ inline int DualEdgeTriangulation::getNumberOfPoints() const
return mPointVector.count();
}
inline Point3D *DualEdgeTriangulation::getPoint( unsigned int i ) const
inline QgsPoint *DualEdgeTriangulation::getPoint( unsigned int i ) const
{
return mPointVector.at( i );
}
@ -233,14 +233,14 @@ inline Point3D *DualEdgeTriangulation::getPoint( unsigned int i ) const
inline bool DualEdgeTriangulation::halfEdgeBBoxTest( int edge, double xlowleft, double ylowleft, double xupright, double yupright ) const
{
return (
( getPoint( mHalfEdge[edge]->getPoint() )->getX() >= xlowleft &&
getPoint( mHalfEdge[edge]->getPoint() )->getX() <= xupright &&
getPoint( mHalfEdge[edge]->getPoint() )->getY() >= ylowleft &&
getPoint( mHalfEdge[edge]->getPoint() )->getY() <= yupright ) ||
( getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->getX() >= xlowleft &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->getX() <= xupright &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->getY() >= ylowleft &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->getY() <= yupright )
( getPoint( mHalfEdge[edge]->getPoint() )->x() >= xlowleft &&
getPoint( mHalfEdge[edge]->getPoint() )->x() <= xupright &&
getPoint( mHalfEdge[edge]->getPoint() )->y() >= ylowleft &&
getPoint( mHalfEdge[edge]->getPoint() )->y() <= yupright ) ||
( getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->x() >= xlowleft &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->x() <= xupright &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->y() >= ylowleft &&
getPoint( mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() )->y() <= yupright )
);
}

30
src/analysis/interpolation/LinTriangleInterpolator.cc
View File

@ -22,9 +22,9 @@ bool LinTriangleInterpolator::calcFirstDerX( double x, double y, Vector3D *vec )
if ( vec && mTIN )
{
Point3D pt1( 0, 0, 0 );
Point3D pt2( 0, 0, 0 );
Point3D pt3( 0, 0, 0 );
QgsPoint pt1( 0, 0, 0 );
QgsPoint pt2( 0, 0, 0 );
QgsPoint pt3( 0, 0, 0 );
if ( !mTIN->getTriangle( x, y, &pt1, &pt2, &pt3 ) )
{
@ -33,7 +33,7 @@ bool LinTriangleInterpolator::calcFirstDerX( double x, double y, Vector3D *vec )
vec->setX( 1.0 );
vec->setY( 0.0 );
vec->setZ( ( pt1.getZ() * ( pt2.getY() - pt3.getY() ) + pt2.getZ() * ( pt3.getY() - pt1.getY() ) + pt3.getZ() * ( pt1.getY() - pt2.getY() ) ) / ( ( pt1.getX() - pt2.getX() ) * ( pt2.getY() - pt3.getY() ) - ( pt2.getX() - pt3.getX() ) * ( pt1.getY() - pt2.getY() ) ) );
vec->setZ( ( pt1.z() * ( pt2.y() - pt3.y() ) + pt2.z() * ( pt3.y() - pt1.y() ) + pt3.z() * ( pt1.y() - pt2.y() ) ) / ( ( pt1.x() - pt2.x() ) * ( pt2.y() - pt3.y() ) - ( pt2.x() - pt3.x() ) * ( pt1.y() - pt2.y() ) ) );
return true;
}
@ -48,9 +48,9 @@ bool LinTriangleInterpolator::calcFirstDerY( double x, double y, Vector3D *vec )
{
if ( vec && mTIN )
{
Point3D pt1( 0, 0, 0 );
Point3D pt2( 0, 0, 0 );
Point3D pt3( 0, 0, 0 );
QgsPoint pt1( 0, 0, 0 );
QgsPoint pt2( 0, 0, 0 );
QgsPoint pt3( 0, 0, 0 );
if ( !mTIN->getTriangle( x, y, &pt1, &pt2, &pt3 ) )
{
@ -59,7 +59,7 @@ bool LinTriangleInterpolator::calcFirstDerY( double x, double y, Vector3D *vec )
vec->setX( 0 );
vec->setY( 1.0 );
vec->setZ( ( pt1.getZ() * ( pt2.getX() - pt3.getX() ) + pt2.getZ() * ( pt3.getX() - pt1.getX() ) + pt3.getZ() * ( pt1.getX() - pt2.getX() ) ) / ( ( pt1.getY() - pt2.getY() ) * ( pt2.getX() - pt3.getX() ) - ( pt2.getY() - pt3.getY() ) * ( pt1.getX() - pt2.getX() ) ) );
vec->setZ( ( pt1.z() * ( pt2.x() - pt3.x() ) + pt2.z() * ( pt3.x() - pt1.x() ) + pt3.z() * ( pt1.x() - pt2.x() ) ) / ( ( pt1.y() - pt2.y() ) * ( pt2.x() - pt3.x() ) - ( pt2.y() - pt3.y() ) * ( pt1.x() - pt2.x() ) ) );
return true;
}
@ -98,22 +98,22 @@ bool LinTriangleInterpolator::calcNormVec( double x, double y, Vector3D *vec )
}
bool LinTriangleInterpolator::calcPoint( double x, double y, Point3D *point )
bool LinTriangleInterpolator::calcPoint( double x, double y, QgsPoint *point )
{
if ( point && mTIN )
{
Point3D pt1( 0, 0, 0 );
Point3D pt2( 0, 0, 0 );
Point3D pt3( 0, 0, 0 );
QgsPoint pt1( 0, 0, 0 );
QgsPoint pt2( 0, 0, 0 );
QgsPoint pt3( 0, 0, 0 );
if ( !mTIN->getTriangle( x, y, &pt1, &pt2, &pt3 ) )
{
return false;//point is outside the convex hull or numerical problems
}
double a = ( pt1.getZ() * ( pt2.getY() - pt3.getY() ) + pt2.getZ() * ( pt3.getY() - pt1.getY() ) + pt3.getZ() * ( pt1.getY() - pt2.getY() ) ) / ( ( pt1.getX() - pt2.getX() ) * ( pt2.getY() - pt3.getY() ) - ( pt2.getX() - pt3.getX() ) * ( pt1.getY() - pt2.getY() ) );
double b = ( pt1.getZ() * ( pt2.getX() - pt3.getX() ) + pt2.getZ() * ( pt3.getX() - pt1.getX() ) + pt3.getZ() * ( pt1.getX() - pt2.getX() ) ) / ( ( pt1.getY() - pt2.getY() ) * ( pt2.getX() - pt3.getX() ) - ( pt2.getY() - pt3.getY() ) * ( pt1.getX() - pt2.getX() ) );
double c = pt1.getZ() - a * pt1.getX() - b * pt1.getY();
double a = ( pt1.z() * ( pt2.y() - pt3.y() ) + pt2.z() * ( pt3.y() - pt1.y() ) + pt3.z() * ( pt1.y() - pt2.y() ) ) / ( ( pt1.x() - pt2.x() ) * ( pt2.y() - pt3.y() ) - ( pt2.x() - pt3.x() ) * ( pt1.y() - pt2.y() ) );
double b = ( pt1.z() * ( pt2.x() - pt3.x() ) + pt2.z() * ( pt3.x() - pt1.x() ) + pt3.z() * ( pt1.x() - pt2.x() ) ) / ( ( pt1.y() - pt2.y() ) * ( pt2.x() - pt3.x() ) - ( pt2.y() - pt3.y() ) * ( pt1.x() - pt2.x() ) );
double c = pt1.z() - a * pt1.x() - b * pt1.y();
point->setX( x );
point->setY( y );

2
src/analysis/interpolation/LinTriangleInterpolator.h
View File

@ -34,7 +34,7 @@ class ANALYSIS_EXPORT LinTriangleInterpolator : public TriangleInterpolator
//! Calculates the normal vector and assigns it to vec
virtual bool calcNormVec( double x, double y, Vector3D *result ) override;
//! Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point
virtual bool calcPoint( double x, double y, Point3D *result ) override;
virtual bool calcPoint( double x, double y, QgsPoint *result ) override;
//! Returns a pointer to the current Triangulation object
virtual DualEdgeTriangulation *getTriangulation() const;
//! Sets a Triangulation

2
src/analysis/interpolation/Line3D.cc
View File

@ -46,7 +46,7 @@ bool Line3D::empty() const
return ( head->getNext() == z );
}
void Line3D::insertPoint( Point3D *p )
void Line3D::insertPoint( QgsPoint *p )
{
if ( currentNode != z )//we can't insert a node behind the z-node
{

10
src/analysis/interpolation/Line3D.h
View File

@ -21,7 +21,7 @@
#include "qgis_analysis.h"
/** \ingroup analysis
* This class represents a line. It is implemented as a single directed linked list of nodes (with related Point3D objects). Attention: the points inserted in a line are not deleted from Line3D*/
* This class represents a line. It is implemented as a single directed linked list of nodes (with related QgsPoint objects). Attention: the points inserted in a line are not deleted from Line3D*/
class ANALYSIS_EXPORT Line3D
{
private:
@ -39,14 +39,14 @@ class ANALYSIS_EXPORT Line3D
public:
Line3D();
~Line3D();
//! Returns true, if the Line contains no Point3D, otherwise false
//! Returns true, if the Line contains no QgsPoint, otherwise false
bool empty() const;
//! Inserts a node behind the current position and sets the current position to this new node
void insertPoint( Point3D *p );
void insertPoint( QgsPoint *p );
//! Removes the point behind the current position
void removePoint();
//! Gets the point at the current position
Point3D *getPoint() const;
QgsPoint *getPoint() const;
//! Returns the current position
unsigned int getCurrent() const;
//! Returns the size of the line (the numbero of inserted Nodes without 'head' and 'z'
@ -72,7 +72,7 @@ inline unsigned int Line3D::getCurrent() const
return current;
}
inline Point3D *Line3D::getPoint() const
inline QgsPoint *Line3D::getPoint() const
{
return ( currentNode->getPoint() );
}

182
src/analysis/interpolation/MathUtils.cc
View File

@ -16,16 +16,16 @@
#include "MathUtils.h"
#include "qgslogger.h"
#include "Point3D.h"
#include "qgspoint.h"
#include "Vector3D.h"
#include <qmath.h>
bool MathUtils::calcBarycentricCoordinates( double x, double y, Point3D *p1, Point3D *p2, Point3D *p3, Point3D *result )
bool MathUtils::calcBarycentricCoordinates( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result )
{
if ( p1 && p2 && p3 && result )
{
Point3D p( x, y, 0 );
QgsPoint p( x, y, 0 );
double area = triArea( p1, p2, p3 );
if ( area == 0 )//p1, p2, p3 are in a line
{
@ -50,7 +50,7 @@ bool MathUtils::calcBarycentricCoordinates( double x, double y, Point3D *p1, Poi
}
}
bool MathUtils::BarycentricToXY( double u, double v, double w, Point3D *p1, Point3D *p2, Point3D *p3, Point3D *result )//this is wrong at the moment. Furthermore, the case, where the denominators are 0 have to be treated (two ways of calculating px and py)
bool MathUtils::BarycentricToXY( double u, double v, double w, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result )//this is wrong at the moment. Furthermore, the case, where the denominators are 0 have to be treated (two ways of calculating px and py)
{
Q_UNUSED( w );
@ -65,29 +65,29 @@ bool MathUtils::BarycentricToXY( double u, double v, double w, Point3D *p1, Poin
return false;
}
double denominator = ( ( p2->getY() - p3->getY() ) * ( p1->getX() - p3->getX() ) - ( p3->getY() - p1->getY() ) * ( p3->getX() - p2->getX() ) );
double denominator = ( ( p2->y() - p3->y() ) * ( p1->x() - p3->x() ) - ( p3->y() - p1->y() ) * ( p3->x() - p2->x() ) );
if ( denominator != 0 )//drop out py in the two equations
{
px = ( 2 * u * area * ( p1->getX() - p3->getX() ) - 2 * v * area * ( p3->getX() - p2->getX() ) - p2->getX() * p3->getY() * ( p1->getX() - p3->getX() ) + p3->getX() * p1->getY() * ( p3->getX() - p2->getX() ) + p3->getX() * p2->getY() * ( p1->getX() - p3->getX() ) - p1->getX() * p3->getY() * ( p3->getX() - p2->getX() ) ) / denominator;
if ( ( p3->getX() - p2->getX() ) != 0 )
px = ( 2 * u * area * ( p1->x() - p3->x() ) - 2 * v * area * ( p3->x() - p2->x() ) - p2->x() * p3->y() * ( p1->x() - p3->x() ) + p3->x() * p1->y() * ( p3->x() - p2->x() ) + p3->x() * p2->y() * ( p1->x() - p3->x() ) - p1->x() * p3->y() * ( p3->x() - p2->x() ) ) / denominator;
if ( ( p3->x() - p2->x() ) != 0 )
{
py = ( 2 * u * area - px * ( p2->getY() - p3->getY() ) - p2->getX() * p3->getY() + p3->getX() * p2->getY() ) / ( p3->getX() - p2->getX() );
py = ( 2 * u * area - px * ( p2->y() - p3->y() ) - p2->x() * p3->y() + p3->x() * p2->y() ) / ( p3->x() - p2->x() );
}
else
{
py = ( 2 * v * area - px * ( p3->getY() - p1->getY() ) - p3->getX() * p1->getY() + p1->getX() * p3->getY() ) / ( p1->getX() - p3->getX() );
py = ( 2 * v * area - px * ( p3->y() - p1->y() ) - p3->x() * p1->y() + p1->x() * p3->y() ) / ( p1->x() - p3->x() );
}
}
else//dorp out px in the two equations(maybe this possibility occurs only, if p1, p2 and p3 are coplanar
{
py = ( 2 * u * area * ( p3->getY() - p1->getY() ) - 2 * v * area * ( p2->getY() - p3->getY() ) - p2->getX() * p3->getY() * ( p3->getY() - p1->getY() ) + p3->getX() * p1->getY() * ( p2->getY() - p3->getY() ) + p3->getX() * p2->getY() * ( p3->getY() - p1->getY() ) - p1->getX() * p3->getY() * ( p2->getY() - p3->getY() ) ) / ( ( p3->getX() - p2->getX() ) * ( p3->getY() - p1->getY() ) - ( p1->getX() - p3->getX() ) * ( p2->getY() - p3->getY() ) );
if ( ( p2->getY() - p3->getY() ) != 0 )
py = ( 2 * u * area * ( p3->y() - p1->y() ) - 2 * v * area * ( p2->y() - p3->y() ) - p2->x() * p3->y() * ( p3->y() - p1->y() ) + p3->x() * p1->y() * ( p2->y() - p3->y() ) + p3->x() * p2->y() * ( p3->y() - p1->y() ) - p1->x() * p3->y() * ( p2->y() - p3->y() ) ) / ( ( p3->x() - p2->x() ) * ( p3->y() - p1->y() ) - ( p1->x() - p3->x() ) * ( p2->y() - p3->y() ) );
if ( ( p2->y() - p3->y() ) != 0 )
{
px = ( 2 * u * area - py * ( p3->getX() - p2->getX() ) - p2->getX() * p3->getY() + p3->getX() * p2->getY() ) / ( p2->getY() - p3->getY() );
px = ( 2 * u * area - py * ( p3->x() - p2->x() ) - p2->x() * p3->y() + p3->x() * p2->y() ) / ( p2->y() - p3->y() );
}
else
{
px = ( 2 * v * area - py * ( p1->getX() - p3->getX() ) - p3->getX() * p1->getY() + p1->getX() * p3->getY() ) / ( p3->getY() - p1->getY() );
px = ( 2 * v * area - py * ( p1->x() - p3->x() ) - p3->x() * p1->y() + p1->x() * p3->y() ) / ( p3->y() - p1->y() );
}
}
result->setX( px );
@ -116,20 +116,20 @@ double MathUtils::cFDerBernsteinPoly( int n, int i, double t )
return n * ( calcBernsteinPoly( n - 1, i - 1, t ) - calcBernsteinPoly( n - 1, i, t ) );
}
bool MathUtils::circumcenter( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *result )//version using the property that the distances from p1, p2, p3 to the circumcenter have to be equal. Possibly there is a bug
bool MathUtils::circumcenter( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result )//version using the property that the distances from p1, p2, p3 to the circumcenter have to be equal. Possibly there is a bug
{
if ( p1 && p2 && p3 && result )
{
double distp1p2 = sqrt( ( p1->getX() - p2->getX() ) * ( p1->getX() - p2->getX() ) + ( p1->getY() - p2->getY() ) * ( p1->getY() - p2->getY() ) );
double distp2p3 = sqrt( ( p2->getX() - p3->getX() ) * ( p2->getX() - p3->getX() ) + ( p2->getY() - p3->getY() ) * ( p2->getY() - p3->getY() ) );
double distp1p2 = sqrt( ( p1->x() - p2->x() ) * ( p1->x() - p2->x() ) + ( p1->y() - p2->y() ) * ( p1->y() - p2->y() ) );
double distp2p3 = sqrt( ( p2->x() - p3->x() ) * ( p2->x() - p3->x() ) + ( p2->y() - p3->y() ) * ( p2->y() - p3->y() ) );
if ( distp1p2 > distp2p3 )
{
//swap p1 and p3 to avoid round-off errors
Point3D *temp = p1;
QgsPoint *temp = p1;
p1 = p3;
p3 = temp;
}
double denominator = -p3->getX() * p2->getY() + p3->getX() * p1->getY() + p1->getX() * p2->getY() + p2->getX() * p3->getY() - p2->getX() * p1->getY() - p1->getX() * p3->getY();
double denominator = -p3->x() * p2->y() + p3->x() * p1->y() + p1->x() * p2->y() + p2->x() * p3->y() - p2->x() * p1->y() - p1->x() * p3->y();
//if one of the denominator is zero we will have problems
if ( denominator == 0 )
{
@ -138,8 +138,8 @@ bool MathUtils::circumcenter( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *re
}
else
{
result->setX( 0.5 * ( p1->getX()*p1->getX()*p2->getY() - p1->getX()*p1->getX()*p3->getY() - p3->getX()*p3->getX()*p2->getY() - p1->getY()*p2->getX()*p2->getX() - p1->getY()*p1->getY()*p3->getY() - p3->getY()*p3->getY()*p2->getY() + p1->getY()*p1->getY()*p2->getY() + p3->getY()*p2->getX()*p2->getX() - p1->getY()*p2->getY()*p2->getY() + p1->getY()*p3->getY()*p3->getY() + p1->getY()*p3->getX()*p3->getX() + p3->getY()*p2->getY()*p2->getY() ) / denominator );
result->setY( -0.5 * ( p3->getX()*p2->getX()*p2->getX() + p2->getX()*p1->getY()*p1->getY() + p3->getX()*p2->getY()*p2->getY() - p3->getX()*p1->getX()*p1->getX() + p1->getX()*p3->getY()*p3->getY() - p3->getX()*p1->getY()*p1->getY() - p1->getX()*p2->getX()*p2->getX() - p2->getX()*p3->getY()*p3->getY() - p1->getX()*p2->getY()*p2->getY() - p2->getX()*p3->getX()*p3->getX() + p1->getX()*p3->getX()*p3->getX() + p2->getX()*p1->getX()*p1->getX() ) / denominator );
result->setX( 0.5 * ( p1->x()*p1->x()*p2->y() - p1->x()*p1->x()*p3->y() - p3->x()*p3->x()*p2->y() - p1->y()*p2->x()*p2->x() - p1->y()*p1->y()*p3->y() - p3->y()*p3->y()*p2->y() + p1->y()*p1->y()*p2->y() + p3->y()*p2->x()*p2->x() - p1->y()*p2->y()*p2->y() + p1->y()*p3->y()*p3->y() + p1->y()*p3->x()*p3->x() + p3->y()*p2->y()*p2->y() ) / denominator );
result->setY( -0.5 * ( p3->x()*p2->x()*p2->x() + p2->x()*p1->y()*p1->y() + p3->x()*p2->y()*p2->y() - p3->x()*p1->x()*p1->x() + p1->x()*p3->y()*p3->y() - p3->x()*p1->y()*p1->y() - p1->x()*p2->x()*p2->x() - p2->x()*p3->y()*p3->y() - p1->x()*p2->y()*p2->y() - p2->x()*p3->x()*p3->x() + p1->x()*p3->x()*p3->x() + p2->x()*p1->x()*p1->x() ) / denominator );
#if 0
//debugging: test, if the distances from p1, p2, p3 to result are equal
@ -164,20 +164,20 @@ bool MathUtils::circumcenter( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *re
}
#if 0
bool MathUtils::circumcenter( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *result )//version imitating the geometric construction
bool MathUtils::circumcenter( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result )//version imitating the geometric construction
{
if ( p1 && p2 && p3 && result )
{
Point3D midpoint12( ( p1->getX() + p2->getX() ) / 2, ( p1->getY() + p2->getY() ) / 2, 0 );
Point3D midpoint23( ( p2->getX() + p3->getX() ) / 2, ( p2->getY() + p3->getY() ) / 2, 0 );
QgsPoint midpoint12( ( p1->getX() + p2->getX() ) / 2, ( p1->getY() + p2->getY() ) / 2, 0 );
QgsPoint midpoint23( ( p2->getX() + p3->getX() ) / 2, ( p2->getY() + p3->getY() ) / 2, 0 );
Vector3D v12( p2->getX() - p1->getX(), p2->getY() - p1->getY(), 0 );
Vector3D v23( p3->getX() - p2->getX(), p3->getY() - p2->getY(), 0 );
Vector3D n12;
MathUtils::normalLeft( &v12, &n12, 10000 );
Vector3D n23;
MathUtils::normalLeft( &v23, &n23, 10000 );
Point3D helppoint1( midpoint12.getX() + n12.getX(), midpoint12.getY() + n12.getY(), 0 );
Point3D helppoint2( midpoint23.getX() + n23.getX(), midpoint23.getY() + n23.getY(), 0 );
QgsPoint helppoint1( midpoint12.getX() + n12.getX(), midpoint12.getY() + n12.getY(), 0 );
QgsPoint helppoint2( midpoint23.getX() + n23.getX(), midpoint23.getY() + n23.getY(), 0 );
MathUtils::lineIntersection( &midpoint12, &helppoint1, &midpoint23, &helppoint2, result );
//debugging: test, if the distances from p1, p2, p3 to result are equal
@ -201,17 +201,17 @@ bool MathUtils::circumcenter( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *re
}
#endif // 0
double MathUtils::distPointFromLine( Point3D *thepoint, Point3D *p1, Point3D *p2 )
double MathUtils::distPointFromLine( QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2 )
{
if ( thepoint && p1 && p2 )
{
Vector3D normal( 0, 0, 0 );
Vector3D line( p2->getX() - p1->getX(), p2->getY() - p1->getY(), 0 );
Vector3D line( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
MathUtils::normalLeft( &line, &normal, 1 );
double a = normal.getX();
double b = normal.getY();
double c = -( normal.getX() * p2->getX() + normal.getY() * p2->getY() );
double distance = qAbs( ( a * thepoint->getX() + b * thepoint->getY() + c ) / ( sqrt( a * a + b * b ) ) );
double c = -( normal.getX() * p2->x() + normal.getY() * p2->y() );
double distance = qAbs( ( a * thepoint->x() + b * thepoint->y() + c ) / ( sqrt( a * a + b * b ) ) );
return distance;
}
else
@ -241,20 +241,20 @@ int MathUtils::faculty( int n )
return result;
}
bool MathUtils::inCircle( Point3D *testp, Point3D *p1, Point3D *p2, Point3D *p3 )
bool MathUtils::inCircle( QgsPoint *testp, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 )
{
double tolerance = 0.0001;//if the amount of aValue is below this, testp is approximately on the circle and we have to define another criterion to tell, if it is inside or outside.
if ( testp && p1 && p2 && p3 )
{
double ax = p1->getX();
double ay = p1->getY();
double bx = p2->getX();
double by = p2->getY();
double cx = p3->getX();
double cy = p3->getY();
double px = testp->getX();
double py = testp->getY();
double ax = p1->x();
double ay = p1->y();
double bx = p2->x();
double by = p2->y();
double cx = p3->x();
double cy = p3->y();
double px = testp->x();
double py = testp->y();
double xmin = qMin( qMin( ax, px ), qMin( bx, cx ) );
double ymin = qMin( qMin( ay, py ), qMin( by, cy ) );
@ -281,22 +281,22 @@ bool MathUtils::inCircle( Point3D *testp, Point3D *p1, Point3D *p2, Point3D *p3
}
}
bool MathUtils::inDiametral( Point3D *p1, Point3D *p2, Point3D *point )
bool MathUtils::inDiametral( QgsPoint *p1, QgsPoint *p2, QgsPoint *point )
{
return angle( p1, point, p2, point ) > 90;
}
#if 0
bool MathUtils::inDiametral( Point3D *p1, Point3D *p2, Point3D *point )
bool MathUtils::inDiametral( QgsPoint *p1, QgsPoint *p2, QgsPoint *point )
{
if ( p1 && p2 && point )
{
Vector3D p1p2( p2->getX() - p1->getX(), p2->getY() - p1->getY(), 0 );
Vector3D orthogonalvec;//vector orthogonal to p1p2 (length radius)
Point3D midpoint( ( p1->getX() + p2->getX() ) / 2, ( p1->getY() + p2->getY() ) / 2, 0 );
QgsPoint midpoint( ( p1->getX() + p2->getX() ) / 2, ( p1->getY() + p2->getY() ) / 2, 0 );
double radius = p1p2.getLength() / 2;
normalLeft( &p1p2, &orthogonalvec, radius );
Point3D p3( midpoint.getX() + orthogonalvec.getX(), midpoint.getY() + orthogonalvec.getY(), 0 );
QgsPoint p3( midpoint.getX() + orthogonalvec.getX(), midpoint.getY() + orthogonalvec.getY(), 0 );
return inCircle( point, p1, p2, &p3 );
}
else
@ -307,15 +307,15 @@ bool MathUtils::inDiametral( Point3D *p1, Point3D *p2, Point3D *point )
}
#endif // 0
double MathUtils::leftOf( Point3D *thepoint, Point3D *p1, Point3D *p2 )
double MathUtils::leftOf( QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2 )
{
if ( thepoint && p1 && p2 )
{
//just for debugging
double f1 = thepoint->getX() - p1->getX();
double f2 = p2->getY() - p1->getY();
double f3 = thepoint->getY() - p1->getY();
double f4 = p2->getX() - p1->getX();
double f1 = thepoint->x() - p1->x();
double f2 = p2->y() - p1->y();
double f3 = thepoint->y() - p1->y();
double f4 = p2->x() - p1->x();
return f1 * f2 - f3 * f4;
//return thepoint->getX()-p1->getX())*(p2->getY()-p1->getY())-(thepoint->getY()-p1->getY())*(p2->getX()-p1->getX();//calculating the vectorproduct
}
@ -326,22 +326,22 @@ double MathUtils::leftOf( Point3D *thepoint, Point3D *p1, Point3D *p2 )
}
}
bool MathUtils::lineIntersection( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *p4 )
bool MathUtils::lineIntersection( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4 )
{
if ( p1 && p2 && p3 && p4 )
{
double t1, t2;
Vector3D p1p2( p2->getX() - p1->getX(), p2->getY() - p1->getY(), 0 );
Vector3D p3p4( p4->getX() - p3->getX(), p4->getY() - p3->getY(), 0 );
Vector3D p1p2( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
Vector3D p3p4( p4->x() - p3->x(), p4->y() - p3->y(), 0 );
if ( ( p3p4.getX()*p1p2.getY() - p3p4.getY()*p1p2.getX() ) != 0 && p1p2.getX() != 0 ) //avoid division through zero
{
t2 = ( p1->getX() * p1p2.getY() - p1->getY() * p1p2.getX() + p3->getY() * p1p2.getX() - p3->getX() * p1p2.getY() ) / ( p3p4.getX() * p1p2.getY() - p3p4.getY() * p1p2.getX() );
t1 = ( p3->getX() - p1->getX() + t2 * p3p4.getX() ) / p1p2.getX();
t2 = ( p1->x() * p1p2.getY() - p1->y() * p1p2.getX() + p3->y() * p1p2.getX() - p3->x() * p1p2.getY() ) / ( p3p4.getX() * p1p2.getY() - p3p4.getY() * p1p2.getX() );
t1 = ( p3->x() - p1->x() + t2 * p3p4.getX() ) / p1p2.getX();
}
else if ( ( p1p2.getX()*p3p4.getY() - p1p2.getY()*p3p4.getX() ) != 0 && p3p4.getX() != 0 )
{
t1 = ( p3->getX() * p3p4.getY() - p3->getY() * p3p4.getX() - p1->getX() * p3p4.getY() + p1->getY() * p3p4.getX() ) / ( p1p2.getX() * p3p4.getY() - p1p2.getY() * p3p4.getX() );
t2 = ( p1->getX() + t1 * p1p2.getX() - p3->getX() ) / p3p4.getX();
t1 = ( p3->x() * p3p4.getY() - p3->y() * p3p4.getX() - p1->x() * p3p4.getY() + p1->y() * p3p4.getX() ) / ( p1p2.getX() * p3p4.getY() - p1p2.getY() * p3p4.getX() );
t2 = ( p1->x() + t1 * p1p2.getX() - p3->x() ) / p3p4.getX();
}
else//the lines are parallel
{
@ -369,22 +369,22 @@ bool MathUtils::lineIntersection( Point3D *p1, Point3D *p2, Point3D *p3, Point3D
}
}
bool MathUtils::lineIntersection( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *p4, Point3D *intersection_point )
bool MathUtils::lineIntersection( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4, QgsPoint *intersection_point )
{
if ( p1 && p2 && p3 && p4 )
{
double t1, t2;
Vector3D p1p2( p2->getX() - p1->getX(), p2->getY() - p1->getY(), 0 );
Vector3D p3p4( p4->getX() - p3->getX(), p4->getY() - p3->getY(), 0 );
Vector3D p1p2( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
Vector3D p3p4( p4->x() - p3->x(), p4->y() - p3->y(), 0 );
if ( ( p3p4.getX()*p1p2.getY() - p3p4.getY()*p1p2.getX() ) != 0 && p1p2.getX() != 0 ) //avoid division through zero
{
t2 = ( p1->getX() * p1p2.getY() - p1->getY() * p1p2.getX() + p3->getY() * p1p2.getX() - p3->getX() * p1p2.getY() ) / ( p3p4.getX() * p1p2.getY() - p3p4.getY() * p1p2.getX() );
t1 = ( p3->getX() - p1->getX() + t2 * p3p4.getX() ) / p1p2.getX();
t2 = ( p1->x() * p1p2.getY() - p1->y() * p1p2.getX() + p3->y() * p1p2.getX() - p3->x() * p1p2.getY() ) / ( p3p4.getX() * p1p2.getY() - p3p4.getY() * p1p2.getX() );
t1 = ( p3->x() - p1->x() + t2 * p3p4.getX() ) / p1p2.getX();
}
else if ( ( p1p2.getX()*p3p4.getY() - p1p2.getY()*p3p4.getX() ) != 0 && p3p4.getX() != 0 )
{
t1 = ( p3->getX() * p3p4.getY() - p3->getY() * p3p4.getX() - p1->getX() * p3p4.getY() + p1->getY() * p3p4.getX() ) / ( p1p2.getX() * p3p4.getY() - p1p2.getY() * p3p4.getX() );
t2 = ( p1->getX() + t1 * p1p2.getX() - p3->getX() ) / p3p4.getX();
t1 = ( p3->x() * p3p4.getY() - p3->y() * p3p4.getX() - p1->x() * p3p4.getY() + p1->y() * p3p4.getX() ) / ( p1p2.getX() * p3p4.getY() - p1p2.getY() * p3p4.getX() );
t2 = ( p1->x() + t1 * p1p2.getX() - p3->x() ) / p3p4.getX();
}
else//the lines are parallel
{
@ -404,8 +404,8 @@ bool MathUtils::lineIntersection( Point3D *p1, Point3D *p2, Point3D *p3, Point3D
return false;
}
//calculate the intersection point
intersection_point->setX( p1->getX() * ( 1 - t1 ) + p2->getX()*t1 );
intersection_point->setY( p1->getY() * ( 1 - t1 ) + p2->getY()*t1 );
intersection_point->setX( p1->x() * ( 1 - t1 ) + p2->x()*t1 );
intersection_point->setY( p1->y() * ( 1 - t1 ) + p2->y()*t1 );
intersection_point->setZ( 0 );
return true;
}
@ -487,11 +487,11 @@ double MathUtils::power( double a, int b )
}
}
double MathUtils::triArea( Point3D *pa, Point3D *pb, Point3D *pc )
double MathUtils::triArea( QgsPoint *pa, QgsPoint *pb, QgsPoint *pc )
{
if ( pa && pb && pc )
{
double deter = ( pa->getX() * pb->getY() + pb->getX() * pc->getY() + pc->getX() * pa->getY() - pa->getX() * pc->getY() - pb->getX() * pa->getY() - pc->getX() * pb->getY() );
double deter = ( pa->x() * pb->y() + pb->x() * pc->y() + pc->x() * pa->y() - pa->x() * pc->y() - pb->x() * pa->y() - pc->x() * pb->y() );
return 0.5 * deter;
}
else//null pointer
@ -623,9 +623,9 @@ bool MathUtils::normalLeft( Vector3D *v1, Vector3D *result, double length )
result->setX( ( -b + sqrt( d ) ) / ( 2 * a ) ); //take one of the two solutions of the quadratic equation
result->setY( ( -result->getX()*v1->getX() ) / v1->getY() );
Point3D point1( 0, 0, 0 );
Point3D point2( v1->getX(), v1->getY(), 0 );
Point3D point3( result->getX(), result->getY(), 0 );
QgsPoint point1( 0, 0, 0 );
QgsPoint point2( v1->getX(), v1->getY(), 0 );
QgsPoint point3( result->getX(), result->getY(), 0 );
if ( !( leftOf( &point1, &point2, &point3 ) < 0 ) )//if we took the solution on the right side, change the sign of the components of the result
{
@ -680,9 +680,9 @@ bool MathUtils::normalRight( Vector3D *v1, Vector3D *result, double length )
result->setX( ( -b + sqrt( d ) ) / ( 2 * a ) ); //take one of the two solutions of the quadratic equation
result->setY( ( -result->getX()*v1->getX() ) / v1->getY() );
Point3D point1( 0, 0, 0 );
Point3D point2( v1->getX(), v1->getY(), 0 );
Point3D point3( result->getX(), result->getY(), 0 );
QgsPoint point1( 0, 0, 0 );
QgsPoint point2( v1->getX(), v1->getY(), 0 );
QgsPoint point3( result->getX(), result->getY(), 0 );
if ( ( leftOf( &point1, &point2, &point3 ) < 0 ) ) //if we took the solution on the right side, change the sign of the components of the result
{
@ -700,16 +700,16 @@ bool MathUtils::normalRight( Vector3D *v1, Vector3D *result, double length )
}
void MathUtils::normalFromPoints( Point3D *p1, Point3D *p2, Point3D *p3, Vector3D *vec )
void MathUtils::normalFromPoints( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, Vector3D *vec )
{
if ( p1 && p2 && p3 && vec )//no null pointers
{
double ax = p2->getX() - p1->getX();
double ay = p2->getY() - p1->getY();
double az = p2->getZ() - p1->getZ();
double bx = p3->getX() - p1->getX();
double by = p3->getY() - p1->getY();
double bz = p3->getZ() - p1->getZ();
double ax = p2->x() - p1->x();
double ay = p2->y() - p1->y();
double az = p2->z() - p1->z();
double bx = p3->x() - p1->x();
double by = p3->y() - p1->y();
double bz = p3->z() - p1->z();
vec->setX( ay * bz - az * by );
vec->setY( az * bx - ax * bz );
@ -718,7 +718,7 @@ void MathUtils::normalFromPoints( Point3D *p1, Point3D *p2, Point3D *p3, Vector3
}
double MathUtils::crossVec( Point3D *first, Vector3D *vec1, Point3D *second, Vector3D *vec2 )
double MathUtils::crossVec( QgsPoint *first, Vector3D *vec1, QgsPoint *second, Vector3D *vec2 )
{
if ( first && vec1 && second && vec2 )
{
@ -730,7 +730,7 @@ double MathUtils::crossVec( Point3D *first, Vector3D *vec1, Point3D *second, Vec
cout << "vec2: " << vec2->getX() << "//" << vec2->getY() << "//" << vec2->getZ() << endl << flush;
cout << "t2: " << ((first->getX()*vec1->getY()-first->getY()*vec1->getX()-second->getX()*vec1->getY()+second->getY()*vec1->getX())/(vec2->getX()*vec1->getY()-vec2->getY()*vec1->getX())) << endl << flush;*/
return ( ( first->getX() * vec1->getY() - first->getY() * vec1->getX() - second->getX() * vec1->getY() + second->getY() * vec1->getX() ) / ( vec2->getX() * vec1->getY() - vec2->getY() * vec1->getX() ) );
return ( ( first->x() * vec1->getY() - first->y() * vec1->getX() - second->x() * vec1->getY() + second->y() * vec1->getX() ) / ( vec2->getX() * vec1->getY() - vec2->getY() * vec1->getX() ) );
}
else//if a division by zero would occur
@ -748,9 +748,9 @@ double MathUtils::crossVec( Point3D *first, Vector3D *vec1, Point3D *second, Vec
}
}
bool MathUtils::pointInsideTriangle( double x, double y, Point3D *p1, Point3D *p2, Point3D *p3 )
bool MathUtils::pointInsideTriangle( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 )
{
Point3D thepoint( x, y, 0 );
QgsPoint thepoint( x, y, 0 );
if ( MathUtils::leftOf( &thepoint, p1, p2 ) > 0 )
{
return false;
@ -837,15 +837,15 @@ bool MathUtils::normalMinDistance( Vector3D *tangent, Vector3D *target, Vector3D
}
double MathUtils::planeTest( Point3D *test, Point3D *pt1, Point3D *pt2, Point3D *pt3 )
double MathUtils::planeTest( QgsPoint *test, QgsPoint *pt1, QgsPoint *pt2, QgsPoint *pt3 )
{
if ( test && pt1 && pt2 && pt3 )
{
double a = ( pt1->getZ() * ( pt2->getY() - pt3->getY() ) + pt2->getZ() * ( pt3->getY() - pt1->getY() ) + pt3->getZ() * ( pt1->getY() - pt2->getY() ) ) / ( ( pt1->getX() - pt2->getX() ) * ( pt2->getY() - pt3->getY() ) - ( pt2->getX() - pt3->getX() ) * ( pt1->getY() - pt2->getY() ) );
double b = ( pt1->getZ() * ( pt2->getX() - pt3->getX() ) + pt2->getZ() * ( pt3->getX() - pt1->getX() ) + pt3->getZ() * ( pt1->getX() - pt2->getX() ) ) / ( ( pt1->getY() - pt2->getY() ) * ( pt2->getX() - pt3->getX() ) - ( pt2->getY() - pt3->getY() ) * ( pt1->getX() - pt2->getX() ) );
double c = pt1->getZ() - a * pt1->getX() - b * pt1->getY();
double zpredicted = test->getX() * a + test->getY() * b + c;
return ( test->getZ() - zpredicted );
double a = ( pt1->z() * ( pt2->y() - pt3->y() ) + pt2->z() * ( pt3->y() - pt1->y() ) + pt3->z() * ( pt1->y() - pt2->y() ) ) / ( ( pt1->x() - pt2->x() ) * ( pt2->y() - pt3->y() ) - ( pt2->x() - pt3->x() ) * ( pt1->y() - pt2->y() ) );
double b = ( pt1->z() * ( pt2->x() - pt3->x() ) + pt2->z() * ( pt3->x() - pt1->x() ) + pt3->z() * ( pt1->x() - pt2->x() ) ) / ( ( pt1->y() - pt2->y() ) * ( pt2->x() - pt3->x() ) - ( pt2->y() - pt3->y() ) * ( pt1->x() - pt2->x() ) );
double c = pt1->z() - a * pt1->x() - b * pt1->y();
double zpredicted = test->x() * a + test->y() * b + c;
return ( test->z() - zpredicted );
}
else
{
@ -854,12 +854,12 @@ double MathUtils::planeTest( Point3D *test, Point3D *pt1, Point3D *pt2, Point3D
}
}
double MathUtils::angle( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *p4 )
double MathUtils::angle( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4 )
{
if ( p1 && p2 && p3 && p4 )
{
Vector3D v1( p2->getX() - p1->getX(), p2->getY() - p1->getY(), 0 );
Vector3D v2( p4->getX() - p3->getX(), p4->getY() - p3->getY(), 0 );
Vector3D v1( p2->x() - p1->x(), p2->y() - p1->y(), 0 );
Vector3D v2( p4->x() - p3->x(), p4->y() - p3->y(), 0 );
double value = acos( ( v1.getX() * v2.getX() + v1.getY() * v2.getY() ) / ( v1.getLength() * v2.getLength() ) ) * 180 / M_PI;
return value;
}

32
src/analysis/interpolation/MathUtils.h
View File

@ -20,14 +20,14 @@
#include <cmath>
#include "qgis_analysis.h"
class Point3D;
class QgsPoint;
class Vector3D;
namespace MathUtils
{
//! Calculates the barycentric coordinates of a point (x,y) with respect to p1, p2, p3 and stores the three barycentric coordinates in 'result'. Thus the u-coordinate is stored in result::x, the v-coordinate in result::y and the w-coordinate in result::z. Attention: p1, p2 and p3 have to be ordered counterclockwise
bool ANALYSIS_EXPORT calcBarycentricCoordinates( double x, double y, Point3D *p1, Point3D *p2, Point3D *p3, Point3D *result );
bool ANALYSIS_EXPORT BarycentricToXY( double u, double v, double w, Point3D *p1, Point3D *p2, Point3D *p3, Point3D *result );
bool ANALYSIS_EXPORT calcBarycentricCoordinates( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result );
bool ANALYSIS_EXPORT BarycentricToXY( double u, double v, double w, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result );
//! Calculates the value of a Bernstein polynomial
double ANALYSIS_EXPORT calcBernsteinPoly( int n, int i, double t );
//! Calculates the first derivative of a Bernstein polynomial with respect to the parameter t
@ -37,21 +37,21 @@ namespace MathUtils
//! Calculates the first derivative of a cubic Hermite polynomial with respect to the parameter t
double ANALYSIS_EXPORT cFDerCubicHermitePoly( int n, int i, double t );
//! Calculates the center of the circle passing through p1, p2 and p3. Returns true in case of success and false otherwise (e.g. all three points on a line)
bool ANALYSIS_EXPORT circumcenter( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *result );
bool ANALYSIS_EXPORT circumcenter( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *result );
//! Calculates the (2 dimensional) distance from 'thepoint' to the line defined by p1 and p2
double ANALYSIS_EXPORT distPointFromLine( Point3D *thepoint, Point3D *p1, Point3D *p2 );
double ANALYSIS_EXPORT distPointFromLine( QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2 );
//! Faculty function
int ANALYSIS_EXPORT faculty( int n );
//! Tests, whether 'testp' is inside the circle through 'p1', 'p2' and 'p3'
bool ANALYSIS_EXPORT inCircle( Point3D *testp, Point3D *p1, Point3D *p2, Point3D *p3 );
bool ANALYSIS_EXPORT inCircle( QgsPoint *testp, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 );
//! Tests, whether 'point' is inside the diametral circle through 'p1' and 'p2'
bool ANALYSIS_EXPORT inDiametral( Point3D *p1, Point3D *p2, Point3D *point );
bool ANALYSIS_EXPORT inDiametral( QgsPoint *p1, QgsPoint *p2, QgsPoint *point );
//! Returns whether 'thepoint' is left or right of the line from 'p1' to 'p2'. Negativ values mean left and positiv values right. There may be numerical instabilities, so a threshold may be useful
double ANALYSIS_EXPORT leftOf( Point3D *thepoint, Point3D *p1, Point3D *p2 );
double ANALYSIS_EXPORT leftOf( QgsPoint *thepoint, QgsPoint *p1, QgsPoint *p2 );
//! Returns true, if line1 (p1 to p2) and line2 (p3 to p4) intersect. If the lines have an endpoint in common, false is returned
bool ANALYSIS_EXPORT lineIntersection( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *p4 );
bool ANALYSIS_EXPORT lineIntersection( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4 );
//! Returns true, if line1 (p1 to p2) and line2 (p3 to p4) intersect. If the lines have an endpoint in common, false is returned. The intersecting point is stored in 'intersection_point.
bool ANALYSIS_EXPORT lineIntersection( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *p4, Point3D *intersection_point );
bool ANALYSIS_EXPORT lineIntersection( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4, QgsPoint *intersection_point );
//! Lower function
int ANALYSIS_EXPORT lower( int n, int i );
//! Returns the maximum of two doubles or the first argument if both are equal
@ -61,25 +61,25 @@ namespace MathUtils
//! Power function for integer coefficients
double ANALYSIS_EXPORT power( double a, int b );//calculates a power b
//! Returns the area of a triangle. If the points are ordered counterclockwise, the value will be positiv. If they are ordered clockwise, the value will be negativ
double ANALYSIS_EXPORT triArea( Point3D *pa, Point3D *pb, Point3D *pc );
double ANALYSIS_EXPORT triArea( QgsPoint *pa, QgsPoint *pb, QgsPoint *pc );
//! Calculates the z-component of a vector with coordinates 'x' and 'y'which is in the same tangent plane as the tangent vectors 'v1' and 'v2'. The result is assigned to 'result'
bool ANALYSIS_EXPORT derVec( const Vector3D *v1, const Vector3D *v2, Vector3D *result, double x, double y );
//! Calculates the intersection of the two vectors vec1 and vec2, which start at first(vec1) and second(vec2) end. The return value is t2(multiplication of v2 with t2 and adding to 'second' results the intersection point)
double ANALYSIS_EXPORT crossVec( Point3D *first, Vector3D *vec1, Point3D *second, Vector3D *vec2 );
double ANALYSIS_EXPORT crossVec( QgsPoint *first, Vector3D *vec1, QgsPoint *second, Vector3D *vec2 );
//! Assigns the vector 'result', which is normal to the vector 'v1', on the left side of v1 and has length 'length'. This method works only with two dimensions.
bool ANALYSIS_EXPORT normalLeft( Vector3D *v1, Vector3D *result, double length );
//! Assigns the vector 'result', which is normal to the vector 'v1', on the right side of v1 and has length 'length'. The calculation is only in two dimensions
bool ANALYSIS_EXPORT normalRight( Vector3D *v1, Vector3D *result, double length );
//! Calculates the normal vector of the plane through the points p1, p2 and p3 and assigns the result to vec. If the points are ordered counterclockwise, the normal will have a positive z-coordinate;
void ANALYSIS_EXPORT normalFromPoints( Point3D *p1, Point3D *p2, Point3D *p3, Vector3D *vec );
void ANALYSIS_EXPORT normalFromPoints( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, Vector3D *vec );
//! Returns true, if the point with coordinates x and y is inside (or at the edge) of the triangle p1,p2,p3 and false, if it is outside. p1, p2 and p3 have to be ordered counterclockwise
bool ANALYSIS_EXPORT pointInsideTriangle( double x, double y, Point3D *p1, Point3D *p2, Point3D *p3 );
bool ANALYSIS_EXPORT pointInsideTriangle( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 );
//! Calculates a Vector orthogonal to 'tangent' with length 1 and closest possible to result. Returns true in case of success and false otherwise
bool ANALYSIS_EXPORT normalMinDistance( Vector3D *tangent, Vector3D *target, Vector3D *result );
//! Tests, if 'test' is in the same plane as 'p1', 'p2' and 'p3' and returns the z-difference from the plane to 'test
double ANALYSIS_EXPORT planeTest( Point3D *test, Point3D *pt1, Point3D *pt2, Point3D *pt3 );
double ANALYSIS_EXPORT planeTest( QgsPoint *test, QgsPoint *pt1, QgsPoint *pt2, QgsPoint *pt3 );
//! Calculates the angle between two segments (in 2 dimension, z-values are ignored)
double ANALYSIS_EXPORT angle( Point3D *p1, Point3D *p2, Point3D *p3, Point3D *p4 );
double ANALYSIS_EXPORT angle( QgsPoint *p1, QgsPoint *p2, QgsPoint *p3, QgsPoint *p4 );
}
#endif

12
src/analysis/interpolation/Node.cc
View File

@ -20,7 +20,7 @@ Node::Node( const Node &n )
{
if ( n.getPoint() )
{
Point3D *point = new Point3D( n.getPoint()->getX(), n.getPoint()->getY(), n.getPoint()->getZ() );
QgsPoint *point = new QgsPoint( n.getPoint()->x(), n.getPoint()->y(), n.getPoint()->z() );
mPoint = point;
}
else
@ -34,11 +34,11 @@ Node::Node( const Node &n )
Node &Node::operator=( const Node &n )
{
Point3D *tmp = mPoint;
QgsPoint *tmp = mPoint;
if ( n.getPoint() )//mPoint of n is not a null pointer
{
mPoint = new Point3D( n.getPoint()->getX(), n.getPoint()->getY(), n.getPoint()->getZ() );
mPoint = new QgsPoint( n.getPoint()->x(), n.getPoint()->y(), n.getPoint()->z() );
if ( !mPoint )//no memory
{
mPoint = tmp;
@ -52,10 +52,8 @@ Node &Node::operator=( const Node &n )
mPoint = nullptr;
}
if ( tmp )
{
delete tmp;
}
delete tmp;
mNext = n.getNext();
return ( *this );

20
src/analysis/interpolation/Node.h
View File

@ -17,16 +17,16 @@
#ifndef NODE_H
#define NODE_H
#include "Point3D.h"
#include "qgspoint.h"
#include "qgis_analysis.h"
/** \ingroup analysis
* Node is a class used by Line3D. It represents a node in the single directed linked list. Associated Point3D objects are deleted when the node is deleted.*/
* Node is a class used by Line3D. It represents a node in the single directed linked list. Associated QgsPoint objects are deleted when the node is deleted.*/
class ANALYSIS_EXPORT Node
{
protected:
//! Pointer to the Point3D object associated with the node
Point3D *mPoint = nullptr;
//! Pointer to the QgsPoint object associated with the node
QgsPoint *mPoint = nullptr;
//! Pointer to the next Node in the linked list
Node *mNext = nullptr;
public:
@ -36,12 +36,12 @@ class ANALYSIS_EXPORT Node
Node &operator=( const Node &n );
//! Returns a pointer to the next element in the linked list
Node *getNext() const;
//! Returns a pointer to the Point3D object associated with the node
Point3D *getPoint() const;
//! Returns a pointer to the QgsPoint object associated with the node
QgsPoint *getPoint() const;
//! Sets the pointer to the next node
void setNext( Node *n );
//! Sets a new pointer to an associated Point3D object
void setPoint( Point3D *p );
//! Sets a new pointer to an associated QgsPoint object
void setPoint( QgsPoint *p );
};
inline Node::Node()
@ -61,7 +61,7 @@ inline Node *Node::getNext() const
return mNext;
}
inline Point3D *Node::getPoint() const
inline QgsPoint *Node::getPoint() const
{
return mPoint;
}
@ -71,7 +71,7 @@ inline void Node::setNext( Node *n )
mNext = n;
}
inline void Node::setPoint( Point3D *p )
inline void Node::setPoint( QgsPoint *p )
{
mPoint = p;
}

14
src/analysis/interpolation/NormVecDecorator.cc
View File

@ -39,7 +39,7 @@ NormVecDecorator::~NormVecDecorator()
}
}
int NormVecDecorator::addPoint( Point3D *p )
int NormVecDecorator::addPoint( QgsPoint *p )
{
if ( mTIN )
{
@ -231,7 +231,7 @@ bool NormVecDecorator::calcNormalForPoint( double x, double y, int point, Vector
}
bool NormVecDecorator::calcPoint( double x, double y, Point3D *result )
bool NormVecDecorator::calcPoint( double x, double y, QgsPoint *result )
{
if ( !alreadyestimated )
@ -252,7 +252,7 @@ bool NormVecDecorator::calcPoint( double x, double y, Point3D *result )
}
}
bool NormVecDecorator::getTriangle( double x, double y, Point3D *p1, Vector3D *v1, Point3D *p2, Vector3D *v2, Point3D *p3, Vector3D *v3 )
bool NormVecDecorator::getTriangle( double x, double y, QgsPoint *p1, Vector3D *v1, QgsPoint *p2, Vector3D *v2, QgsPoint *p3, Vector3D *v3 )
{
if ( p1 && p2 && p3 && v1 && v2 && v3 )
{
@ -322,7 +322,7 @@ NormVecDecorator::PointState NormVecDecorator::getState( int pointno ) const
}
bool NormVecDecorator::getTriangle( double x, double y, Point3D *p1, int *ptn1, Vector3D *v1, PointState *state1, Point3D *p2, int *ptn2, Vector3D *v2, PointState *state2, Point3D *p3, int *ptn3, Vector3D *v3, PointState *state3 )
bool NormVecDecorator::getTriangle( double x, double y, QgsPoint *p1, int *ptn1, Vector3D *v1, PointState *state1, QgsPoint *p2, int *ptn2, Vector3D *v2, PointState *state2, QgsPoint *p3, int *ptn3, Vector3D *v3, PointState *state3 )
{
if ( p1 && p2 && p3 && v1 && v2 && v3 && ptn1 && ptn2 && ptn3 && state1 && state2 && state3 )
{
@ -438,10 +438,10 @@ bool NormVecDecorator::estimateFirstDerivative( int pointno )
if ( p1 != -1 && p2 != -1 && p3 != -1 )//don't calculate normal, if a point is a virtual point
{
MathUtils::normalFromPoints( getPoint( p1 ), getPoint( p2 ), getPoint( p3 ), &part );
double dist1 = getPoint( p3 )->dist3D( getPoint( p1 ) );
double dist2 = getPoint( p3 )->dist3D( getPoint( p2 ) );
double dist1 = getPoint( p3 )->distance3D( *getPoint( p1 ) );
double dist2 = getPoint( p3 )->distance3D( *getPoint( p2 ) );
//don't add the normal if the triangle is horizontal
if ( ( getPoint( p1 )->getZ() != getPoint( p2 )->getZ() ) || ( getPoint( p1 )->getZ() != getPoint( p3 )->getZ() ) )
if ( ( getPoint( p1 )->z() != getPoint( p2 )->z() ) || ( getPoint( p1 )->z() != getPoint( p3 )->z() ) )
{
currentweight = 1 / ( dist1 * dist1 * dist2 * dist2 );
total.setX( total.getX() + part.getX()*currentweight );

8
src/analysis/interpolation/NormVecDecorator.h
View File

@ -38,13 +38,13 @@ class ANALYSIS_EXPORT NormVecDecorator: public TriDecorator
NormVecDecorator( Triangulation *tin );
virtual ~NormVecDecorator();
//! Adds a point to the triangulation
int addPoint( Point3D *p ) override;
int addPoint( QgsPoint *p ) override;
//! Calculates the normal at a point on the surface and assigns it to 'result'. Returns true in case of success and false in case of failure
bool calcNormal( double x, double y, Vector3D *result ) override;
//! Calculates the normal of a triangle-point for the point with coordinates x and y. This is needed, if a point is on a break line and there is no unique normal stored in 'mNormVec'. Returns false, it something went wrong and true otherwise
bool calcNormalForPoint( double x, double y, int point, Vector3D *result );
//! Calculates x-, y and z-value of the point on the surface and assigns it to 'result'. Returns true in case of success and flase in case of failure
bool calcPoint( double x, double y, Point3D *result ) override;
bool calcPoint( double x, double y, QgsPoint *result ) override;
//! Eliminates the horizontal triangles by swapping or by insertion of new points. If alreadyestimated is true, a re-estimation of the normals will be done
virtual void eliminateHorizontalTriangles() override;
//! Estimates the first derivative a point. Return true in case of succes and false otherwise
@ -54,12 +54,12 @@ class ANALYSIS_EXPORT NormVecDecorator: public TriDecorator
//! Returns a pointer to the normal vector for the point with the number n
Vector3D *getNormal( int n ) const;
//! Finds out, in which triangle a point with coordinates x and y is and assigns the triangle points to p1, p2, p3 and the estimated normals to v1, v2, v3. The vectors are normally taken from 'mNormVec', except if p1, p2 or p3 is a point on a breakline. In this case, the normal is calculated on-the-fly. Returns false, if something went wrong and true otherwise
bool getTriangle( double x, double y, Point3D *p1 SIP_OUT, Vector3D *v1 SIP_OUT, Point3D *p2 SIP_OUT, Vector3D *v2 SIP_OUT, Point3D *p3 SIP_OUT, Vector3D *v3 SIP_OUT );
bool getTriangle( double x, double y, QgsPoint *p1 SIP_OUT, Vector3D *v1 SIP_OUT, QgsPoint *p2 SIP_OUT, Vector3D *v2 SIP_OUT, QgsPoint *p3 SIP_OUT, Vector3D *v3 SIP_OUT );
/** This function behaves similar to the one above. Additionally, the numbers of the points are returned (ptn1, ptn2, ptn3) as well as the PointStates of the triangle points (state1, state2, state3)
* \note not available in Python bindings
*/
bool getTriangle( double x, double y, Point3D *p1, int *ptn1, Vector3D *v1, PointState *state1, Point3D *p2, int *ptn2, Vector3D *v2, PointState *state2, Point3D *p3, int *ptn3, Vector3D *v3, PointState *state3 ) SIP_SKIP;
bool getTriangle( double x, double y, QgsPoint *p1, int *ptn1, Vector3D *v1, PointState *state1, QgsPoint *p2, int *ptn2, Vector3D *v2, PointState *state2, QgsPoint *p3, int *ptn3, Vector3D *v3, PointState *state3 ) SIP_SKIP;
//! Returns the state of the point with the number 'pointno'
PointState getState( int pointno ) const;
//! Sets an interpolator

8
src/analysis/interpolation/ParametricLine.cc
View File

@ -37,7 +37,7 @@ void ParametricLine::calcSecDer( float t, Vector3D *v )
QgsDebugMsg( "warning, derive a class from ParametricLine" );
}
void ParametricLine::calcPoint( float t, Point3D *p )
void ParametricLine::calcPoint( float t, QgsPoint *p )
{
Q_UNUSED( t );
Q_UNUSED( p );
@ -56,7 +56,7 @@ void ParametricLine::remove( int i )
QgsDebugMsg( "warning, derive a class from ParametricLine" );
}
void ParametricLine::setControlPoly( QVector<Point3D *> *cp )
void ParametricLine::setControlPoly( QVector<QgsPoint *> *cp )
{
Q_UNUSED( cp );
QgsDebugMsg( "warning, derive a class from ParametricLine" );
@ -74,14 +74,14 @@ int ParametricLine::getDegree() const
return mDegree;
}
const Point3D *ParametricLine::getControlPoint( int number ) const
const QgsPoint *ParametricLine::getControlPoint( int number ) const
{
Q_UNUSED( number );
QgsDebugMsg( "warning, derive a class from ParametricLine" );
return nullptr;
}
const QVector<Point3D *> *ParametricLine::getControlPoly() const
const QVector<QgsPoint *> *ParametricLine::getControlPoly() const
{
QgsDebugMsg( "warning, derive a class from ParametricLine" );
return nullptr;

18
src/analysis/interpolation/ParametricLine.h
View File

@ -17,7 +17,7 @@
#ifndef PARAMETRICLINE_H
#define PARAMETRICLINE_H
#include "Point3D.h"
#include "qgspoint.h"
#include <QVector>
#include "qgis_analysis.h"
@ -34,29 +34,29 @@ class ANALYSIS_EXPORT ParametricLine
//! Pointer to the parent object. If there isn't one, mParent is 0
ParametricLine *mParent = nullptr;
//! MControlPoly stores the points of the control polygon
QVector<Point3D *> *mControlPoly;
QVector<QgsPoint *> *mControlPoly;
public:
//! Default constructor
ParametricLine();
/** Constructor, par is a pointer to the parent object, controlpoly the controlpolygon
*/
ParametricLine( ParametricLine *par, QVector<Point3D *> *controlpoly );
ParametricLine( ParametricLine *par, QVector<QgsPoint *> *controlpoly );
virtual ~ParametricLine();
virtual void add( ParametricLine *pl ) = 0;
virtual void calcFirstDer( float t, Vector3D *v ) = 0;
virtual void calcSecDer( float t, Vector3D *v ) = 0;
//virtual Point3D calcPoint(float t);
virtual void calcPoint( float t, Point3D * ) = 0;
//virtual QgsPoint calcPoint(float t);
virtual void calcPoint( float t, QgsPoint * ) = 0;
virtual void changeDirection() = 0;
//virtual void draw(QPainter* p);
virtual const Point3D *getControlPoint( int number ) const = 0;
virtual const QVector<Point3D *> *getControlPoly() const = 0;
virtual const QgsPoint *getControlPoint( int number ) const = 0;
virtual const QVector<QgsPoint *> *getControlPoly() const = 0;
virtual int getDegree() const = 0;
virtual ParametricLine *getParent() const = 0;
//virtual bool intersects(ParametricLine* pal);
virtual void remove( int i ) = 0;
virtual void setControlPoly( QVector<Point3D *> *cp ) = 0;
virtual void setControlPoly( QVector<QgsPoint *> *cp ) = 0;
virtual void setParent( ParametricLine *paral ) = 0;
};
@ -70,7 +70,7 @@ inline ParametricLine::ParametricLine()
}
inline ParametricLine::ParametricLine( ParametricLine *par, QVector<Point3D *> *controlpoly )
inline ParametricLine::ParametricLine( ParametricLine *par, QVector<QgsPoint *> *controlpoly )
: mDegree( 0 )
, mParent( par )
, mControlPoly( controlpoly )

45
src/analysis/interpolation/Point3D.cc
View File

@ -1,45 +0,0 @@
/***************************************************************************
Point3D.cc - description
-------------------
copyright : (C) 2004 by Marco Hugentobler
email : mhugent@geo.unizh.ch
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#include "Point3D.h"
#include "qgslogger.h"
bool Point3D::operator==( const Point3D &p ) const
{
return ( mX == p.getX() && mY == p.getY() && mZ == p.getZ() );
}
bool Point3D::operator!=( const Point3D &p ) const
{
return ( !( ( *this ) == p ) );
}
double Point3D::dist3D( Point3D *p ) const
{
if ( p )
{
double dx, dy, dz;
dx = p->getX() - getX();
dy = p->getY() - getY();
dz = p->getZ() - getZ();
return sqrt( dx * dx + dy * dy + dz * dz );
}
else
{
QgsDebugMsg( "" );
return 0;
}
}

108
src/analysis/interpolation/Point3D.h
View File

@ -1,108 +0,0 @@
/***************************************************************************
Point3D.h - description
-------------------
copyright : (C) 2004 by Marco Hugentobler
email : mhugent@geo.unizh.ch
***************************************************************************/
/***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************/
#ifndef POINT3D_H
#define POINT3D_H
#include <cmath>
#include "qgis_analysis.h"
/** \ingroup analysis
* Point3D is a class to represent a three dimensional point*/
class ANALYSIS_EXPORT Point3D
{
protected:
//! X-coordinate
double mX;
//! Y-coordinate
double mY;
//! Z-coordinate
double mZ;
public:
Point3D();
//! Constructor with the x-, y- and z-coordinate as arguments
Point3D( double x, double y, double z );
Point3D( const Point3D &p );
bool operator==( const Point3D &p ) const;
bool operator!=( const Point3D &p ) const;
//! Calculates the three-dimensional distance to another point
double dist3D( Point3D *p ) const;
//! Returns the x-coordinate of the point
double getX() const;
//! Returns the y-coordinate of the point
double getY() const;
//! Returns the z-coordinate of the point
double getZ() const;
//! Sets the x-coordinate of the point
void setX( double x );
//! Sets the y-coordinate of the point
void setY( double y );
//! Sets the z-coordinate of the point
void setZ( double z );
};
inline Point3D::Point3D()
: mX( 0 )
, mY( 0 )
, mZ( 0 )
{
}
inline Point3D::Point3D( double x, double y, double z )
: mX( x )
, mY( y )
, mZ( z )
{
}
inline Point3D::Point3D( const Point3D &p ): mX( p.mX ), mY( p.mY ), mZ( p.mZ )
{
}
inline double Point3D::getX() const
{
return mX;
}
inline double Point3D::getY() const
{
return mY;
}
inline double Point3D::getZ() const
{
return mZ;
}
inline void Point3D::setX( double x )
{
mX = x;
}
inline void Point3D::setY( double y )
{
mY = y;
}
inline void Point3D::setZ( double z )
{
mZ = z;
}
#endif

12
src/analysis/interpolation/TriDecorator.cc
View File

@ -29,7 +29,7 @@ void TriDecorator::addLine( Line3D *line, bool breakline )
}
}
int TriDecorator::addPoint( Point3D *p )
int TriDecorator::addPoint( QgsPoint *p )
{
if ( mTIN )
{
@ -69,7 +69,7 @@ bool TriDecorator::calcNormal( double x, double y, Vector3D *result )
}
}
bool TriDecorator::calcPoint( double x, double y, Point3D *result )
bool TriDecorator::calcPoint( double x, double y, QgsPoint *result )
{
if ( mTIN )
{
@ -83,11 +83,11 @@ bool TriDecorator::calcPoint( double x, double y, Point3D *result )
}
}
Point3D *TriDecorator::getPoint( unsigned int i ) const
QgsPoint *TriDecorator::getPoint( unsigned int i ) const
{
if ( mTIN )
{
Point3D *p = mTIN->getPoint( i );
QgsPoint *p = mTIN->getPoint( i );
return p;
}
else
@ -97,7 +97,7 @@ Point3D *TriDecorator::getPoint( unsigned int i ) const
}
}
bool TriDecorator::getTriangle( double x, double y, Point3D *p1, int *n1, Point3D *p2, int *n2, Point3D *p3, int *n3 )
bool TriDecorator::getTriangle( double x, double y, QgsPoint *p1, int *n1, QgsPoint *p2, int *n2, QgsPoint *p3, int *n3 )
{
if ( mTIN )
{
@ -111,7 +111,7 @@ bool TriDecorator::getTriangle( double x, double y, Point3D *p1, int *n1, Point3
}
}
bool TriDecorator::getTriangle( double x, double y, Point3D *p1, Point3D *p2, Point3D *p3 )
bool TriDecorator::getTriangle( double x, double y, QgsPoint *p1, QgsPoint *p2, QgsPoint *p3 )
{
if ( mTIN )
{

10
src/analysis/interpolation/TriDecorator.h
View File

@ -30,17 +30,17 @@ class ANALYSIS_EXPORT TriDecorator : public Triangulation
explicit TriDecorator( Triangulation *t );
virtual ~TriDecorator();
virtual void addLine( Line3D *line, bool breakline ) override;
virtual int addPoint( Point3D *p ) override;
virtual int addPoint( QgsPoint *p ) override;
//! Adds an association to a triangulation
virtual void addTriangulation( Triangulation *t );
//! Performs a consistency check, remove this later
virtual void performConsistencyTest() override;
virtual bool calcNormal( double x, double y, Vector3D *result ) override;
virtual bool calcPoint( double x, double y, Point3D *result ) override;
virtual Point3D *getPoint( unsigned int i ) const override;
virtual bool calcPoint( double x, double y, QgsPoint *result ) override;
virtual QgsPoint *getPoint( unsigned int i ) const override;
virtual int getNumberOfPoints() const override;
bool getTriangle( double x, double y, Point3D *p1, int *n1 SIP_OUT, Point3D *p2 SIP_OUT, int *n2 SIP_OUT, Point3D *p3 SIP_OUT, int *n3 SIP_OUT ) override;
bool getTriangle( double x, double y, Point3D *p1 SIP_OUT, Point3D *p2 SIP_OUT, Point3D *p3 SIP_OUT ) override;
bool getTriangle( double x, double y, QgsPoint *p1, int *n1 SIP_OUT, QgsPoint *p2 SIP_OUT, int *n2 SIP_OUT, QgsPoint *p3 SIP_OUT, int *n3 SIP_OUT ) override;
bool getTriangle( double x, double y, QgsPoint *p1 SIP_OUT, QgsPoint *p2 SIP_OUT, QgsPoint *p3 SIP_OUT ) override;
virtual int getOppositePoint( int p1, int p2 ) override;
virtual QList<int> *getSurroundingTriangles( int pointno ) override;
virtual double getXMax() const override;

4
src/analysis/interpolation/TriangleInterpolator.h
View File

@ -17,7 +17,7 @@
#ifndef TINTERPOLATOR_H
#define TINTERPOLATOR_H
#include <Point3D.h>
#include <qgspoint.h>
#include <Vector3D.h>
#include "qgis_analysis.h"
@ -30,7 +30,7 @@ class ANALYSIS_EXPORT TriangleInterpolator
//! Calculates the normal vector and assigns it to vec
virtual bool calcNormVec( double x, double y, Vector3D *result ) = 0;
//! Performs a linear interpolation in a triangle and assigns the x-,y- and z-coordinates to point
virtual bool calcPoint( double x, double y, Point3D *result ) = 0;
virtual bool calcPoint( double x, double y, QgsPoint *result ) = 0;
};
#endif

10
src/analysis/interpolation/Triangulation.h
View File

@ -49,7 +49,7 @@ class ANALYSIS_EXPORT Triangulation
* Adds a point to the triangulation
* Ownership is transferred to this class
*/
virtual int addPoint( Point3D *p ) = 0;
virtual int addPoint( QgsPoint *p ) = 0;
/**
* Calculates the normal at a point on the surface and assigns it to 'result'.
@ -64,18 +64,18 @@ class ANALYSIS_EXPORT Triangulation
* Calculates x-, y and z-value of the point on the surface and assigns it to 'result'.
* Returns true in case of success and flase in case of failure
*/
virtual bool calcPoint( double x, double y, Point3D *result ) = 0;
virtual bool calcPoint( double x, double y, QgsPoint *result ) = 0;
//! Returns a pointer to the point with number i. Any virtual points must have the number -1
virtual Point3D *getPoint( unsigned int i ) const = 0;
virtual QgsPoint *getPoint( unsigned int i ) const = 0;
/** Finds out in which triangle the point with coordinates x and y is and
* assigns the numbers of the vertices to 'n1', 'n2' and 'n3' and the vertices to 'p1', 'p2' and 'p3'
*/
virtual bool getTriangle( double x, double y, Point3D *p1, int *n1 SIP_OUT, Point3D *p2 SIP_OUT, int *n2 SIP_OUT, Point3D *p3 SIP_OUT, int *n3 SIP_OUT ) = 0;
virtual bool getTriangle( double x, double y, QgsPoint *p1, int *n1 SIP_OUT, QgsPoint *p2 SIP_OUT, int *n2 SIP_OUT, QgsPoint *p3 SIP_OUT, int *n3 SIP_OUT ) = 0;
//! Finds out, in which triangle the point with coordinates x and y is and assigns the points at the vertices to 'p1', 'p2' and 'p3
virtual bool getTriangle( double x, double y, Point3D *p1 SIP_OUT, Point3D *p2 SIP_OUT, Point3D *p3 SIP_OUT ) = 0;
virtual bool getTriangle( double x, double y, QgsPoint *p1 SIP_OUT, QgsPoint *p2 SIP_OUT, QgsPoint *p3 SIP_OUT ) = 0;
//! Returns the number of the point opposite to the triangle points p1, p2 (which have to be on a halfedge)
virtual int getOppositePoint( int p1, int p2 ) = 0;

2
src/analysis/interpolation/Vector3D.h
View File

@ -22,7 +22,7 @@
/** \ingroup analysis
* Class Vector3D represents a 3D-Vector, capable to store x-,y- and
* z-coordinates in double values. In fact, the class is the same as Point3D.
* z-coordinates in double values. In fact, the class is the same as QgsPoint.
* The name 'vector' makes it easier to understand the programs.
*/

24
src/analysis/interpolation/qgstininterpolator.cpp
View File

@ -22,7 +22,7 @@
#include "NormVecDecorator.h"
#include "LinTriangleInterpolator.h"
#include "Line3D.h"
#include "Point3D.h"
#include "qgspoint.h"
#include "qgsfeature.h"
#include "qgsgeometry.h"
#include "qgsvectorlayer.h"
@ -58,12 +58,12 @@ int QgsTINInterpolator::interpolatePoint( double x, double y, double &result )
return 1;
}
Point3D r;
QgsPoint r;
if ( !mTriangleInterpolator->calcPoint( x, y, &r ) )
{
return 2;
}
result = r.getZ();
result = r.z();
return 0;
}
@ -224,7 +224,7 @@ int QgsTINInterpolator::insertData( QgsFeature *f, bool zCoord, int attr, InputT
{
z = attributeValue;
}
Point3D *point = new Point3D( x, y, z );
QgsPoint *point = new QgsPoint( x, y, z );
if ( mTriangulation->addPoint( point ) == -100 )
{
return -1;
@ -279,11 +279,11 @@ int QgsTINInterpolator::insertData( QgsFeature *f, bool zCoord, int attr, InputT
if ( type == POINTS )
{
//todo: handle error code -100
mTriangulation->addPoint( new Point3D( x, y, z ) );
mTriangulation->addPoint( new QgsPoint( x, y, z ) );
}
else
{
line->insertPoint( new Point3D( x, y, z ) );
line->insertPoint( new QgsPoint( x, y, z ) );
}
}
@ -323,11 +323,11 @@ int QgsTINInterpolator::insertData( QgsFeature *f, bool zCoord, int attr, InputT
if ( type == POINTS )
{
//todo: handle error code -100
mTriangulation->addPoint( new Point3D( x, y, z ) );
mTriangulation->addPoint( new QgsPoint( x, y, z ) );
}
else
{
line->insertPoint( new Point3D( x, y, z ) );
line->insertPoint( new QgsPoint( x, y, z ) );
}
}
if ( type != POINTS )
@ -367,11 +367,11 @@ int QgsTINInterpolator::insertData( QgsFeature *f, bool zCoord, int attr, InputT
if ( type == POINTS )
{
//todo: handle error code -100
mTriangulation->addPoint( new Point3D( x, y, z ) );
mTriangulation->addPoint( new QgsPoint( x, y, z ) );
}
else
{
line->insertPoint( new Point3D( x, y, z ) );
line->insertPoint( new QgsPoint( x, y, z ) );
}
}
@ -417,11 +417,11 @@ int QgsTINInterpolator::insertData( QgsFeature *f, bool zCoord, int attr, InputT
if ( type == POINTS )
{
//todo: handle error code -100
mTriangulation->addPoint( new Point3D( x, y, z ) );
mTriangulation->addPoint( new QgsPoint( x, y, z ) );
}
else
{
line->insertPoint( new Point3D( x, y, z ) );
line->insertPoint( new QgsPoint( x, y, z ) );
}
}
if ( type != POINTS )

6
src/core/geometry/qgspoint.cpp
View File

@ -31,12 +31,12 @@
* See details in QEP #17
****************************************************************************/
QgsPoint::QgsPoint( double x, double y )
QgsPoint::QgsPoint( double x, double y, double z, double m )
: QgsAbstractGeometry()
, mX( x )
, mY( y )
, mZ( 0.0 )
, mM( 0.0 )
, mZ( z )
, mM( m )
{
mWkbType = QgsWkbTypes::Point;
}

7
src/core/geometry/qgspoint.h
View File

@ -44,11 +44,10 @@ class CORE_EXPORT QgsPoint: public QgsAbstractGeometry
public:
/** Construct a 2 dimensional point with an initial x and y coordinate.
* \param x x-coordinate of point
* \param y y-coordinate of point
/**
* Construct a point with the provided initial coordinate values.
*/
QgsPoint( double x = 0.0, double y = 0.0 );
QgsPoint( double x = 0.0, double y = 0.0, double z = 0.0, double m = 0.0 );
/** Construct a QgsPoint from a QgsPointXY object
*/

8
tests/src/python/acceptable_missing_doc.py
View File

@ -208,7 +208,7 @@ ACCEPTABLE_MISSING_DOCS = {
"QgsColorBrewerPalette": ["listSchemes()", "listSchemeVariants(const QString &schemeName)", "listSchemeColors(const QString &schemeName, int colors)"],
"QgsGraduatedSymbolRendererWidget": ["updateGraduatedSymbolIcon()", "classifyGraduated()", "selectedRanges()", "sizeScaleFieldChanged(const QString &fldName)", "rowsMoved()", "rangesClicked(const QModelIndex &idx)", "labelFormatChanged()", "disconnectUpdateHandlers()", "refreshRanges(bool reset=false)", "create(QgsVectorLayer *layer, QgsStyle *style, QgsFeatureRenderer *renderer)", "reapplyColorRamp()", "scaleMethodChanged(QgsSymbol::ScaleMethod scaleMethod)", "changeRange(int rangeIdx)", "findSymbolForRange(double lowerBound, double upperBound, const QgsRangeList &ranges) const ", "changeCurrentValue(QStandardItem *item)", "changeRangeSymbol(int rangeIdx)", "connectUpdateHandlers()", "modelDataChanged()", "changeSelectedSymbols()", "showSymbolLevels()", "rowsOrdered()", "graduatedColumnChanged(const QString &field)", "reapplySizes()", "QgsGraduatedSymbolRendererWidget(QgsVectorLayer *layer, QgsStyle *style, QgsFeatureRenderer *renderer)", "changeGraduatedSymbol()", "rangesDoubleClicked(const QModelIndex &idx)", "updateUiFromRenderer(bool updateCount=true)"],
"QgsGpsdConnection": ["QgsGpsdConnection(const QString &host, qint16 port, const QString &device)"],
"Point3D": ["Point3D(const Point3D &p)"],
"QgsPoint": ["QgsPoint(const QgsPoint &p)"],
"QgsVectorLayerUndoCommand": ["cache()", "QgsVectorLayerUndoCommand(QgsVectorLayerEditBuffer *buffer)", "layer()"],
"QgsDataDefinedRotationDialog": ["QgsDataDefinedRotationDialog(const QList< QgsSymbol * > &symbolList, QgsVectorLayer *layer)"],
"QgsException": ["QgsException(QString const &what)"],
@ -256,7 +256,7 @@ ACCEPTABLE_MISSING_DOCS = {
"QgsSymbolRenderContext": ["renderContext() const ", "setMapUnitScale(const QgsMapUnitScale &scale)", "renderContext()", "outputUnit() const ", "outputLineWidth(double width) const ", "outputPixelSize(double size) const ", "selected() const ", "setRenderHints(int hints)", "QgsSymbolRenderContext(QgsRenderContext &c, QgsSymbol::OutputUnit u, qreal alpha=1.0, bool selected=false, int renderHints=0, const QgsFeature *f=nullptr, const QgsFields *fields=nullptr, const QgsMapUnitScale &mapUnitScale=QgsMapUnitScale())", "setSelected(bool selected)", "setOutputUnit(QgsSymbol::OutputUnit u)", "setFeature(const QgsFeature *f)", "renderHints() const ", "mapUnitScale() const "],
"QgsIdentifyMenu": ["MenuLevel", "execWithSingleResult()", "allowMultipleReturn()", "resultsIfExternalAction()", "showFeatureActions()", "maxFeatureDisplay()", "maxLayerDisplay()"],
"QgsPointLocator_VisitorNearestEdge": ["QgsPointLocator_VisitorNearestEdge(QgsPointLocator *pl, QgsPointLocator::Match &m, const QgsPointXY &srcPoint, QgsPointLocator::MatchFilter *filter=nullptr)"],
"ParametricLine": ["getControlPoint(int number) const =0", "getDegree() const =0", "calcPoint(float t, Point3D *)=0", "calcFirstDer(float t, Vector3D *v)=0", "remove(int i)=0", "setParent(ParametricLine *paral)=0", "calcSecDer(float t, Vector3D *v)=0", "getParent() const =0", "changeDirection()=0", "add(ParametricLine *pl)=0", "getControlPoly() const =0", "setControlPoly(QVector< Point3D * > *cp)=0"],
"ParametricLine": ["getControlPoint(int number) const =0", "getDegree() const =0", "calcPoint(float t, QgsPoint *)=0", "calcFirstDer(float t, Vector3D *v)=0", "remove(int i)=0", "setParent(ParametricLine *paral)=0", "calcSecDer(float t, Vector3D *v)=0", "getParent() const =0", "changeDirection()=0", "add(ParametricLine *pl)=0", "getControlPoly() const =0", "setControlPoly(QVector< QgsPoint * > *cp)=0"],
"QgsComposerAttributeTable": ["QgsComposerAttributeTable(QgsComposition *composition)"],
"QgisInterface": ["helpToolBar()=0", "actionToggleFullScreen()=0", "newLayerMenu()=0", "actionNodeTool()=0", "actionCopyFeatures()=0", "actionNewProject()=0", "actionNewVectorLayer()=0", "actionCancelEdits()=0", "actionSplitParts()=0", "actionShowComposerManager()=0", "actionHelpContents()=0", "actionOpenFieldCalculator()=0", "actionSaveProjectAs()=0", "actionLayerProperties()=0", "helpMenu()=0", "actionSimplifyFeature()=0", "actionAddAllToOverview()=0", "databaseMenu()=0", "actionCheckQgisVersion()=0", "actionSaveProject()=0", "actionCancelAllEdits()=0", "layerTreeView()=0", "advancedDigitizeToolBar()=0", "actionAbout()=0", "webMenu()=0", "actionDuplicateLayer()=0", "actionShowPythonDialog()=0", "vectorToolBar()=0", "actionPasteFeatures()=0", "actionPrintComposer()=0", "actionSplitFeatures()=0", "actionPluginListSeparator()=0", "pluginManagerInterface()=0", "actionRemoveAllFromOverview()=0", "digitizeToolBar()=0", "actionSaveActiveLayerEdits()=0", "firstRightStandardMenu()=0", "projectMenu()=0", "actionProjectProperties()=0", "actionCopyLayerStyle()=0", "actionAllEdits()=0", "actionAddOgrLayer()=0", "actionDeleteRing()=0", "actionAddWmsLayer()=0", "settingsMenu()=0", "actionAddRing()=0", "fileToolBar()=0", "editMenu()=0", "pluginMenu()=0", "mapNavToolToolBar()=0", "vectorMenu()=0", "actionCutFeatures()=0", "actionSaveMapAsImage()=0", "actionMoveFeature()=0", "windowMenu()=0", "actionShowSelectedLayers()=0", "actionDeleteSelected()=0", "actionOptions()=0", "layerToolBar()=0", "actionDeletePart()=0", "actionAddPgLayer()=0", "databaseToolBar()=0", "actionRollbackAllEdits()=0", "actionShowAllLayers()=0", "viewMenu()=0", "actionRollbackEdits()=0", "actionSaveAllEdits()=0", "actionPasteLayerStyle()=0", "actionQgisHomePage()=0", "attributesToolBar()=0", "actionExit()=0", "rasterToolBar()=0", "actionHideAllLayers()=0", "actionAddToOverview()=0", "layerMenu()=0", "actionToggleEditing()=0", "actionLayerSaveAs()=0", "actionManagePlugins()=0", "actionCustomProjection()=0", "rasterMenu()=0", "actionAddPart()=0", "pluginToolBar()=0", "actionRemoveLayer()=0", "actionOpenTable()=0", "actionSaveEdits()=0", "actionOpenProject()=0", "webToolBar()=0", "actionAddFeature()=0", "actionHideSelectedLayers()=0", "actionAddRasterLayer()=0"],
"QgsCollapsibleGroupBoxBasic": ["QgsCollapsibleGroupBoxBasic(const QString &title, QWidget *parent=nullptr)", "checkToggled(bool ckd)", "QgsCollapsibleGroupBoxBasic(QWidget *parent=nullptr)", "init()", "clearModifiers()", "titleRect() const ", "checkClicked(bool ckd)", "updateStyle()", "toggleCollapsed()"],
@ -751,7 +751,7 @@ ACCEPTABLE_MISSING_ADDED_NOTE = [
"QgsColorWidgetAction",
"Bezier3D",
"QgsSQLStatement::NodeJoin",
"Point3D",
"QgsPoint",
"QgsSLConnect",
"QgsConditionalStyle",
"QgsExpression::NodeList",
@ -1743,7 +1743,7 @@ ACCEPTABLE_MISSING_BRIEF = ['QgsBrushStyleComboBox',
'QgsTrackedVectorLayerTools',
'QgsSymbolLayerUtils',
'QgsComposerRasterSymbolItem',
'Point3D',
'QgsPoint',
'QgsGeometryGeneratorSymbolLayerWidget',
'QgsVectorLayerFeatureIterator',
'QgsFeatureRenderer',