分形几何

开发平台：
Visual C++

utility.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 "utility.h"
// C++headers
#include <cmath>
// User headers
#include "dcel.h"
#include "base.h"
#include "line2d.h"
Edge* makeEdge(void)
{
    DCEL* dcel = new DCEL();
    return dcel->e;
}
void splice(Edge* a, Edge* b)
{
    a->prev()->eNext = b;
    b->prev()->eNext = a;
    Edge* t1 = a->prev();
    Edge* t2 = b->prev();
    a->ePrev = t2;
    b->ePrev = t1;
}
void deleteEdge(Edge* e)
{
    if (e != NULL)
    {
        splice(e, e->oPrev());
        splice(e->sym(), e->sym()->oPrev());
        delete e->qEdge();
        e = NULL;
    }
}
Edge* connect(Edge* a, Edge* b)
{
    Edge* e = makeEdge();
    splice(e, a->lNext());
    splice(e->sym(), b);
    e->endPoints(a->dest(), b->org());
    return e;
}
void swap(Edge* e)
{
    Edge* a = e->oPrev();
    Edge* b = e->sym()->oPrev();
    splice(e, a);
    splice(e->sym(), b);
    splice(e, a->lNext());
    splice(e->sym(), b->lNext());
    e->endPoints(a->dest(), b->dest());
}
double triarea(const Point2d& a, const Point2d& b, const Point2d& c)
{
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
bool inCircle(const Point2d& a, const Point2d& b,
              const Point2d& c, const Point2d& d)
{
    if (&d == &a || &d == &b || &d == &c)
        return false;
    return (a.x*a.x + a.y*a.y) * triarea(b, c, d) -
           (b.x*b.x + b.y*b.y) * triarea(a, c, d) +
           (c.x*c.x + c.y*c.y) * triarea(a, b, d) -
           (d.x*d.x + d.y*d.y) * triarea(a, b, c) > 0;
}
bool ccw(const Point2d& a, const Point2d& b, const Point2d& c)
{
    return (triarea(a, b, c) > 0);
}
bool rightOf(const Point2d& x, Edge* e)
{
    return ccw(x, e->dest2d(), e->org2d());
}
bool leftOf(const Point2d& x, Edge* e)
{
    return ccw(x, e->org2d(), e->dest2d());
}
bool onEdge(const Point2d& x, Edge* e)
{
    double t1, t2, t3;
    t1 = (x - e->org2d()).squareNorm();
    t2 = (x - e->dest2d()).squareNorm();
    if (t1 < TOLERANCE || t2 < TOLERANCE)
        return true;
    t3 = (e->org2d() - e->dest2d()).squareNorm();
    if (t1 > t3 || t2 > t3)
        return false;
    Line2d line(e->org2d(), e->dest2d());
    return (fabs(line.eval(x)) < TOLERANCE);
}
bool valid(Edge* e, Edge* basel)
{
	return rightOf(*(e->dest()), basel);
}
bool testInside(const Point2d &x, Edge* startEdge)
{
	Edge* e = startEdge->twin();
	do {
		if (leftOf(x, e))
			return false;
		e = e->lNext();
	} while (e != startEdge->twin());
	return true;
}
Edge* locate(const Point2d& x, Edge* startEdge)
{
	if (!testInside(x, startEdge))
		return NULL;
    Edge* e = startEdge;
	int lastState = 0;
    while (true) 
	{
        if (&x == e->org() || &x == e->dest())
            return e;
		else if (onEdge(x, e))
		    return e;
        else if (!rightOf(x, e->oNext()))
		{
			if (lastState != 1)
			{
				lastState = 1;
				startEdge = e;
			}
            e = e->oNext();
			if (e == startEdge)
				return NULL;
		}
        else if (!rightOf(x, e->dPrev()))
		{
			if (lastState != 2)
			{
				lastState = 2;
				startEdge = e;
			}
            e = e->dPrev();
			if (e == startEdge)
				return NULL;
		}
        else
            return e;
    }
}
bool lineIntersected(const Point2d& p1, const Point2d& p2,
                     const Point2d& p3, const Point2d& p4)
{
    if((ccw(p3, p4, p1) != ccw(p3, p4, p2)) &&
		ccw(p1, p2, p3) != ccw(p1, p2, p4))
	    return true;
	else
		return false;
}
Edge* searchIntersectedEdge(Point2d& p1, Point2d& p2, Edge* startEdge)
{
	Edge *e = startEdge->twin()->lNext();
	bool isIntersect = false;
	while (e != startEdge->twin())
	{
		if (lineIntersected(p1, p2, e->org2d(), e->dest2d()) && leftOf(p1, e))
		{
			isIntersect = true;
			return e->twin();
		}
		e = e->lNext();
	}
	if (lineIntersected(p1, p2, e->org2d(), e->dest2d()) && leftOf(p1, e))
		return e->twin();
	else
		return NULL;
}
void intersectedTris(Point2d& p1, Point2d& p2,
                     vector<Edge *>& edges, Edge* startEdge)
{
    // find the first intersected triangle
	Edge *initEdge = locate(p1, startEdge);
	bool outerCase1 = false;
	bool outerCase2 = false;
	if (initEdge == NULL)
	{
		outerCase1 = true;
		initEdge = searchIntersectedEdge(p1, p2, startEdge);  // search outer intersect edge.
		if (initEdge == NULL)  // no intersection.
			return;
	}
	if (locate(p2, startEdge) == NULL)
		outerCase2 = true;
	Edge *iterEdge;
	Edge *firstEdge = initEdge;
	bool isIntersected = false;
	bool isFirst = false;
	// decide whether p1 is on the initEdge
	if(!onEdge(p1, initEdge))
	{
		while(1)
		{
			isIntersected = false;
			iterEdge = initEdge->lNext();
			while(iterEdge != initEdge)
			{
				if (iterEdge == firstEdge || iterEdge == firstEdge->twin())
					return; //break
				if(lineIntersected(p1, p2, iterEdge->org2d(), iterEdge->dest2d()))
				{
					isIntersected = true;
					break;
				}
				iterEdge = iterEdge->lNext();
			}
			if(!isFirst && !isIntersected && !outerCase1)
			{
			    if(lineIntersected(p1, p2, iterEdge->org2d(), iterEdge->dest2d()))
				{
					isIntersected = true;
				}
			}
			// p1, p2 is not surrounded by tri define by initEdge
			if(isIntersected)
			{
				if (!isFirst)
					isFirst = true;
				edges.push_back(iterEdge);
				initEdge = iterEdge->twin();
			}
			else
			{
				if(!isFirst)
				{
					edges.push_back(initEdge);
					return;
				}
				else if (isFirst && !outerCase2)
				{
					if (!outerCase2)
						edges.push_back(initEdge);
					return;
				}
				else
				{
					return;
				}
			}
		}
	}
	else
	{
	}
}
void circumCircle(const Point2d& p1, const Point2d& p2, const Point2d& p3,
                  Point2d* pCenter, double* radius)
{
	double u1 = 0.5 * (pow(p2.x, 2.0) - pow(p1.x, 2.0) + pow(p2.y, 2.0) - pow(p1.y, 2.0));
	double u2 = 0.5 * (pow(p3.x, 2.0) - pow(p1.x, 2.0) + pow(p3.y, 2.0) - pow(p1.y, 2.0));
	double d11 = p2.x - p1.x;
	double d12 = p2.y - p1.y;
	double d21 = p3.x - p1.x;
	double d22 = p3.y - p1.y;
	pCenter->x = (u1 * d22 - u2 * d12) / (d11 * d22 - d21 * d12);
	pCenter->y = (u2 * d11 - u1 * d21) / (d11 * d22 - d21 * d12);
	*radius = sqrt((pCenter->x - p1.x) * (pCenter->x - p1.x)
        + (pCenter->y - p1.y) * (pCenter->y - p1.y));
}

						