+//
// geometry.cpp: Algebraic geometry helper functions
//
// Part of the Architektonas Project
#include "geometry.h"
#include <math.h>
-#include "circle.h"
-#include "dimension.h"
-#include "line.h"
+#include <stdio.h>
+#include "global.h"
#include "mathconstants.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);
-}
-
-
// 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.
// Note: lp1 is the tail, lp2 is the head of the line (vector).
double magnitude = lineSegment.Magnitude();
Vector pointSegment = point - tail;
double t = lineSegment.Dot(pointSegment) / (magnitude * magnitude);
+
return t;
}
+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)
{
return mirroredPoint;
}
-
+//
+// 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);
+ 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:
*/
-#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)
{
- Vector r(l1->position, l1->endpoint);
- Vector s(l2->position, l2->endpoint);
- Vector v1 = l2->position - l1->position; // q - p
-// Vector v2 = l1->position - l2->position; // p - q
-//printf("l1: (%lf, %lf) (%lf, %lf), l2: (%lf, %lf) (%lf, %lf)\n", l1->position.x, l1->position.y, l1->endpoint.x, l1->endpoint.y, l2->position.x, l2->position.y, l2->endpoint.x, l2->endpoint.y);
+ 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;
+ // The angle is zero, so the lines probably overlap, or are parallel. Either way, there is either INFINITE intersections, or zero. Either way, the result is useless to us. It does present a bit of an inconsistency though: lines connected at their endpoints that aren't at a zero angle will show as overlapping (the only case this could even REMOTELY be useful is when the angle between is 180°--not zero, as zero is the degenerate case). :-P
if (rxs == 0)
{
- double qpxr = (v1.x * r.y) - (r.x * v1.y);
-
-//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));
+// WHY would you do this??? The lines overlap at an INFINITE number of points!
+// The assumption is that the angle is 180°, not 0°.
+ // If we want to detect the degenerate case, we can check the parameter of the endpoints of the one line to the other. If any of the parameters are in (0, 1), then it's the degenerate case (we would check the endpoints of the shorter segment against the longer).
+/* double qpxr = (v1.x * r.y) - (r.x * v1.y);
- // Lines are parallel, so no intersection...
+ // Line segments are parallel, so no intersection...
if (qpxr != 0)
- return 0;
+ return;
+
+ // Otherwise, the segments are colinear. Need to check for the 0° (degenerate) vs the 180° (OK) case.
+ Object * larger = l1;
+ Object * smaller = l2;
+
+ if (r->Magnitude() < s->Magnitude())
+ {
+ larger = l2;
+ smaller = l1;
+ }
+
+ double param1 = ParameterOfLineAndPoint(larger->p[0], larger->p[1], smaller->p[0]);
+ double param2 = ParameterOfLineAndPoint(larger->p[0], larger->p[1], smaller->p[1]);
+
+ // Check for the degenerate case, and return if found
+ if ((param1 > 0 && param1 < 1.0) || (param2 > 0 && param2 < 1.0))
+ return;
+
+//// or just use AngleBetween: Nah, won't work...
-#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->position == l2->position)
+ if (l1->p[0] == l2->p[0])
t = 0, u = 0;
- else if (l1->position == l2->endpoint)
+ else if (l1->p[0] == l2->p[1])
t = 0, u = 1.0;
- else if (l1->endpoint == l2->position)
+ else if (l1->p[1] == l2->p[0])
t = 1.0, u = 0;
- else if (l1->endpoint == l2->endpoint)
+ else if (l1->p[1] == l2->p[1])
t = 1.0, u = 1.0;
- else
- return 0;
-#endif
+ else*/
+ return;
}
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 (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.
-
-2. If r × s = 0 and (q − p) × r = 0, but neither 0 ≤ (q − p) · r ≤ r · r nor 0 ≤ (p − q) · s ≤ s · s, then the two lines are collinear but disjoint.
-3. If r × s = 0 and (q − p) × r ≠ 0, then the two lines are parallel and non-intersecting.
+ // Check that the parameters are above the epsilon, otherwise clamp them to
+ // zero or one, as the case may be
+ if (fabs(t) < EPSILON)
+ t = 0;
+ else if (fabs(1.0 - t) < EPSILON)
+ t = 1.0;
-4. If r × s ≠ 0 and 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, the two line segments meet at the point p + t r = q + u s.
+ if (fabs(u) < EPSILON)
+ u = 0;
+ else if (fabs(1.0 - u) < EPSILON)
+ u = 1.0;
-5. Otherwise, the two line segments are not parallel but do not intersect.
-*/
- // Return parameter values, if user passed in valid pointers
- if (tp)
- *tp = t;
-
- if (up)
- *up = u;
+ 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))
- return 1;
-
- return 0;
+ Global::numIntersectParams = 2;
}
-
-// 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 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*/)
+void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2)
{
-#if 0
- Vector center = c->position;
- Vector v1 = l->position - center;
- Vector v2 = l->endpoint - center;
- Vector d = v2 - v1;
- double dr = d.Magnitude();
- double determinant = (v1.x * v2.y) - (v1.y * v2.x);
+ // Set up global vars
+ Global::numIntersectPoints = Global::numIntersectParams = 0;
- double discriminant = ((c->radius * c->radius) * (dr * dr)) - (determinant * determinant);
-
- if (discriminant < 0)
- return false;
-
-
-
- return true;
-#else
-/*
-I'm thinking a better approach to this might be as follows:
-
--- Get the distance of the circle's center from the line segment. If it's
- > the radius, it doesn't intersect.
--- If the parameter is off the line segment, check distance to endpoints. (Not sure
- how to proceed from here, it's different than the following.)
- [Actually, you can use the following for all of it. You only know if you have
- an intersection at the last step, which is OK.]
--- If the radius == distance, we have a tangent line.
--- If radius > distance, use Pythagorus to find the length on either side of the
- normal to the spots where the hypotenuse (== radius' length) contact the line.
--- Use those points to find the parameter on the line segment; if they're not on
- the line segment, no intersection.
-*/
- double t = ParameterOfLineAndPoint(l->position, l->endpoint, c->position);
-//small problem here: it clamps the result to the line segment. NOT what we want
-//here! !!! FIX !!! [DONE]
- Vector p = l->GetPointAtParameter(t);
- double distance = Vector::Magnitude(c->position, p);
-
- // If the center of the circle is farther from the line than the radius, fail.
- if (distance > c->radius)
- return 0;
-
- // 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();
+ // 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();
- return 1;
- }
+ // 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;
- // The line intersects the circle in two points (possibly). Use Pythagorus
- // to find them for testing.
- double offset = sqrt((c->radius * c->radius) - (distance * distance));
-//need to convert distance to paramter value... :-/
-//t = position on line / length of line segment, so if we divide the offset by length,
-//that should give us what we want.
- double length = Vector::Magnitude(l->position, l->endpoint);
- double t1 = t + (offset / length);
- double t2 = t - (offset / length);
-
-//need to find angles for the circle...
- Vector cp1 = l->position + (Vector(l->position, l->endpoint) * (length * t1));
- Vector cp2 = l->position + (Vector(l->position, l->endpoint) * (length * t2));
- double a1 = Vector(c->position, cp1).Angle();
- double a2 = Vector(c->position, cp2).Angle();
-
-//instead of this, return a # which is the # of intersections. [DONE]
- int intersections = 0;
-
- // Now check for if the parameters are in range
- if ((t1 >= 0) && (t1 <= 1.0))
+ // 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]))
{
- intersections++;
+ 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;
}
-
- if ((t2 >= 0) && (t2 <= 1.0))
+ else if (d == fabs(c1->radius[0] - c2->radius[0]))
{
- intersections++;
+ 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;
}
- return intersections;
-#endif
+/*
+ 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;
}
-
-// 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*/)
+//
+// N.B.: l is the line, c is the circle
+//
+void Geometry::CheckLineToCircleIntersection(Object * l, Object * c)
{
- // 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);
+ // Set up global vars
+ Global::numIntersectPoints = Global::numIntersectParams = 0;
- // 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;
+ // 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();
- // 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();
+ // Step 2: See if we have 0, 1, or 2 intersection points
- if (up)
- *up = centerLine.Angle() + PI;
+ // 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 1;
+ return;
}
- // Find the distance from the center of c1 to the perpendicular chord
- // (which contains the points of intersection)
- 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;
+ // Case #3: Two intersection points (possibly--secant)
- if (vp)
- *vp = centerLine.Angle() + angleC1;
+ // 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));
- if (wp)
- *wp = (centerLine.Angle() + PI) + angleC2;
+ // 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);
- if (p1)
- *p1 = c1->position + (centerLine.Unit() * x) + (Vector::Normal(Vector(), centerLine) * (a / 2.0));
+ // 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 (p2)
- *p2 = c1->position + (centerLine.Unit() * x) - (Vector::Normal(Vector(), centerLine) * (a / 2.0));
+ if ((u >= 0.0) && (u <= 1.0))
+ {
+ Global::intersectPoint[Global::numIntersectPoints] = i1;
+ Global::numIntersectPoints++;
+ }
- return 2;
+ 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^2 = a^2 + b^2 -2ab*cos(C)
+// c² = a² + b² - 2ab * cos(C)
// Solving for C:
-// cos(C) = (c^2 - a^2 - b^2) / -2ab = (a^2 + b^2 - c^2) / 2ab
+// cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab
// Law of Sines:
// a / sin A = b / sin B = c / sin C
// 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));
*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 of the point and the center of the circle
+ • Get the length of the line segment from point 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
+
+Another way:
+
+ • Find angle between radius and line between point and center
+ • Angle +/- from line (point to center) are the tangent points
+
+*/
+void Geometry::FindTangents(Object * c, Point p)
+{
+ // Set up global vars
+ Global::numIntersectPoints = Global::numIntersectParams = 0;
+
+ // Get the distance between the point and the center of the circle, the
+ // length, and the angle.
+ Vector centerLine(c->p[0], p);
+ double d = centerLine.Magnitude();
+ double clAngle = centerLine.Angle();
+
+ // If the point is on or inside the circle, there are no tangents.
+ if (d <= c->radius[0])
+ return;
+
+ // Use 'cos(a) = adjacent / hypotenuse' to get the tangent angle
+ double a = acos(c->radius[0] / d);
+
+ // Finally, find the points of intersection by using +/- the angle found
+ // from the centerline's angle
+ Global::intersectPoint[0].x = c->p[0].x + (cos(clAngle + a) * c->radius[0]);
+ Global::intersectPoint[0].y = c->p[0].y + (sin(clAngle + a) * c->radius[0]);
+ Global::intersectPoint[1].x = c->p[0].x + (cos(clAngle - a) * c->radius[0]);
+ Global::intersectPoint[1].y = c->p[0].y + (sin(clAngle - a) * c->radius[0]);
+ Global::numIntersectPoints = 2;
+}
+
+/*
+The extangents can be found by collapsing the circle of smaller radius to a point, and diminishing the larger circle by the radius of the smaller, then treating it like the point-circle case (you have to translate the tangent back out by the length of the radius of the smaller circle once found).
+
+Not sure what the analogous method for finding the intangents are... :-/
+Looks like the intangent can be found by augmenting the larger circle by the smaller radius, and doing the center of the smaller as the point-circle case again, only translating the line back by the radius of the smaller.
+*/
+void Geometry::FindTangents(Object * c1, Object * c2)
+{
+ // Set up global vars
+ Global::numIntersectPoints = Global::numIntersectParams = 0;
+
+ // Find the larger and smaller of the two:
+ Object * cLg = c1, * cSm = c2;
+
+ if (c2->radius[0] > c1->radius[0])
+ cLg = c2, cSm = c1;
+
+ // Get the distance between the point and the center of the circle, the
+ // length, and the angle.
+ Vector centerLine(cLg->p[0], cSm->p[0]);
+ double d = centerLine.Magnitude();
+ double clAngle = centerLine.Angle();
+
+ // If one circle is completely inside the other, there are no tangents.
+ if ((d + cSm->radius[0]) <= cLg->radius[0])
+ return;
+
+ // Subtract the radius of the smaller from the larger, and use the point-circle method to find the extangent:
+ double a = acos((cLg->radius[0] - cSm->radius[0]) / d);
+
+ // Finally, find the points of intersection by using +/- the angle found
+ // from the centerline's angle
+ Global::intersectPoint[0].x = cLg->p[0].x + (cos(clAngle + a) * cLg->radius[0]);
+ Global::intersectPoint[0].y = cLg->p[0].y + (sin(clAngle + a) * cLg->radius[0]);
+ Global::intersectPoint[1].x = cSm->p[0].x + (cos(clAngle + a) * cSm->radius[0]);
+ Global::intersectPoint[1].y = cSm->p[0].y + (sin(clAngle + a) * cSm->radius[0]);
+
+ Global::intersectPoint[2].x = cLg->p[0].x + (cos(clAngle - a) * cLg->radius[0]);
+ Global::intersectPoint[2].y = cLg->p[0].y + (sin(clAngle - a) * cLg->radius[0]);
+ Global::intersectPoint[3].x = cSm->p[0].x + (cos(clAngle - a) * cSm->radius[0]);
+ Global::intersectPoint[3].y = cSm->p[0].y + (sin(clAngle - a) * cSm->radius[0]);
+ Global::numIntersectPoints = 4;
+
+ // If the circles overlap, there are no intangents.
+ if (d <= (cLg->radius[0] + cSm->radius[0]))
+ return;
+
+ // Add the radius of the smaller from the larger, and use the point-circle method to find the intangent:
+ a = acos((cLg->radius[0] + cSm->radius[0]) / d);
+
+ // Finally, find the points of intersection by using +/- the angle found
+ // from the centerline's angle
+ Global::intersectPoint[4].x = cLg->p[0].x + (cos(clAngle + a) * cLg->radius[0]);
+ Global::intersectPoint[4].y = cLg->p[0].y + (sin(clAngle + a) * cLg->radius[0]);
+ Global::intersectPoint[5].x = cSm->p[0].x + (cos(clAngle - a) * cSm->radius[0]);
+ Global::intersectPoint[5].y = cSm->p[0].y + (sin(clAngle - a) * cSm->radius[0]);
+
+ Global::intersectPoint[6].x = cLg->p[0].x + (cos(clAngle - a) * cLg->radius[0]);
+ Global::intersectPoint[6].y = cLg->p[0].y + (sin(clAngle - a) * cLg->radius[0]);
+ Global::intersectPoint[7].x = cSm->p[0].x + (cos(clAngle + a) * cSm->radius[0]);
+ Global::intersectPoint[7].y = cSm->p[0].y + (sin(clAngle + a) * cSm->radius[0]);
+ Global::numIntersectPoints = 8;
+}
+
+//
+// Parameter 1: point in question
+// Parameter 2, 3: points we are comparing to
+//
+Point Geometry::NearestTo(Point point, Point p1, Point p2)
+{
+ double l1 = Vector::Magnitude(point, p1);
+ double l2 = Vector::Magnitude(point, p2);
+
+ return (l1 < l2 ? p1 : p2);
+}
+
+Vector Geometry::GetNormalOfPointAndLine(Point p, Line * l)
+{
+ Vector normal = Vector::Normal(l->p[0], l->p[1]);
+ normal.Unit();
+ double dir = Vector::Dot(p - l->p[0], normal);
+
+ if (dir < 0)
+ normal = -normal;
+
+ return normal;
+}
+
+Circle Geometry::FindCircleForThreePoints(Point p1, Point p2, Point p3)
+{
+ // We use matrices and determinants to find the center and radius of the
+ // circle, given the three points passed in:
+ //
+ // [(x^2 + y^2) x y 1 ]
+ // [(x1^2 + y1^2) x1 y1 1 ]
+ // [(x2^2 + y2^2) x2 y2 1 ]
+ // [(x3^2 + y3^2) x3 y3 1 ]
+ //
+ // Then, center <x0, y0> is found by:
+ //
+ // x0 = 1/2 • M12/M11, y0 = -1/2 • M13/M11
+ //
+ // The radius can be found with the center and any of the points via
+ // Pythagoras, or with:
+ //
+ // r^2 = x0^2 + y0^2 + M14/M11
+ //
+ // If M11 = 0, then the three points are colinear and there is no circle.
+ // (M## is the minor determinant of the 4x4 matrix, where a 3x3 matrix is
+ // created by deleting the row and column of the 4x4 with the indices
+ // given.)
+
+ Circle c;
+ double m[3][4];
+
+ m[0][0] = (p1.x * p1.x) + (p1.y * p1.y);
+ m[0][1] = p1.x;
+ m[0][2] = p1.y;
+ m[0][3] = 1.0;
+ m[1][0] = (p2.x * p2.x) + (p2.y * p2.y);
+ m[1][1] = p2.x;
+ m[1][2] = p2.y;
+ m[1][3] = 1.0;
+ m[2][0] = (p3.x * p3.x) + (p3.y * p3.y);
+ m[2][1] = p3.x;
+ m[2][2] = p3.y;
+ m[2][3] = 1.0;
+
+ double minor11 = Determinant3x3(m[0][1], m[0][2], m[0][3], m[1][1], m[1][2], m[1][3], m[2][1], m[2][2], m[2][3]);
+
+ // Sanity check: if M11 is zero, the points are colinear.
+ if (minor11 == 0)
+ return c;
+
+ double minor12 = Determinant3x3(m[0][0], m[0][2], m[0][3], m[1][0], m[1][2], m[1][3], m[2][0], m[2][2], m[2][3]);
+ double minor13 = Determinant3x3(m[0][0], m[0][1], m[0][3], m[1][0], m[1][1], m[1][3], m[2][0], m[2][1], m[2][3]);
+ double minor14 = Determinant3x3(m[0][0], m[0][1], m[0][2], m[1][0], m[1][1], m[1][2], m[2][0], m[2][1], m[2][2]);
+
+ c.p[0].x = 0.5 * (minor12 / minor11);
+ c.p[0].y = -0.5 * (minor13 / minor11);
+ c.radius[0] = sqrt((c.p[0].x * c.p[0].x) + (c.p[0].y * c.p[0].y) + (minor14 / minor11));
+
+ return c;
+}
+
+double Geometry::Determinant3x3(double a11, double a12, double a13, double a21, double a22, double a23, double a31, double a32, double a33)
+{
+ return (a11 * ((a22 * a33) - (a32 * a23)))
+ - (a12 * ((a21 * a33) - (a31 * a23)))
+ + (a13 * ((a21 * a32) - (a31 * a22)));
+}
+
+Arc Geometry::Unpack(Point tail, Point head, double bump)
+{
+ double length = Vector::Magnitude(tail, head) / 2.0;
+ double bumpLen = length * fabs(bump);
+ Point midpoint = Vector::Midpoint(tail, head);
+ Vector mpNormal = Vector::Normal(tail, head); // Normal points to the left
+
+ // Flip the normal if the bump is pointing left
+ if (bump < 0)
+ mpNormal = -mpNormal;
+
+ // N.B.: The radius can also be found with r = (a² + h²) / 2h where a is
+ // the 1/2 length of the line segment and h is the bump length.
+// double radius = (length + (1.0 / length)) / 2.0;
+ double radius = 0.5 * (((length * length) + (bumpLen * bumpLen)) / bumpLen);
+ Vector ctrVec = mpNormal * (radius - bumpLen);
+ Point center = midpoint + ctrVec;
+
+ // Find arc endpoints...
+ double angle1 = Vector::Angle(center, tail);
+ double angle2 = Vector::Angle(center, head);
+ double span = angle2 - angle1;
+
+ // Fix span depending on which way the arc is being drawn...
+ if (bump > 0 && span < 0)
+ span += TAU;
+
+ if (bump < 0 && span > 0)
+ span -= TAU;
+
+ return Arc(center, radius, angle1, span);
+}
projects / architektonas / blobdiff