]> Shamusworld >> Repos - architektonas/blobdiff - src/geometry.cpp
Added preliminary circle-to-circle intersection code.
[architektonas] / src / geometry.cpp
index eacaec7feabc20ee4b8cddce4a15e538f2905bd4..6ed6fab1020c70449e64955f4bc9c6c0548f28b8 100644 (file)
 
 #include "geometry.h"
 #include <math.h>
-#include "line.h"
-#include "circle.h"
+#include <stdio.h>
+#include "global.h"
+#include "mathconstants.h"
 
 
+// This is unused
+#if 0
 Point Geometry::IntersectionOfLineAndLine(Point p1, Point p2, Point p3, Point p4)
 {
        // Find the intersection of the lines by formula:
@@ -40,6 +43,7 @@ Point Geometry::IntersectionOfLineAndLine(Point p1, Point p2, Point p3, Point p4
 
        return Point(px / d, py / d, 0);
 }
+#endif
 
 
 // Returns the parameter of a point in space to this vector. If the parameter
@@ -61,15 +65,15 @@ double Geometry::ParameterOfLineAndPoint(Point tail, Point head, Point point)
 }
 
 
-Point Geometry::MirrorPointAroundLine(Point point, Point p1, Point p2)
+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);
@@ -81,9 +85,13 @@ 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));
 
@@ -97,6 +105,17 @@ double Geometry::Determinant(Point p1, Point p2)
 }
 
 
+void Geometry::Intersects(Object * obj1, Object * obj2)//, double * tp/*= 0*/, double * up/*= 0*/, double * vp/*= 0*/, double * wp/*= 0*/)
+{
+       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);
+}
+
+
 /*
 Intersecting line segments:
 An easier way:
@@ -116,30 +135,68 @@ So check if the above two numbers are both >=0 and <=1.
 */
 
 
-#if 0
-// Finds the intesection between two objects (if any)
-bool Geometry::Intersects(Object * obj1, Object * obj2, double * t, double * s)
-{
-}
-#endif
-
 // Finds the intersection between two lines (if any)
-int Geometry::Intersects(Line * l1, Line * l2, double * tp/*= 0*/, double * up/*= 0*/)
+void Geometry::CheckLineToLineIntersection(Object * l1, Object * l2)//, double * tp, double * up)
 {
-       Vector r(l1->position, l1->endpoint);
-       Vector s(l2->position, l2->endpoint);
-       Vector v1 = l2->position - l1->position;
-//     Vector v1 = l1->position - l2->position;
+       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
+#if 0
+       Vector v2 = l1->p[0] - l2->p[0];        // p - q
+printf("l1: (%lf, %lf) (%lf, %lf), l2: (%lf, %lf) (%lf, %lf)\n", l1->p[0].x, l1->p[0].y, l1->p[1].x, l1->p[1].y, l2->p[0].x, l2->p[0].y, l2->p[1].x, l2->p[1].y);
+#endif
        double rxs = (r.x * s.y) - (s.x * r.y);
+       double t, u;
 
        if (rxs == 0)
-               return 0;
+       {
+               double qpxr = (v1.x * r.y) - (r.x * v1.y);
 
-       double t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
-       double u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
+#if 0
+printf("  --> R x S = 0! (q - p) x r = %lf\n", qpxr);
+printf("  -->(q - p) . r = %lf, r . r = %lf\n", v1.Dot(r), r.Dot(r));
+printf("  -->(p - q) . s = %lf, s . s = %lf\n", v2.Dot(s), s.Dot(s));
+printf("  -->(q - p) . s = %lf, (p - q) . r = %lf\n", v1.Dot(s), v2.Dot(r));
+#endif
+
+               // Lines are parallel, so no intersection...
+               if (qpxr != 0)
+                       return;
+
+#if 0
+//this works IFF the vectors are pointing in the same direction. everything else
+//is fucked!
+               // If (q - p) . r == r . r, t = 1, u = 0
+               if (v1.Dot(r) == r.Dot(r))
+                       t = 1.0, u = 0;
+               // If (p - q) . s == s . s, t = 0, u = 1
+               else if (v2.Dot(s) == s.Dot(s))
+                       t = 0, u = 1.0;
+               else
+                       return 0;
+#else
+               // 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;
+#endif
+       }
+       else
+       {
+               t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
+               u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
+       }
 /*
-Now there are five cases:
+Now there are five cases (NOTE: only valid if vectors face the same way!):
 
 1. If r × s = 0 and (q − p) × r = 0, then the two lines are collinear. If in addition, either 0 ≤ (q − p) · r ≤ r · r or 0 ≤ (p − q) · s ≤ s · s, then the two lines are overlapping.
 
@@ -152,6 +209,7 @@ Now there are five cases:
 5. Otherwise, the two line segments are not parallel but do not intersect.
 */
        // Return parameter values, if user passed in valid pointers
+#if 0
        if (tp)
                *tp = t;
 
@@ -163,10 +221,69 @@ Now there are five cases:
                return 1;
 
        return 0;
+#else
+       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;
+#endif
+}
+
+
+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]))
+               || (d == fabs(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;
+       }
+
+       // 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;
+}
+
+
+#if 0
+// Finds the intersection between two lines (if any)
+int Geometry::Intersects(Line * l1, Dimension * d1, double * tp/*= 0*/, double * up/*= 0*/)
+{
+       Line l2(d1->position, d1->endpoint);
+       return Intersects(l1, &l2, tp, up);
 }
 
 
-// Finds the intesection(s) between a line and a circle (if any)
+// Finds the intersection(s) between a line and a circle (if any)
 int Geometry::Intersects(Line * l, Circle * c, double * tp/*= 0*/, double * up/*= 0*/, double * vp/*= 0*/, double * wp/*= 0*/)
 {
 #if 0
@@ -214,7 +331,16 @@ I'm thinking a better approach to this might be as follows:
        // Now we have to check for intersection points.
        // Tangent case: (needs to return something)
        if ((distance == c->radius) && (t >= 0.0) && (t <= 1.0))
+       {
+               // Need to set tp & up to something... !!! FIX !!!
+               if (tp)
+                       *tp = t;
+
+               if (up)
+                       *up = Vector(c->position, p).Angle();
+
                return 1;
+       }
 
        // The line intersects the circle in two points (possibly). Use Pythagorus
        // to find them for testing.
@@ -251,3 +377,144 @@ I'm thinking a better approach to this might be as follows:
 }
 
 
+// Finds the intersection(s) between a circle and a circle (if any)
+// There can be 0, 1, or 2 intersections.
+// Returns the angles of the points of intersection in tp thru wp, with the
+// angles returned as c1, c2, c1, c2 (if applicable--in the 1 intersection case,
+// only the first two angles are returned: c1, c2).
+int Geometry::Intersects(Circle * c1, Circle * c2, double * tp/*= 0*/, double * up/*= 0*/, double * vp/*= 0*/, double * wp/*= 0*/, Point * p1/*= 0*/, Point * p2/*= 0*/)
+{
+       // Get the distance between centers. If the distance plus the radius of the
+       // smaller circle is less than the radius of the larger circle, there is no
+       // intersection. If the distance is greater than the sum of the radii,
+       // there is no intersection. If the distance is equal to the sum of the
+       // radii, they are tangent and intersect at one point. Otherwise, they
+       // intersect at two points.
+       Vector centerLine(c1->position, c2->position);
+       double d = centerLine.Magnitude();
+//printf("Circle #1: pos=<%lf, %lf>, r=%lf\n", c1->position.x, c1->position.y, c1->radius);
+//printf("Circle #2: pos=<%lf, %lf>, r=%lf\n", c2->position.x, c2->position.y, c2->radius);
+//printf("Distance between #1 & #2: %lf\n", d);
+
+       // Check to see if we actually have an intersection, and return failure if not
+       if ((fabs(c1->radius - c2->radius) > d) || ((c1->radius + c2->radius) < d))
+               return 0;
+
+       // There are *two* tangent cases!
+       if (((c1->radius + c2->radius) == d) || (fabs(c1->radius - c2->radius) == d))
+       {
+               // Need to return something in tp & up!! !!! FIX !!! [DONE]
+               if (tp)
+                       *tp = centerLine.Angle();
+
+               if (up)
+                       *up = centerLine.Angle() + PI;
+
+               return 1;
+       }
+
+       // Find the distance from the center of c1 to the perpendicular chord
+       // (which contains the points of intersection)
+       // [N.B.: This is derived from Pythagorus by using the unknown distance
+       //        from the center line to the point where the two radii coincide as
+       //        a common unknown to two instances of the formula.]
+       double x = ((d * d) - (c2->radius * c2->radius) + (c1->radius * c1->radius))
+               / (2.0 * d);
+       // Find the the length of the perpendicular chord
+// Not needed...!
+       double a = sqrt((-d + c2->radius - c1->radius) * (-d - c2->radius + c1->radius) * (-d + c2->radius + c1->radius) * (d + c2->radius + c1->radius)) / d;
+
+       // Now, you can use pythagorus to find the length of the hypotenuse, but we
+       // already know that length, it's the radius! :-P
+       // What's needed is the angle of the center line and the radial line. Since
+       // there's two intersection points, there's also four angles (two for each
+       // circle)!
+       // We can use the arccos to find the angle using just the radius and the
+       // distance to the perpendicular chord...!
+       double angleC1 = acos(x / c1->radius);
+       double angleC2 = acos((d - x) / c2->radius);
+
+       if (tp)
+               *tp = centerLine.Angle() - angleC1;
+
+       if (up)
+               *up = (centerLine.Angle() + PI) - angleC2;
+
+       if (vp)
+               *vp =  centerLine.Angle() + angleC1;
+
+       if (wp)
+               *wp = (centerLine.Angle() + PI) + angleC2;
+
+       if (p1)
+               *p1 = c1->position + (centerLine.Unit() * x) + (Vector::Normal(Vector(), centerLine) * (a / 2.0));
+
+       if (p2)
+               *p2 = c1->position + (centerLine.Unit() * x) - (Vector::Normal(Vector(), centerLine) * (a / 2.0));
+
+       return 2;
+}
+#endif
+
+// 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);
+}
+