Skip to content
Permalink
Browse files
Fix broken TIN interpolation (loss of point z coordinates)
  • Loading branch information
nyalldawson committed Aug 18, 2017
1 parent 356588a commit 4bba95f
Show file tree
Hide file tree
Showing 5 changed files with 48 additions and 39 deletions.
@@ -164,7 +164,6 @@ bool CloughTocherInterpolator::calcNormVec( double x, double y, Vector3D *result
return true;
}

QgsDebugMsg( "warning, point outside the triangle" );
result->setX( 0 );//return a vertical normal if failed
result->setY( 0 );
result->setZ( 1 );
@@ -240,7 +239,6 @@ bool CloughTocherInterpolator::calcPoint( double x, double y, QgsPoint *result )
}
else
{
QgsDebugMsg( "warning, point outside the triangle" );
result->setX( x );
result->setY( y );
result->setZ( 0 );
@@ -33,27 +33,43 @@ 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
QgsPoint point1;
QgsPoint point1 = QgsPoint( 0, 0, 0 );
//! Second point of the triangle in x-,y-,z-coordinates
QgsPoint point2;
QgsPoint point2 = QgsPoint( 0, 0, 0 );
//! Third point of the triangle in x-,y-,z-coordinates
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;
QgsPoint point3 = QgsPoint( 0, 0, 0 );
//! Control point 1
QgsPoint cp1 = QgsPoint( 0, 0, 0 );
//! Control point 2
QgsPoint cp2 = QgsPoint( 0, 0, 0 );
//! Control point 3
QgsPoint cp3 = QgsPoint( 0, 0, 0 );
//! Control point 4
QgsPoint cp4 = QgsPoint( 0, 0, 0 );
//! Control point 5
QgsPoint cp5 = QgsPoint( 0, 0, 0 );
//! Control point 6
QgsPoint cp6 = QgsPoint( 0, 0, 0 );
//! Control point 7
QgsPoint cp7 = QgsPoint( 0, 0, 0 );
//! Control point 8
QgsPoint cp8 = QgsPoint( 0, 0, 0 );
//! Control point 9
QgsPoint cp9 = QgsPoint( 0, 0, 0 );
//! Control point 10
QgsPoint cp10 = QgsPoint( 0, 0, 0 );
//! Control point 11
QgsPoint cp11 = QgsPoint( 0, 0, 0 );
//! Control point 12
QgsPoint cp12 = QgsPoint( 0, 0, 0 );
//! Control point 13
QgsPoint cp13 = QgsPoint( 0, 0, 0 );
//! Control point 14
QgsPoint cp14 = QgsPoint( 0, 0, 0 );
//! Control point 15
QgsPoint cp15 = QgsPoint( 0, 0, 0 );
//! Control point 16
QgsPoint cp16 = QgsPoint( 0, 0, 0 );
//! Derivative in x-direction at point1
double der1X;
//! Derivative in y-direction at point1
@@ -67,11 +83,11 @@ class ANALYSIS_EXPORT CloughTocherInterpolator : public TriangleInterpolator
//! Derivative in y-direction at point3
double der3Y;
//! Stores point1 of the last run
QgsPoint lpoint1;
QgsPoint lpoint1 = QgsPoint( 0, 0, 0 );
//! Stores point2 of the last run
QgsPoint lpoint2;
QgsPoint lpoint2 = QgsPoint( 0, 0, 0 );
//! Stores point3 of the last run
QgsPoint lpoint3;
QgsPoint lpoint3 = QgsPoint( 0, 0, 0 );
//! 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
@@ -976,7 +976,6 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, QgsPoint *p1, int *
int edge = baseEdgeOfTriangle( &point );
if ( edge == -10 )//the point is outside the convex hull
{
QgsDebugMsg( "edge outside the convex hull" );
return false;
}

@@ -1094,7 +1093,6 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, QgsPoint *p1, QgsPo
int edge = baseEdgeOfTriangle( &point );
if ( edge == -10 )//the point is outside the convex hull
{
QgsDebugMsg( "edge outside the convex hull" );
return false;
}
else if ( edge >= 0 )//the point is inside the convex hull
@@ -1178,14 +1176,12 @@ bool DualEdgeTriangulation::getTriangle( double x, double y, QgsPoint *p1, QgsPo
}
else//problems
{
QgsDebugMsg( QString( "problems: the edge is: %1" ).arg( edge ) );
return false;
}
}

else
{
QgsDebugMsg( "warning, null pointer" );
return false;
}
}
@@ -1270,7 +1266,7 @@ 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
{
QgsPoint crosspoint;
QgsPoint crosspoint( 0, 0, 0 );
int p3, p4;
p3 = mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getDual()]->getPoint();
@@ -1292,7 +1288,7 @@ 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
{
QgsPoint crosspoint;
QgsPoint crosspoint( 0, 0, 0 );
int p3, p4;
p3 = mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getDual()]->getPoint();
@@ -1358,7 +1354,7 @@ 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
{
QgsPoint crosspoint;
QgsPoint crosspoint( 0, 0, 0 );
int p3, p4;
p3 = mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
@@ -1380,7 +1376,7 @@ 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
{
QgsPoint crosspoint;
QgsPoint crosspoint( 0, 0, 0 );
int p3, p4;
p3 = mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
@@ -1405,7 +1401,7 @@ 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
{
QgsPoint crosspoint;
QgsPoint crosspoint( 0, 0, 0 );
int p3, p4;
p3 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint();
@@ -1427,7 +1423,7 @@ 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
{
QgsPoint crosspoint;
QgsPoint crosspoint( 0, 0, 0 );
int p3, p4;
p3 = mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint();
p4 = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint();
@@ -1769,7 +1765,7 @@ void DualEdgeTriangulation::ruppertRefinement()
int minedgenext;
int minedgenextnext;

QgsPoint circumcenter;
QgsPoint circumcenter( 0, 0, 0 );

while ( !edge_angle.empty() )
{
@@ -2205,7 +2201,7 @@ void DualEdgeTriangulation::ruppertRefinement()

/*******otherwise, try to add the circumcenter to the triangulation************************************************************************************************/

QgsPoint *p = new QgsPoint();
QgsPoint *p = new QgsPoint( 0, 0, 0 );
mDecorator->calcPoint( circumcenter.x(), circumcenter.y(), p );
int pointno = mDecorator->addPoint( p );

@@ -3218,7 +3214,7 @@ int DualEdgeTriangulation::splitHalfEdge( int edge, float position )
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
QgsPoint zvaluepoint;
QgsPoint zvaluepoint( 0, 0, 0 );
mDecorator->calcPoint( p->x(), p->y(), &zvaluepoint );
p->setZ( zvaluepoint.z() );

@@ -350,7 +350,6 @@ bool NormVecDecorator::getTriangle( double x, double y, QgsPoint *p1, int *ptn1,
}
else
{
QgsDebugMsg( "warning, getTriangle returned false" );
return false;
}

@@ -58,7 +58,7 @@ int QgsTINInterpolator::interpolatePoint( double x, double y, double &result )
return 1;
}

QgsPoint r;
QgsPoint r( 0, 0, 0 );
if ( !mTriangleInterpolator->calcPoint( x, y, &r ) )
{
return 2;

0 comments on commit 4bba95f

Please sign in to comment.