mirror of
https://github.com/qgis/QGIS.git
synced 2025-04-15 00:04:00 -04:00
aceddae65a
Deprecate some QgsFeature methods which take or return pointers, update other classes as required
3292 lines
120 KiB
C++
3292 lines
120 KiB
C++
/***************************************************************************
|
|
DualEdgeTriangulation.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 "DualEdgeTriangulation.h"
|
|
#include <map>
|
|
#include "Line3D.h"
|
|
#include "MathUtils.h"
|
|
#include "qgsgeometry.h"
|
|
#include "qgslogger.h"
|
|
#include "qgsvectorfilewriter.h"
|
|
|
|
double leftOfTresh = 0.00000001;
|
|
|
|
DualEdgeTriangulation::~DualEdgeTriangulation()
|
|
{
|
|
//remove all the points
|
|
if ( !mPointVector.isEmpty() )
|
|
{
|
|
for ( int i = 0; i < mPointVector.count(); i++ )
|
|
{
|
|
delete mPointVector[i];
|
|
}
|
|
}
|
|
|
|
//remove all the HalfEdge
|
|
if ( !mHalfEdge.isEmpty() )
|
|
{
|
|
for ( int i = 0; i < mHalfEdge.count(); i++ )
|
|
{
|
|
delete mHalfEdge[i];
|
|
}
|
|
}
|
|
}
|
|
|
|
void DualEdgeTriangulation::performConsistencyTest()
|
|
{
|
|
QgsDebugMsg( "performing consistency test" );
|
|
|
|
for ( int i = 0; i < mHalfEdge.count(); i++ )
|
|
{
|
|
int a = mHalfEdge[mHalfEdge[i]->getDual()]->getDual();
|
|
int b = mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getNext();
|
|
if ( i != a )
|
|
{
|
|
QgsDebugMsg( "warning, first test failed" );
|
|
}
|
|
if ( i != b )
|
|
{
|
|
QgsDebugMsg( "warning, second test failed" );
|
|
}
|
|
}
|
|
QgsDebugMsg( "consistency test finished" );
|
|
}
|
|
|
|
void DualEdgeTriangulation::addLine( Line3D* line, bool breakline )
|
|
{
|
|
int actpoint = -10;//number of the last point, which has been inserted from the line
|
|
int currentpoint = -10;//number of the point, which is currently inserted from the line
|
|
if ( line )
|
|
{
|
|
//first, find the first point
|
|
unsigned int i;
|
|
line->goToBegin();
|
|
|
|
for ( i = 0; i < line->getSize(); i++ )
|
|
{
|
|
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() ) );
|
|
if ( actpoint != -100 )
|
|
{
|
|
i++;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if ( actpoint == -100 )//no point of the line could be inserted
|
|
{
|
|
delete line;
|
|
return;
|
|
}
|
|
|
|
for ( ; i < line->getSize(); i++ )
|
|
{
|
|
line->goToNext();
|
|
currentpoint = mDecorator->addPoint( new Point3D( *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 );
|
|
}
|
|
actpoint = currentpoint;
|
|
}
|
|
}
|
|
delete line;
|
|
}
|
|
|
|
int DualEdgeTriangulation::addPoint( Point3D* p )
|
|
{
|
|
if ( p )
|
|
{
|
|
// QgsDebugMsg( QString("inserting point %1,%2//%3//%4").arg(mPointVector.count()).arg(p->getX()).arg(p->getY()).arg(p->getZ()));
|
|
|
|
//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();
|
|
}
|
|
|
|
else//update bounding box else
|
|
{
|
|
if (( *p ).getX() < xMin )
|
|
{
|
|
xMin = ( *p ).getX();
|
|
}
|
|
else if (( *p ).getX() > xMax )
|
|
{
|
|
xMax = ( *p ).getX();
|
|
}
|
|
|
|
if (( *p ).getY() < yMin )
|
|
{
|
|
yMin = ( *p ).getY();
|
|
}
|
|
else if (( *p ).getY() > yMax )
|
|
{
|
|
yMax = ( *p ).getY();
|
|
}
|
|
}
|
|
|
|
//then update mPointVector
|
|
mPointVector.append( p );
|
|
|
|
//then update the HalfEdgeStructure
|
|
if ( mPointVector.count() == 1 )//insert the first point into the triangulation
|
|
{
|
|
unsigned int zedge = insertEdge( -10, -10, -1, false, false );//edge pointing from p to the virtual point
|
|
unsigned int fedge = insertEdge(( int )zedge, ( int )zedge, 0, false, false );//edge pointing from the virtual point to p
|
|
( mHalfEdge.at( zedge ) )->setDual(( int )fedge );
|
|
( mHalfEdge.at( zedge ) )->setNext(( int )fedge );
|
|
|
|
}
|
|
|
|
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() )
|
|
{
|
|
QgsDebugMsg( "second point is the same as the first point, it thus has not been inserted" );
|
|
Point3D* p = mPointVector[1];
|
|
mPointVector.remove( 1 );
|
|
delete p;
|
|
return -100;
|
|
}
|
|
|
|
unsigned int sedge = insertEdge( -10, -10, 1, false, false );//edge pointing from point 0 to point 1
|
|
unsigned int tedge = insertEdge(( int )sedge, 0, 0, false, false );//edge pointing from point 1 to point 0
|
|
unsigned int foedge = insertEdge( -10, 4, 1, false, false );//edge pointing from the virtual point to point 1
|
|
unsigned int fiedge = insertEdge(( int )foedge, 1, -1, false, false );//edge pointing from point 2 to the virtual point
|
|
mHalfEdge.at( sedge )->setDual(( int )tedge );
|
|
mHalfEdge.at( sedge )->setNext(( int )fiedge );
|
|
mHalfEdge.at( foedge )->setDual(( int )fiedge );
|
|
mHalfEdge.at( foedge )->setNext(( int )tedge );
|
|
mHalfEdge.at( 0 )->setNext(( int )foedge );
|
|
mHalfEdge.at( 1 )->setNext(( int )sedge );
|
|
|
|
mEdgeInside = 3;
|
|
}
|
|
|
|
else if ( mPointVector.count() == 3 )//insert the third point into the triangulation
|
|
{
|
|
//we first have to decide, if the third point is to the left or to the right of the line through p0 and p1
|
|
double number = MathUtils::leftOf( p, mPointVector[0], mPointVector[1] );
|
|
if ( number < -leftOfTresh )//p is on the left side
|
|
{
|
|
//insert six new edges
|
|
unsigned int edgea = insertEdge( -10, -10, 2, false, false );//edge pointing from point1 to point2
|
|
unsigned int edgeb = insertEdge(( int )edgea, 5, 1, false, false );//edge pointing from point2 to point1
|
|
unsigned int edgec = insertEdge( -10, ( int )edgeb, 2, false, false );//edge pointing from the virtual point to p2
|
|
unsigned int edged = insertEdge( -10, 2, 0, false, false );//edge pointing from point2 to point0
|
|
unsigned int edgee = insertEdge(( int )edged, -10, 2, false, false );//edge pointing from point0 to point2
|
|
unsigned int edgef = insertEdge(( int )edgec, 1, -1, false, false );//edge pointing from point2 to the virtual point
|
|
mHalfEdge.at( edgea )->setDual(( int )edgeb );
|
|
mHalfEdge.at( edgea )->setNext(( int )edged );
|
|
mHalfEdge.at( edgec )->setDual(( int )edgef );
|
|
mHalfEdge.at( edged )->setDual(( int )edgee );
|
|
mHalfEdge.at( edgee )->setNext(( int )edgef );
|
|
mHalfEdge.at( 5 )->setNext(( int )edgec );
|
|
mHalfEdge.at( 1 )->setNext(( int )edgee );
|
|
mHalfEdge.at( 2 )->setNext(( int )edgea );
|
|
}
|
|
|
|
else if ( number > leftOfTresh )//p is on the right side
|
|
{
|
|
//insert six new edges
|
|
unsigned int edgea = insertEdge( -10, -10, 2, false, false );//edge pointing from p0 to p2
|
|
unsigned int edgeb = insertEdge(( int )edgea, 0, 0, false, false );//edge pointing from p2 to p0
|
|
unsigned int edgec = insertEdge( -10, ( int )edgeb, 2, false, false );//edge pointing from the virtual point to p2
|
|
unsigned int edged = insertEdge( -10, 3, 1, false, false );//edge pointing from p2 to p1
|
|
unsigned int edgee = insertEdge(( int )edged, -10, 2, false, false );//edge pointing from p1 to p2
|
|
unsigned int edgef = insertEdge(( int )edgec, 4, -1, false, false );//edge pointing from p2 to the virtual point
|
|
mHalfEdge.at( edgea )->setDual(( int )edgeb );
|
|
mHalfEdge.at( edgea )->setNext(( int )edged );
|
|
mHalfEdge.at( edgec )->setDual(( int )edgef );
|
|
mHalfEdge.at( edged )->setDual(( int )edgee );
|
|
mHalfEdge.at( edgee )->setNext(( int )edgef );
|
|
mHalfEdge.at( 0 )->setNext(( int )edgec );
|
|
mHalfEdge.at( 4 )->setNext(( int )edgee );
|
|
mHalfEdge.at( 3 )->setNext(( int )edgea );
|
|
}
|
|
|
|
else//p is in a line with p0 and p1
|
|
{
|
|
mPointVector.remove( mPointVector.count() - 1 );
|
|
QgsDebugMsg( "error: third point is on the same line as the first and the second point. It thus has not been inserted into the triangulation" );
|
|
return -100;
|
|
}
|
|
}
|
|
|
|
else if ( mPointVector.count() >= 4 )//insert an arbitrary point into the triangulation
|
|
{
|
|
int number = baseEdgeOfTriangle( p );
|
|
|
|
//point is outside the convex hull----------------------------------------------------
|
|
if ( number == -10 )
|
|
{
|
|
unsigned int cwedge = mEdgeOutside;//the last visible edge clockwise from mEdgeOutside
|
|
unsigned int ccwedge = mEdgeOutside;//the last visible edge counterclockwise from mEdgeOutside
|
|
|
|
//mEdgeOutside is in each case visible
|
|
mHalfEdge[mHalfEdge[mEdgeOutside]->getNext()]->setPoint( mPointVector.count() - 1 );
|
|
|
|
//find cwedge and replace the virtual point with the new point when necessary
|
|
while ( MathUtils::leftOf( mPointVector[( unsigned int ) mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[cwedge]->getNext()]->getDual()]->getNext()]->getPoint()], p, mPointVector[( unsigned int ) mHalfEdge[cwedge]->getPoint()] ) < ( -leftOfTresh ) )
|
|
{
|
|
//set the point number of the necessary edge to the actual point instead of the virtual point
|
|
mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[cwedge]->getNext()]->getDual()]->getNext()]->getNext()]->setPoint( mPointVector.count() - 1 );
|
|
//advance cwedge one edge further clockwise
|
|
cwedge = ( unsigned int )mHalfEdge[mHalfEdge[mHalfEdge[cwedge]->getNext()]->getDual()]->getNext();
|
|
}
|
|
|
|
//build the necessary connections with the virtual point
|
|
unsigned int edge1 = insertEdge( mHalfEdge[cwedge]->getNext(), -10, mHalfEdge[cwedge]->getPoint(), false, false );//edge pointing from the new point to the last visible point clockwise
|
|
unsigned int edge2 = insertEdge( mHalfEdge[mHalfEdge[cwedge]->getNext()]->getDual(), -10, -1, false, false );//edge pointing from the last visible point to the virtual point
|
|
unsigned int edge3 = insertEdge( -10, edge1, mPointVector.count() - 1, false, false );//edge pointing from the virtual point to new point
|
|
|
|
//adjust the other pointers
|
|
mHalfEdge[mHalfEdge[mHalfEdge[cwedge]->getNext()]->getDual()]->setDual( edge2 );
|
|
mHalfEdge[mHalfEdge[cwedge]->getNext()]->setDual( edge1 );
|
|
mHalfEdge[edge1]->setNext( edge2 );
|
|
mHalfEdge[edge2]->setNext( edge3 );
|
|
|
|
|
|
|
|
//find ccwedge and replace the virtual point with the new point when necessary
|
|
while ( MathUtils::leftOf( mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[ccwedge]->getNext()]->getNext()]->getPoint()], mPointVector[mPointVector.count()-1], mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwedge]->getNext()]->getNext()]->getDual()]->getNext()]->getPoint()] ) < ( -leftOfTresh ) )
|
|
{
|
|
//set the point number of the necessary edge to the actual point instead of the virtual point
|
|
mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwedge]->getNext()]->getNext()]->getDual()]->setPoint( mPointVector.count() - 1 );
|
|
//advance ccwedge one edge further counterclockwise
|
|
ccwedge = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwedge]->getNext()]->getNext()]->getDual()]->getNext()]->getNext();
|
|
}
|
|
|
|
//build the necessary connections with the virtual point
|
|
unsigned int edge4 = insertEdge( mHalfEdge[mHalfEdge[ccwedge]->getNext()]->getNext(), -10, mPointVector.count() - 1, false, false );//points from the last visible point counterclockwise to the new point
|
|
unsigned int edge5 = insertEdge( edge3, -10, -1, false, false );//points from the new point to the virtual point
|
|
unsigned int edge6 = insertEdge( mHalfEdge[mHalfEdge[mHalfEdge[ccwedge]->getNext()]->getNext()]->getDual(), edge4, mHalfEdge[mHalfEdge[ccwedge]->getDual()]->getPoint(), false, false );//points from the virtual point to the last visible point counterclockwise
|
|
|
|
|
|
|
|
//adjust the other pointers
|
|
mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[ccwedge]->getNext()]->getNext()]->getDual()]->setDual( edge6 );
|
|
mHalfEdge[mHalfEdge[mHalfEdge[ccwedge]->getNext()]->getNext()]->setDual( edge4 );
|
|
mHalfEdge[edge4]->setNext( edge5 );
|
|
mHalfEdge[edge5]->setNext( edge6 );
|
|
mHalfEdge[edge3]->setDual( edge5 );
|
|
|
|
//now test the HalfEdge at the former convex hull for swappint
|
|
unsigned int index = ccwedge;
|
|
unsigned int toswap;
|
|
while ( true )
|
|
{
|
|
toswap = index;
|
|
index = mHalfEdge[mHalfEdge[mHalfEdge[index]->getNext()]->getDual()]->getNext();
|
|
checkSwap( toswap, 0 );
|
|
if ( toswap == cwedge )
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
//point is inside the convex hull------------------------------------------------------
|
|
|
|
|
|
else if ( number >= 0 )
|
|
{
|
|
int nextnumber = mHalfEdge[number]->getNext();
|
|
int nextnextnumber = mHalfEdge[mHalfEdge[number]->getNext()]->getNext();
|
|
|
|
//insert 6 new HalfEdges for the connections to the vertices of the triangle
|
|
unsigned int edge1 = insertEdge( -10, nextnumber, mHalfEdge[number]->getPoint(), false, false );
|
|
unsigned int edge2 = insertEdge(( int )edge1, -10, mPointVector.count() - 1, false, false );
|
|
unsigned int edge3 = insertEdge( -10, nextnextnumber, mHalfEdge[nextnumber]->getPoint(), false, false );
|
|
unsigned int edge4 = insertEdge(( int )edge3, ( int )edge1, mPointVector.count() - 1, false, false );
|
|
unsigned int edge5 = insertEdge( -10, number, mHalfEdge[nextnextnumber]->getPoint(), false, false );
|
|
unsigned int edge6 = insertEdge(( int )edge5, ( int )edge3, mPointVector.count() - 1, false, false );
|
|
|
|
|
|
mHalfEdge.at( edge1 )->setDual(( int )edge2 );
|
|
mHalfEdge.at( edge2 )->setNext(( int )edge5 );
|
|
mHalfEdge.at( edge3 )->setDual(( int )edge4 );
|
|
mHalfEdge.at( edge5 )->setDual(( int )edge6 );
|
|
mHalfEdge.at( number )->setNext(( int )edge2 );
|
|
mHalfEdge.at( nextnumber )->setNext(( int )edge4 );
|
|
mHalfEdge.at( nextnextnumber )->setNext(( int )edge6 );
|
|
|
|
//check, if there are swaps necessary
|
|
checkSwap( number, 0 );
|
|
checkSwap( nextnumber, 0 );
|
|
checkSwap( nextnextnumber, 0 );
|
|
}
|
|
|
|
//the point is exactly on an existing edge (the number of the edge is stored in the variable 'mEdgeWithPoint'---------------
|
|
else if ( number == -20 )
|
|
{
|
|
int edgea = mEdgeWithPoint;
|
|
int edgeb = mHalfEdge[mEdgeWithPoint]->getDual();
|
|
int edgec = mHalfEdge[edgea]->getNext();
|
|
int edged = mHalfEdge[edgec]->getNext();
|
|
int edgee = mHalfEdge[edgeb]->getNext();
|
|
int edgef = mHalfEdge[edgee]->getNext();
|
|
|
|
//insert the six new edges
|
|
int nedge1 = insertEdge( -10, mHalfEdge[edgea]->getNext(), mHalfEdge[edgea]->getPoint(), false, false );
|
|
int nedge2 = insertEdge( nedge1, -10, mPointVector.count() - 1, false, false );
|
|
int nedge3 = insertEdge( -10, edged, mHalfEdge[edgec]->getPoint(), false, false );
|
|
int nedge4 = insertEdge( nedge3, nedge1, mPointVector.count() - 1, false, false );
|
|
int nedge5 = insertEdge( -10, edgef, mHalfEdge[edgee]->getPoint(), false, false );
|
|
int nedge6 = insertEdge( nedge5, edgeb, mPointVector.count() - 1, false, false );
|
|
|
|
//adjust the triangular structure
|
|
mHalfEdge[nedge1]->setDual( nedge2 );
|
|
mHalfEdge[nedge2]->setNext( nedge5 );
|
|
mHalfEdge[nedge3]->setDual( nedge4 );
|
|
mHalfEdge[nedge5]->setDual( nedge6 );
|
|
mHalfEdge[edgea]->setPoint( mPointVector.count() - 1 );
|
|
mHalfEdge[edgea]->setNext( nedge3 );
|
|
mHalfEdge[edgec]->setNext( nedge4 );
|
|
mHalfEdge[edgee]->setNext( nedge6 );
|
|
mHalfEdge[edgef]->setNext( nedge2 );
|
|
|
|
//swap edges if necessary
|
|
checkSwap( edgec, 0 );
|
|
checkSwap( edged, 0 );
|
|
checkSwap( edgee, 0 );
|
|
checkSwap( edgef, 0 );
|
|
}
|
|
|
|
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];
|
|
mPointVector.remove( mPointVector.count() - 1 );
|
|
delete p;
|
|
return -100;
|
|
}
|
|
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() ) );
|
|
|
|
mPointVector.remove( mPointVector.count() - 1 );
|
|
delete newPoint;
|
|
return mTwiceInsPoint;
|
|
}
|
|
}
|
|
|
|
return ( mPointVector.count() - 1 );
|
|
}
|
|
else
|
|
{
|
|
QgsDebugMsg( "warning: null pointer" );
|
|
return -100;
|
|
}
|
|
}
|
|
|
|
int DualEdgeTriangulation::baseEdgeOfPoint( int point )
|
|
{
|
|
unsigned int actedge = mEdgeInside;//starting edge
|
|
|
|
if ( mPointVector.count() < 4 || point == -1 )//at the beginning, mEdgeInside is not defined yet
|
|
{
|
|
//first find pointingedge(an edge pointing to p1)
|
|
for ( int i = 0; i < mHalfEdge.count(); i++ )
|
|
{
|
|
if ( mHalfEdge[i]->getPoint() == point )//we found it
|
|
{
|
|
return i;
|
|
}
|
|
}
|
|
}
|
|
|
|
int control = 0;
|
|
|
|
while ( true )//otherwise, start the search
|
|
{
|
|
control += 1;
|
|
if ( control > 1000000 )
|
|
{
|
|
// QgsDebugMsg( "warning, endless loop" );
|
|
|
|
//use the secure and slow method
|
|
//qWarning( "******************warning, using the slow method in baseEdgeOfPoint****************************************" );
|
|
for ( int i = 0; i < mHalfEdge.count(); i++ )
|
|
{
|
|
if ( mHalfEdge[i]->getPoint() == point && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 )//we found it
|
|
{
|
|
return i;
|
|
}
|
|
}
|
|
}
|
|
|
|
int frompoint = mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint();
|
|
int topoint = mHalfEdge[actedge]->getPoint();
|
|
|
|
if ( frompoint == -1 || topoint == -1 )//this would cause a crash. Therefore we use the slow method in this case
|
|
{
|
|
for ( int i = 0; i < mHalfEdge.count(); i++ )
|
|
{
|
|
if ( mHalfEdge[i]->getPoint() == point && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 )//we found it
|
|
{
|
|
mEdgeInside = i;
|
|
return i;
|
|
}
|
|
}
|
|
}
|
|
|
|
double leftofnumber = MathUtils::leftOf( mPointVector[point], mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] );
|
|
|
|
|
|
if ( mHalfEdge[actedge]->getPoint() == point && mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint() != -1 )//we found the edge
|
|
{
|
|
mEdgeInside = actedge;
|
|
return actedge;
|
|
}
|
|
|
|
else if ( leftofnumber <= 0.0 )
|
|
{
|
|
actedge = mHalfEdge[actedge]->getNext();
|
|
}
|
|
|
|
else
|
|
{
|
|
actedge = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getDual()]->getNext()]->getNext()]->getDual();
|
|
}
|
|
}
|
|
}
|
|
|
|
int DualEdgeTriangulation::baseEdgeOfTriangle( Point3D* point )
|
|
{
|
|
unsigned int actedge = mEdgeInside;//start with an edge which does not point to the virtual point (usualy number 3)
|
|
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
|
|
int numinstabs = 0;//number of suspect left-of-tests due to 'leftOfTresh'
|
|
int runs = 0;//counter for the number of iterations in the loop to prevent an endless loop
|
|
int firstendp = 0, secendp = 0, thendp = 0, fouendp = 0; //four numbers of endpoints in cases when two left-of-test are 0
|
|
|
|
while ( true )
|
|
{
|
|
if ( runs > nBaseOfRuns )//prevents endless loops
|
|
{
|
|
// QgsDebugMsg("warning, probable endless loop detected");
|
|
return -100;
|
|
}
|
|
|
|
double leftofvalue = MathUtils::leftOf( point, mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] );
|
|
|
|
if ( leftofvalue < ( -leftOfTresh ) )//point is on the left side
|
|
{
|
|
counter += 1;
|
|
if ( counter == 3 )//three successful passes means that we have found the triangle
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
else if ( leftofvalue == 0 )//point is exactly in the line of the edge
|
|
{
|
|
if ( nulls == 0 )
|
|
{
|
|
//store the numbers of the two endpoints of the line
|
|
firstendp = mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint();
|
|
secendp = mHalfEdge[actedge]->getPoint();
|
|
}
|
|
else if ( nulls == 1 )
|
|
{
|
|
//store the numbers of the two endpoints of the line
|
|
thendp = mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint();
|
|
fouendp = mHalfEdge[actedge]->getPoint();
|
|
}
|
|
counter += 1;
|
|
mEdgeWithPoint = actedge;
|
|
nulls += 1;
|
|
if ( counter == 3 )//three successful passes means that we have found the triangle
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
else if ( leftofvalue < leftOfTresh )//numerical problems
|
|
{
|
|
counter += 1;
|
|
numinstabs += 1;
|
|
if ( counter == 3 )//three successful passes means that we have found the triangle
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
else//point is on the right side
|
|
{
|
|
actedge = mHalfEdge[actedge]->getDual();
|
|
counter = 1;
|
|
nulls = 0;
|
|
numinstabs = 0;
|
|
}
|
|
|
|
actedge = mHalfEdge[actedge]->getNext();
|
|
if ( mHalfEdge[actedge]->getPoint() == -1 )//the half edge points to the virtual point
|
|
{
|
|
if ( nulls == 1 )//point is exactly on the convex hull
|
|
{
|
|
return -20;
|
|
}
|
|
mEdgeOutside = ( unsigned int )mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext();
|
|
mEdgeInside = mHalfEdge[mHalfEdge[mEdgeOutside]->getDual()]->getNext();
|
|
return -10;//the point is outside the convex hull
|
|
}
|
|
runs++;
|
|
}
|
|
|
|
if ( numinstabs > 0 )//we hit an existing point or a numerical instability occurred
|
|
{
|
|
// QgsDebugMsg("numerical instability occurred");
|
|
mUnstableEdge = actedge;
|
|
return -5;
|
|
}
|
|
|
|
if ( nulls == 2 )
|
|
{
|
|
//find out the number of the point, which has already been inserted
|
|
if ( firstendp == thendp || firstendp == fouendp )
|
|
{
|
|
//firstendp is the number of the point which has been inserted twice
|
|
mTwiceInsPoint = firstendp;
|
|
// QgsDebugMsg(QString("point nr %1 already inserted").arg(firstendp));
|
|
}
|
|
else if ( secendp == thendp || secendp == fouendp )
|
|
{
|
|
//secendp is the number of the point which has been inserted twice
|
|
mTwiceInsPoint = secendp;
|
|
// QgsDebugMsg(QString("point nr %1 already inserted").arg(secendp));
|
|
}
|
|
|
|
return -25;//return the code for a point that is already contained in the triangulation
|
|
}
|
|
|
|
if ( nulls == 1 )//point is on an existing edge
|
|
{
|
|
return -20;
|
|
}
|
|
|
|
mEdgeInside = actedge;
|
|
|
|
int nr1, nr2, nr3;
|
|
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();
|
|
|
|
//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
|
|
{
|
|
return actedge;
|
|
}
|
|
else if ( x2 < x1 && x2 < x3 )
|
|
{
|
|
return mHalfEdge[actedge]->getNext();
|
|
}
|
|
else if ( x3 < x1 && x3 < x2 )
|
|
{
|
|
return mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext();
|
|
}
|
|
//in case two x-coordinates are the same, the edge pointing to the point with the lower y-coordinate is returned
|
|
else if ( x1 == x2 )
|
|
{
|
|
if ( y1 < y2 )
|
|
{
|
|
return actedge;
|
|
}
|
|
else if ( y2 < y1 )
|
|
{
|
|
return mHalfEdge[actedge]->getNext();
|
|
}
|
|
}
|
|
else if ( x2 == x3 )
|
|
{
|
|
if ( y2 < y3 )
|
|
{
|
|
return mHalfEdge[actedge]->getNext();
|
|
}
|
|
else if ( y3 < y2 )
|
|
{
|
|
return mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext();
|
|
}
|
|
}
|
|
else if ( x1 == x3 )
|
|
{
|
|
if ( y1 < y3 )
|
|
{
|
|
return actedge;
|
|
}
|
|
else if ( y3 < y1 )
|
|
{
|
|
return mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext();
|
|
}
|
|
}
|
|
return -100;//this means a bug happened
|
|
}
|
|
|
|
bool DualEdgeTriangulation::calcNormal( double x, double y, Vector3D* result )
|
|
{
|
|
if ( result && mTriangleInterpolator )
|
|
{
|
|
if ( !mTriangleInterpolator->calcNormVec( x, y, result ) )
|
|
{return false;}
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
QgsDebugMsg( "warning, null pointer" );
|
|
return false;
|
|
}
|
|
}
|
|
|
|
bool DualEdgeTriangulation::calcPoint( double x, double y, Point3D* result )
|
|
{
|
|
if ( result && mTriangleInterpolator )
|
|
{
|
|
if ( !mTriangleInterpolator->calcPoint( x, y, result ) )
|
|
{return false;}
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
QgsDebugMsg( "warning, null pointer" );
|
|
return false;
|
|
}
|
|
}
|
|
|
|
bool DualEdgeTriangulation::checkSwap( unsigned int edge, unsigned int recursiveDeep )
|
|
{
|
|
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()];
|
|
if ( MathUtils::inCircle( ptd, pta, ptb, ptc ) && recursiveDeep < 100 )//empty circle criterion violated
|
|
{
|
|
doSwap( edge, recursiveDeep );//swap the edge (recursive)
|
|
return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
void DualEdgeTriangulation::doOnlySwap( unsigned int edge )
|
|
{
|
|
unsigned int edge1 = edge;
|
|
unsigned int edge2 = mHalfEdge[edge]->getDual();
|
|
unsigned int edge3 = mHalfEdge[edge]->getNext();
|
|
unsigned int edge4 = mHalfEdge[mHalfEdge[edge]->getNext()]->getNext();
|
|
unsigned int edge5 = mHalfEdge[mHalfEdge[edge]->getDual()]->getNext();
|
|
unsigned int edge6 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getNext();
|
|
mHalfEdge[edge1]->setNext( edge4 );//set the necessary nexts
|
|
mHalfEdge[edge2]->setNext( edge6 );
|
|
mHalfEdge[edge3]->setNext( edge2 );
|
|
mHalfEdge[edge4]->setNext( edge5 );
|
|
mHalfEdge[edge5]->setNext( edge1 );
|
|
mHalfEdge[edge6]->setNext( edge3 );
|
|
mHalfEdge[edge1]->setPoint( mHalfEdge[edge3]->getPoint() );//change the points to which edge1 and edge2 point
|
|
mHalfEdge[edge2]->setPoint( mHalfEdge[edge5]->getPoint() );
|
|
}
|
|
|
|
void DualEdgeTriangulation::doSwap( unsigned int edge, unsigned int recursiveDeep )
|
|
{
|
|
unsigned int edge1 = edge;
|
|
unsigned int edge2 = mHalfEdge[edge]->getDual();
|
|
unsigned int edge3 = mHalfEdge[edge]->getNext();
|
|
unsigned int edge4 = mHalfEdge[mHalfEdge[edge]->getNext()]->getNext();
|
|
unsigned int edge5 = mHalfEdge[mHalfEdge[edge]->getDual()]->getNext();
|
|
unsigned int edge6 = mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getNext();
|
|
mHalfEdge[edge1]->setNext( edge4 );//set the necessary nexts
|
|
mHalfEdge[edge2]->setNext( edge6 );
|
|
mHalfEdge[edge3]->setNext( edge2 );
|
|
mHalfEdge[edge4]->setNext( edge5 );
|
|
mHalfEdge[edge5]->setNext( edge1 );
|
|
mHalfEdge[edge6]->setNext( edge3 );
|
|
mHalfEdge[edge1]->setPoint( mHalfEdge[edge3]->getPoint() );//change the points to which edge1 and edge2 point
|
|
mHalfEdge[edge2]->setPoint( mHalfEdge[edge5]->getPoint() );
|
|
recursiveDeep++;
|
|
checkSwap( edge3, recursiveDeep );
|
|
checkSwap( edge6, recursiveDeep );
|
|
checkSwap( edge4, recursiveDeep );
|
|
checkSwap( edge5, recursiveDeep );
|
|
}
|
|
|
|
#if 0
|
|
void DualEdgeTriangulation::draw( QPainter* p, double xlowleft, double ylowleft, double xupright, double yupright, double width, double height ) const
|
|
{
|
|
//if mPointVector is empty, there is nothing to do
|
|
if ( mPointVector.isEmpty() )
|
|
{
|
|
return;
|
|
}
|
|
|
|
p->setPen( mEdgeColor );
|
|
|
|
bool* control = new bool[mHalfEdge.count()];//controllarray that no edge is painted twice
|
|
bool* control2 = new bool[mHalfEdge.count()];//controllarray for the flat triangles
|
|
|
|
for ( unsigned int i = 0; i <= mHalfEdge.count() - 1; i++ )
|
|
{
|
|
control[i] = false;
|
|
control2[i] = false;
|
|
}
|
|
|
|
if ((( xupright - xlowleft ) / width ) > (( yupright - ylowleft ) / height ) )
|
|
{
|
|
double lowerborder = -( height * ( xupright - xlowleft ) / width - yupright );//real world coordinates of the lower widget border. This is useful to know because of the HalfEdge bounding box test
|
|
for ( unsigned int i = 0; i < mHalfEdge.count() - 1; i++ )
|
|
{
|
|
if ( mHalfEdge[i]->getPoint() == -1 || mHalfEdge[mHalfEdge[i]->getDual()]->getPoint() == -1 )
|
|
{continue;}
|
|
|
|
//check, if the edge belongs to a flat triangle, remove this later
|
|
if ( !control2[i] )
|
|
{
|
|
double p1, p2, p3;
|
|
if ( mHalfEdge[i]->getPoint() != -1 && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 && mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint() != -1 )
|
|
{
|
|
p1 = mPointVector[mHalfEdge[i]->getPoint()]->getZ();
|
|
p2 = mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getZ();
|
|
p3 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getZ();
|
|
if ( p1 == p2 && p2 == p3 && halfEdgeBBoxTest( i, xlowleft, lowerborder, xupright, yupright ) && halfEdgeBBoxTest( mHalfEdge[i]->getNext(), xlowleft, lowerborder, xupright, yupright ) && halfEdgeBBoxTest( mHalfEdge[mHalfEdge[i]->getNext()]->getNext(), xlowleft, lowerborder, xupright, yupright ) )//draw the triangle
|
|
{
|
|
QPointArray pa( 3 );
|
|
pa.setPoint( 0, ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
|
|
pa.setPoint( 1, ( mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
|
|
pa.setPoint( 2, ( mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
|
|
QColor c( 255, 0, 0 );
|
|
p->setBrush( c );
|
|
p->drawPolygon( pa );
|
|
}
|
|
}
|
|
|
|
control2[i] = true;
|
|
control2[mHalfEdge[i]->getNext()] = true;
|
|
control2[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()] = true;
|
|
}//end of the section, which has to be removed later
|
|
|
|
if ( control[i] )//check, if edge has already been drawn
|
|
{continue;}
|
|
|
|
//draw the edge;
|
|
if ( halfEdgeBBoxTest( i, xlowleft, lowerborder, xupright, yupright ) )//only draw the halfedge if its bounding box intersects the painted area
|
|
{
|
|
if ( mHalfEdge[i]->getBreak() )//change the color it the edge is a breakline
|
|
{
|
|
p->setPen( mBreakEdgeColor );
|
|
}
|
|
else if ( mHalfEdge[i]->getForced() )//change the color if the edge is forced
|
|
{
|
|
p->setPen( mForcedEdgeColor );
|
|
}
|
|
|
|
|
|
p->drawLine(( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( xupright - xlowleft )*width, ( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getX() - xlowleft ) / ( xupright - xlowleft )*width, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getY() ) / ( xupright - xlowleft )*width );
|
|
|
|
if ( mHalfEdge[i]->getForced() )
|
|
{
|
|
p->setPen( mEdgeColor );
|
|
}
|
|
|
|
|
|
}
|
|
control[i] = true;
|
|
control[mHalfEdge[i]->getDual()] = true;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
double rightborder = width * ( yupright - ylowleft ) / height + xlowleft;//real world coordinates of the right widget border. This is useful to know because of the HalfEdge bounding box test
|
|
for ( unsigned int i = 0; i < mHalfEdge.count() - 1; i++ )
|
|
{
|
|
if ( mHalfEdge[i]->getPoint() == -1 || mHalfEdge[mHalfEdge[i]->getDual()]->getPoint() == -1 )
|
|
{continue;}
|
|
|
|
//check, if the edge belongs to a flat triangle, remove this section later
|
|
if ( !control2[i] )
|
|
{
|
|
double p1, p2, p3;
|
|
if ( mHalfEdge[i]->getPoint() != -1 && mHalfEdge[mHalfEdge[i]->getNext()]->getPoint() != -1 && mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint() != -1 )
|
|
{
|
|
p1 = mPointVector[mHalfEdge[i]->getPoint()]->getZ();
|
|
p2 = mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getZ();
|
|
p3 = mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getZ();
|
|
if ( p1 == p2 && p2 == p3 && halfEdgeBBoxTest( i, xlowleft, ylowleft, rightborder, yupright ) && halfEdgeBBoxTest( mHalfEdge[i]->getNext(), xlowleft, ylowleft, rightborder, yupright ) && halfEdgeBBoxTest( mHalfEdge[mHalfEdge[i]->getNext()]->getNext(), xlowleft, ylowleft, rightborder, yupright ) )//draw the triangle
|
|
{
|
|
QPointArray pa( 3 );
|
|
pa.setPoint( 0, ( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
|
|
pa.setPoint( 1, ( mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getNext()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
|
|
pa.setPoint( 2, ( mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
|
|
QColor c( 255, 0, 0 );
|
|
p->setBrush( c );
|
|
p->drawPolygon( pa );
|
|
}
|
|
}
|
|
|
|
control2[i] = true;
|
|
control2[mHalfEdge[i]->getNext()] = true;
|
|
control2[mHalfEdge[mHalfEdge[i]->getNext()]->getNext()] = true;
|
|
}//end of the section, which has to be removed later
|
|
|
|
|
|
if ( control[i] )//check, if edge has already been drawn
|
|
{continue;}
|
|
|
|
//draw the edge
|
|
if ( halfEdgeBBoxTest( i, xlowleft, ylowleft, rightborder, yupright ) )//only draw the edge if its bounding box intersects with the painted area
|
|
{
|
|
if ( mHalfEdge[i]->getBreak() )//change the color if the edge is a breakline
|
|
{
|
|
p->setPen( mBreakEdgeColor );
|
|
}
|
|
else if ( mHalfEdge[i]->getForced() )//change the color if the edge is forced
|
|
{
|
|
p->setPen( mForcedEdgeColor );
|
|
}
|
|
|
|
p->drawLine(( mPointVector[mHalfEdge[i]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[i]->getPoint()]->getY() ) / ( yupright - ylowleft )*height, ( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getX() - xlowleft ) / ( yupright - ylowleft )*height, ( yupright - mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()]->getY() ) / ( yupright - ylowleft )*height );
|
|
|
|
if ( mHalfEdge[i]->getForced() )
|
|
{
|
|
p->setPen( mEdgeColor );
|
|
}
|
|
|
|
}
|
|
control[i] = true;
|
|
control[mHalfEdge[i]->getDual()] = true;
|
|
}
|
|
}
|
|
|
|
delete[] control;
|
|
delete[] control2;
|
|
}
|
|
#endif
|
|
|
|
int DualEdgeTriangulation::getOppositePoint( int p1, int p2 )
|
|
{
|
|
|
|
//first find a half edge which points to p2
|
|
int firstedge = baseEdgeOfPoint( p2 );
|
|
|
|
//then find the edge which comes from p1 or print an error message if there isn't any
|
|
int theedge = -10;
|
|
int nextnextedge = firstedge;
|
|
int edge, nextedge;
|
|
do
|
|
{
|
|
edge = mHalfEdge[nextnextedge]->getDual();
|
|
if ( mHalfEdge[edge]->getPoint() == p1 )
|
|
{
|
|
theedge = nextnextedge;
|
|
break;
|
|
}//we found the edge
|
|
nextedge = mHalfEdge[edge]->getNext();
|
|
nextnextedge = mHalfEdge[nextedge]->getNext();
|
|
}
|
|
while ( nextnextedge != firstedge );
|
|
|
|
if ( theedge == -10 )//there is no edge between p1 and p2
|
|
{
|
|
QgsDebugMsg( QString( "warning, error: the points are: %1 and %2" ).arg( p1 ).arg( p2 ) );
|
|
return -10;
|
|
}
|
|
|
|
//finally find the opposite point
|
|
return mHalfEdge[mHalfEdge[mHalfEdge[theedge]->getDual()]->getNext()]->getPoint();
|
|
|
|
}
|
|
|
|
QList<int>* DualEdgeTriangulation::getSurroundingTriangles( int pointno )
|
|
{
|
|
int firstedge = baseEdgeOfPoint( pointno );
|
|
|
|
if ( firstedge == -1 )//an error occurred
|
|
{
|
|
return nullptr;
|
|
}
|
|
|
|
QList<int>* vlist = new QList<int>();//create the value list on the heap
|
|
|
|
int actedge = firstedge;
|
|
int edge, nextedge, nextnextedge;
|
|
do
|
|
{
|
|
edge = mHalfEdge[actedge]->getDual();
|
|
vlist->append( mHalfEdge[edge]->getPoint() );//add the number of the endpoint of the first edge to the value list
|
|
nextedge = mHalfEdge[edge]->getNext();
|
|
vlist->append( mHalfEdge[nextedge]->getPoint() );//add the number of the endpoint of the second edge to the value list
|
|
nextnextedge = mHalfEdge[nextedge]->getNext();
|
|
vlist->append( mHalfEdge[nextnextedge]->getPoint() );//add the number of endpoint of the third edge to the value list
|
|
if ( mHalfEdge[nextnextedge]->getBreak() )//add, whether the third edge is a breakline or not
|
|
{
|
|
vlist->append( -10 );
|
|
}
|
|
else
|
|
{
|
|
vlist->append( -20 );
|
|
}
|
|
actedge = nextnextedge;
|
|
}
|
|
while ( nextnextedge != firstedge );
|
|
|
|
return vlist;
|
|
|
|
}
|
|
|
|
bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D* p1, int* n1, Point3D* p2, int* n2, Point3D* p3, int* n3 )
|
|
{
|
|
if ( mPointVector.size() < 3 )
|
|
{
|
|
return false;
|
|
}
|
|
|
|
if ( p1 && p2 && p3 )
|
|
{
|
|
Point3D point( x, y, 0 );
|
|
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
|
|
{
|
|
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() );
|
|
( *n1 ) = ptnr1;
|
|
( *n2 ) = ptnr2;
|
|
( *n3 ) = ptnr3;
|
|
return true;
|
|
}
|
|
else if ( edge == -20 )//the point is exactly on an edge
|
|
{
|
|
int ptnr1 = mHalfEdge[mEdgeWithPoint]->getPoint();
|
|
int ptnr2 = mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getPoint();
|
|
int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getNext()]->getPoint();
|
|
if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
|
|
{
|
|
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() );
|
|
( *n1 ) = ptnr1;
|
|
( *n2 ) = ptnr2;
|
|
( *n3 ) = ptnr3;
|
|
return true;
|
|
}
|
|
else if ( edge == -25 )//x and y are the coordinates of an existing point
|
|
{
|
|
int edge1 = baseEdgeOfPoint( mTwiceInsPoint );
|
|
int edge2 = mHalfEdge[edge1]->getNext();
|
|
int edge3 = mHalfEdge[edge2]->getNext();
|
|
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() );
|
|
( *n1 ) = ptnr1;
|
|
( *n2 ) = ptnr2;
|
|
( *n3 ) = ptnr3;
|
|
return true;
|
|
}
|
|
else if ( edge == -5 )//numerical problems in 'baseEdgeOfTriangle'
|
|
{
|
|
int ptnr1 = mHalfEdge[mUnstableEdge]->getPoint();
|
|
int ptnr2 = mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getPoint();
|
|
int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getNext()]->getPoint();
|
|
if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
|
|
{
|
|
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() );
|
|
( *n1 ) = ptnr1;
|
|
( *n2 ) = ptnr2;
|
|
( *n3 ) = ptnr3;
|
|
return true;
|
|
}
|
|
else//problems
|
|
{
|
|
QgsDebugMsg( QString( "problem: the edge is: %1" ).arg( edge ) );
|
|
return false;
|
|
}
|
|
}
|
|
|
|
else
|
|
{
|
|
QgsDebugMsg( "warning, null pointer" );
|
|
return false;
|
|
}
|
|
}
|
|
|
|
bool DualEdgeTriangulation::getTriangle( double x, double y, Point3D* p1, Point3D* p2, Point3D* p3 )
|
|
{
|
|
if ( mPointVector.size() < 3 )
|
|
{
|
|
return false;
|
|
}
|
|
|
|
if ( p1 && p2 && p3 )
|
|
{
|
|
Point3D point( x, y, 0 );
|
|
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
|
|
{
|
|
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() );
|
|
return true;
|
|
}
|
|
else if ( edge == -20 )//the point is exactly on an edge
|
|
{
|
|
int ptnr1 = mHalfEdge[mEdgeWithPoint]->getPoint();
|
|
int ptnr2 = mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getPoint();
|
|
int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mEdgeWithPoint]->getNext()]->getNext()]->getPoint();
|
|
if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
|
|
{
|
|
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() );
|
|
return true;
|
|
}
|
|
else if ( edge == -25 )//x and y are the coordinates of an existing point
|
|
{
|
|
int edge1 = baseEdgeOfPoint( mTwiceInsPoint );
|
|
int edge2 = mHalfEdge[edge1]->getNext();
|
|
int edge3 = mHalfEdge[edge2]->getNext();
|
|
int ptnr1 = mHalfEdge[edge1]->getPoint();
|
|
int ptnr2 = mHalfEdge[edge2]->getPoint();
|
|
int ptnr3 = mHalfEdge[edge3]->getPoint();
|
|
if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
|
|
{
|
|
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() );
|
|
return true;
|
|
}
|
|
else if ( edge == -5 )//numerical problems in 'baseEdgeOfTriangle'
|
|
{
|
|
int ptnr1 = mHalfEdge[mUnstableEdge]->getPoint();
|
|
int ptnr2 = mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getPoint();
|
|
int ptnr3 = mHalfEdge[mHalfEdge[mHalfEdge[mUnstableEdge]->getNext()]->getNext()]->getPoint();
|
|
if ( ptnr1 == -1 || ptnr2 == -1 || ptnr3 == -1 )
|
|
{
|
|
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() );
|
|
return true;
|
|
}
|
|
else//problems
|
|
{
|
|
QgsDebugMsg( QString( "problems: the edge is: %1" ).arg( edge ) );
|
|
return false;
|
|
}
|
|
}
|
|
|
|
else
|
|
{
|
|
QgsDebugMsg( "warning, null pointer" );
|
|
return false;
|
|
}
|
|
}
|
|
|
|
unsigned int DualEdgeTriangulation::insertEdge( int dual, int next, int point, bool mbreak, bool forced )
|
|
{
|
|
HalfEdge* edge = new HalfEdge( dual, next, point, mbreak, forced );
|
|
mHalfEdge.append( edge );
|
|
return mHalfEdge.count() - 1;
|
|
|
|
}
|
|
|
|
int DualEdgeTriangulation::insertForcedSegment( int p1, int p2, bool breakline )
|
|
{
|
|
if ( p1 == p2 )
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
//list with (half of) the crossed edges
|
|
QList<int> crossedEdges;
|
|
|
|
//an edge pointing to p1
|
|
int pointingedge = 0;
|
|
|
|
pointingedge = baseEdgeOfPoint( p1 );
|
|
|
|
if ( pointingedge == -1 )
|
|
{
|
|
return -100;//return an error code
|
|
}
|
|
|
|
//go around p1 and find out, if the segment already exists and if not, which is the first cutted edge
|
|
int actedge = mHalfEdge[pointingedge]->getDual();
|
|
//number to prevent endless loops
|
|
int control = 0;
|
|
|
|
while ( true )//if it's an endless loop, something went wrong
|
|
{
|
|
control += 1;
|
|
if ( control > 17000 )
|
|
{
|
|
QgsDebugMsg( "warning, endless loop" );
|
|
return -100;//return an error code
|
|
}
|
|
|
|
if ( mHalfEdge[actedge]->getPoint() == -1 )//actedge points to the virtual point
|
|
{
|
|
actedge = mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext()]->getDual();
|
|
continue;
|
|
}
|
|
|
|
//test, if actedge is already the forced edge
|
|
if ( mHalfEdge[actedge]->getPoint() == p2 )
|
|
{
|
|
mHalfEdge[actedge]->setForced( true );
|
|
mHalfEdge[actedge]->setBreak( breakline );
|
|
mHalfEdge[mHalfEdge[actedge]->getDual()]->setForced( true );
|
|
mHalfEdge[mHalfEdge[actedge]->getDual()]->setBreak( breakline );
|
|
return actedge;
|
|
}
|
|
|
|
//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 ) )
|
|
{
|
|
//mark actedge and Dual(actedge) as forced, reset p1 and start the method from the beginning
|
|
mHalfEdge[actedge]->setForced( true );
|
|
mHalfEdge[actedge]->setBreak( breakline );
|
|
mHalfEdge[mHalfEdge[actedge]->getDual()]->setForced( true );
|
|
mHalfEdge[mHalfEdge[actedge]->getDual()]->setBreak( breakline );
|
|
int a = insertForcedSegment( mHalfEdge[actedge]->getPoint(), p2, breakline );
|
|
return a;
|
|
}
|
|
|
|
//test, if the forced segment intersects Next(actedge)
|
|
if ( mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint() == -1 )//intersection with line to the virtual point makes no sense
|
|
{
|
|
actedge = mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext()]->getDual();
|
|
continue;
|
|
}
|
|
else if ( MathUtils::lineIntersection( mPointVector[p1], mPointVector[p2], mPointVector[mHalfEdge[mHalfEdge[actedge]->getNext()]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getDual()]->getPoint()] ) )
|
|
{
|
|
if ( mHalfEdge[mHalfEdge[actedge]->getNext()]->getForced() && mForcedCrossBehaviour == Triangulation::SnappingType_VERTICE )//if the crossed edge is a forced edge, we have to snap the forced line to the next node
|
|
{
|
|
Point3D 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() ) );
|
|
if ( dista <= distb )
|
|
{
|
|
insertForcedSegment( p1, p3, breakline );
|
|
int e = insertForcedSegment( p3, p2, breakline );
|
|
return e;
|
|
}
|
|
else if ( distb <= dista )
|
|
{
|
|
insertForcedSegment( p1, p4, breakline );
|
|
int e = insertForcedSegment( p4, p2, breakline );
|
|
return e;
|
|
}
|
|
}
|
|
else if ( mHalfEdge[mHalfEdge[actedge]->getNext()]->getForced() && mForcedCrossBehaviour == Triangulation::INSERT_VERTICE )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
|
|
{
|
|
Point3D 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() ) );
|
|
float frac = distpart / disttot;
|
|
|
|
if ( frac == 0 || frac == 1 )//just in case MathUtils::lineIntersection does not work as expected
|
|
{
|
|
if ( frac == 0 )
|
|
{
|
|
//mark actedge and Dual(actedge) as forced, reset p1 and start the method from the beginning
|
|
mHalfEdge[actedge]->setForced( true );
|
|
mHalfEdge[actedge]->setBreak( breakline );
|
|
mHalfEdge[mHalfEdge[actedge]->getDual()]->setForced( true );
|
|
mHalfEdge[mHalfEdge[actedge]->getDual()]->setBreak( breakline );
|
|
int a = insertForcedSegment( p4, p2, breakline );
|
|
return a;
|
|
}
|
|
else if ( frac == 1 )
|
|
{
|
|
//mark actedge and Dual(actedge) as forced, reset p1 and start the method from the beginning
|
|
mHalfEdge[actedge]->setForced( true );
|
|
mHalfEdge[actedge]->setBreak( breakline );
|
|
mHalfEdge[mHalfEdge[actedge]->getDual()]->setForced( true );
|
|
mHalfEdge[mHalfEdge[actedge]->getDual()]->setBreak( breakline );
|
|
if ( p3 != p2 )
|
|
{
|
|
int a = insertForcedSegment( p3, p2, breakline );
|
|
return a;
|
|
}
|
|
else
|
|
{
|
|
return actedge;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
|
|
else
|
|
{
|
|
int newpoint = splitHalfEdge( mHalfEdge[actedge]->getNext(), frac );
|
|
insertForcedSegment( p1, newpoint, breakline );
|
|
int e = insertForcedSegment( newpoint, p2, breakline );
|
|
return e;
|
|
}
|
|
}
|
|
|
|
//add the first HalfEdge to the list of crossed edges
|
|
crossedEdges.append( mHalfEdge[actedge]->getNext() );
|
|
break;
|
|
}
|
|
actedge = mHalfEdge[mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext()]->getDual();
|
|
}
|
|
|
|
//we found the first edge, terminated the method or called the method with other points. Lets search for all the other crossed edges
|
|
|
|
while ( true )//if its an endless loop, something went wrong.
|
|
{
|
|
if ( MathUtils::lineIntersection( mPointVector[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint()], mPointVector[p1], mPointVector[p2] ) )
|
|
{
|
|
if ( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getForced() && mForcedCrossBehaviour == Triangulation::SnappingType_VERTICE )//if the crossed edge is a forced edge and mForcedCrossBehaviour is SnappingType_VERTICE, we have to snap the forced line to the next node
|
|
{
|
|
Point3D 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() ) );
|
|
if ( dista <= distb )
|
|
{
|
|
insertForcedSegment( p1, p3, breakline );
|
|
int e = insertForcedSegment( p3, p2, breakline );
|
|
return e;
|
|
}
|
|
else if ( distb <= dista )
|
|
{
|
|
insertForcedSegment( p1, p4, breakline );
|
|
int e = insertForcedSegment( p4, p2, breakline );
|
|
return e;
|
|
}
|
|
}
|
|
else if ( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getForced() && mForcedCrossBehaviour == Triangulation::INSERT_VERTICE )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
|
|
{
|
|
Point3D 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() ) );
|
|
float frac = distpart / disttot;
|
|
if ( frac == 0 || frac == 1 )
|
|
{
|
|
break;//seems that a roundoff error occurred. We found the endpoint
|
|
}
|
|
int newpoint = splitHalfEdge( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext(), frac );
|
|
insertForcedSegment( p1, newpoint, breakline );
|
|
int e = insertForcedSegment( newpoint, p2, breakline );
|
|
return e;
|
|
}
|
|
|
|
crossedEdges.append( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext() );
|
|
continue;
|
|
}
|
|
else if ( MathUtils::lineIntersection( mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getPoint()], mPointVector[p1], mPointVector[p2] ) )
|
|
{
|
|
if ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getForced() && mForcedCrossBehaviour == Triangulation::SnappingType_VERTICE )//if the crossed edge is a forced edge and mForcedCrossBehaviour is SnappingType_VERTICE, we have to snap the forced line to the next node
|
|
{
|
|
Point3D 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() ) );
|
|
if ( dista <= distb )
|
|
{
|
|
insertForcedSegment( p1, p3, breakline );
|
|
int e = insertForcedSegment( p3, p2, breakline );
|
|
return e;
|
|
}
|
|
else if ( distb <= dista )
|
|
{
|
|
insertForcedSegment( p1, p4, breakline );
|
|
int e = insertForcedSegment( p4, p2, breakline );
|
|
return e;
|
|
}
|
|
}
|
|
else if ( mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext()]->getForced() && mForcedCrossBehaviour == Triangulation::INSERT_VERTICE )//if the crossed edge is a forced edge, we have to insert a new vertice on this edge
|
|
{
|
|
Point3D 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() ) );
|
|
float frac = distpart / disttot;
|
|
if ( frac == 0 || frac == 1 )
|
|
{
|
|
break;//seems that a roundoff error occurred. We found the endpoint
|
|
}
|
|
int newpoint = splitHalfEdge( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext(), frac );
|
|
insertForcedSegment( p1, newpoint, breakline );
|
|
int e = insertForcedSegment( newpoint, p2, breakline );
|
|
return e;
|
|
}
|
|
|
|
crossedEdges.append( mHalfEdge[mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext()]->getNext() );
|
|
continue;
|
|
}
|
|
else//forced edge terminates
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
//set the flags 'forced' and 'break' to false for every edge and dualedge of 'crossEdges'
|
|
QList<int>::const_iterator iter;
|
|
for ( iter = crossedEdges.constBegin(); iter != crossedEdges.constEnd(); ++iter )
|
|
{
|
|
mHalfEdge[( *( iter ) )]->setForced( false );
|
|
mHalfEdge[( *( iter ) )]->setBreak( false );
|
|
mHalfEdge[mHalfEdge[( *( iter ) )]->getDual()]->setForced( false );
|
|
mHalfEdge[mHalfEdge[( *( iter ) )]->getDual()]->setBreak( false );
|
|
}
|
|
|
|
//crossed edges is filled, now the two polygons to be retriangulated can be build
|
|
|
|
|
|
QList<int> freelist = crossedEdges;//copy the list with the crossed edges to remove the edges already reused
|
|
|
|
//create the left polygon as a list of the numbers of the halfedges
|
|
QList<int> leftPolygon;
|
|
QList<int> rightPolygon;
|
|
|
|
//insert the forced edge and enter the corresponding halfedges as the first edges in the left and right polygons. The nexts and points are set later because of the algorithm to build two polygons from 'crossedEdges'
|
|
int firstedge = freelist.first();//edge pointing from p1 to p2
|
|
mHalfEdge[firstedge]->setForced( true );
|
|
mHalfEdge[firstedge]->setBreak( breakline );
|
|
leftPolygon.append( firstedge );
|
|
int dualfirstedge = mHalfEdge[freelist.first()]->getDual();//edge pointing from p2 to p1
|
|
mHalfEdge[dualfirstedge]->setForced( true );
|
|
mHalfEdge[dualfirstedge]->setBreak( breakline );
|
|
rightPolygon.append( dualfirstedge );
|
|
freelist.pop_front();//delete the first entry from the freelist
|
|
|
|
//finish the polygon on the left side
|
|
int actpointl = p2;
|
|
QList<int>::const_iterator leftiter; //todo: is there a better way to set an iterator to the last list element?
|
|
leftiter = crossedEdges.constEnd();
|
|
--leftiter;
|
|
while ( true )
|
|
{
|
|
int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext()]->getPoint();
|
|
if ( newpoint != actpointl )
|
|
{
|
|
//insert the edge into the leftPolygon
|
|
actpointl = newpoint;
|
|
int theedge = mHalfEdge[mHalfEdge[mHalfEdge[( *leftiter )]->getDual()]->getNext()]->getNext();
|
|
leftPolygon.append( theedge );
|
|
}
|
|
if ( leftiter == crossedEdges.constBegin() )
|
|
{break;}
|
|
--leftiter;
|
|
}
|
|
|
|
//insert the last element into leftPolygon
|
|
leftPolygon.append( mHalfEdge[crossedEdges.first()]->getNext() );
|
|
|
|
//finish the polygon on the right side
|
|
QList<int>::const_iterator rightiter;
|
|
int actpointr = p1;
|
|
for ( rightiter = crossedEdges.constBegin(); rightiter != crossedEdges.constEnd(); ++rightiter )
|
|
{
|
|
int newpoint = mHalfEdge[mHalfEdge[mHalfEdge[( *rightiter )]->getNext()]->getNext()]->getPoint();
|
|
if ( newpoint != actpointr )
|
|
{
|
|
//insert the edge into the right polygon
|
|
actpointr = newpoint;
|
|
int theedge = mHalfEdge[mHalfEdge[( *rightiter )]->getNext()]->getNext();
|
|
rightPolygon.append( theedge );
|
|
}
|
|
}
|
|
|
|
|
|
//insert the last element into rightPolygon
|
|
rightPolygon.append( mHalfEdge[mHalfEdge[crossedEdges.last()]->getDual()]->getNext() );
|
|
mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );//set 'Next' of the last edge to dualfirstedge
|
|
|
|
//set the necessary nexts of leftPolygon(exept the first)
|
|
int actedgel = leftPolygon[1];
|
|
leftiter = leftPolygon.constBegin();
|
|
leftiter += 2;
|
|
for ( ; leftiter != leftPolygon.constEnd(); ++leftiter )
|
|
{
|
|
mHalfEdge[actedgel]->setNext(( *leftiter ) );
|
|
actedgel = ( *leftiter );
|
|
}
|
|
|
|
//set all the necessary nexts of rightPolygon
|
|
int actedger = rightPolygon[1];
|
|
rightiter = rightPolygon.constBegin();
|
|
rightiter += 2;
|
|
for ( ; rightiter != rightPolygon.constEnd(); ++rightiter )
|
|
{
|
|
mHalfEdge[actedger]->setNext(( *rightiter ) );
|
|
actedger = ( *( rightiter ) );
|
|
}
|
|
|
|
|
|
//setNext and setPoint for the forced edge because this would disturb the building of 'leftpoly' and 'rightpoly' otherwise
|
|
mHalfEdge[leftPolygon.first()]->setNext(( *( ++( leftiter = leftPolygon.begin() ) ) ) );
|
|
mHalfEdge[leftPolygon.first()]->setPoint( p2 );
|
|
mHalfEdge[leftPolygon.last()]->setNext( firstedge );
|
|
mHalfEdge[rightPolygon.first()]->setNext(( *( ++( rightiter = rightPolygon.begin() ) ) ) );
|
|
mHalfEdge[rightPolygon.first()]->setPoint( p1 );
|
|
mHalfEdge[rightPolygon.last()]->setNext( dualfirstedge );
|
|
|
|
triangulatePolygon( &leftPolygon, &freelist, firstedge );
|
|
triangulatePolygon( &rightPolygon, &freelist, dualfirstedge );
|
|
|
|
//optimisation of the new edges
|
|
for ( iter = crossedEdges.begin(); iter != crossedEdges.end(); ++iter )
|
|
{
|
|
checkSwap(( *( iter ) ), 0 );
|
|
}
|
|
|
|
return leftPolygon.first();
|
|
}
|
|
|
|
void DualEdgeTriangulation::setForcedCrossBehaviour( Triangulation::forcedCrossBehaviour b )
|
|
{
|
|
mForcedCrossBehaviour = b;
|
|
}
|
|
|
|
void DualEdgeTriangulation::setEdgeColor( int r, int g, int b )
|
|
{
|
|
mEdgeColor.setRgb( r, g, b );
|
|
}
|
|
|
|
void DualEdgeTriangulation::setForcedEdgeColor( int r, int g, int b )
|
|
{
|
|
mForcedEdgeColor.setRgb( r, g, b );
|
|
}
|
|
|
|
void DualEdgeTriangulation::setBreakEdgeColor( int r, int g, int b )
|
|
{
|
|
mBreakEdgeColor.setRgb( r, g, b );
|
|
}
|
|
|
|
void DualEdgeTriangulation::setTriangleInterpolator( TriangleInterpolator* interpolator )
|
|
{
|
|
mTriangleInterpolator = interpolator;
|
|
}
|
|
|
|
void DualEdgeTriangulation::eliminateHorizontalTriangles()
|
|
{
|
|
QgsDebugMsg( "am in eliminateHorizontalTriangles" );
|
|
double minangle = 0;//minimum angle for swapped triangles. If triangles generated by a swap would have a minimum angle (in degrees) below that value, the swap will not be done.
|
|
|
|
while ( true )
|
|
{
|
|
bool swapped = false;//flag which allows exiting the loop
|
|
bool* control = new bool[mHalfEdge.count()];//controlarray
|
|
|
|
for ( int i = 0; i <= mHalfEdge.count() - 1; i++ )
|
|
{
|
|
control[i] = false;
|
|
}
|
|
|
|
|
|
for ( int i = 0; i <= mHalfEdge.count() - 1; i++ )
|
|
{
|
|
if ( control[i] )//edge has already been examined
|
|
{
|
|
continue;
|
|
}
|
|
|
|
int e1, e2, e3;//numbers of the three edges
|
|
e1 = i;
|
|
e2 = mHalfEdge[e1]->getNext();
|
|
e3 = mHalfEdge[e2]->getNext();
|
|
|
|
int p1, p2, p3;//numbers of the three points
|
|
p1 = mHalfEdge[e1]->getPoint();
|
|
p2 = mHalfEdge[e2]->getPoint();
|
|
p3 = mHalfEdge[e3]->getPoint();
|
|
|
|
//skip the iteration, if one point is the virtual point
|
|
if ( p1 == -1 || p2 == -1 || p3 == -1 )
|
|
{
|
|
control[e1] = true;
|
|
control[e2] = true;
|
|
control[e3] = true;
|
|
continue;
|
|
}
|
|
|
|
double el1, el2, el3;//elevations of the points
|
|
el1 = mPointVector[p1]->getZ();
|
|
el2 = mPointVector[p2]->getZ();
|
|
el3 = mPointVector[p3]->getZ();
|
|
|
|
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 )
|
|
{
|
|
doOnlySwap(( uint )e1 );
|
|
swapped = true;
|
|
}
|
|
else if ( swapPossible(( uint )e2 ) && mPointVector[mHalfEdge[mHalfEdge[mHalfEdge[e2]->getDual()]->getNext()]->getPoint()]->getZ() != 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 )
|
|
{
|
|
doOnlySwap(( uint )e3 );
|
|
swapped = true;
|
|
}
|
|
control[e1] = true;
|
|
control[e2] = true;
|
|
control[e3] = true;
|
|
continue;
|
|
}
|
|
else//no horizontal triangle, go to the next one
|
|
{
|
|
control[e1] = true;
|
|
control[e2] = true;
|
|
control[e3] = true;
|
|
continue;
|
|
}
|
|
}
|
|
if ( !swapped )
|
|
{
|
|
delete[] control;
|
|
break;
|
|
}
|
|
delete[] control;
|
|
}
|
|
|
|
QgsDebugMsg( "end of method" );
|
|
}
|
|
|
|
void DualEdgeTriangulation::ruppertRefinement()
|
|
{
|
|
//minimum angle
|
|
double mintol = 17;//refinement stops after the minimum angle reached this tolerance
|
|
|
|
//data structures
|
|
std::map<int, double> edge_angle;//search tree with the edge number as key
|
|
std::multimap<double, int> angle_edge;//multimap (map with not unique keys) with angle as key
|
|
QSet<int> dontexamine;//search tree containing the edges which do not have to be examined (because of numerical problems)
|
|
|
|
|
|
//first, go through all the forced edges and subdivide if they are encroached by a point
|
|
bool stop = false;//flag to ensure that the for-loop is repeated until no half edge is split any more
|
|
|
|
while ( !stop )
|
|
{
|
|
stop = true;
|
|
int nhalfedges = mHalfEdge.count();
|
|
|
|
for ( int i = 0; i < nhalfedges - 1; i++ )
|
|
{
|
|
int next = mHalfEdge[i]->getNext();
|
|
int nextnext = mHalfEdge[next]->getNext();
|
|
|
|
if ( mHalfEdge[next]->getPoint() != -1 && ( mHalfEdge[i]->getForced() || mHalfEdge[mHalfEdge[mHalfEdge[i]->getDual()]->getNext()]->getPoint() == -1 ) )//check for encroached points on forced segments and segments on the inner side of the convex hull, but don't consider edges on the outer side of the convex hull
|
|
{
|
|
if ( !(( mHalfEdge[next]->getForced() || edgeOnConvexHull( next ) ) || ( mHalfEdge[nextnext]->getForced() || edgeOnConvexHull( nextnext ) ) ) )//don't consider triangles where all three edges are forced edges or hull edges
|
|
{
|
|
//test for encroachment
|
|
while ( MathUtils::inDiametral( mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()], mPointVector[mHalfEdge[i]->getPoint()], mPointVector[mHalfEdge[next]->getPoint()] ) )
|
|
{
|
|
//split segment
|
|
int pointno = splitHalfEdge( i, 0.5 );
|
|
Q_UNUSED( pointno );
|
|
stop = false;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
//examine the triangulation for angles below the minimum and insert the edges into angle_edge and edge_angle, except the small angle is between forced segments or convex hull edges
|
|
double angle;//angle between edge i and the consecutive edge
|
|
int p1, p2, p3;//numbers of the triangle points
|
|
for ( int i = 0; i < mHalfEdge.count() - 1; i++ )
|
|
{
|
|
p1 = mHalfEdge[mHalfEdge[i]->getDual()]->getPoint();
|
|
p2 = mHalfEdge[i]->getPoint();
|
|
p3 = mHalfEdge[mHalfEdge[i]->getNext()]->getPoint();
|
|
|
|
if ( p1 == -1 || p2 == -1 || p3 == -1 )//don't consider triangles with the virtual point
|
|
{
|
|
continue;
|
|
}
|
|
angle = MathUtils::angle( mPointVector[p1], mPointVector[p2], mPointVector[p3], mPointVector[p2] );
|
|
|
|
bool twoforcededges;//flag to decide, if edges should be added to the maps. Do not add them if true
|
|
|
|
|
|
if (( mHalfEdge[i]->getForced() || edgeOnConvexHull( i ) ) && ( mHalfEdge[mHalfEdge[i]->getNext()]->getForced() || edgeOnConvexHull( mHalfEdge[i]->getNext() ) ) )
|
|
{
|
|
twoforcededges = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges = false;
|
|
}
|
|
|
|
if ( angle < mintol && !twoforcededges )
|
|
{
|
|
edge_angle.insert( std::make_pair( i, angle ) );
|
|
angle_edge.insert( std::make_pair( angle, i ) );
|
|
}
|
|
}
|
|
|
|
//debugging: print out all the angles below the minimum for a test
|
|
for ( std::multimap<double, int>::const_iterator it = angle_edge.begin(); it != angle_edge.end(); ++it )
|
|
{
|
|
QgsDebugMsg( QString( "angle: %1" ).arg( it->first ) );
|
|
}
|
|
|
|
|
|
double minangle = 0;//actual minimum angle
|
|
int minedge;//first edge adjacent to the minimum angle
|
|
int minedgenext;
|
|
int minedgenextnext;
|
|
|
|
Point3D circumcenter;
|
|
|
|
while ( !edge_angle.empty() )
|
|
{
|
|
minangle = angle_edge.begin()->first;
|
|
QgsDebugMsg( QString( "minangle: %1" ).arg( minangle ) );
|
|
minedge = angle_edge.begin()->second;
|
|
minedgenext = mHalfEdge[minedge]->getNext();
|
|
minedgenextnext = mHalfEdge[minedgenext]->getNext();
|
|
|
|
//calculate the circumcenter
|
|
if ( !MathUtils::circumcenter( mPointVector[mHalfEdge[minedge]->getPoint()], mPointVector[mHalfEdge[minedgenext]->getPoint()], mPointVector[mHalfEdge[minedgenextnext]->getPoint()], &circumcenter ) )
|
|
{
|
|
QgsDebugMsg( "warning, calculation of circumcenter failed" );
|
|
//put all three edges to dontexamine and remove them from the other maps
|
|
dontexamine.insert( minedge );
|
|
edge_angle.erase( minedge );
|
|
std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
|
|
while ( minedgeiter->second != minedge )
|
|
{
|
|
++minedgeiter;
|
|
}
|
|
angle_edge.erase( minedgeiter );
|
|
continue;
|
|
}
|
|
|
|
if ( !pointInside( circumcenter.getX(), circumcenter.getY() ) )
|
|
{
|
|
//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() ) );
|
|
dontexamine.insert( minedge );
|
|
edge_angle.erase( minedge );
|
|
std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
|
|
while ( minedgeiter->second != minedge )
|
|
{
|
|
++minedgeiter;
|
|
}
|
|
angle_edge.erase( minedgeiter );
|
|
continue;
|
|
}
|
|
|
|
/*******find out, if any forced segment or convex hull segment is in the influence region of the circumcenter. In this case, split the segment********************/
|
|
|
|
bool encroached = false;
|
|
|
|
#if 0 //slow version
|
|
int numhalfedges = mHalfEdge.count();//begin slow version
|
|
for ( int i = 0; i < numhalfedges; i++ )
|
|
{
|
|
if ( mHalfEdge[i]->getForced() || edgeOnConvexHull( i ) )
|
|
{
|
|
if ( MathUtils::inDiametral( mPointVector[mHalfEdge[i]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[i]->getDual()]->getPoint()], &circumcenter ) )
|
|
{
|
|
encroached = true;
|
|
//split segment
|
|
QgsDebugMsg( "segment split" );
|
|
int pointno = splitHalfEdge( i, 0.5 );
|
|
|
|
//update dontexmine, angle_edge, edge_angle
|
|
int pointingedge = baseEdgeOfPoint( pointno );
|
|
|
|
int actedge = pointingedge;
|
|
int ed1, ed2, ed3;//numbers of the three edges
|
|
int pt1, pt2, pt3;//numbers of the three points
|
|
double angle1, angle2, angle3;
|
|
|
|
do
|
|
{
|
|
ed1 = mHalfEdge[actedge]->getDual();
|
|
pt1 = mHalfEdge[ed1]->getPoint();
|
|
ed2 = mHalfEdge[ed1]->getNext();
|
|
pt2 = mHalfEdge[ed2]->getPoint();
|
|
ed3 = mHalfEdge[ed2]->getNext();
|
|
pt3 = mHalfEdge[ed3]->getPoint();
|
|
actedge = ed3;
|
|
|
|
if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )//don't consider triangles with the virtual point
|
|
{
|
|
continue;
|
|
}
|
|
|
|
angle1 = MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
|
|
angle2 = MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
|
|
angle3 = MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
|
|
|
|
//don't put the edges on the maps if two segments are forced or on a hull
|
|
bool twoforcededges1, twoforcededges2, twoforcededges3;//flag to indicate, if angle1, angle2 and angle3 are between forced edges or hull edges
|
|
|
|
if (( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) )
|
|
{
|
|
twoforcededges1 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges1 = false;
|
|
}
|
|
|
|
if (( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) )
|
|
{
|
|
twoforcededges2 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges2 = false;
|
|
}
|
|
|
|
if (( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) )
|
|
{
|
|
twoforcededges3 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges3 = false;
|
|
}
|
|
|
|
//update the settings related to ed1
|
|
QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
|
|
if ( ed1iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed1iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed1 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed1 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed1 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle1 < mintol && !twoforcededges1 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed1, angle1 ) );
|
|
angle_edge.insert( std::make_pair( angle1, ed1 ) );
|
|
}
|
|
|
|
//update the settings related to ed2
|
|
QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
|
|
if ( ed2iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed2iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed2 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed2 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed2 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle2 < mintol && !twoforcededges2 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed2, angle2 ) );
|
|
angle_edge.insert( std::make_pair( angle2, ed2 ) );
|
|
}
|
|
|
|
//update the settings related to ed3
|
|
QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
|
|
if ( ed3iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed3iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed3 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed3 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed3 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle3 < mintol && !twoforcededges3 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed3, angle3 ) );
|
|
angle_edge.insert( std::make_pair( angle3, ed3 ) );
|
|
}
|
|
|
|
|
|
}
|
|
while ( actedge != pointingedge );
|
|
}
|
|
}
|
|
}
|
|
#endif //end slow version
|
|
|
|
|
|
//fast version. Maybe this does not work
|
|
QSet<int> influenceedges;//begin fast method
|
|
int baseedge = baseEdgeOfTriangle( &circumcenter );
|
|
if ( baseedge == -5 )//a numerical instability occurred or the circumcenter already exists in the triangulation
|
|
{
|
|
//delete minedge from edge_angle and minangle from angle_edge
|
|
edge_angle.erase( minedge );
|
|
std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
|
|
while ( minedgeiter->second != minedge )
|
|
{
|
|
++minedgeiter;
|
|
}
|
|
angle_edge.erase( minedgeiter );
|
|
continue;
|
|
}
|
|
else if ( baseedge == -25 )//circumcenter already exists in the triangulation
|
|
{
|
|
//delete minedge from edge_angle and minangle from angle_edge
|
|
edge_angle.erase( minedge );
|
|
std::multimap<double, int>::iterator minedgeiter = angle_edge.find( minangle );
|
|
while ( minedgeiter->second != minedge )
|
|
{
|
|
++minedgeiter;
|
|
}
|
|
angle_edge.erase( minedgeiter );
|
|
continue;
|
|
}
|
|
else if ( baseedge == -20 )
|
|
{
|
|
baseedge = mEdgeWithPoint;
|
|
}
|
|
|
|
evaluateInfluenceRegion( &circumcenter, baseedge, influenceedges );
|
|
evaluateInfluenceRegion( &circumcenter, mHalfEdge[baseedge]->getNext(), influenceedges );
|
|
evaluateInfluenceRegion( &circumcenter, mHalfEdge[mHalfEdge[baseedge]->getNext()]->getNext(), influenceedges );
|
|
|
|
for ( QSet<int>::iterator it = influenceedges.begin(); it != influenceedges.end(); ++it )
|
|
{
|
|
if (( mHalfEdge[*it]->getForced() || edgeOnConvexHull( *it ) ) && MathUtils::inDiametral( mPointVector[mHalfEdge[*it]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[*it]->getDual()]->getPoint()], &circumcenter ) )
|
|
{
|
|
//split segment
|
|
QgsDebugMsg( "segment split" );
|
|
int pointno = splitHalfEdge( *it, 0.5 );
|
|
encroached = true;
|
|
|
|
//update dontexmine, angle_edge, edge_angle
|
|
int pointingedge = baseEdgeOfPoint( pointno );
|
|
|
|
int actedge = pointingedge;
|
|
int ed1, ed2, ed3;//numbers of the three edges
|
|
int pt1, pt2, pt3;//numbers of the three points
|
|
double angle1, angle2, angle3;
|
|
|
|
do
|
|
{
|
|
ed1 = mHalfEdge[actedge]->getDual();
|
|
pt1 = mHalfEdge[ed1]->getPoint();
|
|
ed2 = mHalfEdge[ed1]->getNext();
|
|
pt2 = mHalfEdge[ed2]->getPoint();
|
|
ed3 = mHalfEdge[ed2]->getNext();
|
|
pt3 = mHalfEdge[ed3]->getPoint();
|
|
actedge = ed3;
|
|
|
|
if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )//don't consider triangles with the virtual point
|
|
{
|
|
continue;
|
|
}
|
|
|
|
angle1 = MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
|
|
angle2 = MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
|
|
angle3 = MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
|
|
|
|
//don't put the edges on the maps if two segments are forced or on a hull
|
|
bool twoforcededges1, twoforcededges2, twoforcededges3;//flag to decide, if edges should be added to the maps. Do not add them if true
|
|
|
|
|
|
|
|
if (( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) )
|
|
{
|
|
twoforcededges1 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges1 = false;
|
|
}
|
|
|
|
if (( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) )
|
|
{
|
|
twoforcededges2 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges2 = false;
|
|
}
|
|
|
|
if (( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) )
|
|
{
|
|
twoforcededges3 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges3 = false;
|
|
}
|
|
|
|
|
|
//update the settings related to ed1
|
|
QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
|
|
if ( ed1iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed1iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed1 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed1 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed1 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle1 < mintol && !twoforcededges1 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed1, angle1 ) );
|
|
angle_edge.insert( std::make_pair( angle1, ed1 ) );
|
|
}
|
|
|
|
//update the settings related to ed2
|
|
QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
|
|
if ( ed2iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed2iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed2 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed2 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed2 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle2 < mintol && !twoforcededges2 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed2, angle2 ) );
|
|
angle_edge.insert( std::make_pair( angle2, ed2 ) );
|
|
}
|
|
|
|
//update the settings related to ed3
|
|
QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
|
|
if ( ed3iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed3iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed3 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed3 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed3 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle3 < mintol && !twoforcededges3 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed3, angle3 ) );
|
|
angle_edge.insert( std::make_pair( angle3, ed3 ) );
|
|
}
|
|
|
|
|
|
}
|
|
while ( actedge != pointingedge );
|
|
}
|
|
} //end fast method
|
|
|
|
|
|
if ( encroached )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
/*******otherwise, try to add the circumcenter to the triangulation************************************************************************************************/
|
|
|
|
Point3D* p = new Point3D();
|
|
mDecorator->calcPoint( circumcenter.getX(), circumcenter.getY(), p );
|
|
int pointno = mDecorator->addPoint( p );
|
|
|
|
if ( pointno == -100 || pointno == mTwiceInsPoint )
|
|
{
|
|
if ( pointno == -100 )
|
|
{
|
|
QgsDebugMsg( QString( "put circumcenter %1//%2 on dontexamine list because of numerical instabilities" )
|
|
.arg( circumcenter.getX() ).arg( circumcenter.getY() ) );
|
|
}
|
|
else if ( pointno == mTwiceInsPoint )
|
|
{
|
|
QgsDebugMsg( QString( "put circumcenter %1//%2 on dontexamine list because it is already inserted" )
|
|
.arg( circumcenter.getX() ).arg( circumcenter.getY() ) );
|
|
//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() )
|
|
{
|
|
flag = true;
|
|
}
|
|
}
|
|
if ( !flag )
|
|
{
|
|
QgsDebugMsg( "point is not present in the triangulation" );
|
|
}
|
|
}
|
|
//put all three edges to dontexamine and remove them from the other maps
|
|
dontexamine.insert( minedge );
|
|
edge_angle.erase( minedge );
|
|
std::multimap<double, int>::iterator minedgeiter = angle_edge.lower_bound( minangle );
|
|
while ( minedgeiter->second != minedge )
|
|
{
|
|
++minedgeiter;
|
|
}
|
|
angle_edge.erase( minedgeiter );
|
|
continue;
|
|
}
|
|
else//insertion successful
|
|
{
|
|
QgsDebugMsg( "circumcenter added" );
|
|
|
|
//update the maps
|
|
//go around the inserted point and make changes for every half edge
|
|
int pointingedge = baseEdgeOfPoint( pointno );
|
|
|
|
int actedge = pointingedge;
|
|
int ed1, ed2, ed3;//numbers of the three edges
|
|
int pt1, pt2, pt3;//numbers of the three points
|
|
double angle1, angle2, angle3;
|
|
|
|
do
|
|
{
|
|
ed1 = mHalfEdge[actedge]->getDual();
|
|
pt1 = mHalfEdge[ed1]->getPoint();
|
|
ed2 = mHalfEdge[ed1]->getNext();
|
|
pt2 = mHalfEdge[ed2]->getPoint();
|
|
ed3 = mHalfEdge[ed2]->getNext();
|
|
pt3 = mHalfEdge[ed3]->getPoint();
|
|
actedge = ed3;
|
|
|
|
if ( pt1 == -1 || pt2 == -1 || pt3 == -1 )//don't consider triangles with the virtual point
|
|
{
|
|
continue;
|
|
}
|
|
|
|
angle1 = MathUtils::angle( mPointVector[pt3], mPointVector[pt1], mPointVector[pt2], mPointVector[pt1] );
|
|
angle2 = MathUtils::angle( mPointVector[pt1], mPointVector[pt2], mPointVector[pt3], mPointVector[pt2] );
|
|
angle3 = MathUtils::angle( mPointVector[pt2], mPointVector[pt3], mPointVector[pt1], mPointVector[pt3] );
|
|
|
|
//todo: put all three edges on the dontexamine list if two edges are forced or convex hull edges
|
|
bool twoforcededges1, twoforcededges2, twoforcededges3;
|
|
|
|
if (( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) && ( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) )
|
|
{
|
|
twoforcededges1 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges1 = false;
|
|
}
|
|
|
|
if (( mHalfEdge[ed2]->getForced() || edgeOnConvexHull( ed2 ) ) && ( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) )
|
|
{
|
|
twoforcededges2 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges2 = false;
|
|
}
|
|
|
|
if (( mHalfEdge[ed3]->getForced() || edgeOnConvexHull( ed3 ) ) && ( mHalfEdge[ed1]->getForced() || edgeOnConvexHull( ed1 ) ) )
|
|
{
|
|
twoforcededges3 = true;
|
|
}
|
|
else
|
|
{
|
|
twoforcededges3 = false;
|
|
}
|
|
|
|
|
|
//update the settings related to ed1
|
|
QSet<int>::iterator ed1iter = dontexamine.find( ed1 );
|
|
if ( ed1iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed1iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed1 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed1 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed1 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle1 < mintol && !twoforcededges1 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed1, angle1 ) );
|
|
angle_edge.insert( std::make_pair( angle1, ed1 ) );
|
|
}
|
|
|
|
|
|
//update the settings related to ed2
|
|
QSet<int>::iterator ed2iter = dontexamine.find( ed2 );
|
|
if ( ed2iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed2iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed2 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed2 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed2 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle2 < mintol && !twoforcededges2 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed2, angle2 ) );
|
|
angle_edge.insert( std::make_pair( angle2, ed2 ) );
|
|
}
|
|
|
|
//update the settings related to ed3
|
|
QSet<int>::iterator ed3iter = dontexamine.find( ed3 );
|
|
if ( ed3iter != dontexamine.end() )
|
|
{
|
|
//edge number is on dontexamine list
|
|
dontexamine.erase( ed3iter );
|
|
}
|
|
else
|
|
{
|
|
//test, if it is on edge_angle and angle_edge and erase them if yes
|
|
std::map<int, double>::iterator tempit1;
|
|
tempit1 = edge_angle.find( ed3 );
|
|
if ( tempit1 != edge_angle.end() )
|
|
{
|
|
//erase the entries
|
|
double angle = tempit1->second;
|
|
edge_angle.erase( ed3 );
|
|
std::multimap<double, int>::iterator tempit2 = angle_edge.lower_bound( angle );
|
|
while ( tempit2->second != ed3 )
|
|
{
|
|
++tempit2;
|
|
}
|
|
angle_edge.erase( tempit2 );
|
|
}
|
|
}
|
|
|
|
if ( angle3 < mintol && !twoforcededges3 )
|
|
{
|
|
edge_angle.insert( std::make_pair( ed3, angle3 ) );
|
|
angle_edge.insert( std::make_pair( angle3, ed3 ) );
|
|
}
|
|
|
|
|
|
}
|
|
while ( actedge != pointingedge );
|
|
}
|
|
}
|
|
|
|
#if 0
|
|
//debugging: print out all edge of dontexamine
|
|
for ( QSet<int>::iterator it = dontexamine.begin(); it != dontexamine.end(); ++it )
|
|
{
|
|
QgsDebugMsg( QString( "edge nr. %1 is in dontexamine" ).arg( *it ) );
|
|
}
|
|
#endif
|
|
}
|
|
|
|
|
|
bool DualEdgeTriangulation::swapPossible( unsigned int edge )
|
|
{
|
|
//test, if edge belongs to a forced edge
|
|
if ( mHalfEdge[edge]->getForced() )
|
|
{
|
|
return false;
|
|
}
|
|
|
|
//test, if the edge is on the convex hull or is connected to the virtual point
|
|
if ( mHalfEdge[edge]->getPoint() == -1 || mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getPoint() == -1 || mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint() == -1 )
|
|
{
|
|
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()];
|
|
if ( MathUtils::leftOf( ptc, pta, ptb ) > leftOfTresh )
|
|
{
|
|
return false;
|
|
}
|
|
else if ( MathUtils::leftOf( ptd, ptb, ptc ) > leftOfTresh )
|
|
{
|
|
return false;
|
|
}
|
|
else if ( MathUtils::leftOf( pta, ptc, ptd ) > leftOfTresh )
|
|
{
|
|
return false;
|
|
}
|
|
else if ( MathUtils::leftOf( ptb, ptd, pta ) > leftOfTresh )
|
|
{
|
|
return false;
|
|
}
|
|
return true;
|
|
}
|
|
|
|
void DualEdgeTriangulation::triangulatePolygon( QList<int>* poly, QList<int>* free, int mainedge )
|
|
{
|
|
if ( poly && free )
|
|
{
|
|
if ( poly->count() == 3 )//polygon is already a triangle
|
|
{
|
|
return;
|
|
}
|
|
|
|
//search for the edge pointing on the closest point(distedge) and for the next(nextdistedge)
|
|
QList<int>::const_iterator iterator = ++( poly->constBegin() );//go to the second edge
|
|
double distance = MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] );
|
|
int distedge = ( *iterator );
|
|
int nextdistedge = mHalfEdge[( *iterator )]->getNext();
|
|
++iterator;
|
|
|
|
while ( iterator != --( poly->constEnd() ) )
|
|
{
|
|
if ( MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] ) < distance )
|
|
{
|
|
distedge = ( *iterator );
|
|
nextdistedge = mHalfEdge[( *iterator )]->getNext();
|
|
distance = MathUtils::distPointFromLine( mPointVector[mHalfEdge[( *iterator )]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[mainedge]->getPoint()] );
|
|
}
|
|
++iterator;
|
|
}
|
|
|
|
if ( nextdistedge == ( *( --poly->end() ) ) )//the nearest point is connected to the endpoint of mainedge
|
|
{
|
|
int inserta = free->first();//take an edge from the freelist
|
|
int insertb = mHalfEdge[inserta]->getDual();
|
|
free->pop_front();
|
|
|
|
mHalfEdge[inserta]->setNext(( poly->at( 1 ) ) );
|
|
mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
|
|
mHalfEdge[insertb]->setNext( nextdistedge );
|
|
mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
|
|
mHalfEdge[distedge]->setNext( inserta );
|
|
mHalfEdge[mainedge]->setNext( insertb );
|
|
|
|
QList<int> polya;
|
|
for ( iterator = ( ++( poly->constBegin() ) ); ( *iterator ) != nextdistedge; ++iterator )
|
|
{
|
|
polya.append(( *iterator ) );
|
|
}
|
|
polya.prepend( inserta );
|
|
|
|
#if 0
|
|
//print out all the elements of polya for a test
|
|
for ( iterator = polya.begin(); iterator != polya.end(); ++iterator )
|
|
{
|
|
QgsDebugMsg( *iterator );
|
|
}
|
|
#endif
|
|
|
|
triangulatePolygon( &polya, free, inserta );
|
|
}
|
|
|
|
else if ( distedge == ( *( ++poly->begin() ) ) )//the nearest point is connected to the beginpoint of mainedge
|
|
{
|
|
int inserta = free->first();//take an edge from the freelist
|
|
int insertb = mHalfEdge[inserta]->getDual();
|
|
free->pop_front();
|
|
|
|
mHalfEdge[inserta]->setNext(( poly->at( 2 ) ) );
|
|
mHalfEdge[inserta]->setPoint( mHalfEdge[distedge]->getPoint() );
|
|
mHalfEdge[insertb]->setNext( mainedge );
|
|
mHalfEdge[insertb]->setPoint( mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint() );
|
|
mHalfEdge[distedge]->setNext( insertb );
|
|
mHalfEdge[( *( --poly->end() ) )]->setNext( inserta );
|
|
|
|
QList<int> polya;
|
|
iterator = poly->constBegin();
|
|
iterator += 2;
|
|
while ( iterator != poly->constEnd() )
|
|
{
|
|
polya.append(( *iterator ) );
|
|
++iterator;
|
|
}
|
|
polya.prepend( inserta );
|
|
|
|
triangulatePolygon( &polya, free, inserta );
|
|
}
|
|
|
|
else//the nearest point is not connected to an endpoint of mainedge
|
|
{
|
|
int inserta = free->first();//take an edge from the freelist
|
|
int insertb = mHalfEdge[inserta]->getDual();
|
|
free->pop_front();
|
|
|
|
int insertc = free->first();
|
|
int insertd = mHalfEdge[insertc]->getDual();
|
|
free->pop_front();
|
|
|
|
mHalfEdge[inserta]->setNext(( poly->at( 1 ) ) );
|
|
mHalfEdge[inserta]->setPoint( mHalfEdge[mainedge]->getPoint() );
|
|
mHalfEdge[insertb]->setNext( insertd );
|
|
mHalfEdge[insertb]->setPoint( mHalfEdge[distedge]->getPoint() );
|
|
mHalfEdge[insertc]->setNext( nextdistedge );
|
|
mHalfEdge[insertc]->setPoint( mHalfEdge[distedge]->getPoint() );
|
|
mHalfEdge[insertd]->setNext( mainedge );
|
|
mHalfEdge[insertd]->setPoint( mHalfEdge[mHalfEdge[mainedge]->getDual()]->getPoint() );
|
|
|
|
mHalfEdge[distedge]->setNext( inserta );
|
|
mHalfEdge[mainedge]->setNext( insertb );
|
|
mHalfEdge[( *( --poly->end() ) )]->setNext( insertc );
|
|
|
|
//build two new polygons for recursive triangulation
|
|
QList<int> polya;
|
|
QList<int> polyb;
|
|
|
|
for ( iterator = ++( poly->constBegin() ); ( *iterator ) != nextdistedge; ++iterator )
|
|
{
|
|
polya.append(( *iterator ) );
|
|
}
|
|
polya.prepend( inserta );
|
|
|
|
|
|
while ( iterator != poly->constEnd() )
|
|
{
|
|
polyb.append(( *iterator ) );
|
|
++iterator;
|
|
}
|
|
polyb.prepend( insertc );
|
|
|
|
triangulatePolygon( &polya, free, inserta );
|
|
triangulatePolygon( &polyb, free, insertc );
|
|
}
|
|
}
|
|
|
|
else
|
|
{
|
|
QgsDebugMsg( "warning, null pointer" );
|
|
}
|
|
|
|
}
|
|
|
|
bool DualEdgeTriangulation::pointInside( double x, double y )
|
|
{
|
|
Point3D 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
|
|
int numinstabs = 0;//number of suspect left-of-tests due to 'leftOfTresh'
|
|
int runs = 0;//counter for the number of iterations in the loop to prevent an endless loop
|
|
|
|
while ( true )
|
|
{
|
|
if ( runs > nBaseOfRuns )//prevents endless loops
|
|
{
|
|
QgsDebugMsg( QString( "warning, instability detected: Point coordinates: %1//%2" ).arg( x ).arg( y ) );
|
|
return false;
|
|
}
|
|
|
|
if ( MathUtils::leftOf( &point, mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] ) < ( -leftOfTresh ) )//point is on the left side
|
|
{
|
|
counter += 1;
|
|
if ( counter == 3 )//three successful passes means that we have found the triangle
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
else if ( MathUtils::leftOf( &point, mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] ) == 0 )//point is exactly in the line of the edge
|
|
{
|
|
counter += 1;
|
|
mEdgeWithPoint = actedge;
|
|
nulls += 1;
|
|
if ( counter == 3 )//three successful passes means that we have found the triangle
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
else if ( MathUtils::leftOf( &point, mPointVector[mHalfEdge[mHalfEdge[actedge]->getDual()]->getPoint()], mPointVector[mHalfEdge[actedge]->getPoint()] ) < leftOfTresh )//numerical problems
|
|
{
|
|
counter += 1;
|
|
numinstabs += 1;
|
|
if ( counter == 3 )//three successful passes means that we have found the triangle
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
else//point is on the right side
|
|
{
|
|
actedge = mHalfEdge[actedge]->getDual();
|
|
counter = 1;
|
|
nulls = 0;
|
|
numinstabs = 0;
|
|
}
|
|
|
|
actedge = mHalfEdge[actedge]->getNext();
|
|
if ( mHalfEdge[actedge]->getPoint() == -1 )//the half edge points to the virtual point
|
|
{
|
|
if ( nulls == 1 )//point is exactly on the convex hull
|
|
{
|
|
return true;
|
|
}
|
|
mEdgeOutside = ( unsigned int )mHalfEdge[mHalfEdge[actedge]->getNext()]->getNext();
|
|
return false;//the point is outside the convex hull
|
|
}
|
|
runs++;
|
|
}
|
|
|
|
if ( nulls == 2 )//we hit an existing point
|
|
{
|
|
return true;
|
|
}
|
|
if ( numinstabs > 0 )//a numerical instability occurred
|
|
{
|
|
QgsDebugMsg( "numerical instabilities" );
|
|
return true;
|
|
}
|
|
|
|
if ( nulls == 1 )//point is on an existing edge
|
|
{
|
|
return true;
|
|
}
|
|
mEdgeInside = actedge;
|
|
return true;
|
|
}
|
|
|
|
#if 0
|
|
bool DualEdgeTriangulation::readFromTAFF( QString filename )
|
|
{
|
|
QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );//change the cursor
|
|
|
|
QFile file( filename );//the file to be read
|
|
file.open( IO_Raw | IO_ReadOnly );
|
|
QBuffer buffer( file.readAll() );//the buffer to copy the file to
|
|
file.close();
|
|
|
|
QTextStream textstream( &buffer );
|
|
buffer.open( IO_ReadOnly );
|
|
|
|
QString buff;
|
|
int numberofhalfedges;
|
|
int numberofpoints;
|
|
|
|
//edge section
|
|
while ( buff.mid( 0, 4 ) != "TRIA" )//search for the TRIA-section
|
|
{
|
|
buff = textstream.readLine();
|
|
}
|
|
while ( buff.mid( 0, 4 ) != "NEDG" )
|
|
{
|
|
buff = textstream.readLine();
|
|
}
|
|
numberofhalfedges = buff.section( ' ', 1, 1 ).toInt();
|
|
mHalfEdge.resize( numberofhalfedges );
|
|
|
|
|
|
while ( buff.mid( 0, 4 ) != "DATA" )//search for the data section
|
|
{
|
|
textstream >> buff;
|
|
}
|
|
|
|
int nr1, nr2, dual1, dual2, point1, point2, next1, next2, fo1, fo2, b1, b2;
|
|
bool forced1, forced2, break1, break2;
|
|
|
|
|
|
//import all the DualEdges
|
|
QProgressBar* edgebar = new QProgressBar();//use a progress dialog so that it is not boring for the user
|
|
edgebar->setCaption( "Reading edges..." );
|
|
edgebar->setTotalSteps( numberofhalfedges / 2 );
|
|
edgebar->setMinimumWidth( 400 );
|
|
edgebar->move( 500, 500 );
|
|
edgebar->show();
|
|
|
|
for ( int i = 0; i < numberofhalfedges / 2; i++ )
|
|
{
|
|
if ( i % 1000 == 0 )
|
|
{
|
|
edgebar->setProgress( i );
|
|
}
|
|
|
|
textstream >> nr1;
|
|
textstream >> point1;
|
|
textstream >> next1;
|
|
textstream >> fo1;
|
|
|
|
if ( fo1 == 0 )
|
|
{
|
|
forced1 = false;
|
|
}
|
|
else
|
|
{
|
|
forced1 = true;
|
|
}
|
|
|
|
textstream >> b1;
|
|
|
|
if ( b1 == 0 )
|
|
{
|
|
break1 = false;
|
|
}
|
|
else
|
|
{
|
|
break1 = true;
|
|
}
|
|
|
|
textstream >> nr2;
|
|
textstream >> point2;
|
|
textstream >> next2;
|
|
textstream >> fo2;
|
|
|
|
if ( fo2 == 0 )
|
|
{
|
|
forced2 = false;
|
|
}
|
|
else
|
|
{
|
|
forced2 = true;
|
|
}
|
|
|
|
textstream >> b2;
|
|
|
|
if ( b2 == 0 )
|
|
{
|
|
break2 = false;
|
|
}
|
|
else
|
|
{
|
|
break2 = true;
|
|
}
|
|
|
|
HalfEdge* hf1 = new HalfEdge();
|
|
hf1->setDual( nr2 );
|
|
hf1->setNext( next1 );
|
|
hf1->setPoint( point1 );
|
|
hf1->setBreak( break1 );
|
|
hf1->setForced( forced1 );
|
|
|
|
HalfEdge* hf2 = new HalfEdge();
|
|
hf2->setDual( nr1 );
|
|
hf2->setNext( next2 );
|
|
hf2->setPoint( point2 );
|
|
hf2->setBreak( break2 );
|
|
hf2->setForced( forced2 );
|
|
|
|
// QgsDebugMsg( QString( "inserting half edge pair %1" ).arg( i ) );
|
|
mHalfEdge.insert( nr1, hf1 );
|
|
mHalfEdge.insert( nr2, hf2 );
|
|
|
|
}
|
|
|
|
edgebar->setProgress( numberofhalfedges / 2 );
|
|
delete edgebar;
|
|
|
|
//set mEdgeInside to a reasonable value
|
|
for ( int i = 0; i < numberofhalfedges; i++ )
|
|
{
|
|
int a, b, c, d;
|
|
a = mHalfEdge[i]->getPoint();
|
|
b = mHalfEdge[mHalfEdge[i]->getDual()]->getPoint();
|
|
c = mHalfEdge[mHalfEdge[i]->getNext()]->getPoint();
|
|
d = mHalfEdge[mHalfEdge[mHalfEdge[i]->getDual()]->getNext()]->getPoint();
|
|
if ( a != -1 && b != -1 && c != -1 && d != -1 )
|
|
{
|
|
mEdgeInside = i;
|
|
break;
|
|
}
|
|
}
|
|
|
|
//point section
|
|
while ( buff.mid( 0, 4 ) != "POIN" )
|
|
{
|
|
buff = textstream.readLine();
|
|
QgsDebugMsg( buff );
|
|
}
|
|
while ( buff.mid( 0, 4 ) != "NPTS" )
|
|
{
|
|
buff = textstream.readLine();
|
|
QgsDebugMsg( buff );
|
|
}
|
|
numberofpoints = buff.section( ' ', 1, 1 ).toInt();
|
|
mPointVector.resize( numberofpoints );
|
|
|
|
while ( buff.mid( 0, 4 ) != "DATA" )
|
|
{
|
|
textstream >> buff;
|
|
}
|
|
|
|
QProgressBar* pointbar = new QProgressBar();
|
|
pointbar->setCaption( "Reading points..." );
|
|
pointbar->setTotalSteps( numberofpoints );
|
|
pointbar->setMinimumWidth( 400 );
|
|
pointbar->move( 500, 500 );
|
|
pointbar->show();
|
|
|
|
|
|
double x, y, z;
|
|
for ( int i = 0; i < numberofpoints; i++ )
|
|
{
|
|
if ( i % 1000 == 0 )
|
|
{
|
|
pointbar->setProgress( i );
|
|
}
|
|
|
|
textstream >> x;
|
|
textstream >> y;
|
|
textstream >> z;
|
|
|
|
Point3D* p = new Point3D( x, y, z );
|
|
|
|
// QgsDebugMsg( QString( "inserting point %1" ).arg( i ) );
|
|
mPointVector.insert( i, p );
|
|
|
|
if ( i == 0 )
|
|
{
|
|
xMin = x;
|
|
xMax = x;
|
|
yMin = y;
|
|
yMax = y;
|
|
}
|
|
else
|
|
{
|
|
//update the bounding box
|
|
if ( x < xMin )
|
|
{
|
|
xMin = x;
|
|
}
|
|
else if ( x > xMax )
|
|
{
|
|
xMax = x;
|
|
}
|
|
|
|
if ( y < yMin )
|
|
{
|
|
yMin = y;
|
|
}
|
|
else if ( y > yMax )
|
|
{
|
|
yMax = y;
|
|
}
|
|
}
|
|
}
|
|
|
|
pointbar->setProgress( numberofpoints );
|
|
delete pointbar;
|
|
QApplication::restoreOverrideCursor();
|
|
|
|
}
|
|
|
|
bool DualEdgeTriangulation::saveToTAFF( QString filename ) const
|
|
{
|
|
QFile outputfile( filename );
|
|
if ( !outputfile.open( IO_WriteOnly ) )
|
|
{
|
|
QMessageBox::warning( 0, "warning", "File could not be written", QMessageBox::Ok, QMessageBox::NoButton, QMessageBox::NoButton );
|
|
return false;
|
|
}
|
|
|
|
QTextStream outstream( &outputfile );
|
|
outstream.precision( 9 );
|
|
|
|
//export the edges. Attention, dual edges must be adjacent in the TAFF-file
|
|
outstream << "TRIA" << std::endl << std::flush;
|
|
outstream << "NEDG " << mHalfEdge.count() << std::endl << std::flush;
|
|
outstream << "PANO 1" << std::endl << std::flush;
|
|
outstream << "DATA ";
|
|
|
|
bool* cont = new bool[mHalfEdge.count()];
|
|
for ( unsigned int i = 0; i <= mHalfEdge.count() - 1; i++ )
|
|
{
|
|
cont[i] = false;
|
|
}
|
|
|
|
for ( unsigned int i = 0; i < mHalfEdge.count(); i++ )
|
|
{
|
|
if ( cont[i] )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
int dual = mHalfEdge[i]->getDual();
|
|
outstream << i << " " << mHalfEdge[i]->getPoint() << " " << mHalfEdge[i]->getNext() << " " << mHalfEdge[i]->getForced() << " " << mHalfEdge[i]->getBreak() << " ";
|
|
outstream << dual << " " << mHalfEdge[dual]->getPoint() << " " << mHalfEdge[dual]->getNext() << " " << mHalfEdge[dual]->getForced() << " " << mHalfEdge[dual]->getBreak() << " ";
|
|
cont[i] = true;
|
|
cont[dual] = true;
|
|
}
|
|
outstream << std::endl << std::flush;
|
|
outstream << std::endl << std::flush;
|
|
|
|
delete[] cont;
|
|
|
|
//export the points to the file
|
|
outstream << "POIN" << std::endl << std::flush;
|
|
outstream << "NPTS " << getNumberOfPoints() << std::endl << std::flush;
|
|
outstream << "PATT 3" << std::endl << std::flush;
|
|
outstream << "DATA ";
|
|
|
|
for ( int i = 0; i < getNumberOfPoints(); i++ )
|
|
{
|
|
Point3D* p = mPointVector[i];
|
|
outstream << p->getX() << " " << p->getY() << " " << p->getZ() << " ";
|
|
}
|
|
outstream << std::endl << std::flush;
|
|
outstream << std::endl << std::flush;
|
|
|
|
return true;
|
|
}
|
|
#endif //0
|
|
|
|
bool DualEdgeTriangulation::swapEdge( double x, double y )
|
|
{
|
|
Point3D p( x, y, 0 );
|
|
int edge1 = baseEdgeOfTriangle( &p );
|
|
if ( edge1 >= 0 )
|
|
{
|
|
int edge2, edge3;
|
|
Point3D* point1;
|
|
Point3D* point2;
|
|
Point3D* point3;
|
|
edge2 = mHalfEdge[edge1]->getNext();
|
|
edge3 = mHalfEdge[edge2]->getNext();
|
|
point1 = getPoint( mHalfEdge[edge1]->getPoint() );
|
|
point2 = getPoint( mHalfEdge[edge2]->getPoint() );
|
|
point3 = getPoint( mHalfEdge[edge3]->getPoint() );
|
|
if ( point1 && point2 && point3 )
|
|
{
|
|
//find out the closest edge to the point and swap this edge
|
|
double dist1, dist2, dist3;
|
|
dist1 = MathUtils::distPointFromLine( &p, point3, point1 );
|
|
dist2 = MathUtils::distPointFromLine( &p, point1, point2 );
|
|
dist3 = MathUtils::distPointFromLine( &p, point2, point3 );
|
|
if ( dist1 <= dist2 && dist1 <= dist3 )
|
|
{
|
|
//qWarning("edge "+QString::number(edge1)+" is closest");
|
|
if ( swapPossible( edge1 ) )
|
|
{
|
|
doOnlySwap( edge1 );
|
|
}
|
|
}
|
|
else if ( dist2 <= dist1 && dist2 <= dist3 )
|
|
{
|
|
//qWarning("edge "+QString::number(edge2)+" is closest");
|
|
if ( swapPossible( edge2 ) )
|
|
{
|
|
doOnlySwap( edge2 );
|
|
}
|
|
}
|
|
else if ( dist3 <= dist1 && dist3 <= dist2 )
|
|
{
|
|
//qWarning("edge "+QString::number(edge3)+" is closest");
|
|
if ( swapPossible( edge3 ) )
|
|
{
|
|
doOnlySwap( edge3 );
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
// QgsDebugMsg("warning: null pointer");
|
|
return false;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// QgsDebugMsg("Edge number negative");
|
|
return false;
|
|
}
|
|
}
|
|
|
|
QList<int>* DualEdgeTriangulation::getPointsAroundEdge( double x, double y )
|
|
{
|
|
Point3D p( x, y, 0 );
|
|
int p1, p2, p3, p4;
|
|
int edge1 = baseEdgeOfTriangle( &p );
|
|
if ( edge1 >= 0 )
|
|
{
|
|
int edge2, edge3;
|
|
Point3D* point1;
|
|
Point3D* point2;
|
|
Point3D* point3;
|
|
edge2 = mHalfEdge[edge1]->getNext();
|
|
edge3 = mHalfEdge[edge2]->getNext();
|
|
point1 = getPoint( mHalfEdge[edge1]->getPoint() );
|
|
point2 = getPoint( mHalfEdge[edge2]->getPoint() );
|
|
point3 = getPoint( mHalfEdge[edge3]->getPoint() );
|
|
if ( point1 && point2 && point3 )
|
|
{
|
|
double dist1, dist2, dist3;
|
|
dist1 = MathUtils::distPointFromLine( &p, point3, point1 );
|
|
dist2 = MathUtils::distPointFromLine( &p, point1, point2 );
|
|
dist3 = MathUtils::distPointFromLine( &p, point2, point3 );
|
|
if ( dist1 <= dist2 && dist1 <= dist3 )
|
|
{
|
|
p1 = mHalfEdge[edge1]->getPoint();
|
|
p2 = mHalfEdge[mHalfEdge[edge1]->getNext()]->getPoint();
|
|
p3 = mHalfEdge[mHalfEdge[edge1]->getDual()]->getPoint();
|
|
p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge1]->getDual()]->getNext()]->getPoint();
|
|
}
|
|
else if ( dist2 <= dist1 && dist2 <= dist3 )
|
|
{
|
|
p1 = mHalfEdge[edge2]->getPoint();
|
|
p2 = mHalfEdge[mHalfEdge[edge2]->getNext()]->getPoint();
|
|
p3 = mHalfEdge[mHalfEdge[edge2]->getDual()]->getPoint();
|
|
p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge2]->getDual()]->getNext()]->getPoint();
|
|
}
|
|
else if ( dist3 <= dist1 && dist3 <= dist2 )
|
|
{
|
|
p1 = mHalfEdge[edge3]->getPoint();
|
|
p2 = mHalfEdge[mHalfEdge[edge3]->getNext()]->getPoint();
|
|
p3 = mHalfEdge[mHalfEdge[edge3]->getDual()]->getPoint();
|
|
p4 = mHalfEdge[mHalfEdge[mHalfEdge[edge3]->getDual()]->getNext()]->getPoint();
|
|
}
|
|
QList<int>* list = new QList<int>();
|
|
list->append( p1 );
|
|
list->append( p2 );
|
|
list->append( p3 );
|
|
list->append( p4 );
|
|
return list;
|
|
}
|
|
else
|
|
{
|
|
QgsDebugMsg( "warning: null pointer" );
|
|
return nullptr;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
QgsDebugMsg( "Edge number negative" );
|
|
return nullptr;
|
|
}
|
|
}
|
|
|
|
bool DualEdgeTriangulation::saveAsShapefile( const QString& fileName ) const
|
|
{
|
|
QString shapeFileName = fileName;
|
|
|
|
QgsFields fields;
|
|
fields.append( QgsField( "type", QVariant::String, "String" ) );
|
|
|
|
// add the extension if not present
|
|
if ( shapeFileName.indexOf( ".shp" ) == -1 )
|
|
{
|
|
shapeFileName += ".shp";
|
|
}
|
|
|
|
//delete already existing files
|
|
if ( QFile::exists( shapeFileName ) )
|
|
{
|
|
if ( !QgsVectorFileWriter::deleteShapeFile( shapeFileName ) )
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
|
|
QgsVectorFileWriter writer( shapeFileName, "Utf-8", fields, Qgis::WKBLineString );
|
|
if ( writer.hasError() != QgsVectorFileWriter::NoError )
|
|
{
|
|
return false;
|
|
}
|
|
|
|
bool *alreadyVisitedEdges = new bool[mHalfEdge.size()];
|
|
if ( !alreadyVisitedEdges )
|
|
{
|
|
QgsDebugMsg( "out of memory" );
|
|
return false;
|
|
}
|
|
|
|
for ( int i = 0; i < mHalfEdge.size(); ++i )
|
|
{
|
|
alreadyVisitedEdges[i] = false;
|
|
}
|
|
|
|
for ( int i = 0; i < mHalfEdge.size(); ++i )
|
|
{
|
|
HalfEdge* currentEdge = mHalfEdge[i];
|
|
if ( currentEdge->getPoint() != -1 && mHalfEdge[currentEdge->getDual()]->getPoint() != -1 && !alreadyVisitedEdges[currentEdge->getDual()] )
|
|
{
|
|
QgsFeature edgeLineFeature;
|
|
|
|
//geometry
|
|
Point3D* p1 = mPointVector[currentEdge->getPoint()];
|
|
Point3D* p2 = mPointVector[mHalfEdge[currentEdge->getDual()]->getPoint()];
|
|
QgsPolyline lineGeom;
|
|
lineGeom.push_back( QgsPoint( p1->getX(), p1->getY() ) );
|
|
lineGeom.push_back( QgsPoint( p2->getX(), p2->getY() ) );
|
|
QgsGeometry* geom = QgsGeometry::fromPolyline( lineGeom );
|
|
edgeLineFeature.setGeometry( *geom );
|
|
delete geom;
|
|
edgeLineFeature.initAttributes( 1 );
|
|
|
|
//attributes
|
|
QString attributeString;
|
|
if ( currentEdge->getForced() )
|
|
{
|
|
if ( currentEdge->getBreak() )
|
|
{
|
|
attributeString = "break line";
|
|
}
|
|
else
|
|
{
|
|
attributeString = "structure line";
|
|
}
|
|
}
|
|
edgeLineFeature.setAttribute( 0, attributeString );
|
|
|
|
writer.addFeature( edgeLineFeature );
|
|
}
|
|
alreadyVisitedEdges[i] = true;
|
|
}
|
|
|
|
delete [] alreadyVisitedEdges;
|
|
|
|
return true;
|
|
}
|
|
|
|
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() );
|
|
|
|
//search for the minimum angle (it is important, which directions the lines have!)
|
|
double minangle;
|
|
double angle1 = MathUtils::angle( p1, p2, p4, p2 );
|
|
minangle = angle1;
|
|
double angle2 = MathUtils::angle( p3, p2, p4, p2 );
|
|
if ( angle2 < minangle )
|
|
{
|
|
minangle = angle2;
|
|
}
|
|
double angle3 = MathUtils::angle( p2, p3, p4, p3 );
|
|
if ( angle3 < minangle )
|
|
{
|
|
minangle = angle3;
|
|
}
|
|
double angle4 = MathUtils::angle( p3, p4, p2, p4 );
|
|
if ( angle4 < minangle )
|
|
{
|
|
minangle = angle4;
|
|
}
|
|
double angle5 = MathUtils::angle( p2, p4, p1, p4 );
|
|
if ( angle5 < minangle )
|
|
{
|
|
minangle = angle5;
|
|
}
|
|
double angle6 = MathUtils::angle( p4, p1, p2, p1 );
|
|
if ( angle6 < minangle )
|
|
{
|
|
minangle = angle6;
|
|
}
|
|
|
|
return minangle;
|
|
}
|
|
|
|
int DualEdgeTriangulation::splitHalfEdge( int edge, float position )
|
|
{
|
|
//just a short test if position is between 0 and 1
|
|
if ( position < 0 || position > 1 )
|
|
{
|
|
QgsDebugMsg( "warning, position is not between 0 and 1" );
|
|
}
|
|
|
|
//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 );
|
|
|
|
//calculate the z-value of the point to insert
|
|
Point3D zvaluepoint;
|
|
mDecorator->calcPoint( p->getX(), p->getY(), &zvaluepoint );
|
|
p->setZ( zvaluepoint.getZ() );
|
|
|
|
//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() ) );
|
|
mPointVector.insert( mPointVector.count(), p );
|
|
|
|
//insert the six new halfedges
|
|
int dualedge = mHalfEdge[edge]->getDual();
|
|
int edge1 = insertEdge( -10, -10, mPointVector.count() - 1, false, false );
|
|
int edge2 = insertEdge( edge1, mHalfEdge[mHalfEdge[edge]->getNext()]->getNext(), mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint(), false, false );
|
|
int edge3 = insertEdge( -10, mHalfEdge[mHalfEdge[dualedge]->getNext()]->getNext(), mHalfEdge[mHalfEdge[dualedge]->getNext()]->getPoint(), false, false );
|
|
int edge4 = insertEdge( edge3, dualedge, mPointVector.count() - 1, false, false );
|
|
int edge5 = insertEdge( -10, mHalfEdge[edge]->getNext(), mHalfEdge[edge]->getPoint(), mHalfEdge[edge]->getBreak(), mHalfEdge[edge]->getForced() );
|
|
int edge6 = insertEdge( edge5, edge3, mPointVector.count() - 1, mHalfEdge[dualedge]->getBreak(), mHalfEdge[dualedge]->getForced() );
|
|
mHalfEdge[edge1]->setDual( edge2 );
|
|
mHalfEdge[edge1]->setNext( edge5 );
|
|
mHalfEdge[edge3]->setDual( edge4 );
|
|
mHalfEdge[edge5]->setDual( edge6 );
|
|
|
|
//adjust the already existing halfedges
|
|
mHalfEdge[mHalfEdge[edge]->getNext()]->setNext( edge1 );
|
|
mHalfEdge[mHalfEdge[dualedge]->getNext()]->setNext( edge4 );
|
|
mHalfEdge[edge]->setNext( edge2 );
|
|
mHalfEdge[edge]->setPoint( mPointVector.count() - 1 );
|
|
mHalfEdge[mHalfEdge[edge3]->getNext()]->setNext( edge6 );
|
|
|
|
//test four times recursively for swapping
|
|
checkSwap( mHalfEdge[edge5]->getNext(), 0 );
|
|
checkSwap( mHalfEdge[edge2]->getNext(), 0 );
|
|
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
|
|
|
|
return mPointVector.count() - 1;
|
|
}
|
|
|
|
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 )
|
|
{
|
|
if ( set.find( edge ) == set.end() )
|
|
{
|
|
set.insert( edge );
|
|
}
|
|
else//prevent endless loops
|
|
{
|
|
return;
|
|
}
|
|
|
|
if ( !mHalfEdge[edge]->getForced() && !edgeOnConvexHull( edge ) )
|
|
{
|
|
//test, if point is in the circle through both endpoints of edge and the endpoint of edge->dual->next->point
|
|
if ( MathUtils::inCircle( point, mPointVector[mHalfEdge[mHalfEdge[edge]->getDual()]->getPoint()], mPointVector[mHalfEdge[edge]->getPoint()], mPointVector[mHalfEdge[mHalfEdge[edge]->getNext()]->getPoint()] ) )
|
|
{
|
|
evaluateInfluenceRegion( point, mHalfEdge[mHalfEdge[edge]->getDual()]->getNext(), set );
|
|
evaluateInfluenceRegion( point, mHalfEdge[mHalfEdge[mHalfEdge[edge]->getDual()]->getNext()]->getNext(), set );
|
|
}
|
|
}
|
|
}
|