X-Git-Url: http://shamusworld.gotdns.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fgeometry.cpp;h=95d12c65068b16715086191e1d70c7e01483dcbd;hb=bf5a50feb0f84a4627a65c5b82c3ca2d2eefe54b;hp=5f9a30b964f07348db24af44482508c52bacf668;hpb=921bf050ffe5fc81a9ab377e634180e659ee5d5d;p=architektonas

diff --git a/src/geometry.cpp b/src/geometry.cpp
index 5f9a30b..95d12c6 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -15,33 +15,15 @@
 
 #include "geometry.h"
 #include <math.h>
-
-Point Geometry::IntersectionOfLineAndLine(Point p1, Point p2, Point p3, Point p4)
-{
-	// Find the intersection of the lines by formula:
-	// px = (x1y2 - y1x2)(x3 - x4) - (x1 - x2)(x3y4 - y3x4)
-	// py = (x1y2 - y1x2)(y3 - y4) - (y1 - y2)(x3y4 - y3x4)
-	// d = (x1 - x2)(y3 - y4) - (y1 - y2)(x3 - x4) = 0 if lines are parallel
-	// Intersection is (px / d, py / d)
-
-	double d = ((p1.x - p2.x) * (p3.y - p4.y)) - ((p1.y - p2.y) * (p3.x - p4.x));
-
-	// Check for parallel lines, and return sentinel if so
-	if (d == 0)
-		return Point(0, 0, -1);
-
-	double px = (((p1.x * p2.y) - (p1.y * p2.x)) * (p3.x - p4.x))
-		- ((p1.x - p2.x) * ((p3.x * p4.y) - (p3.y * p4.x)));
-	double py = (((p1.x * p2.y) - (p1.y * p2.x)) * (p3.y - p4.y))
-		- ((p1.y - p2.y) * ((p3.x * p4.y) - (p3.y * p4.x)));
-
-	return Point(px / d, py / d, 0);
-}
+#include <stdio.h>
+#include "global.h"
+#include "mathconstants.h"
 
 
 // Returns the parameter of a point in space to this vector. If the parameter
 // is between 0 and 1, the normal of the vector to the point is on the vector.
-double Geometry::ParameterOfLineAndPoint(Point lp1, Point lp2, Point point)
+// Note: lp1 is the tail, lp2 is the head of the line (vector).
+double Geometry::ParameterOfLineAndPoint(Point tail, Point head, Point point)
 {
 	// Geometric interpretation:
 	// The parameterized point on the vector lineSegment is where the normal of
@@ -49,23 +31,46 @@ double Geometry::ParameterOfLineAndPoint(Point lp1, Point lp2, Point point)
 	// the perpendicular lies beyond the 1st endpoint. If pp > 1, then the
 	// perpendicular lies beyond the 2nd endpoint.
 
-	Vector lineSegment = lp1 - lp2;
+	Vector lineSegment = head - tail;
 	double magnitude = lineSegment.Magnitude();
-	Vector pointSegment = point - lp2;
+	Vector pointSegment = point - tail;
 	double t = lineSegment.Dot(pointSegment) / (magnitude * magnitude);
 	return t;
 }
 
 
-Point Geometry::MirrorPointAroundLine(Point point, Point p1, Point p2)
+double Geometry::DistanceToLineFromPoint(Point tail, Point head, Point point)
+{
+	// Interpretation: given a line in the form x = a + tu, where u is the
+	// unit vector of the line, a is the tail and t is a parameter which
+	// describes the line, the distance of a point p to the line is given by:
+	// || (a - p) - ((a - p) . u) u ||
+	// We go an extra step: we set the sign to reflect which side of the line
+	// it's on (+ == to the left if head points away from you, - == to the
+	// right)
+	Vector line(tail, head);
+	Vector u = line.Unit();
+	Vector a_p = tail - point;
+	Vector dist = a_p - (u * (a_p).Dot(u));
+
+	double angle = Vector::Angle(tail, point) - line.Angle();
+
+	if (angle < 0)
+		angle += TAU;
+
+	return dist.Magnitude() * (angle < HALF_TAU ? +1.0 : -1.0);
+}
+
+
+Point Geometry::MirrorPointAroundLine(Point point, Point tail, Point head)
 {
 	// Get the vector of the intersection of the line and the normal on the
 	// line to the point in question.
-	double t = ParameterOfLineAndPoint(p1, p2, point);
-	Vector v = Vector(p1, p2) * t;
+	double t = ParameterOfLineAndPoint(tail, head, point);
+	Vector v = Vector(tail, head) * t;
 
-	// Get the point normal to point to the line passed in (p2 is the tail)
-	Point normalOnLine = p2 + v;
+	// Get the point normal to point to the line passed in
+	Point normalOnLine = tail + v;
 
 	// Make our mirrored vector (head - tail)
 	Vector mirror = -(point - normalOnLine);
@@ -77,12 +82,303 @@ Point Geometry::MirrorPointAroundLine(Point point, Point p1, Point p2)
 }
 
 
+//
+// point: The point we're rotating
+// rotationPoint: The point we're rotating around
+//
 Point Geometry::RotatePointAroundPoint(Point point, Point rotationPoint, double angle)
 {
-	Vector v = Vector(point, rotationPoint);
+	Vector v = Vector(rotationPoint, point);
 	double px = (v.x * cos(angle)) - (v.y * sin(angle));
 	double py = (v.x * sin(angle)) + (v.y * cos(angle));
 
 	return Vector(rotationPoint.x + px, rotationPoint.y + py, 0);
 }
 
+
+double Geometry::Determinant(Point p1, Point p2)
+{
+	return (p1.x * p2.y) - (p2.x * p1.y);
+}
+
+
+void Geometry::Intersects(Object * obj1, Object * obj2)
+{
+	Global::numIntersectPoints = Global::numIntersectParams = 0;
+
+	if ((obj1->type == OTLine) && (obj2->type == OTLine))
+		CheckLineToLineIntersection(obj1, obj2);
+	else if ((obj1->type == OTCircle) && (obj2->type == OTCircle))
+		CheckCircleToCircleIntersection(obj1, obj2);
+	else if ((obj1->type == OTLine) && (obj2->type == OTCircle))
+		CheckLineToCircleIntersection(obj1, obj2);
+	else if ((obj1->type == OTCircle) && (obj2->type == OTLine))
+		CheckLineToCircleIntersection(obj2, obj1);
+}
+
+
+/*
+Intersecting line segments:
+An easier way:
+Segment L1 has edges A=(a1,a2), A'=(a1',a2').
+Segment L2 has edges B=(b1,b2), B'=(b1',b2').
+Segment L1 is the set of points tA'+(1-t)A, where 0<=t<=1.
+Segment L2 is the set of points sB'+(1-s)B, where 0<=s<=1.
+Segment L1 meet segment L2 if and only if for some t and s we have
+tA'+(1-t)A=sB'+(1-s)B
+The solution of this with respect to t and s is
+
+t=((-b?'a?+b?'b?+b?a?+a?b?'-a?b?-b?b?')/(b?'a?'-b?'a?-b?a?'+b?a?-a?'b?'+a?'b?+a?b?'-a?b?))
+
+s=((-a?b?+a?'b?-a?a?'+b?a?+a?'a?-b?a?')/(b?'a?'-b?'a?-b?a?'+b?a?-a?'b??+a?'b?+a?b?'-a?b?))
+
+So check if the above two numbers are both >=0 and <=1.
+*/
+
+
+// Finds the intersection between two lines (if any)
+void Geometry::CheckLineToLineIntersection(Object * l1, Object * l2)
+{
+	Global::numIntersectPoints = Global::numIntersectParams = 0;
+
+	Vector r(l1->p[0], l1->p[1]);
+	Vector s(l2->p[0], l2->p[1]);
+	Vector v1 = l2->p[0] - l1->p[0];	// q - p
+
+	double rxs = (r.x * s.y) - (s.x * r.y);
+	double t, u;
+
+	if (rxs == 0)
+	{
+		double qpxr = (v1.x * r.y) - (r.x * v1.y);
+
+		// Lines are parallel, so no intersection...
+		if (qpxr != 0)
+			return;
+
+		// Check to see which endpoints are connected... Four possibilities:
+		if (l1->p[0] == l2->p[0])
+			t = 0, u = 0;
+		else if (l1->p[0] == l2->p[1])
+			t = 0, u = 1.0;
+		else if (l1->p[1] == l2->p[0])
+			t = 1.0, u = 0;
+		else if (l1->p[1] == l2->p[1])
+			t = 1.0, u = 1.0;
+		else
+			return;
+	}
+	else
+	{
+		t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
+		u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
+	}
+
+	Global::intersectParam[0] = t;
+	Global::intersectParam[1] = u;
+
+	// If the parameters are in range, we have overlap!
+	if ((t >= 0) && (t <= 1.0) && (u >= 0) && (u <= 1.0))
+		Global::numIntersectParams = 1;
+}
+
+
+void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2)
+{
+	// Set up global vars
+	Global::numIntersectPoints = Global::numIntersectParams = 0;
+
+	// Get the distance between the centers of the circles
+	Vector centerLine(c1->p[0], c2->p[0]);
+	double d = centerLine.Magnitude();
+	double clAngle = centerLine.Angle();
+
+	// If the distance between centers is greater than the sum of the radii or
+	// less than the difference between the radii, there is NO intersection
+	if ((d > (c1->radius[0] + c2->radius[0]))
+		|| (d < fabs(c1->radius[0] - c2->radius[0])))
+		return;
+
+	// If the distance between centers is equal to the sum of the radii or
+	// equal to the difference between the radii, the intersection is tangent
+	// to both circles.
+	if (d == (c1->radius[0] + c2->radius[0]))
+	{
+		Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0]);
+		Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0]);
+		Global::numIntersectPoints = 1;
+		return;
+	}
+	else if (d == fabs(c1->radius[0] - c2->radius[0]))
+	{
+		double sign = (c1->radius[0] > c2->radius[0] ? +1 : -1);
+		Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0] * sign);
+		Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0] * sign);
+		Global::numIntersectPoints = 1;
+		return;
+	}
+
+/*
+       c² = a² + b² - 2ab·cos µ
+2ab·cos µ = a² + b² - c²
+    cos µ = (a² + b² - c²) / 2ab
+*/
+	// Use the Law of Cosines to find the angle between the centerline and the
+	// radial line on Circle #1
+	double a = acos(((c1->radius[0] * c1->radius[0]) + (d * d) - (c2->radius[0] * c2->radius[0])) / (2.0 * c1->radius[0] * d));
+
+	// Finally, find the points of intersection by using +/- the angle found
+	// from the centerline's angle
+	Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle + a) * c1->radius[0]);
+	Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle + a) * c1->radius[0]);
+	Global::intersectPoint[1].x = c1->p[0].x + (cos(clAngle - a) * c1->radius[0]);
+	Global::intersectPoint[1].y = c1->p[0].y + (sin(clAngle - a) * c1->radius[0]);
+	Global::numIntersectPoints = 2;
+}
+
+
+//
+// N.B.: l is the line, c is the circle
+//
+void Geometry::CheckLineToCircleIntersection(Object * l, Object * c)
+{
+	// Set up global vars
+	Global::numIntersectPoints = Global::numIntersectParams = 0;
+
+	// Step 1: Find shortest distance from center of circle to the infinite line
+	double t = ParameterOfLineAndPoint(l->p[0], l->p[1], c->p[0]);
+	Point p = l->p[0] + (Vector(l->p[0], l->p[1]) * t);
+	Vector radial = Vector(c->p[0], p);
+	double distance = radial.Magnitude();
+
+	// Step 2: See if we have 0, 1, or 2 intersection points
+
+	// Case #1: No intersection points
+	if (distance > c->radius[0])
+		return;
+	// Case #2: One intersection point (possibly--tangent)
+	else if (distance == c->radius[0])
+	{
+		// Only intersects if the parameter is on the line segment!
+		if ((t >= 0.0) && (t <= 1.0))
+		{
+			Global::intersectPoint[0] = c->p[0] + radial;
+			Global::numIntersectPoints = 1;
+		}
+
+		return;
+	}
+
+	// Case #3: Two intersection points (possibly--secant)
+
+	// So, we have the line, and the perpendicular from the center of the
+	// circle to the line. Now figure out where the intersection points are.
+	// This is a right triangle, though do we really know all the sides?
+	// Don't need to, 2 is enough for Pythagoras :-)
+	// Radius is the hypotenuse, so we have to use c² = a² + b² => a² = c² - b²
+	double perpendicularLength = sqrt((c->radius[0] * c->radius[0]) - (distance * distance));
+
+	// Now, find the intersection points using the length...
+	Vector lineUnit = Vector(l->p[0], l->p[1]).Unit();
+	Point i1 = p + (lineUnit * perpendicularLength);
+	Point i2 = p - (lineUnit * perpendicularLength);
+
+	// Next we need to see if they are on the line segment...
+	double u = ParameterOfLineAndPoint(l->p[0], l->p[1], i1);
+	double v = ParameterOfLineAndPoint(l->p[0], l->p[1], i2);
+
+	if ((u >= 0.0) && (u <= 1.0))
+	{
+		Global::intersectPoint[Global::numIntersectPoints] = i1;
+		Global::numIntersectPoints++;
+	}
+
+	if ((v >= 0.0) && (v <= 1.0))
+	{
+		Global::intersectPoint[Global::numIntersectPoints] = i2;
+		Global::numIntersectPoints++;
+	}
+}
+
+
+// should we just do common trig solves, like AAS, ASA, SAS, SSA?
+// Law of Cosines:
+// c² = a² + b² - 2ab * cos(C)
+// Solving for C:
+// cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab
+// Law of Sines:
+// a / sin A = b / sin B = c / sin C
+
+// Solve the angles of the triangle given the sides. Angles returned are
+// opposite of the given sides (so a1 consists of sides s2 & s3, and so on).
+void Geometry::FindAnglesForSides(double s1, double s2, double s3, double * a1, double * a2, double * a3)
+{
+	// Use law of cosines to find 1st angle
+	double cosine1 = ((s2 * s2) + (s3 * s3) - (s1 * s1)) / (2.0 * s2 * s3);
+
+	// Check for a valid triangle
+	if ((cosine1 < -1.0) || (cosine1 > 1.0))
+		return;
+
+	double angle1 = acos(cosine1);
+
+	// Use law of sines to find 2nd & 3rd angles
+// sin A / a = sin B / b
+// sin B = (sin A / a) * b
+// B = arcsin( sin A * (b / a))
+// ??? ==> B = A * arcsin(b / a)
+/*
+Well, look here:
+sin B = sin A * (b / a)
+sin B / sin A = b / a
+arcsin( sin B / sin A ) = arcsin( b / a )
+
+hmm... dunno...
+*/
+
+	double angle2 = asin(s2 * (sin(angle1) / s1));
+	double angle3 = asin(s3 * (sin(angle1) / s1));
+
+	if (a1)
+		*a1 = angle1;
+
+	if (a2)
+		*a2 = angle2;
+
+	if (a3)
+		*a3 = angle3;
+}
+
+
+Point Geometry::GetPointForParameter(Object * obj, double t)
+{
+	if (obj->type == OTLine)
+	{
+		// Translate line vector to the origin, then add the scaled vector to
+		// initial point of the line.
+		Vector v = obj->p[1] - obj->p[0];
+		return obj->p[0] + (v * t);
+	}
+
+	return Point(0, 0);
+}
+
+
+Point Geometry::Midpoint(Line * line)
+{
+	return Point((line->p[0].x + line->p[1].x) / 2.0,
+		(line->p[0].y + line->p[1].y) / 2.0);
+}
+
+
+/*
+How to find the tangent of a point off a circle:
+
+ •  Calculate the midpoint on the point and the center of the circle
+ •  Get the length of the line segment from the and the center divided by two
+ •  Use that length to construct a circle with the point at the center and the
+    radius equal to that length
+ •  The intersection of the two circles are the tangent points
+
+*/
+