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"
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 Point Geometry::MirrorPointAroundLine(Point point, Point tail, Point head)
44 // Get the vector of the intersection of the line and the normal on the
45 // line to the point in question.
46 double t = ParameterOfLineAndPoint(tail, head, point);
47 Vector v = Vector(tail, head) * t;
49 // Get the point normal to point to the line passed in
50 Point normalOnLine = tail + v;
52 // Make our mirrored vector (head - tail)
53 Vector mirror = -(point - normalOnLine);
55 // Find the mirrored point
56 Point mirroredPoint = normalOnLine + mirror;
63 // point: The point we're rotating
64 // rotationPoint: The point we're rotating around
66 Point Geometry::RotatePointAroundPoint(Point point, Point rotationPoint, double angle)
68 Vector v = Vector(rotationPoint, point);
69 double px = (v.x * cos(angle)) - (v.y * sin(angle));
70 double py = (v.x * sin(angle)) + (v.y * cos(angle));
72 return Vector(rotationPoint.x + px, rotationPoint.y + py, 0);
76 double Geometry::Determinant(Point p1, Point p2)
78 return (p1.x * p2.y) - (p2.x * p1.y);
82 void Geometry::Intersects(Object * obj1, Object * obj2)
84 Global::numIntersectPoints = Global::numIntersectParams = 0;
86 if ((obj1->type == OTLine) && (obj2->type == OTLine))
87 CheckLineToLineIntersection(obj1, obj2);
88 else if ((obj1->type == OTCircle) && (obj2->type == OTCircle))
89 CheckCircleToCircleIntersection(obj1, obj2);
90 else if ((obj1->type == OTLine) && (obj2->type == OTCircle))
91 CheckLineToCircleIntersection(obj1, obj2);
92 else if ((obj1->type == OTCircle) && (obj2->type == OTLine))
93 CheckLineToCircleIntersection(obj2, obj1);
98 Intersecting line segments:
100 Segment L1 has edges A=(a1,a2), A'=(a1',a2').
101 Segment L2 has edges B=(b1,b2), B'=(b1',b2').
102 Segment L1 is the set of points tA'+(1-t)A, where 0<=t<=1.
103 Segment L2 is the set of points sB'+(1-s)B, where 0<=s<=1.
104 Segment L1 meet segment L2 if and only if for some t and s we have
105 tA'+(1-t)A=sB'+(1-s)B
106 The solution of this with respect to t and s is
108 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?))
110 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?))
112 So check if the above two numbers are both >=0 and <=1.
116 // Finds the intersection between two lines (if any)
117 void Geometry::CheckLineToLineIntersection(Object * l1, Object * l2)
119 Global::numIntersectPoints = Global::numIntersectParams = 0;
121 Vector r(l1->p[0], l1->p[1]);
122 Vector s(l2->p[0], l2->p[1]);
123 Vector v1 = l2->p[0] - l1->p[0]; // q - p
125 double rxs = (r.x * s.y) - (s.x * r.y);
130 double qpxr = (v1.x * r.y) - (r.x * v1.y);
132 // Lines are parallel, so no intersection...
136 // Check to see which endpoints are connected... Four possibilities:
137 if (l1->p[0] == l2->p[0])
139 else if (l1->p[0] == l2->p[1])
141 else if (l1->p[1] == l2->p[0])
143 else if (l1->p[1] == l2->p[1])
150 t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
151 u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
154 Global::intersectParam[0] = t;
155 Global::intersectParam[1] = u;
157 // If the parameters are in range, we have overlap!
158 if ((t >= 0) && (t <= 1.0) && (u >= 0) && (u <= 1.0))
159 Global::numIntersectParams = 1;
163 void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2)
165 // Set up global vars
166 Global::numIntersectPoints = Global::numIntersectParams = 0;
168 // Get the distance between the centers of the circles
169 Vector centerLine(c1->p[0], c2->p[0]);
170 double d = centerLine.Magnitude();
171 double clAngle = centerLine.Angle();
173 // If the distance between centers is greater than the sum of the radii or
174 // less than the difference between the radii, there is NO intersection
175 if ((d > (c1->radius[0] + c2->radius[0]))
176 || (d < fabs(c1->radius[0] - c2->radius[0])))
179 // If the distance between centers is equal to the sum of the radii or
180 // equal to the difference between the radii, the intersection is tangent
182 if (d == (c1->radius[0] + c2->radius[0]))
184 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0]);
185 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0]);
186 Global::numIntersectPoints = 1;
189 else if (d == fabs(c1->radius[0] - c2->radius[0]))
191 double sign = (c1->radius[0] > c2->radius[0] ? +1 : -1);
192 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0] * sign);
193 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0] * sign);
194 Global::numIntersectPoints = 1;
199 c² = a² + b² - 2ab·cos µ
200 2ab·cos µ = a² + b² - c²
201 cos µ = (a² + b² - c²) / 2ab
203 // Use the Law of Cosines to find the angle between the centerline and the
204 // radial line on Circle #1
205 double a = acos(((c1->radius[0] * c1->radius[0]) + (d * d) - (c2->radius[0] * c2->radius[0])) / (2.0 * c1->radius[0] * d));
207 // Finally, find the points of intersection by using +/- the angle found
208 // from the centerline's angle
209 Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle + a) * c1->radius[0]);
210 Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle + a) * c1->radius[0]);
211 Global::intersectPoint[1].x = c1->p[0].x + (cos(clAngle - a) * c1->radius[0]);
212 Global::intersectPoint[1].y = c1->p[0].y + (sin(clAngle - a) * c1->radius[0]);
213 Global::numIntersectPoints = 2;
218 // N.B.: l is the line, c is the circle
220 void Geometry::CheckLineToCircleIntersection(Object * l, Object * c)
222 // Set up global vars
223 Global::numIntersectPoints = Global::numIntersectParams = 0;
225 // Step 1: Find shortest distance from center of circle to the infinite line
226 double t = ParameterOfLineAndPoint(l->p[0], l->p[1], c->p[0]);
227 Point p = l->p[0] + (Vector(l->p[0], l->p[1]) * t);
228 Vector radial = Vector(c->p[0], p);
229 double distance = radial.Magnitude();
231 // Step 2: See if we have 0, 1, or 2 intersection points
233 // Case #1: No intersection points
234 if (distance > c->radius[0])
236 // Case #2: One intersection point (possibly--tangent)
237 else if (distance == c->radius[0])
239 // Only intersects if the parameter is on the line segment!
240 if ((t >= 0.0) && (t <= 1.0))
242 Global::intersectPoint[0] = c->p[0] + radial;
243 Global::numIntersectPoints = 1;
249 // Case #3: Two intersection points (possibly--secant)
251 // So, we have the line, and the perpendicular from the center of the
252 // circle to the line. Now figure out where the intersection points are.
253 // This is a right triangle, though do we really know all the sides?
254 // Don't need to, 2 is enough for Pythagoras :-)
255 // Radius is the hypotenuse, so we have to use c² = a² + b² => a² = c² - b²
256 double perpendicularLength = sqrt((c->radius[0] * c->radius[0]) - (distance * distance));
258 // Now, find the points using the length, then check to see if they are on
260 Vector lineUnit = Vector(l->p[0], l->p[1]).Unit();
261 Point i1 = p + (lineUnit * perpendicularLength);
262 Point i2 = p - (lineUnit * perpendicularLength);
264 // Now we have our intersection points, next we need to see if they are on
265 // the line segment...
266 double u = ParameterOfLineAndPoint(l->p[0], l->p[1], i1);
267 double v = ParameterOfLineAndPoint(l->p[0], l->p[1], i2);
269 if ((u >= 0.0) && (u <= 1.0))
271 Global::intersectPoint[Global::numIntersectPoints] = i1;
272 Global::numIntersectPoints++;
275 if ((v >= 0.0) && (v <= 1.0))
277 Global::intersectPoint[Global::numIntersectPoints] = i2;
278 Global::numIntersectPoints++;
283 // should we just do common trig solves, like AAS, ASA, SAS, SSA?
285 // c² = a² + b² - 2ab * cos(C)
287 // cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab
289 // a / sin A = b / sin B = c / sin C
291 // Solve the angles of the triangle given the sides. Angles returned are
292 // opposite of the given sides (so a1 consists of sides s2 & s3, and so on).
293 void Geometry::FindAnglesForSides(double s1, double s2, double s3, double * a1, double * a2, double * a3)
295 // Use law of cosines to find 1st angle
296 double cosine1 = ((s2 * s2) + (s3 * s3) - (s1 * s1)) / (2.0 * s2 * s3);
298 // Check for a valid triangle
299 if ((cosine1 < -1.0) || (cosine1 > 1.0))
302 double angle1 = acos(cosine1);
304 // Use law of sines to find 2nd & 3rd angles
305 // sin A / a = sin B / b
306 // sin B = (sin A / a) * b
307 // B = arcsin( sin A * (b / a))
308 // ??? ==> B = A * arcsin(b / a)
311 sin B = sin A * (b / a)
312 sin B / sin A = b / a
313 arcsin( sin B / sin A ) = arcsin( b / a )
318 double angle2 = asin(s2 * (sin(angle1) / s1));
319 double angle3 = asin(s3 * (sin(angle1) / s1));
332 Point Geometry::GetPointForParameter(Object * obj, double t)
334 if (obj->type == OTLine)
336 // Translate line vector to the origin, then add the scaled vector to
337 // initial point of the line.
338 Vector v = obj->p[1] - obj->p[0];
339 return obj->p[0] + (v * t);