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.
22 Point Geometry::IntersectionOfLineAndLine(Point p1, Point p2, Point p3, Point p4)
24 // Find the intersection of the lines by formula:
25 // px = (x1y2 - y1x2)(x3 - x4) - (x1 - x2)(x3y4 - y3x4)
26 // py = (x1y2 - y1x2)(y3 - y4) - (y1 - y2)(x3y4 - y3x4)
27 // d = (x1 - x2)(y3 - y4) - (y1 - y2)(x3 - x4) = 0 if lines are parallel
28 // Intersection is (px / d, py / d)
30 double d = ((p1.x - p2.x) * (p3.y - p4.y)) - ((p1.y - p2.y) * (p3.x - p4.x));
32 // Check for parallel lines, and return sentinel if so
34 return Point(0, 0, -1);
36 double px = (((p1.x * p2.y) - (p1.y * p2.x)) * (p3.x - p4.x))
37 - ((p1.x - p2.x) * ((p3.x * p4.y) - (p3.y * p4.x)));
38 double py = (((p1.x * p2.y) - (p1.y * p2.x)) * (p3.y - p4.y))
39 - ((p1.y - p2.y) * ((p3.x * p4.y) - (p3.y * p4.x)));
41 return Point(px / d, py / d, 0);
45 // Returns the parameter of a point in space to this vector. If the parameter
46 // is between 0 and 1, the normal of the vector to the point is on the vector.
47 // Note: lp1 is the tail, lp2 is the head of the line (vector).
48 double Geometry::ParameterOfLineAndPoint(Point tail, Point head, Point point)
50 // Geometric interpretation:
51 // The parameterized point on the vector lineSegment is where the normal of
52 // the lineSegment to the point intersects lineSegment. If the pp < 0, then
53 // the perpendicular lies beyond the 1st endpoint. If pp > 1, then the
54 // perpendicular lies beyond the 2nd endpoint.
56 Vector lineSegment = head - tail;
57 double magnitude = lineSegment.Magnitude();
58 Vector pointSegment = point - tail;
59 double t = lineSegment.Dot(pointSegment) / (magnitude * magnitude);
64 Point Geometry::MirrorPointAroundLine(Point point, Point p1, Point p2)
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(p1, p2, point);
69 Vector v = Vector(p1, p2) * t;
71 // Get the point normal to point to the line passed in (p2 is the tail)
72 Point normalOnLine = p2 + 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 Geometry::RotatePointAroundPoint(Point point, Point rotationPoint, double angle)
86 Vector v = Vector(point, rotationPoint);
87 double px = (v.x * cos(angle)) - (v.y * sin(angle));
88 double py = (v.x * sin(angle)) + (v.y * cos(angle));
90 return Vector(rotationPoint.x + px, rotationPoint.y + py, 0);
94 double Geometry::Determinant(Point p1, Point p2)
96 return (p1.x * p2.y) - (p2.x * p1.y);
101 Intersecting line segments:
103 Segment L1 has edges A=(a1,a2), A'=(a1',a2').
104 Segment L2 has edges B=(b1,b2), B'=(b1',b2').
105 Segment L1 is the set of points tA'+(1-t)A, where 0<=t<=1.
106 Segment L2 is the set of points sB'+(1-s)B, where 0<=s<=1.
107 Segment L1 meet segment L2 if and only if for some t and s we have
108 tA'+(1-t)A=sB'+(1-s)B
109 The solution of this with respect to t and s is
111 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?))
113 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?))
115 So check if the above two numbers are both >=0 and <=1.
120 // Finds the intesection between two objects (if any)
121 bool Geometry::Intersects(Object * obj1, Object * obj2, double * t, double * s)
126 // Finds the intersection between two lines (if any)
127 int Geometry::Intersects(Line * l1, Line * l2, double * tp/*= 0*/, double * up/*= 0*/)
129 Vector r(l1->position, l1->endpoint);
130 Vector s(l2->position, l2->endpoint);
131 Vector v1 = l2->position - l1->position;
132 // Vector v1 = l1->position - l2->position;
134 double rxs = (r.x * s.y) - (s.x * r.y);
139 double t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
140 double u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
142 Now there are five cases:
144 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.
146 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.
148 3. If r × s = 0 and (q − p) × r ≠ 0, then the two lines are parallel and non-intersecting.
150 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.
152 5. Otherwise, the two line segments are not parallel but do not intersect.
154 // Return parameter values, if user passed in valid pointers
161 // If the parameters are in range, we have overlap!
162 if ((t >= 0) && (t <= 1.0) && (u >= 0) && (u <= 1.0))
169 // Finds the intesection(s) between a line and a circle (if any)
170 int Geometry::Intersects(Line * l, Circle * c, double * tp/*= 0*/, double * up/*= 0*/, double * vp/*= 0*/, double * wp/*= 0*/)
173 Vector center = c->position;
174 Vector v1 = l->position - center;
175 Vector v2 = l->endpoint - center;
177 double dr = d.Magnitude();
178 double determinant = (v1.x * v2.y) - (v1.y * v2.x);
180 double discriminant = ((c->radius * c->radius) * (dr * dr)) - (determinant * determinant);
182 if (discriminant < 0)
190 I'm thinking a better approach to this might be as follows:
192 -- Get the distance of the circle's center from the line segment. If it's
193 > the radius, it doesn't intersect.
194 -- If the parameter is off the line segment, check distance to endpoints. (Not sure
195 how to proceed from here, it's different than the following.)
196 [Actually, you can use the following for all of it. You only know if you have
197 an intersection at the last step, which is OK.]
198 -- If the radius == distance, we have a tangent line.
199 -- If radius > distance, use Pythagorus to find the length on either side of the
200 normal to the spots where the hypotenuse (== radius' length) contact the line.
201 -- Use those points to find the parameter on the line segment; if they're not on
202 the line segment, no intersection.
204 double t = ParameterOfLineAndPoint(l->position, l->endpoint, c->position);
205 //small problem here: it clamps the result to the line segment. NOT what we want
206 //here! !!! FIX !!! [DONE]
207 Vector p = l->GetPointAtParameter(t);
208 double distance = Vector::Magnitude(c->position, p);
210 // If the center of the circle is farther from the line than the radius, fail.
211 if (distance > c->radius)
214 // Now we have to check for intersection points.
215 // Tangent case: (needs to return something)
216 if ((distance == c->radius) && (t >= 0.0) && (t <= 1.0))
219 // The line intersects the circle in two points (possibly). Use Pythagorus
220 // to find them for testing.
221 double offset = sqrt((c->radius * c->radius) - (distance * distance));
222 //need to convert distance to paramter value... :-/
223 //t = position on line / length of line segment, so if we divide the offset by length,
224 //that should give us what we want.
225 double length = Vector::Magnitude(l->position, l->endpoint);
226 double t1 = t + (offset / length);
227 double t2 = t - (offset / length);
229 //need to find angles for the circle...
230 Vector cp1 = l->position + (Vector(l->position, l->endpoint) * (length * t1));
231 Vector cp2 = l->position + (Vector(l->position, l->endpoint) * (length * t2));
232 double a1 = Vector(c->position, cp1).Angle();
233 double a2 = Vector(c->position, cp2).Angle();
235 //instead of this, return a # which is the # of intersections. [DONE]
236 int intersections = 0;
238 // Now check for if the parameters are in range
239 if ((t1 >= 0) && (t1 <= 1.0))
244 if ((t2 >= 0) && (t2 <= 1.0))
249 return intersections;