分形几何

开发平台：
Visual C++

delaunaytriangulation.cpp：源码内容
							//
// Delaunay Triangulation
//
// Homework of CG lesson (Fall 2009) in Tsinghua University.
// All rights reserved.
//
// temp include, delete later if necessary
#include "stdafx.h"
#include "delaunaytriangulation.h"
// C++ headers
#include <cassert>
#include <algorithm>
#include <stack>
// User headers
#include "utility.h"
#include "dcel.h"
using namespace std;
DelaunayTriangulation::DelaunayTriangulation(void)
{
    isFinished = false;
    nPoints = 0;
    points = NULL;
}
Point2d* DelaunayTriangulation::point(int pos) const
{
    assert(pos >= 0 && pos < nPoints);
    return (points + pos);
}
void DelaunayTriangulation::initial(int* count, Point2d* ps)
{
    nPoints = *count;
    points = ps;
    std::sort(points, points + nPoints);
    // Deletes the repeated elements.
    int j = 1;
    for (int i = 1; i < nPoints; ++i)
    {
        if (points[i].x != points[i-1].x || points[i].y != points[i-1].y)
        {
            points[j++] = points[i];
        }
    }
    nPoints = *count = j;
    isFinished = false;
}
void DelaunayTriangulation::destroy(void)
{
	if(!maxEdge.le)
	{
	    return;
	}
	// collect dcels
	vector<DCEL *>dcels;
	dcels.clear();
	collectDcel(maxEdge.le, dcels);
	collectDcel(maxEdge.re, dcels);
	// delete the DCEL.
	vector<DCEL *>::iterator iter;
	for(iter = dcels.begin(); iter != dcels.end(); iter++)
	    delete (*iter);
    maxEdge.le = NULL;
	maxEdge.re = NULL;
	nPoints = 0;
    isFinished = false;
    points = NULL;
}
void DelaunayTriangulation::collectDcel(Edge* e, vector<DCEL*> &dcels)
{
    DCEL *pDcel;
	Edge *edge;
	
	std::stack<Edge *>dcelStack;
	dcelStack.push(e);
	while(!dcelStack.empty())
	{
		edge = dcelStack.top();
		pDcel = (DCEL *)(edge->qEdge());
		dcelStack.pop();
	    if(!pDcel->visited)
		{
			dcels.push_back(pDcel);
			pDcel->visited = true;
			dcelStack.push(edge->oNext());
			dcelStack.push(edge->oPrev());
			dcelStack.push(edge->dNext());
			dcelStack.push(edge->dPrev());
		}
	}
	
}
void DelaunayTriangulation::doDelaunayTriangulation(void)
{
    if (nPoints >= 2)
    {
        maxEdge = delaunay(0, nPoints - 1);
        isFinished = true;
    }
}
MaxEdge DelaunayTriangulation::delaunay(int begin, int end)
{
    // the number of points
    int size = end - begin + 1;
    // max edge
    MaxEdge ret;
    // delaunay triangulation
    if(size == 2) // only two points
    {
        // let s1, s2 be two sites in sorted order.
        // create an edge a from s1 to s2
        Edge* a = makeEdge();
        a->endPoints(point(begin), point(end));
        ret.le = a;
        ret.re = a->twin();
        return ret;
    }
    else if(size == 3) // only three points
    {
        // let s1, s2, s3 be the three sites, in sorted order,
        // create edges a connecting s1 to s2 and b connecting s2 to s3
        Edge* a = makeEdge();
        Edge* b = makeEdge();
        splice(a->twin(), b);
        a->endPoints(point(begin), point(begin + 1));
        b->endPoints(point(begin + 1), point(end));
        // close the triangle
        Edge* c;
        if(ccw(*point(begin), *point(begin + 1), *point(end)))
        {
            c = connect(b, a);
            ret.le = a;
            ret.re = b->twin();
            return ret;
        }
        else if(ccw(*point(begin), *point(end), *point(begin + 1)))
        {
            c = connect(b, a);
            ret.le = c->twin();
            ret.re = c;
            return ret;
        }
        else  // the three points are collinear
        {
            ret.le = a;
            ret.re = b->twin();
            return ret;
        }
    }
    else // |sites| >= 4
    {
        Edge *ldo, *ldi; // left half result
        Edge *rdo, *rdi; // right half result
        // recursively delaunay triangulation L and R halves
        MaxEdge leftRet = delaunay(begin, begin + (size / 2) - 1);
        MaxEdge rightRet = delaunay(begin + (size / 2), end);
        ldo = leftRet.le;
        ldi = leftRet.re;
        rdi = rightRet.le;
        rdo = rightRet.re;
        // Compute the lower commom tangent of L and R
        while(1)
        {
            if(leftOf(*(rdi->org()), ldi))
            {
                ldi = ldi->lNext();
            }
            else if(rightOf(*(ldi->org()), rdi))
            {
                rdi = rdi->rPrev();
            }
            else
            {
                break;
            }
        }
        // create a first edge basel from rdi.org to ldi.org
        Edge* basel = connect(rdi->twin(), ldi);
        if((ldi->org()) == (ldo->org()))
        {
            ldo = basel->twin();
        }
        if((rdi->org()) == (rdo->org()))
        {
            rdo = basel;
        }
        // merge
        while(1)
        {
            // locate the first L point (lcand.Dest) to be encountered by the rising bubble
            // and delete L edges out of basel.Dest that fail the circle test
            Edge* lcand = basel->twin()->oNext();
            Edge* t;
            if(valid(lcand, basel))
            {
                while(inCircle(*(basel->dest()), *(basel->org()), *(lcand->dest()), 
                               *(lcand->oNext()->dest())))
                {
                    t = lcand->oNext();
                    deleteEdge(lcand);
                    lcand = t;
                }
            }
            // symmetrically, locate the first R point to be hit, and delete R edges
            Edge* rcand = basel->oPrev();
            if(valid(rcand, basel))
            {
                while(inCircle(*(basel->dest()), *(basel->org()), *(rcand->dest()), 
                               *(rcand->oPrev()->dest())))
                {
                    t = rcand->oPrev();
                    deleteEdge(rcand);
                    rcand = t;
                }
            }
            // if both lcand and rcand are invalid, then basel is the upper common tangent
            if((!valid(lcand, basel)) && (!valid(rcand, basel)))
            {
                break;
            }
            // the next cross edge is to be connected to either lcand.Dest or rcand.Dest
            // if both are valid, then choose the appropriate one using the inCircle test
            if((!valid(lcand, basel)) || (valid(rcand, basel)
                    && inCircle(*(lcand->dest()), *(lcand->org()),
                                *(rcand->org()), *(rcand->dest()))))
            {
                // add cross edge basel from rcand.Dest to basel.Dest
                basel = connect(rcand, basel->twin());
            }
            else
            {
                // add cross edge basel from basel.org to lcand.Dest
                basel = connect(basel->twin(), lcand->twin());
            }
        }
        ret.le = ldo;
        ret.re = rdo;
        return ret;
    }
}

						