2 // geometry.cpp: Algebraic geometry helper functions
4 // Part of the Architektonas Project
5 // (C) 2011 Underground Software
6 // See the README and GPLv3 files for licensing and warranty information
8 // JLH = James Hammons <jlhamm@acm.org>
11 // --- ---------- ------------------------------------------------------------
12 // JLH 08/31/2011 Created this file
14 // NOTE: All methods in this class are static.
21 #include "mathconstants.h"
23 // Returns the parameter of a point in space to this vector. If the parameter
24 // is between 0 and 1, the normal of the vector to the point is on the vector.
25 // Note: lp1 is the tail, lp2 is the head of the line (vector).
26 double Geometry::ParameterOfLineAndPoint(Point tail, Point head, Point point)
28 // Geometric interpretation:
29 // The parameterized point on the vector lineSegment is where the normal of
30 // the lineSegment to the point intersects lineSegment. If the pp < 0, then
31 // the perpendicular lies beyond the 1st endpoint. If pp > 1, then the
32 // perpendicular lies beyond the 2nd endpoint.
34 Vector lineSegment = head - tail;
35 double magnitude = lineSegment.Magnitude();
36 Vector pointSegment = point - tail;
37 double t = lineSegment.Dot(pointSegment) / (magnitude * magnitude);
42 double Geometry::DistanceToLineFromPoint(Point tail, Point head, Point point)
44 // Interpretation: given a line in the form x = a + tu, where u is the
45 // unit vector of the line, a is the tail and t is a parameter which
46 // describes the line, the distance of a point p to the line is given by:
47 // || (a - p) - ((a - p) • u) u ||
48 // We go an extra step: we set the sign to reflect which side of the line
49 // it's on (+ == to the left if head points away from you, - == to the
51 Vector line(tail, head);
52 Vector u = line.Unit();
53 Vector a_p = tail - point;
54 Vector dist = a_p - (u * (a_p).Dot(u));
56 double angle = Vector::Angle(tail, point) - line.Angle();
61 return dist.Magnitude() * (angle < HALF_TAU ? +1.0 : -1.0);
64 Point Geometry::MirrorPointAroundLine(Point point, Point tail, Point head)
66 // Get the vector of the intersection of the line and the normal on the
67 // line to the point in question.
68 double t = ParameterOfLineAndPoint(tail, head, point);
69 Vector v = Vector(tail, head) * t;
71 // Get the point normal to point to the line passed in
72 Point normalOnLine = tail + v;
74 // Make our mirrored vector (head - tail)
75 Vector mirror = -(point - normalOnLine);
77 // Find the mirrored point
78 Point mirroredPoint = normalOnLine + mirror;
84 // point: The point we're rotating
85 // rotationPoint: The point we're rotating around
87 Point Geometry::RotatePointAroundPoint(Point point, Point rotationPoint, double angle)
89 Vector v = Vector(rotationPoint, point);
90 double px = (v.x * cos(angle)) - (v.y * sin(angle));
91 double py = (v.x * sin(angle)) + (v.y * cos(angle));
93 return Vector(rotationPoint.x + px, rotationPoint.y + py, 0);
96 double Geometry::Determinant(Point p1, Point p2)
98 return (p1.x * p2.y) - (p2.x * p1.y);
101 void Geometry::Intersects(Object * obj1, Object * obj2)
103 Global::numIntersectPoints = Global::numIntersectParams = 0;
105 if ((obj1->type == OTLine) && (obj2->type == OTLine))
106 CheckLineToLineIntersection(obj1, obj2);
107 else if ((obj1->type == OTCircle) && (obj2->type == OTCircle))
108 CheckCircleToCircleIntersection(obj1, obj2);
109 else if ((obj1->type == OTLine) && (obj2->type == OTCircle))
110 CheckLineToCircleIntersection(obj1, obj2);
111 else if ((obj1->type == OTCircle) && (obj2->type == OTLine))
112 CheckLineToCircleIntersection(obj2, obj1);
116 Intersecting line segments:
118 Segment L1 has edges A=(a1,a2), A'=(a1',a2').
119 Segment L2 has edges B=(b1,b2), B'=(b1',b2').
120 Segment L1 is the set of points tA'+(1-t)A, where 0<=t<=1.
121 Segment L2 is the set of points sB'+(1-s)B, where 0<=s<=1.
122 Segment L1 meet segment L2 if and only if for some t and s we have
123 tA'+(1-t)A=sB'+(1-s)B
124 The solution of this with respect to t and s is
126 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?))
128 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?))
130 So check if the above two numbers are both >=0 and <=1.
134 // Finds the intersection between two lines (if any)
135 void Geometry::CheckLineToLineIntersection(Object * l1, Object * l2)
137 Global::numIntersectPoints = Global::numIntersectParams = 0;
139 Vector r(l1->p[0], l1->p[1]);
140 Vector s(l2->p[0], l2->p[1]);
141 Vector v1 = l2->p[0] - l1->p[0]; // q - p
143 double rxs = (r.x * s.y) - (s.x * r.y);
146 // 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
149 // WHY would you do this??? The lines overlap at an INFINITE number of points!
150 // The assumption is that the angle is 180°, not 0°.
151 // 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).
152 /* double qpxr = (v1.x * r.y) - (r.x * v1.y);
154 // Line segments are parallel, so no intersection...
158 // Otherwise, the segments are colinear. Need to check for the 0° (degenerate) vs the 180° (OK) case.
159 Object * larger = l1;
160 Object * smaller = l2;
162 if (r->Magnitude() < s->Magnitude())
168 double param1 = ParameterOfLineAndPoint(larger->p[0], larger->p[1], smaller->p[0]);
169 double param2 = ParameterOfLineAndPoint(larger->p[0], larger->p[1], smaller->p[1]);
171 // Check for the degenerate case, and return if found
172 if ((param1 > 0 && param1 < 1.0) || (param2 > 0 && param2 < 1.0))
175 //// or just use AngleBetween: Nah, won't work...
177 // Check to see which endpoints are connected... Four possibilities:
178 if (l1->p[0] == l2->p[0])
180 else if (l1->p[0] == l2->p[1])
182 else if (l1->p[1] == l2->p[0])
184 else if (l1->p[1] == l2->p[1])
191 t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
192 u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
195 // Check that the parameters are above the epsilon, otherwise clamp them to
196 // zero or one, as the case may be
197 if (fabs(t) < EPSILON)
199 else if (fabs(1.0 - t) < EPSILON)
202 if (fabs(u) < EPSILON)
204 else if (fabs(1.0 - u) < EPSILON)
207 Global::intersectParam[0] = t;
208 Global::intersectParam[1] = u;
210 // If the parameters are in range, we have overlap!
211 if ((t >= 0) && (t <= 1.0) && (u >= 0) && (u <= 1.0))
212 Global::numIntersectParams = 2;
215 void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2)
217 // Set up global vars
218 Global::numIntersectPoints = Global::numIntersectParams = 0;
220 // Get the distance between the centers of the circles
221 Vector centerLine(c1->p[0], c2->p[0]);
222 double d = centerLine.Magnitude();
223 double clAngle = centerLine.Angle();
225 // If the distance between centers is greater than the sum of the radii or
226 // less than the difference between the radii, there is NO intersection
227 if ((d > (c1->radius[0] + c2->radius[0]))
228 || (d < fabs(c1->radius[0] - c2->radius[0])))
231 // If the distance between centers is equal to the sum of the radii or
232 // equal to the difference between the radii, the intersection is tangent
234 if (d == (c1->radius[0] + c2->radius[0]))
236 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0]);
237 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0]);
238 Global::numIntersectPoints = 1;
241 else if (d == fabs(c1->radius[0] - c2->radius[0]))
243 double sign = (c1->radius[0] > c2->radius[0] ? +1 : -1);
244 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0] * sign);
245 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0] * sign);
246 Global::numIntersectPoints = 1;
251 c² = a² + b² - 2ab·cos µ
252 2ab·cos µ = a² + b² - c²
253 cos µ = (a² + b² - c²) / 2ab
255 // Use the Law of Cosines to find the angle between the centerline and the
256 // radial line on Circle #1
257 double a = acos(((c1->radius[0] * c1->radius[0]) + (d * d) - (c2->radius[0] * c2->radius[0])) / (2.0 * c1->radius[0] * d));
259 // Finally, find the points of intersection by using +/- the angle found
260 // from the centerline's angle
261 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle + a) * c1->radius[0]);
262 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle + a) * c1->radius[0]);
263 Global::intersectPoint[1].x = c1->p[0].x + (cos(clAngle - a) * c1->radius[0]);
264 Global::intersectPoint[1].y = c1->p[0].y + (sin(clAngle - a) * c1->radius[0]);
265 Global::numIntersectPoints = 2;
269 // N.B.: l is the line, c is the circle
271 void Geometry::CheckLineToCircleIntersection(Object * l, Object * c)
273 // Set up global vars
274 Global::numIntersectPoints = Global::numIntersectParams = 0;
276 // Step 1: Find shortest distance from center of circle to the infinite line
277 double t = ParameterOfLineAndPoint(l->p[0], l->p[1], c->p[0]);
278 Point p = l->p[0] + (Vector(l->p[0], l->p[1]) * t);
279 Vector radial = Vector(c->p[0], p);
280 double distance = radial.Magnitude();
282 // Step 2: See if we have 0, 1, or 2 intersection points
284 // Case #1: No intersection points
285 if (distance > c->radius[0])
287 // Case #2: One intersection point (possibly--tangent)
288 else if (distance == c->radius[0])
290 // Only intersects if the parameter is on the line segment!
291 if ((t >= 0.0) && (t <= 1.0))
293 Global::intersectPoint[0] = c->p[0] + radial;
294 Global::numIntersectPoints = 1;
300 // Case #3: Two intersection points (possibly--secant)
302 // So, we have the line, and the perpendicular from the center of the
303 // circle to the line. Now figure out where the intersection points are.
304 // This is a right triangle, though do we really know all the sides?
305 // Don't need to, 2 is enough for Pythagoras :-)
306 // Radius is the hypotenuse, so we have to use c² = a² + b² => a² = c² - b²
307 double perpendicularLength = sqrt((c->radius[0] * c->radius[0]) - (distance * distance));
309 // Now, find the intersection points using the length...
310 Vector lineUnit = Vector(l->p[0], l->p[1]).Unit();
311 Point i1 = p + (lineUnit * perpendicularLength);
312 Point i2 = p - (lineUnit * perpendicularLength);
314 // Next we need to see if they are on the line segment...
315 double u = ParameterOfLineAndPoint(l->p[0], l->p[1], i1);
316 double v = ParameterOfLineAndPoint(l->p[0], l->p[1], i2);
318 if ((u >= 0.0) && (u <= 1.0))
320 Global::intersectPoint[Global::numIntersectPoints] = i1;
321 Global::numIntersectPoints++;
324 if ((v >= 0.0) && (v <= 1.0))
326 Global::intersectPoint[Global::numIntersectPoints] = i2;
327 Global::numIntersectPoints++;
331 // should we just do common trig solves, like AAS, ASA, SAS, SSA?
333 // c² = a² + b² - 2ab * cos(C)
335 // cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab
337 // a / sin A = b / sin B = c / sin C
339 // Solve the angles of the triangle given the sides. Angles returned are
340 // opposite of the given sides (so a1 consists of sides s2 & s3, and so on).
341 void Geometry::FindAnglesForSides(double s1, double s2, double s3, double * a1, double * a2, double * a3)
343 // Use law of cosines to find 1st angle
344 double cosine1 = ((s2 * s2) + (s3 * s3) - (s1 * s1)) / (2.0 * s2 * s3);
346 // Check for a valid triangle
347 if ((cosine1 < -1.0) || (cosine1 > 1.0))
350 double angle1 = acos(cosine1);
352 // Use law of sines to find 2nd & 3rd angles
353 // sin A / a = sin B / b
354 // sin B = (sin A / a) * b
355 // B = arcsin( sin A * (b / a))
356 // ??? ==> B = A * arcsin(b / a)
359 sin B = sin A * (b / a)
360 sin B / sin A = b / a
361 arcsin( sin B / sin A ) = arcsin( b / a )
366 double angle2 = asin(s2 * (sin(angle1) / s1));
367 double angle3 = asin(s3 * (sin(angle1) / s1));
379 Point Geometry::GetPointForParameter(Object * obj, double t)
381 if (obj->type == OTLine)
383 // Translate line vector to the origin, then add the scaled vector to
384 // initial point of the line.
385 Vector v = obj->p[1] - obj->p[0];
386 return obj->p[0] + (v * t);
392 Point Geometry::Midpoint(Line * line)
394 return Point((line->p[0].x + line->p[1].x) / 2.0,
395 (line->p[0].y + line->p[1].y) / 2.0);
399 How to find the tangent of a point off a circle:
401 • Calculate the midpoint of the point and the center of the circle
402 • Get the length of the line segment from point and the center divided by two
403 • Use that length to construct a circle with the point at the center and the
404 radius equal to that length
405 • The intersection of the two circles are the tangent points
409 • Find angle between radius and line between point and center
410 • Angle +/- from line (point to center) are the tangent points
413 void Geometry::FindTangents(Object * c, Point p)
415 // Set up global vars
416 Global::numIntersectPoints = Global::numIntersectParams = 0;
418 // Get the distance between the point and the center of the circle, the
419 // length, and the angle.
420 Vector centerLine(c->p[0], p);
421 double d = centerLine.Magnitude();
422 double clAngle = centerLine.Angle();
424 // If the point is on or inside the circle, there are no tangents.
425 if (d <= c->radius[0])
428 // Use 'cos(a) = adjacent / hypotenuse' to get the tangent angle
429 double a = acos(c->radius[0] / d);
431 // Finally, find the points of intersection by using +/- the angle found
432 // from the centerline's angle
433 Global::intersectPoint[0].x = c->p[0].x + (cos(clAngle + a) * c->radius[0]);
434 Global::intersectPoint[0].y = c->p[0].y + (sin(clAngle + a) * c->radius[0]);
435 Global::intersectPoint[1].x = c->p[0].x + (cos(clAngle - a) * c->radius[0]);
436 Global::intersectPoint[1].y = c->p[0].y + (sin(clAngle - a) * c->radius[0]);
437 Global::numIntersectPoints = 2;
441 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).
443 Not sure what the analogous method for finding the intangents are... :-/
444 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.
446 void Geometry::FindTangents(Object * c1, Object * c2)
448 // Set up global vars
449 Global::numIntersectPoints = Global::numIntersectParams = 0;
451 // Find the larger and smaller of the two:
452 Object * cLg = c1, * cSm = c2;
454 if (c2->radius[0] > c1->radius[0])
457 // Get the distance between the point and the center of the circle, the
458 // length, and the angle.
459 Vector centerLine(cLg->p[0], cSm->p[0]);
460 double d = centerLine.Magnitude();
461 double clAngle = centerLine.Angle();
463 // If one circle is completely inside the other, there are no tangents.
464 if ((d + cSm->radius[0]) <= cLg->radius[0])
467 // Subtract the radius of the smaller from the larger, and use the point-circle method to find the extangent:
468 double a = acos((cLg->radius[0] - cSm->radius[0]) / d);
470 // Finally, find the points of intersection by using +/- the angle found
471 // from the centerline's angle
472 Global::intersectPoint[0].x = cLg->p[0].x + (cos(clAngle + a) * cLg->radius[0]);
473 Global::intersectPoint[0].y = cLg->p[0].y + (sin(clAngle + a) * cLg->radius[0]);
474 Global::intersectPoint[1].x = cSm->p[0].x + (cos(clAngle + a) * cSm->radius[0]);
475 Global::intersectPoint[1].y = cSm->p[0].y + (sin(clAngle + a) * cSm->radius[0]);
477 Global::intersectPoint[2].x = cLg->p[0].x + (cos(clAngle - a) * cLg->radius[0]);
478 Global::intersectPoint[2].y = cLg->p[0].y + (sin(clAngle - a) * cLg->radius[0]);
479 Global::intersectPoint[3].x = cSm->p[0].x + (cos(clAngle - a) * cSm->radius[0]);
480 Global::intersectPoint[3].y = cSm->p[0].y + (sin(clAngle - a) * cSm->radius[0]);
481 Global::numIntersectPoints = 4;
483 // If the circles overlap, there are no intangents.
484 if (d <= (cLg->radius[0] + cSm->radius[0]))
487 // Add the radius of the smaller from the larger, and use the point-circle method to find the intangent:
488 a = acos((cLg->radius[0] + cSm->radius[0]) / d);
490 // Finally, find the points of intersection by using +/- the angle found
491 // from the centerline's angle
492 Global::intersectPoint[4].x = cLg->p[0].x + (cos(clAngle + a) * cLg->radius[0]);
493 Global::intersectPoint[4].y = cLg->p[0].y + (sin(clAngle + a) * cLg->radius[0]);
494 Global::intersectPoint[5].x = cSm->p[0].x + (cos(clAngle - a) * cSm->radius[0]);
495 Global::intersectPoint[5].y = cSm->p[0].y + (sin(clAngle - a) * cSm->radius[0]);
497 Global::intersectPoint[6].x = cLg->p[0].x + (cos(clAngle - a) * cLg->radius[0]);
498 Global::intersectPoint[6].y = cLg->p[0].y + (sin(clAngle - a) * cLg->radius[0]);
499 Global::intersectPoint[7].x = cSm->p[0].x + (cos(clAngle + a) * cSm->radius[0]);
500 Global::intersectPoint[7].y = cSm->p[0].y + (sin(clAngle + a) * cSm->radius[0]);
501 Global::numIntersectPoints = 8;
505 // Parameter 1: point in question
506 // Parameter 2, 3: points we are comparing to
508 Point Geometry::NearestTo(Point point, Point p1, Point p2)
510 double l1 = Vector::Magnitude(point, p1);
511 double l2 = Vector::Magnitude(point, p2);
513 return (l1 < l2 ? p1 : p2);
516 Vector Geometry::GetNormalOfPointAndLine(Point p, Line * l)
518 Vector normal = Vector::Normal(l->p[0], l->p[1]);
520 double dir = Vector::Dot(p - l->p[0], normal);
528 Circle Geometry::FindCircleForThreePoints(Point p1, Point p2, Point p3)
530 // We use matrices and determinants to find the center and radius of the
531 // circle, given the three points passed in:
533 // [(x^2 + y^2) x y 1 ]
534 // [(x1^2 + y1^2) x1 y1 1 ]
535 // [(x2^2 + y2^2) x2 y2 1 ]
536 // [(x3^2 + y3^2) x3 y3 1 ]
538 // Then, center <x0, y0> is found by:
540 // x0 = 1/2 • M12/M11, y0 = -1/2 • M13/M11
542 // The radius can be found with the center and any of the points via
543 // Pythagoras, or with:
545 // r^2 = x0^2 + y0^2 + M14/M11
547 // If M11 = 0, then the three points are colinear and there is no circle.
548 // (M## is the minor determinant of the 4x4 matrix, where a 3x3 matrix is
549 // created by deleting the row and column of the 4x4 with the indices
555 m[0][0] = (p1.x * p1.x) + (p1.y * p1.y);
559 m[1][0] = (p2.x * p2.x) + (p2.y * p2.y);
563 m[2][0] = (p3.x * p3.x) + (p3.y * p3.y);
568 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]);
570 // Sanity check: if M11 is zero, the points are colinear.
574 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]);
575 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]);
576 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]);
578 c.p[0].x = 0.5 * (minor12 / minor11);
579 c.p[0].y = -0.5 * (minor13 / minor11);
580 c.radius[0] = sqrt((c.p[0].x * c.p[0].x) + (c.p[0].y * c.p[0].y) + (minor14 / minor11));
585 double Geometry::Determinant3x3(double a11, double a12, double a13, double a21, double a22, double a23, double a31, double a32, double a33)
587 return (a11 * ((a22 * a33) - (a32 * a23)))
588 - (a12 * ((a21 * a33) - (a31 * a23)))
589 + (a13 * ((a21 * a32) - (a31 * a22)));
592 Arc Geometry::Unpack(Point tail, Point head, double bump)
594 double length = Vector::Magnitude(tail, head) / 2.0;
595 double bumpLen = length * fabs(bump);
596 Point midpoint = Vector::Midpoint(tail, head);
597 Vector mpNormal = Vector::Normal(tail, head); // Normal points to the left
599 // Flip the normal if the bump is pointing left
601 mpNormal = -mpNormal;
603 // N.B.: The radius can also be found with r = (a² + h²) / 2h where a is
604 // the 1/2 length of the line segment and h is the bump length.
605 // double radius = (length + (1.0 / length)) / 2.0;
606 double radius = 0.5 * (((length * length) + (bumpLen * bumpLen)) / bumpLen);
607 Vector ctrVec = mpNormal * (radius - bumpLen);
608 Point center = midpoint + ctrVec;
610 // Find arc endpoints...
611 double angle1 = Vector::Angle(center, tail);
612 double angle2 = Vector::Angle(center, head);
613 double span = angle2 - angle1;
615 // Fix span depending on which way the arc is being drawn...
616 if (bump > 0 && span < 0)
619 if (bump < 0 && span > 0)
622 return Arc(center, radius, angle1, span);