www.pudn.com > cad3d.zip > DelaunayTriangulation.cpp
// DelaunayTriangulation.cpp: implementation of the CDelaunayTriangulation class.
//
//////////////////////////////////////////////////////////////////////
#include "Topology\stdafx.h"
#include "Topology\DelaunayTriangulation.h"
#include "3DMath\3DSegment.h"
#include
/*
class CDtVertex : public CVertex
{
public:
CDtVertex(double x = 0,double y = 0,double z = 0,double w = 1)
: CVertex (x,y,z,w) {}
};
*/
#define CLIP_CONCAVE
class CDtEdge : public CSameFacesEdge
{
public:
CDtEdge(CDtVertex* pV0 = NULL, CDtVertex* pV1 = NULL)
: CSameFacesEdge(pV0, pV1) {}
C3DSegment GetSegment() const
{
return C3DSegment(GetVertexRef(0), GetVertexRef(1));
}
C3DSegment GetNormalSeg() const
{
C3DPoint mPt( GetSegment().Evaluate(0.5) );
C3DVector v( GetSegment().GetVector());
C3DVector n( -v(1), v(0), 0);
C3DPoint oPt( mPt + n);
return C3DSegment(mPt, oPt);
}
};
static C3DSegment sTestSeg;
static C3DSegment sEdgeSeg;
static C3DSegment sInnerSeg;
class CDtTri : public CLinkedPoly
{
public:
CDtTri( CDtVertex* pVer0 = NULL, CDtVertex* pVer1 = NULL, CDtVertex* pVer2 = NULL,
CDtEdge* pEdge0 = NULL,CDtEdge* pEdge1 = NULL,CDtEdge* pEdge2 = NULL)
: CLinkedPoly()
{
SetVertex(0, pVer0);
SetVertex(1, pVer1);
SetVertex(2, pVer2);
SetLink(0, pEdge0);
SetLink(1, pEdge1);
SetLink(2, pEdge2);
}
bool IsIn(const C3DPoint& rPt, CDtEdge*& pOnEdge)
{
C3DSegment testSeg( (*this)(0), rPt);
C3DSegment edgeSeg( GetLinkRef(0).GetSegment() );
sTestSeg = testSeg;
sEdgeSeg = edgeSeg;
CDoubleRange testRange, edgeRange, inRange(0,1);
if ( testSeg.IsPoint() || edgeSeg.IsPoint() )
return false;
testRange = testSeg.Intersect(edgeSeg, edgeRange);
if ( testRange.IsEmpty() )
{
return false;
}
CHECK( testRange.IsPoint() );
CHECK( edgeRange.IsPoint() );
if ( inRange.IsIn( edgeRange.GetMax() ) )
{//the line (*this)(0) intersects the edge(0)
C3DPoint intPt = edgeSeg.Evaluate(edgeRange.GetMax());
C3DSegment innerSeg( (*this)(0), intPt);
sInnerSeg = innerSeg;
double dParam = innerSeg.GetEvalParamFast(rPt);
if ( inRange.IsIn( (math_real)dParam ) )
{
if ( IsZero(edgeRange.GetMax()) )
{// point lies on edge #2 (see SetupEdges())
pOnEdge = GetLinkPtr(2);
}
else if ( Equal((math_real)edgeRange.GetMax(), (math_real)1.0) )
{// point lies on edge #1 (see SetupEdges())
pOnEdge = GetLinkPtr(1);
}
else if ( Equal((math_real)dParam, (math_real)1.0) )
{// point lies on edge #0 (see SetupEdges())
pOnEdge = GetLinkPtr((top_int)0);
}
else
{
pOnEdge = NULL;
}
return true;
}
}
return false;
}
void SetupEdges()
{
CHECK(GetLinkPtr(0) != NULL);
CHECK(GetLinkPtr(1) != NULL);
CHECK(GetLinkPtr(2) != NULL);
GetLinkRef(2).SetVertex(0, GetVertexPtr((top_int)0) );
GetLinkRef(2).SetVertex(1, GetVertexPtr(1) );
CHECK_RES( GetLinkRef(2).SetOnNull(this) );
GetLinkRef(0).SetVertex(0, GetVertexPtr(1) );
GetLinkRef(0).SetVertex(1, GetVertexPtr(2) );
CHECK_RES( GetLinkRef(0).SetOnNull(this) );
GetLinkRef(1).SetVertex(0, GetVertexPtr(2) );
GetLinkRef(1).SetVertex(1, GetVertexPtr((top_int)0) );
CHECK_RES( GetLinkRef(1).SetOnNull(this) );
}
C3DPoint GetOutCircleCenter() const
{
C3DSegment n0( GetLinkRef(0).GetNormalSeg() );
C3DSegment n1( GetLinkRef(1).GetNormalSeg() );
CDoubleRange n0Range, n1Range;
n0Range = n0.Intersect(n1, n1Range);
if ( !n0Range.IsPoint() || !n1Range.IsPoint() )
return C3DVector();
CHECK( n0Range.IsPoint() );
CHECK( n1Range.IsPoint() );
return n0.Evaluate( n0Range.GetMax() );
}
double GetOutCircleRadius() const
{
C3DPoint center( GetOutCircleCenter() );
if ( center.IsVector() )
return 0;
C3DVector v( center - (*this)(0) );
return v.Norm();
}
//doesn't check if rPt is on the plane(faster)
bool IsInOutCircle(const C3DPoint& rPt) const
{
C3DPoint center( GetOutCircleCenter() );
if ( center.IsVector() )
return false;
C3DPoint pt(rPt);
pt.Normalize();
C3DVector v( pt - center );
if ( IsLt((math_real)v.Norm(), (math_real)GetOutCircleRadius() ) )
return true;
return false;
}
void DrawOutCircle()
{
C3DPoint center( GetOutCircleCenter() );
if ( center.IsVector() )
return;
/*
double dR = GetOutCircleRadius();
glBegin(GL_LINE_STRIP );
for (double dFi = 0; dFi < 2*M_PI; dFi+= 0.01)
{
double x = center(0) + cos(dFi)*dR;
double y = center(1) + sin(dFi)*dR;
glVertex3d(x,y,0);
}
glEnd();
*/
}
void DrawVoroniDiagram()
{
C3DPoint center( GetOutCircleCenter() );
if ( center.IsVector() )
return;
/*
glBegin(GL_LINES);
for (int i = 0; i < 3; ++i)
{
C3DPoint middle( GetLinkRef(i).GetSegment().Evaluate(0.5) );
glVertex3dv(middle);
glVertex3dv(center);
}
glEnd();
*/
}
bool IsToClip(CDtVertex* pVertices, int nVertCount, double dLimit)
{
C3DPoint center( GetOutCircleCenter() );
if ( center.IsVector() )
return true;
//checks if the center is in/out of the contour
/* CDtEdge* pOnEdge = NULL;
if ( !IsIn(center, pOnEdge) || pOnEdge != NULL)
return true;
*/
C3DPoint limitPt(center(0), (math_real)dLimit);
C3DSegment ray(center, limitPt);
bool bIsIn = false;
CDoubleRange rayRange, testRange, inRange(0,1);
for (int i = 1; i < nVertCount; ++i)
{
C3DSegment testSeg(pVertices[i], pVertices[i-1]);
rayRange = ray.Intersect(testSeg, testRange);
if ( rayRange.IsPoint() && testRange.IsPoint() )
{
if ( inRange.IsIn(testRange.GetMax()) && inRange.IsIn(rayRange.GetMax()) )
{
bIsIn = !bIsIn;
}
}
}
return !bIsIn;
}
// void SetupVertices();
};
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CDelaunayTriangulation::CDelaunayTriangulation (C3DPoint* pVertices, int nVertsCount)
: m_nVertsCount(nVertsCount), m_pVertices(pVertices), m_nTrisCount(0), m_pTris(NULL),
m_nEdgesCount(0), m_pEdges(NULL)
{
Init();
}
CDelaunayTriangulation::~CDelaunayTriangulation()
{
delete [] m_pTris;
delete [] m_pEdges;
delete [] m_pVertices;
}
void CDelaunayTriangulation::Init()
{
m_nVertsCount += 3;
CDtVertex* pTempVertices = m_pVertices;
m_pVertices = new CDtVertex[m_nVertsCount];
memcpy(m_pVertices, pTempVertices, sizeof(CDtVertex)*(m_nVertsCount - 3) );
int nMaxEdges = 3*m_nVertsCount - 3;
int nMaxTris = 2*m_nVertsCount /*- 2*/;
m_pEdges = new CDtEdge[nMaxEdges];
m_pTris = new CDtTri[nMaxTris];
AddLimitTri();
}
void CDelaunayTriangulation::AddLimitTri()
{
double dMaxWidth = 0;
for (int i = 0; i < m_nVertsCount - 3; ++i)
{
for (int j = 0; j < 3; ++j)
{
if ( dMaxWidth < fabs(m_pVertices[i](j)) )
{
dMaxWidth = fabs(m_pVertices[i](j));
}
}
}
//the last three vers...
double dLimit = 3.5*dMaxWidth;
m_pVertices[m_nVertsCount - 1] = CDtVertex((math_real)dLimit,0,0);
m_pVertices[m_nVertsCount - 2] = CDtVertex(0,(math_real)dLimit,0);
m_pVertices[m_nVertsCount - 3] = CDtVertex((math_real)-dLimit,(math_real)-dLimit,0);
//the first tri is the limit tri
m_pTris[0].SetVertex(0,&m_pVertices[m_nVertsCount - 1]);
m_pTris[0].SetVertex(1,&m_pVertices[m_nVertsCount - 2]);
m_pTris[0].SetVertex(2,&m_pVertices[m_nVertsCount - 3]);
m_nTrisCount++;
m_pTris[0].SetLink(0, &m_pEdges[0]);
m_pTris[0].SetLink(1, &m_pEdges[1]);
m_pTris[0].SetLink(2, &m_pEdges[2]);
m_pTris[0].SetupEdges();
m_nEdgesCount += 3;
m_dLimit = dLimit;
}
CDtTri* CDelaunayTriangulation::FindContainingTri(const CDtVertex& rPt, CDtEdge*& pOnEdge)
{
//skipps the limit tri
for (int i = 0; i < m_nTrisCount; ++i)
{
if ( m_pTris[i].IsIn(rPt, pOnEdge) )
{
return (m_pTris + i);
}
}
//returns the limit tri
return NULL;
}
static CDtTri* spTri0;
static CDtTri* spTri1;
static CDtTri* spTri2;
void CDelaunayTriangulation::SplitTri(CDtVertex& rPt, CDtTri* pOnTri)
{
CDtTri tempTri = *pOnTri;
//pEdgeX is the new edge containing the vertex(X) of the triangle pTri
CDtEdge* pEdge0 = m_pEdges + m_nEdgesCount++;
CDtEdge* pEdge1 = m_pEdges + m_nEdgesCount++;
CDtEdge* pEdge2 = m_pEdges + m_nEdgesCount++;
pEdge0->SetVertex(0,&rPt);
pEdge0->SetVertex(1,tempTri.GetVertexPtr((top_int)0));
pEdge1->SetVertex(0,&rPt);
pEdge1->SetVertex(1,tempTri.GetVertexPtr(1));
pEdge2->SetVertex(0,&rPt);
pEdge2->SetVertex(1,tempTri.GetVertexPtr(2));
CDtTri* pTri0 = /*(pOnTri == m_pLimitTri) ? (m_pTris + m_nTrisCount++) :*/ pOnTri;
CDtTri* pTri1 = m_pTris + m_nTrisCount++;
CDtTri* pTri2 = m_pTris + m_nTrisCount++;
*pTri0 = CDtTri( tempTri.GetVertexPtr(1), tempTri.GetVertexPtr(2), &rPt,
pEdge2, pEdge1, tempTri.GetLinkPtr((top_int)0));
*pTri1 = CDtTri( tempTri.GetVertexPtr(2), tempTri.GetVertexPtr((top_int)0), &rPt,
pEdge0, pEdge2, tempTri.GetLinkPtr(1));
*pTri2 = CDtTri( tempTri.GetVertexPtr((top_int)0), tempTri.GetVertexPtr(1), &rPt,
pEdge1, pEdge0, tempTri.GetLinkPtr(2));
tempTri.GetLinkPtr((top_int)0)->ChangeFace(pOnTri, pTri0);
tempTri.GetLinkPtr(1)->ChangeFace(pOnTri, pTri1);
tempTri.GetLinkPtr(2)->ChangeFace(pOnTri, pTri2);
pEdge0->SetLeftFace(pTri1);
pEdge0->SetRightFace(pTri2);
pEdge1->SetLeftFace(pTri0);
pEdge1->SetRightFace(pTri2);
pEdge2->SetLeftFace(pTri0);
pEdge2->SetRightFace(pTri1);
spTri0 = pTri0;
spTri1 = pTri1;
spTri2 = pTri2;
LegalizeTri(&rPt, pTri0);
LegalizeTri(&rPt, pTri1);
LegalizeTri(&rPt, pTri2);
}
void CDelaunayTriangulation::SplitEdge(CDtVertex& rPt, CDtEdge* pOnEdge)
{
CDtEdge tempEdge = *pOnEdge;
CDtTri* pLeftTri = tempEdge.GetLeftFacePtr();
CDtTri* pRightTri = tempEdge.GetRightFacePtr();
if ( pLeftTri == NULL || pRightTri == NULL) return;
CDtTri tempLeftTri(*pLeftTri);
CDtTri tempRightTri(*pRightTri);
//we have 3 new edges and we use one old - this is the pOnEdge, which is splitted
CDtEdge* pEdge0 = pOnEdge;
CDtEdge* pEdge1 = m_pEdges + m_nEdgesCount++;
CDtEdge* pEdge2 = m_pEdges + m_nEdgesCount++;
CDtEdge* pEdge3 = m_pEdges + m_nEdgesCount++;
//setting up the middle edges
//pEdge0 & pEdge1 are the edges colinear with pOnEdge
pEdge0->SetVertex(0, tempEdge.GetVertexPtr(0) );
pEdge0->SetVertex(1, &rPt );
pEdge1->SetVertex(0, &rPt );
pEdge1->SetVertex(1, tempEdge.GetVertexPtr(1) );
//pEdge2 splits left triangle in two
pEdge2->SetVertex(0, tempLeftTri.GetVertexPtr(pOnEdge) );
pEdge2->SetVertex(1, &rPt);
//pEdge3 splits rigth triangle in two
pEdge3->SetVertex(0, &rPt );
pEdge3->SetVertex(1, tempRightTri.GetVertexPtr(pOnEdge) );
//two old tris will be used, as two new too
{
//the corresponding edges opposite to the splitting edge's vertices
//for the 'left' tri
CDtEdge* pLtEdge0 = tempLeftTri.GetLinkPtr( tempEdge.GetVertexPtr(0) );
CDtEdge* pLtEdge1 = tempLeftTri.GetLinkPtr( tempEdge.GetVertexPtr(1) );
CDtTri* pTri0 = pLeftTri;
CDtTri* pTri1 = pRightTri;
//pTri0 & pTri1 are the remains of pLeftTri
//pLtEdge1 is the edge opposite of pOnEdge's #1 vertex, so now it'll be opposite for rPt
*pTri0 = CDtTri(pEdge0->GetVertexPtr(0), &rPt, pEdge2->GetVertexPtr(0),
pEdge2, pLtEdge1 , pEdge0);
pLtEdge1->ChangeFace(pLeftTri, pTri0);//this call is unnessasaty, but just in case...
//pLtEdge0 is the edge opposite of pOnEdge's #0 vertex, so now it'll be opposite for rPt
*pTri1 = CDtTri(pEdge1->GetVertexPtr(1), &rPt, pEdge2->GetVertexPtr(0),
pEdge2, pLtEdge0 , pEdge1);
pLtEdge0->ChangeFace(pLeftTri, pTri1);
pEdge2->SetLeftFace(pTri0);
pEdge2->SetRightFace(pTri1);
pEdge0->SetLeftFace(pTri0);
pEdge1->SetLeftFace(pTri1);
}
{
//the corresponding edges opposite to the splitting edge's vertices
//for the 'right' tri
CDtEdge* pRtEdge0 = tempRightTri.GetLinkPtr( tempEdge.GetVertexPtr(0) );
CDtEdge* pRtEdge1 = tempRightTri.GetLinkPtr( tempEdge.GetVertexPtr(1) );
//pTri2 & pTri3 are the remains of pRightTri
CDtTri* pTri2 = m_pTris + m_nTrisCount++;
CDtTri* pTri3 = m_pTris + m_nTrisCount++;
//pRtEdge1 is the edge opposite of pOnEdge's #1 vertex, so now it'll be opposite for rPt
*pTri2 = CDtTri(pEdge0->GetVertexPtr(0), &rPt, pEdge3->GetVertexPtr(1),
pEdge3, pRtEdge1 , pEdge0);
pRtEdge1->ChangeFace(pRightTri, pTri2);
//pRtEdge0 is the edge opposite of pOnEdge's #0 vertex, so now it'll be opposite for rPt
*pTri3 = CDtTri(pEdge1->GetVertexPtr(1), &rPt, pEdge3->GetVertexPtr(1),
pEdge3, pRtEdge0 , pEdge1);
pRtEdge0->ChangeFace(pRightTri, pTri3);
pEdge3->SetLeftFace(pTri2);
pEdge3->SetRightFace(pTri3);
pEdge0->SetRightFace(pTri2);
pEdge1->SetRightFace(pTri3);
}
}
void CDelaunayTriangulation::LegalizeTri(CDtVertex* pPt, CDtTri* pTri)
{
CDtEdge* pTestEdge = pTri->GetLinkPtr(pPt);
CDtTri* pTestTri = pTestEdge->GetOtherFacePtr(pTri);
if ( pTestTri == NULL) return;
CDtVertex *pTestPt = pTestTri->GetVertexPtr(pTestEdge);
// bool bLegal = true;
CDtVertex* pVer0 = pPt;
CDtVertex* pVer1 = pTestPt;
CDtVertex* pVer2 = pTestEdge->GetVertexPtr(0);
CDtVertex* pVer3 = pTestEdge->GetVertexPtr(1);
if ( pTri->IsInOutCircle(*pTestPt) )
{//Flip pTestEdge
CDtTri* pTri0 = pTri;
CDtTri* pTri1 = pTestTri;
//flips edge
pTestEdge->SetVertex(0, pVer0);
pTestEdge->SetVertex(1, pVer1);
CDtEdge* pEdge0 = pTri0->GetLinkPtr(pVer3);
CDtEdge* pEdge1 = pTri1->GetLinkPtr(pVer3);
CDtEdge* pEdge2 = pTri0->GetLinkPtr(pVer2);
CDtEdge* pEdge3 = pTri1->GetLinkPtr(pVer2);
*pTri0 = CDtTri(pVer0, pVer1, pVer2, pEdge1, pEdge0, pTestEdge);
pEdge1->ChangeFace(pTri1, pTri0);
//pEdge0->ChangeFace(pTri0, pTri0);
*pTri1 = CDtTri(pVer0, pVer1, pVer3, pEdge3, pEdge2, pTestEdge);
pEdge2->ChangeFace(pTri0, pTri1);
//pEdge3->ChangeFace(pTri1, pTri1);
LegalizeTri(pPt, pTri0);
LegalizeTri(pPt, pTri1);
}
}
/*
void CDelaunayTriangulation::LegalizeEdge(CDtVertex& rPt, CDtEdge* pEdge)
{
}
*/
bool CDelaunayTriangulation::IsLimitTri(const CDtTri& rTri)
{
for (int i = 0; i < 3; ++i)
{
if ( IsLimitVertex( rTri.GetVertexPtr(i) ) )
return true;
}
return false;
}
bool CDelaunayTriangulation::IsLimitEdge(const CDtEdge& rEdge)
{
for (int i = 0; i < 2; ++i)
{
if ( IsLimitVertex( rEdge.GetVertexPtr(i) ) )
return true;
}
return false;
}
bool CDelaunayTriangulation::IsLimitVertex(const CDtVertex* pVertex)
{
return ((pVertex - m_pVertices) >= (m_nVertsCount - 3));
}
void CDelaunayTriangulation::Perform(const COglRenderer& rRc)
{
int nOrigVertsCount = m_nVertsCount - 3;
for (int i = 0; i < nOrigVertsCount; ++i)
{
CDtVertex& rPt = m_pVertices[i];
CDtEdge* pOnEdge = NULL;
CDtTri* pOnTri = FindContainingTri(rPt, pOnEdge);
// CHECK(pOnTri != NULL);
if ( pOnTri == NULL)
continue;
if (pOnEdge == NULL)
{
SplitTri(rPt, pOnTri);
}
else
{
SplitEdge(rPt, pOnEdge);
}
}
DrawTris(rRc);
}
void CDelaunayTriangulation::DrawEdges(const COglRenderer& /*rRc*/)
{
/*
rRc.SetDrawColor(1,1,0);
glBegin(GL_LINES);
for (int i = 0; i < m_nEdgesCount; ++i)
{
if (IsLimitEdge(m_pEdges[i]) )
continue;
rRc.AddVertex(m_pEdges[i].GetVertexRef(0));
rRc.AddVertex(m_pEdges[i].GetVertexRef(1));
}
glEnd();
rRc.SetDrawColor(0,0,1);
*/
}
void CDelaunayTriangulation::DrawTris(const COglRenderer& /*rRc*/)
{
#if 0
/*
CDtVertex testv0(-50,0,0);
CDtVertex testv1( 50,0,0);
CDtVertex testv2( 0,50,0);
CDtEdge e0( &testv0, &testv1);
CDtEdge e1( &testv1, &testv2);
CDtEdge e2( &testv2, &testv0);
CDtTri testt(&testv0, &testv1, &testv2, &e1, &e2, &e0);
rRc.SetDrawColor(0,0,1);
testt.DrawOutCircle();
return;
*/
rRc.SetDrawColor(0.4f, 0.3f, 0.7f);
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
for (int i = 0; i < m_nTrisCount; ++i)
{
/*
if ( (m_pTris + i) == spTri0 )
rRc.SetDrawColor(0,0,1);
else if ( (m_pTris + i) == spTri1 )
rRc.SetDrawColor(0,1,0);
else if ( (m_pTris + i) == spTri2 )
rRc.SetDrawColor(1,0,0);
else
{
// rRc.SetDrawColor(0.2,0.2,0.2);
float fI = (float)(i+1)/(float)m_nTrisCount;
rRc.SetDrawColor(fI, fI, fI);
}
*/
if ( IsLimitTri(m_pTris[i]) )
{
continue;
rRc.SetDrawColor(1,1,0);
}
#ifdef CLIP_CONCAVE
else if ( m_pTris[i].IsToClip( m_pVertices, m_nVertsCount - 3, m_dLimit ) )
{
continue;
}
#endif
else
{
rRc.SetDrawColor(1,1,1);
}
// rRc.SetDrawColor((float)i/(float)m_nTrisCount, 0.3f, 0.7f);
/*
rRc.BeginVertices(GL_TRIANGLES);
rRc.AddVertex((m_pTris[i])(0));
rRc.AddVertex((m_pTris[i])(1));
rRc.AddVertex((m_pTris[i])(2));
rRc.EndVertices();
*/
/*
rRc.SetDrawColor(0.3f,0.3f,0.3f);
m_pTris[i].DrawOutCircle();
*/
rRc.SetDrawColor(1.0f,0.1f,0.7f);
m_pTris[i].DrawVoroniDiagram();
/* if ( (m_pTris + i) == spTri0 )
{
rRc.SetDrawColor(0,0,1);
m_pTris[i].DrawOutCircle();
glLineWidth(1);
}
else if ( (m_pTris + i) == spTri1 )
{
rRc.SetDrawColor(0,1,0);
m_pTris[i].DrawOutCircle();
glLineWidth(1);
}
else if ( (m_pTris + i) == spTri2 )
{
rRc.SetDrawColor(1,0,0);
m_pTris[i].DrawOutCircle();
glLineWidth(1);
}
*/
}
rRc.SetDrawColor(0,0,1);
/*
rRc.BeginVertices(GL_LINES);
glLineWidth(10);
rRc.SetDrawColor(1,1,0);
rRc.AddVertex( sTestSeg.GetPtRef(0) );
rRc.AddVertex( sTestSeg.GetPtRef(1) );
rRc.SetDrawColor(1,1,1);
rRc.AddVertex( sEdgeSeg.GetPtRef(0) );
rRc.AddVertex( sEdgeSeg.GetPtRef(1) );
rRc.SetDrawColor(0,1,1);
rRc.AddVertex( sInnerSeg.GetPtRef(0) );
rRc.AddVertex( sInnerSeg.GetPtRef(1) );
rRc.EndVertices();
glLineWidth(1);
*/
#endif
}