/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2012 Excensus LLC.
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Licence as published
 * by the Free Software Foundation. 
 * See the COPYING file for more information.
 *
 **********************************************************************
 *
 * Last port: triangulate/Vertex.java r524
 *
 **********************************************************************/

#include <geos/triangulate/quadedge/Vertex.h>

#include <geos/triangulate/quadedge/TrianglePredicate.h>
#include <geos/triangulate/quadedge/QuadEdge.h>
#include <geos/algorithm/NotRepresentableException.h>


namespace geos {
namespace triangulate { //geos.triangulate
namespace quadedge { //geos.triangulate.quadedge

using namespace algorithm;
using namespace geom;

Vertex::Vertex(double _x, double _y) : p(_x, _y)
{
}

Vertex::Vertex(double _x, double _y, double _z): p( _x, _y, _z)
{
}

Vertex::Vertex(const Coordinate &_p) : p(_p)
{
}

Vertex::Vertex() : p() {
}

int Vertex::classify(const Vertex &p0, const Vertex &p1)
{
	Vertex &p2 = *this;
	std::auto_ptr<Vertex> a = p1.sub(p0);
	std::auto_ptr<Vertex> b = p2.sub(p0);
	double sa = a->crossProduct(*b);
	int ret;

	if (sa > 0.0)
		ret =  LEFT;
	if (sa < 0.0)
		ret =  RIGHT;
	if ((a->getX() * b->getX() < 0.0) || (a->getY() * b->getY() < 0.0))
		ret =  BEHIND;
	if (a->magn() < b->magn())
		ret =  BEYOND;
	if (p0.equals(p2))
		ret =  ORIGIN;
	if (p1.equals(p2))
		ret =  DESTINATION;
	else
		ret =  BETWEEN;

	return ret;
}

bool Vertex::isInCircle(const Vertex &a, const Vertex &b, const Vertex &c) const
{
	return TrianglePredicate::isInCircleRobust(a.p, b.p, c.p, this->p);
	// non-robust - best to not use
	//return TrianglePredicate.isInCircle(a.p, b.p, c.p, this->p);
}

bool Vertex::rightOf(const QuadEdge &e) const {
	return isCCW(e.dest(), e.orig());
}

bool Vertex::leftOf(const QuadEdge &e) const {
	return isCCW(e.orig(), e.dest());
}

std::auto_ptr<HCoordinate> Vertex::bisector(const Vertex &a, const Vertex &b)
{
	// returns the perpendicular bisector of the line segment ab
	double dx = b.getX() - a.getX();
	double dy = b.getY() - a.getY();
	HCoordinate l1 = HCoordinate(a.getX() + dx / 2.0, a.getY() + dy / 2.0, 1.0);
	HCoordinate l2 = HCoordinate(a.getX() - dy + dx / 2.0, a.getY() + dx + dy / 2.0, 1.0);

	return std::auto_ptr<HCoordinate>(new HCoordinate(l1, l2));
}

double Vertex::circumRadiusRatio(const Vertex &b, const Vertex &c)
{
	std::auto_ptr<Vertex> x(circleCenter(b, c));
	double radius = distance(*x, b);
	double edgeLength = distance(*this, b);
	double el = distance(b, c);
	if (el < edgeLength)
	{
		edgeLength = el;
	}
	el = distance(c, *this);
	if (el < edgeLength)
	{
		edgeLength = el;
	}

	return radius / edgeLength;
}

std::auto_ptr<Vertex> Vertex::midPoint(const Vertex &a)
{
	double xm = (p.x + a.getX()) / 2.0;
	double ym = (p.y + a.getY()) / 2.0;
	double zm = (p.z + a.getZ()) / 2.0;
	return std::auto_ptr<Vertex>(new Vertex(xm, ym, zm));
}

std::auto_ptr<Vertex> Vertex::circleCenter(const Vertex &b, const Vertex &c) const
{
	std::auto_ptr<Vertex> a(new Vertex(getX(), getY()));
	// compute the perpendicular bisector of cord ab
	std::auto_ptr<HCoordinate> cab = bisector(*a, b);
	// compute the perpendicular bisector of cord bc
	std::auto_ptr<HCoordinate> cbc = bisector(b, c);
	// compute the intersection of the bisectors (circle radii)
	std::auto_ptr<HCoordinate> hcc(new HCoordinate(*cab, *cbc));
	std::auto_ptr<Vertex> cc;

	try
	{
		cc.reset(new Vertex(hcc->getX(), hcc->getY()));
	} catch (NotRepresentableException nre) {
	}

	return cc;
}

double Vertex::interpolateZValue(const Vertex &v0, const Vertex &v1,
		const Vertex &v2) const
{
	double x0 = v0.getX();
	double y0 = v0.getY();
	double a = v1.getX() - x0;
	double b = v2.getX() - x0;
	double c = v1.getY() - y0;
	double d = v2.getY() - y0;
	double det = a * d - b * c;
	double dx = this->getX() - x0;
	double dy = this->getY() - y0;
	double t = (d * dx - b * dy) / det;
	double u = (-c * dx + a * dy) / det;
	double z = v0.getZ() + t * (v1.getZ() - v0.getZ()) + u * (v2.getZ() - v0.getZ());
	return z;
}

double Vertex::interpolateZ(const Coordinate &p, const Coordinate &v0, 
		const Coordinate &v1, const Coordinate &v2)
{
	double x0 = v0.x;
	double y0 = v0.y;
	double a = v1.x - x0;
	double b = v2.x - x0;
	double c = v1.y - y0;
	double d = v2.y - y0;
	double det = a * d - b * c;
	double dx = p.x - x0;
	double dy = p.y - y0;
	double t = (d * dx - b * dy) / det;
	double u = (-c * dx + a * dy) / det;
	double z = v0.z + t * (v1.z - v0.z) + u * (v2.z - v0.z);
	return z;
}

double Vertex::interpolateZ(const Coordinate &p, const Coordinate &p0, 
		const Coordinate &p1)
{
	double segLen = p0.distance(p1);
	double ptLen = p.distance(p0);
	double dz = p1.z - p0.z;
	double pz = p0.z + dz * (ptLen / segLen);
	return pz;
}

} //namespace geos.triangulate.quadedge
} //namespace geos.triangulate
} //namespace geos