1 // geometry.cpp: Algebraic geometry helper functions
3 // Part of the Architektonas Project
4 // (C) 2011 Underground Software
5 // See the README and GPLv3 files for licensing and warranty information
7 // JLH = James Hammons <jlhamm@acm.org>
10 // --- ---------- ------------------------------------------------------------
11 // JLH 08/31/2011 Created this file
13 // NOTE: All methods in this class are static.
20 #include "mathconstants.h"
25 Point Geometry::IntersectionOfLineAndLine(Point p1, Point p2, Point p3, Point p4)
27 // Find the intersection of the lines by formula:
28 // px = (x1y2 - y1x2)(x3 - x4) - (x1 - x2)(x3y4 - y3x4)
29 // py = (x1y2 - y1x2)(y3 - y4) - (y1 - y2)(x3y4 - y3x4)
30 // d = (x1 - x2)(y3 - y4) - (y1 - y2)(x3 - x4) = 0 if lines are parallel
31 // Intersection is (px / d, py / d)
33 double d = ((p1.x - p2.x) * (p3.y - p4.y)) - ((p1.y - p2.y) * (p3.x - p4.x));
35 // Check for parallel lines, and return sentinel if so
37 return Point(0, 0, -1);
39 double px = (((p1.x * p2.y) - (p1.y * p2.x)) * (p3.x - p4.x))
40 - ((p1.x - p2.x) * ((p3.x * p4.y) - (p3.y * p4.x)));
41 double py = (((p1.x * p2.y) - (p1.y * p2.x)) * (p3.y - p4.y))
42 - ((p1.y - p2.y) * ((p3.x * p4.y) - (p3.y * p4.x)));
44 return Point(px / d, py / d, 0);
49 // Returns the parameter of a point in space to this vector. If the parameter
50 // is between 0 and 1, the normal of the vector to the point is on the vector.
51 // Note: lp1 is the tail, lp2 is the head of the line (vector).
52 double Geometry::ParameterOfLineAndPoint(Point tail, Point head, Point point)
54 // Geometric interpretation:
55 // The parameterized point on the vector lineSegment is where the normal of
56 // the lineSegment to the point intersects lineSegment. If the pp < 0, then
57 // the perpendicular lies beyond the 1st endpoint. If pp > 1, then the
58 // perpendicular lies beyond the 2nd endpoint.
60 Vector lineSegment = head - tail;
61 double magnitude = lineSegment.Magnitude();
62 Vector pointSegment = point - tail;
63 double t = lineSegment.Dot(pointSegment) / (magnitude * magnitude);
68 Point Geometry::MirrorPointAroundLine(Point point, Point tail, Point head)
70 // Get the vector of the intersection of the line and the normal on the
71 // line to the point in question.
72 double t = ParameterOfLineAndPoint(tail, head, point);
73 Vector v = Vector(tail, head) * t;
75 // Get the point normal to point to the line passed in
76 Point normalOnLine = tail + v;
78 // Make our mirrored vector (head - tail)
79 Vector mirror = -(point - normalOnLine);
81 // Find the mirrored point
82 Point mirroredPoint = normalOnLine + mirror;
89 // point: The point we're rotating
90 // rotationPoint: The point we're rotating around
92 Point Geometry::RotatePointAroundPoint(Point point, Point rotationPoint, double angle)
94 Vector v = Vector(rotationPoint, point);
95 double px = (v.x * cos(angle)) - (v.y * sin(angle));
96 double py = (v.x * sin(angle)) + (v.y * cos(angle));
98 return Vector(rotationPoint.x + px, rotationPoint.y + py, 0);
102 double Geometry::Determinant(Point p1, Point p2)
104 return (p1.x * p2.y) - (p2.x * p1.y);
108 void Geometry::Intersects(Object * obj1, Object * obj2)//, double * tp/*= 0*/, double * up/*= 0*/, double * vp/*= 0*/, double * wp/*= 0*/)
110 Global::numIntersectPoints = Global::numIntersectParams = 0;
112 if ((obj1->type == OTLine) && (obj2->type == OTLine))
113 CheckLineToLineIntersection(obj1, obj2);
114 else if ((obj1->type == OTCircle) && (obj2->type == OTCircle))
115 CheckCircleToCircleIntersection(obj1, obj2);
120 Intersecting line segments:
122 Segment L1 has edges A=(a1,a2), A'=(a1',a2').
123 Segment L2 has edges B=(b1,b2), B'=(b1',b2').
124 Segment L1 is the set of points tA'+(1-t)A, where 0<=t<=1.
125 Segment L2 is the set of points sB'+(1-s)B, where 0<=s<=1.
126 Segment L1 meet segment L2 if and only if for some t and s we have
127 tA'+(1-t)A=sB'+(1-s)B
128 The solution of this with respect to t and s is
130 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?))
132 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?))
134 So check if the above two numbers are both >=0 and <=1.
138 // Finds the intersection between two lines (if any)
139 void Geometry::CheckLineToLineIntersection(Object * l1, Object * l2)//, double * tp, double * up)
141 Global::numIntersectPoints = Global::numIntersectParams = 0;
143 Vector r(l1->p[0], l1->p[1]);
144 Vector s(l2->p[0], l2->p[1]);
145 Vector v1 = l2->p[0] - l1->p[0]; // q - p
147 Vector v2 = l1->p[0] - l2->p[0]; // p - q
148 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);
150 double rxs = (r.x * s.y) - (s.x * r.y);
155 double qpxr = (v1.x * r.y) - (r.x * v1.y);
158 printf(" --> R x S = 0! (q - p) x r = %lf\n", qpxr);
159 printf(" -->(q - p) . r = %lf, r . r = %lf\n", v1.Dot(r), r.Dot(r));
160 printf(" -->(p - q) . s = %lf, s . s = %lf\n", v2.Dot(s), s.Dot(s));
161 printf(" -->(q - p) . s = %lf, (p - q) . r = %lf\n", v1.Dot(s), v2.Dot(r));
164 // Lines are parallel, so no intersection...
169 //this works IFF the vectors are pointing in the same direction. everything else
171 // If (q - p) . r == r . r, t = 1, u = 0
172 if (v1.Dot(r) == r.Dot(r))
174 // If (p - q) . s == s . s, t = 0, u = 1
175 else if (v2.Dot(s) == s.Dot(s))
180 // Check to see which endpoints are connected... Four possibilities:
181 if (l1->p[0] == l2->p[0])
183 else if (l1->p[0] == l2->p[1])
185 else if (l1->p[1] == l2->p[0])
187 else if (l1->p[1] == l2->p[1])
195 t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
196 u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
199 Now there are five cases (NOTE: only valid if vectors face the same way!):
201 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.
203 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.
205 3. If r × s = 0 and (q − p) × r ≠ 0, then the two lines are parallel and non-intersecting.
207 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.
209 5. Otherwise, the two line segments are not parallel but do not intersect.
211 // Return parameter values, if user passed in valid pointers
219 // If the parameters are in range, we have overlap!
220 if ((t >= 0) && (t <= 1.0) && (u >= 0) && (u <= 1.0))
225 Global::intersectParam[0] = t;
226 Global::intersectParam[1] = u;
228 // If the parameters are in range, we have overlap!
229 if ((t >= 0) && (t <= 1.0) && (u >= 0) && (u <= 1.0))
230 Global::numIntersectParams = 1;
235 void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2)
237 // Set up global vars
238 Global::numIntersectPoints = Global::numIntersectParams = 0;
240 // Get the distance between the centers of the circles
241 Vector centerLine(c1->p[0], c2->p[0]);
242 double d = centerLine.Magnitude();
243 double clAngle = centerLine.Angle();
245 // If the distance between centers is greater than the sum of the radii or
246 // less than the difference between the radii, there is NO intersection
247 if ((d > (c1->radius[0] + c2->radius[0]))
248 || (d < fabs(c1->radius[0] - c2->radius[0])))
251 // If the distance between centers is equal to the sum of the radii or
252 // equal to the difference between the radii, the intersection is tangent
254 if ((d == (c1->radius[0] + c2->radius[0]))
255 || (d == fabs(c1->radius[0] - c2->radius[0])))
257 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0]);
258 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0]);
259 Global::numIntersectPoints = 1;
263 // Use the Law of Cosines to find the angle between the centerline and the
264 // radial line on Circle #1
265 double a = acos(((c1->radius[0] * c1->radius[0]) + (d * d) - (c2->radius[0] * c2->radius[0])) / (2.0 * c1->radius[0] * d));
267 // Finally, find the points of intersection by using +/- the angle found
268 // from the centerline's angle
269 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle + a) * c1->radius[0]);
270 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle + a) * c1->radius[0]);
271 Global::intersectPoint[1].x = c1->p[0].x + (cos(clAngle - a) * c1->radius[0]);
272 Global::intersectPoint[1].y = c1->p[0].y + (sin(clAngle - a) * c1->radius[0]);
273 Global::numIntersectPoints = 2;
278 // Finds the intersection between two lines (if any)
279 int Geometry::Intersects(Line * l1, Dimension * d1, double * tp/*= 0*/, double * up/*= 0*/)
281 Line l2(d1->position, d1->endpoint);
282 return Intersects(l1, &l2, tp, up);
286 // Finds the intersection(s) between a line and a circle (if any)
287 int Geometry::Intersects(Line * l, Circle * c, double * tp/*= 0*/, double * up/*= 0*/, double * vp/*= 0*/, double * wp/*= 0*/)
290 Vector center = c->position;
291 Vector v1 = l->position - center;
292 Vector v2 = l->endpoint - center;
294 double dr = d.Magnitude();
295 double determinant = (v1.x * v2.y) - (v1.y * v2.x);
297 double discriminant = ((c->radius * c->radius) * (dr * dr)) - (determinant * determinant);
299 if (discriminant < 0)
307 I'm thinking a better approach to this might be as follows:
309 -- Get the distance of the circle's center from the line segment. If it's
310 > the radius, it doesn't intersect.
311 -- If the parameter is off the line segment, check distance to endpoints. (Not sure
312 how to proceed from here, it's different than the following.)
313 [Actually, you can use the following for all of it. You only know if you have
314 an intersection at the last step, which is OK.]
315 -- If the radius == distance, we have a tangent line.
316 -- If radius > distance, use Pythagorus to find the length on either side of the
317 normal to the spots where the hypotenuse (== radius' length) contact the line.
318 -- Use those points to find the parameter on the line segment; if they're not on
319 the line segment, no intersection.
321 double t = ParameterOfLineAndPoint(l->position, l->endpoint, c->position);
322 //small problem here: it clamps the result to the line segment. NOT what we want
323 //here! !!! FIX !!! [DONE]
324 Vector p = l->GetPointAtParameter(t);
325 double distance = Vector::Magnitude(c->position, p);
327 // If the center of the circle is farther from the line than the radius, fail.
328 if (distance > c->radius)
331 // Now we have to check for intersection points.
332 // Tangent case: (needs to return something)
333 if ((distance == c->radius) && (t >= 0.0) && (t <= 1.0))
335 // Need to set tp & up to something... !!! FIX !!!
340 *up = Vector(c->position, p).Angle();
345 // The line intersects the circle in two points (possibly). Use Pythagorus
346 // to find them for testing.
347 double offset = sqrt((c->radius * c->radius) - (distance * distance));
348 //need to convert distance to paramter value... :-/
349 //t = position on line / length of line segment, so if we divide the offset by length,
350 //that should give us what we want.
351 double length = Vector::Magnitude(l->position, l->endpoint);
352 double t1 = t + (offset / length);
353 double t2 = t - (offset / length);
355 //need to find angles for the circle...
356 Vector cp1 = l->position + (Vector(l->position, l->endpoint) * (length * t1));
357 Vector cp2 = l->position + (Vector(l->position, l->endpoint) * (length * t2));
358 double a1 = Vector(c->position, cp1).Angle();
359 double a2 = Vector(c->position, cp2).Angle();
361 //instead of this, return a # which is the # of intersections. [DONE]
362 int intersections = 0;
364 // Now check for if the parameters are in range
365 if ((t1 >= 0) && (t1 <= 1.0))
370 if ((t2 >= 0) && (t2 <= 1.0))
375 return intersections;
380 // Finds the intersection(s) between a circle and a circle (if any)
381 // There can be 0, 1, or 2 intersections.
382 // Returns the angles of the points of intersection in tp thru wp, with the
383 // angles returned as c1, c2, c1, c2 (if applicable--in the 1 intersection case,
384 // only the first two angles are returned: c1, c2).
385 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*/)
387 // Get the distance between centers. If the distance plus the radius of the
388 // smaller circle is less than the radius of the larger circle, there is no
389 // intersection. If the distance is greater than the sum of the radii,
390 // there is no intersection. If the distance is equal to the sum of the
391 // radii, they are tangent and intersect at one point. Otherwise, they
392 // intersect at two points.
393 Vector centerLine(c1->position, c2->position);
394 double d = centerLine.Magnitude();
395 //printf("Circle #1: pos=<%lf, %lf>, r=%lf\n", c1->position.x, c1->position.y, c1->radius);
396 //printf("Circle #2: pos=<%lf, %lf>, r=%lf\n", c2->position.x, c2->position.y, c2->radius);
397 //printf("Distance between #1 & #2: %lf\n", d);
399 // Check to see if we actually have an intersection, and return failure if not
400 if ((fabs(c1->radius - c2->radius) > d) || ((c1->radius + c2->radius) < d))
403 // There are *two* tangent cases!
404 if (((c1->radius + c2->radius) == d) || (fabs(c1->radius - c2->radius) == d))
406 // Need to return something in tp & up!! !!! FIX !!! [DONE]
408 *tp = centerLine.Angle();
411 *up = centerLine.Angle() + PI;
416 // Find the distance from the center of c1 to the perpendicular chord
417 // (which contains the points of intersection)
418 // [N.B.: This is derived from Pythagorus by using the unknown distance
419 // from the center line to the point where the two radii coincide as
420 // a common unknown to two instances of the formula.]
421 double x = ((d * d) - (c2->radius * c2->radius) + (c1->radius * c1->radius))
423 // Find the the length of the perpendicular chord
425 double a = sqrt((-d + c2->radius - c1->radius) * (-d - c2->radius + c1->radius) * (-d + c2->radius + c1->radius) * (d + c2->radius + c1->radius)) / d;
427 // Now, you can use pythagorus to find the length of the hypotenuse, but we
428 // already know that length, it's the radius! :-P
429 // What's needed is the angle of the center line and the radial line. Since
430 // there's two intersection points, there's also four angles (two for each
432 // We can use the arccos to find the angle using just the radius and the
433 // distance to the perpendicular chord...!
434 double angleC1 = acos(x / c1->radius);
435 double angleC2 = acos((d - x) / c2->radius);
438 *tp = centerLine.Angle() - angleC1;
441 *up = (centerLine.Angle() + PI) - angleC2;
444 *vp = centerLine.Angle() + angleC1;
447 *wp = (centerLine.Angle() + PI) + angleC2;
450 *p1 = c1->position + (centerLine.Unit() * x) + (Vector::Normal(Vector(), centerLine) * (a / 2.0));
453 *p2 = c1->position + (centerLine.Unit() * x) - (Vector::Normal(Vector(), centerLine) * (a / 2.0));
459 // should we just do common trig solves, like AAS, ASA, SAS, SSA?
461 // c² = a² + b² - 2ab * cos(C)
463 // cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab
465 // a / sin A = b / sin B = c / sin C
467 // Solve the angles of the triangle given the sides. Angles returned are
468 // opposite of the given sides (so a1 consists of sides s2 & s3, and so on).
469 void Geometry::FindAnglesForSides(double s1, double s2, double s3, double * a1, double * a2, double * a3)
471 // Use law of cosines to find 1st angle
472 double cosine1 = ((s2 * s2) + (s3 * s3) - (s1 * s1)) / (2.0 * s2 * s3);
474 // Check for a valid triangle
475 if ((cosine1 < -1.0) || (cosine1 > 1.0))
478 double angle1 = acos(cosine1);
480 // Use law of sines to find 2nd & 3rd angles
481 // sin A / a = sin B / b
482 // sin B = (sin A / a) * b
483 // B = arcsin( sin A * (b / a))
484 // ??? ==> B = A * arcsin(b / a)
487 sin B = sin A * (b / a)
488 sin B / sin A = b / a
489 arcsin( sin B / sin A ) = arcsin( b / a )
494 double angle2 = asin(s2 * (sin(angle1) / s1));
495 double angle3 = asin(s3 * (sin(angle1) / s1));
508 Point Geometry::GetPointForParameter(Object * obj, double t)
510 if (obj->type == OTLine)
512 // Translate line vector to the origin, then add the scaled vector to
513 // initial point of the line.
514 Vector v = obj->p[1] - obj->p[0];
515 return obj->p[0] + (v * t);