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);
148 double qpxr = (v1.x * r.y) - (r.x * v1.y);
150 // Lines are parallel, so no intersection...
154 // Check to see which endpoints are connected... Four possibilities:
155 if (l1->p[0] == l2->p[0])
157 else if (l1->p[0] == l2->p[1])
159 else if (l1->p[1] == l2->p[0])
161 else if (l1->p[1] == l2->p[1])
168 t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
169 u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
172 Global::intersectParam[0] = t;
173 Global::intersectParam[1] = u;
175 // If the parameters are in range, we have overlap!
176 if ((t >= 0) && (t <= 1.0) && (u >= 0) && (u <= 1.0))
177 Global::numIntersectParams = 1;
180 void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2)
182 // Set up global vars
183 Global::numIntersectPoints = Global::numIntersectParams = 0;
185 // Get the distance between the centers of the circles
186 Vector centerLine(c1->p[0], c2->p[0]);
187 double d = centerLine.Magnitude();
188 double clAngle = centerLine.Angle();
190 // If the distance between centers is greater than the sum of the radii or
191 // less than the difference between the radii, there is NO intersection
192 if ((d > (c1->radius[0] + c2->radius[0]))
193 || (d < fabs(c1->radius[0] - c2->radius[0])))
196 // If the distance between centers is equal to the sum of the radii or
197 // equal to the difference between the radii, the intersection is tangent
199 if (d == (c1->radius[0] + c2->radius[0]))
201 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0]);
202 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0]);
203 Global::numIntersectPoints = 1;
206 else if (d == fabs(c1->radius[0] - c2->radius[0]))
208 double sign = (c1->radius[0] > c2->radius[0] ? +1 : -1);
209 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0] * sign);
210 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0] * sign);
211 Global::numIntersectPoints = 1;
216 c² = a² + b² - 2ab·cos µ
217 2ab·cos µ = a² + b² - c²
218 cos µ = (a² + b² - c²) / 2ab
220 // Use the Law of Cosines to find the angle between the centerline and the
221 // radial line on Circle #1
222 double a = acos(((c1->radius[0] * c1->radius[0]) + (d * d) - (c2->radius[0] * c2->radius[0])) / (2.0 * c1->radius[0] * d));
224 // Finally, find the points of intersection by using +/- the angle found
225 // from the centerline's angle
226 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle + a) * c1->radius[0]);
227 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle + a) * c1->radius[0]);
228 Global::intersectPoint[1].x = c1->p[0].x + (cos(clAngle - a) * c1->radius[0]);
229 Global::intersectPoint[1].y = c1->p[0].y + (sin(clAngle - a) * c1->radius[0]);
230 Global::numIntersectPoints = 2;
234 // N.B.: l is the line, c is the circle
236 void Geometry::CheckLineToCircleIntersection(Object * l, Object * c)
238 // Set up global vars
239 Global::numIntersectPoints = Global::numIntersectParams = 0;
241 // Step 1: Find shortest distance from center of circle to the infinite line
242 double t = ParameterOfLineAndPoint(l->p[0], l->p[1], c->p[0]);
243 Point p = l->p[0] + (Vector(l->p[0], l->p[1]) * t);
244 Vector radial = Vector(c->p[0], p);
245 double distance = radial.Magnitude();
247 // Step 2: See if we have 0, 1, or 2 intersection points
249 // Case #1: No intersection points
250 if (distance > c->radius[0])
252 // Case #2: One intersection point (possibly--tangent)
253 else if (distance == c->radius[0])
255 // Only intersects if the parameter is on the line segment!
256 if ((t >= 0.0) && (t <= 1.0))
258 Global::intersectPoint[0] = c->p[0] + radial;
259 Global::numIntersectPoints = 1;
265 // Case #3: Two intersection points (possibly--secant)
267 // So, we have the line, and the perpendicular from the center of the
268 // circle to the line. Now figure out where the intersection points are.
269 // This is a right triangle, though do we really know all the sides?
270 // Don't need to, 2 is enough for Pythagoras :-)
271 // Radius is the hypotenuse, so we have to use c² = a² + b² => a² = c² - b²
272 double perpendicularLength = sqrt((c->radius[0] * c->radius[0]) - (distance * distance));
274 // Now, find the intersection points using the length...
275 Vector lineUnit = Vector(l->p[0], l->p[1]).Unit();
276 Point i1 = p + (lineUnit * perpendicularLength);
277 Point i2 = p - (lineUnit * perpendicularLength);
279 // Next we need to see if they are on the line segment...
280 double u = ParameterOfLineAndPoint(l->p[0], l->p[1], i1);
281 double v = ParameterOfLineAndPoint(l->p[0], l->p[1], i2);
283 if ((u >= 0.0) && (u <= 1.0))
285 Global::intersectPoint[Global::numIntersectPoints] = i1;
286 Global::numIntersectPoints++;
289 if ((v >= 0.0) && (v <= 1.0))
291 Global::intersectPoint[Global::numIntersectPoints] = i2;
292 Global::numIntersectPoints++;
296 // should we just do common trig solves, like AAS, ASA, SAS, SSA?
298 // c² = a² + b² - 2ab * cos(C)
300 // cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab
302 // a / sin A = b / sin B = c / sin C
304 // Solve the angles of the triangle given the sides. Angles returned are
305 // opposite of the given sides (so a1 consists of sides s2 & s3, and so on).
306 void Geometry::FindAnglesForSides(double s1, double s2, double s3, double * a1, double * a2, double * a3)
308 // Use law of cosines to find 1st angle
309 double cosine1 = ((s2 * s2) + (s3 * s3) - (s1 * s1)) / (2.0 * s2 * s3);
311 // Check for a valid triangle
312 if ((cosine1 < -1.0) || (cosine1 > 1.0))
315 double angle1 = acos(cosine1);
317 // Use law of sines to find 2nd & 3rd angles
318 // sin A / a = sin B / b
319 // sin B = (sin A / a) * b
320 // B = arcsin( sin A * (b / a))
321 // ??? ==> B = A * arcsin(b / a)
324 sin B = sin A * (b / a)
325 sin B / sin A = b / a
326 arcsin( sin B / sin A ) = arcsin( b / a )
331 double angle2 = asin(s2 * (sin(angle1) / s1));
332 double angle3 = asin(s3 * (sin(angle1) / s1));
344 Point Geometry::GetPointForParameter(Object * obj, double t)
346 if (obj->type == OTLine)
348 // Translate line vector to the origin, then add the scaled vector to
349 // initial point of the line.
350 Vector v = obj->p[1] - obj->p[0];
351 return obj->p[0] + (v * t);
357 Point Geometry::Midpoint(Line * line)
359 return Point((line->p[0].x + line->p[1].x) / 2.0,
360 (line->p[0].y + line->p[1].y) / 2.0);
364 How to find the tangent of a point off a circle:
366 • Calculate the midpoint of the point and the center of the circle
367 • Get the length of the line segment from point and the center divided by two
368 • Use that length to construct a circle with the point at the center and the
369 radius equal to that length
370 • The intersection of the two circles are the tangent points
374 • Find angle between radius and line between point and center
375 • Angle +/- from line (point to center) are the tangent points
378 void Geometry::FindTangents(Object * c, Point p)
380 // Set up global vars
381 Global::numIntersectPoints = Global::numIntersectParams = 0;
383 // Get the distance between the point and the center of the circle, the
384 // length, and the angle.
385 Vector centerLine(c->p[0], p);
386 double d = centerLine.Magnitude();
387 double clAngle = centerLine.Angle();
389 // If the point is on or inside the circle, there are no tangents.
390 if (d <= c->radius[0])
393 // Use 'cos(a) = adjacent / hypotenuse' to get the tangent angle
394 double a = acos(c->radius[0] / d);
396 // Finally, find the points of intersection by using +/- the angle found
397 // from the centerline's angle
398 Global::intersectPoint[0].x = c->p[0].x + (cos(clAngle + a) * c->radius[0]);
399 Global::intersectPoint[0].y = c->p[0].y + (sin(clAngle + a) * c->radius[0]);
400 Global::intersectPoint[1].x = c->p[0].x + (cos(clAngle - a) * c->radius[0]);
401 Global::intersectPoint[1].y = c->p[0].y + (sin(clAngle - a) * c->radius[0]);
402 Global::numIntersectPoints = 2;
406 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).
408 Not sure what the analogous method for finding the intangents are... :-/
409 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.
411 void Geometry::FindTangents(Object * c1, Object * c2)
413 // Set up global vars
414 Global::numIntersectPoints = Global::numIntersectParams = 0;
416 // Find the larger and smaller of the two:
417 Object * cLg = c1, * cSm = c2;
419 if (c2->radius[0] > c1->radius[0])
422 // Get the distance between the point and the center of the circle, the
423 // length, and the angle.
424 Vector centerLine(cLg->p[0], cSm->p[0]);
425 double d = centerLine.Magnitude();
426 double clAngle = centerLine.Angle();
428 // If one circle is completely inside the other, there are no tangents.
429 if ((d + cSm->radius[0]) <= cLg->radius[0])
432 // Subtract the radius of the smaller from the larger, and use the point-circle method to find the extangent:
433 double a = acos((cLg->radius[0] - cSm->radius[0]) / d);
435 // Finally, find the points of intersection by using +/- the angle found
436 // from the centerline's angle
437 Global::intersectPoint[0].x = cLg->p[0].x + (cos(clAngle + a) * cLg->radius[0]);
438 Global::intersectPoint[0].y = cLg->p[0].y + (sin(clAngle + a) * cLg->radius[0]);
439 Global::intersectPoint[1].x = cSm->p[0].x + (cos(clAngle + a) * cSm->radius[0]);
440 Global::intersectPoint[1].y = cSm->p[0].y + (sin(clAngle + a) * cSm->radius[0]);
442 Global::intersectPoint[2].x = cLg->p[0].x + (cos(clAngle - a) * cLg->radius[0]);
443 Global::intersectPoint[2].y = cLg->p[0].y + (sin(clAngle - a) * cLg->radius[0]);
444 Global::intersectPoint[3].x = cSm->p[0].x + (cos(clAngle - a) * cSm->radius[0]);
445 Global::intersectPoint[3].y = cSm->p[0].y + (sin(clAngle - a) * cSm->radius[0]);
446 Global::numIntersectPoints = 4;
448 // If the circles overlap, there are no intangents.
449 if (d <= (cLg->radius[0] + cSm->radius[0]))
452 // Add the radius of the smaller from the larger, and use the point-circle method to find the intangent:
453 a = acos((cLg->radius[0] + cSm->radius[0]) / d);
455 // Finally, find the points of intersection by using +/- the angle found
456 // from the centerline's angle
457 Global::intersectPoint[4].x = cLg->p[0].x + (cos(clAngle + a) * cLg->radius[0]);
458 Global::intersectPoint[4].y = cLg->p[0].y + (sin(clAngle + a) * cLg->radius[0]);
459 Global::intersectPoint[5].x = cSm->p[0].x + (cos(clAngle - a) * cSm->radius[0]);
460 Global::intersectPoint[5].y = cSm->p[0].y + (sin(clAngle - a) * cSm->radius[0]);
462 Global::intersectPoint[6].x = cLg->p[0].x + (cos(clAngle - a) * cLg->radius[0]);
463 Global::intersectPoint[6].y = cLg->p[0].y + (sin(clAngle - a) * cLg->radius[0]);
464 Global::intersectPoint[7].x = cSm->p[0].x + (cos(clAngle + a) * cSm->radius[0]);
465 Global::intersectPoint[7].y = cSm->p[0].y + (sin(clAngle + a) * cSm->radius[0]);
466 Global::numIntersectPoints = 8;
470 // Parameter 1: point in question
471 // Parameter 2, 3: points we are comparing to
473 Point Geometry::NearestTo(Point point, Point p1, Point p2)
475 double l1 = Vector::Magnitude(point, p1);
476 double l2 = Vector::Magnitude(point, p2);
478 return (l1 < l2 ? p1 : p2);