// K-3D
// Copyright (c) 1995-2004, Timothy M. Shead
//
// Contact: tshead@k-3d.com
//
// This library 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.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

/** \file
		\brief 2D geometry routines implementation
		\author Tim Shead (tshead@k-3d.com)
*/

#include "basic_math.h"
#include "geometry.h"

namespace k3d
{

double Area(const k3d::vector2& Point1, const k3d::vector2& Point2, const k3d::vector2& Point3)
{
	double result = 0;

	result += (Point1[0] * Point2[1]) - (Point2[0] * Point1[1]);
	result += (Point2[0] * Point3[1]) - (Point3[0] * Point2[1]);
	result += (Point3[0] * Point1[1]) - (Point1[0] * Point3[1]);

	// For robustness ...
	//if(std::fabs(result) < FLT_EPSILON)
	//	return 0;

	return result * 0.5;
}

/////////////////////////////////////////////////////////////////////////////
// rectangle implementation

rectangle::rectangle()
{ left = right = top = bottom = 0.0; }

rectangle::rectangle(const double Left, const double Right, const double Top, const double Bottom)
{ left = Left; right = Right; top = Top; bottom = Bottom; }

rectangle::rectangle(const rectangle& v)
{ left = v.left; right = v.right; top = v.top; bottom = v.bottom; }

rectangle& rectangle::operator = (const rectangle& v)
{
	left = v.left; right = v.right; top = v.top; bottom = v.bottom;
	return *this;
}

double rectangle::Left() const { return left; }
double rectangle::Right() const { return right; }
double rectangle::Top() const { return top; }
double rectangle::Bottom() const { return bottom; }
double rectangle::Width() const { return std::fabs(right - left); }
double rectangle::Height() const { return std::fabs(bottom - top); }

void rectangle::Offset(const double Horizontal, const double Vertical)
{
	left += Horizontal;
	right += Horizontal;
	top += Vertical;
	bottom += Vertical;
}

void rectangle::Grow(const double Thickness)
{
	Grow(Thickness, Thickness);
}

void rectangle::Grow(const double Horizontal, const double Vertical)
{
	left -= Horizontal;
	right += Horizontal;
	top -= Vertical;
	bottom += Vertical;
}

void rectangle::Shrink(const double Thickness)
{
	Shrink(Thickness, Thickness);
}

void rectangle::Shrink(const double Horizontal, const double Vertical)
{
	Grow(-Horizontal, -Vertical);
}

void rectangle::Normalize()
{
	if(left > right)
		std::swap(left, right);

	if(top > bottom)
		std::swap(top, bottom);
}

vector2 rectangle::Center() const
{
	return vector2(left+right, top+bottom) * 0.5;
}

bool rectangle::Contains(const vector2& Point) const
{
	return (left <= Point[0] && Point[0] <= right && top <= Point[1] && Point[1] <= bottom) ? true : false;
}

rectangle ConstrainProportions(const rectangle& Source, const rectangle& Destination)
{
	rectangle source(Source);
	rectangle destination(Destination);

	source.Normalize();
	destination.Normalize();

	source.Offset(-source.Left(), -source.Top());

	if(source.Width() != destination.Width())
	{
		source.bottom *= destination.Width() / source.Width();
		source.right *= destination.Width() / source.Width();
	}

	if(source.Height() > destination.Height())
	{
		source.right *= destination.Height() / source.Height();
		source.bottom *= destination.Height() / source.Height();
	}

	source.Offset(destination.Left() + (0.5 * (destination.Width() - source.Width())),
			destination.Top() + (0.5 * (destination.Height() - source.Height())));

	return source;
}

std::ostream& operator<<(std::ostream& Stream, const rectangle& Arg)
{
	Stream << Arg.Left() << " " << Arg.Right() << " " << Arg.Top() << " " << Arg.Bottom();
	return Stream;
}

std::istream& operator>>(std::istream& Stream, rectangle& Arg)
{
	Stream >> Arg.left >> Arg.right >> Arg.top >> Arg.bottom;
	return Stream;
}

/////////////////////////////////////////////////////////////////////////////
// plane implementation

plane::plane(const vector3& Normal, const double Distance)
{
	m_Normal = Normal;
	m_Normal.Normalize();

	m_Distance = Distance;
}

plane::plane(const vector3& Normal, const vector3& ContainedPoint)
{
	m_Normal = Normal;
	m_Normal.Normalize();

	m_Distance = -ContainedPoint * m_Normal;
}

plane::plane(const vector3& ContainedPointA, const vector3& ContainedPointB, const vector3& ContainedPointC)
{
	m_Normal = ((ContainedPointA - ContainedPointB) ^ (ContainedPointC - ContainedPointB)).Normalize();
	m_Distance = -ContainedPointA * m_Normal;
}

vector3 plane::Normal() const
{
	return m_Normal;
}

double plane::Distance() const
{
	return m_Distance;
}

void plane::SetNormal(const vector3& Normal)
{
	m_Normal = Normal;
}

void plane::SetDistance(const double Distance)
{
	m_Distance = Distance;
}

std::ostream& operator<<(std::ostream& Stream, const plane& Arg)
{
	Stream << Arg.Normal() << " " << Arg.Distance();
	return Stream;
}

std::istream& operator>>(std::istream& Stream, plane& Arg)
{
	Stream >> Arg.m_Normal >> Arg.m_Distance;
	return Stream;
}

///
bool PlaneLineIntersection(const plane& Plane, const k3d::vector3 LineOrigin, const k3d::vector3 LineDirection, k3d::vector3& Result)
{
	// Calculate the angle (dot product) between line and plane ...
	const double vd = Plane.Normal() * LineDirection;

	// Make sure the line & plane aren't parallel ...
	if(0.0 == vd)
		return false;

	const double v0 = -(Plane.Normal() * LineOrigin + Plane.Distance());
	const double t = v0 / vd;

	Result = LineOrigin + (t * LineDirection);

	return true;
}

bool LineIntersection(const k3d::vector2& A1, const k3d::vector2& A2, const k3d::vector2& B1, const k3d::vector2& B2, k3d::vector2& Result)
{
	const k3d::vector2 a = (A2 - A1).Normalize();
	const k3d::vector2 b = (B2 - B1).Normalize();
	const k3d::vector2 c = B1 - A1;
	const k3d::vector2 bperp = b.Perpendicular();

	const double denominator = bperp * a;
	if(0 == denominator)
		return false;

	const double t = (bperp * c) / denominator;

	Result = A1 + (t * a);

	return true;
}

line_segment_intersection_type LineSegmentIntersection(const k3d::vector2& A1, const k3d::vector2& A2, const k3d::vector2& B1, const k3d::vector2& B2)
{
	const double sign1 = k3d::sign(k3d::Area(A1, A2, B1));
	const double sign2 = k3d::sign(k3d::Area(A1, A2, B2));

	// The three most likely possibilities ...
	if(sign1 && sign2)
	{
		// No intersection ...
		if(sign1 == sign2)
			return SDP_NOINTERSECTION;

		k3d::vector2 intersection;
		if(!LineIntersection(A1, A2, B1, B2, intersection))
			{
				std::cerr << __PRETTY_FUNCTION__ << " : LineIntersection 1 failed." << std::endl;
				return SDP_ERROR;
			}

		const double sign3 = k3d::sign(k3d::Area(B1, intersection, A1));
		const double sign4 = k3d::sign(k3d::Area(B1, intersection, A2));

		if(sign3 != sign4)
			return SDP_INTERIORINTERSECTION;

		return SDP_EXTERIORINTERSECTION;
	}
	// Two variations with single-point intersections ...
	else if(sign1 || sign2)
	{
		k3d::vector2 intersection;
		if(!LineIntersection(A1, A2, B1, B2, intersection))
			{
				std::cerr << __PRETTY_FUNCTION__ << " : LineIntersection 2 failed." << std::endl;
				return SDP_ERROR;
			}

		const double sign3 = k3d::sign(k3d::Area(B1, intersection, A1));
		const double sign4 = k3d::sign(k3d::Area(B1, intersection, A2));

		if(sign3 != sign4)
			return SDP_INTERIORCOLINEARPOINT;

		return SDP_EXTERIORCOLINEARPOINT;
	}
	// The four colinear possibilities ...
	else
	{
		// This ray is getting crowded!  Find something somewhere else
		k3d::vector2 start = B1 + (B2 - B1).Perpendicular();

		const double sign3 = k3d::sign(k3d::Area(start, B1, A1));
		const double sign4 = k3d::sign(k3d::Area(start, B1, A2));
		const double sign5 = k3d::sign(k3d::Area(start, B2, A1));
		const double sign6 = k3d::sign(k3d::Area(start, B2, A2));

		if((sign3 != sign4) && (sign5 != sign6))
			return SDP_INTERIORCOLINEARPOINTS;

		if((sign3 == sign4) && (sign4 == sign5) && (sign5 == sign6))
			return SDP_EXTERIORCOLINEARPOINTSNOINTERSECTION;

		if((sign3 == sign4) && (sign4 != sign5) && (sign5 == sign6))
			return SDP_EXTERIORCOLINEARPOINTS;

		return SDP_INTERIOREXTERIORCOLINEARPOINTS;
	}
}

double BernsteinBasis(const unsigned long Order, const unsigned long ControlPoint, const double Parameter)
{
	// Sanity checks ...
	assert(Order > 1);
	assert(ControlPoint <= Order);

	const unsigned long n = Order - 1;
	const unsigned long i = ControlPoint;
	const double t = Parameter;

	const double ni = factorial(n) / (factorial(i) * factorial(n - i));

	return ni * pow(t, i) * pow((1 - t), (n - i));
}

bool intersect_lines(const k3d::vector3& P1, const k3d::vector3& T1, const k3d::vector3& P2, const k3d::vector3& T2, k3d::vector3& Result)
{
	// Code originally from Aqsis, http://www.aqsis.com

    k3d::vector3	v, px;

    px = T1 ^ ( P1 - T2 );
    v = px ^ T1;

    double	t = ( P1 - P2 ) * v;
    double vw = v * T2;
    if ( ( vw * vw ) < 1.0e-07 )
        return false;
    t /= vw;
    Result = P2 + ( ( ( P1 - P2 ) * v ) / vw ) * T2 ;
    return true;
}

} // namespace k3d