Shamusworld >> Repos - architektonas/blob - src/geometry.cpp

   1 // geometry.cpp: Algebraic geometry helper functions
   2 //
   3 // Part of the Architektonas Project
   4 // (C) 2011 Underground Software
   5 // See the README and GPLv3 files for licensing and warranty information
   6 //
   7 // JLH = James Hammons <jlhamm@acm.org>
   8 //
   9 // WHO  WHEN        WHAT
  10 // ---  ----------  ------------------------------------------------------------
  11 // JLH  08/31/2011  Created this file
  12 //
  13 // NOTE: All methods in this class are static.
  14 //
  15
  16 #include "geometry.h"
  17 #include <math.h>
  18 #include <stdio.h>
  19 #include "global.h"
  20 #include "mathconstants.h"
  21
  22
  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)
  27 {
  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.
  33
  34         Vector lineSegment = head - tail;
  35         double magnitude = lineSegment.Magnitude();
  36         Vector pointSegment = point - tail;
  37         double t = lineSegment.Dot(pointSegment) / (magnitude * magnitude);
  38         return t;
  39 }
  40
  41
  42 Point Geometry::MirrorPointAroundLine(Point point, Point tail, Point head)
  43 {
  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;
  48
  49         // Get the point normal to point to the line passed in
  50         Point normalOnLine = tail + v;
  51
  52         // Make our mirrored vector (head - tail)
  53         Vector mirror = -(point - normalOnLine);
  54
  55         // Find the mirrored point
  56         Point mirroredPoint = normalOnLine + mirror;
  57
  58         return mirroredPoint;
  59 }
  60
  61
  62 //
  63 // point: The point we're rotating
  64 // rotationPoint: The point we're rotating around
  65 //
  66 Point Geometry::RotatePointAroundPoint(Point point, Point rotationPoint, double angle)
  67 {
  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));
  71
  72         return Vector(rotationPoint.x + px, rotationPoint.y + py, 0);
  73 }
  74
  75
  76 double Geometry::Determinant(Point p1, Point p2)
  77 {
  78         return (p1.x * p2.y) - (p2.x * p1.y);
  79 }
  80
  81
  82 void Geometry::Intersects(Object * obj1, Object * obj2)
  83 {
  84         Global::numIntersectPoints = Global::numIntersectParams = 0;
  85
  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);
  94 }
  95
  96
  97 /*
  98 Intersecting line segments:
  99 An easier way:
 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
 107
 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?))
 109
 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?))
 111
 112 So check if the above two numbers are both >=0 and <=1.
 113 */
 114
 115
 116 // Finds the intersection between two lines (if any)
 117 void Geometry::CheckLineToLineIntersection(Object * l1, Object * l2)
 118 {
 119         Global::numIntersectPoints = Global::numIntersectParams = 0;
 120
 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
 124
 125         double rxs = (r.x * s.y) - (s.x * r.y);
 126         double t, u;
 127
 128         if (rxs == 0)
 129         {
 130                 double qpxr = (v1.x * r.y) - (r.x * v1.y);
 131
 132                 // Lines are parallel, so no intersection...
 133                 if (qpxr != 0)
 134                         return;
 135
 136                 // Check to see which endpoints are connected... Four possibilities:
 137                 if (l1->p[0] == l2->p[0])
 138                         t = 0, u = 0;
 139                 else if (l1->p[0] == l2->p[1])
 140                         t = 0, u = 1.0;
 141                 else if (l1->p[1] == l2->p[0])
 142                         t = 1.0, u = 0;
 143                 else if (l1->p[1] == l2->p[1])
 144                         t = 1.0, u = 1.0;
 145                 else
 146                         return;
 147         }
 148         else
 149         {
 150                 t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
 151                 u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
 152         }
 153
 154         Global::intersectParam[0] = t;
 155         Global::intersectParam[1] = u;
 156
 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;
 160 }
 161
 162
 163 void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2)
 164 {
 165         // Set up global vars
 166         Global::numIntersectPoints = Global::numIntersectParams = 0;
 167
 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();
 172
 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])))
 177                 return;
 178
 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
 181         // to both circles.
 182         if (d == (c1->radius[0] + c2->radius[0]))
 183         {
 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;
 187                 return;
 188         }
 189         else if (d == fabs(c1->radius[0] - c2->radius[0]))
 190         {
 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;
 195                 return;
 196         }
 197
 198 /*
 199        c² = a² + b² - 2ab·cos µ
 200 2ab·cos µ = a² + b² - c²
 201     cos µ = (a² + b² - c²) / 2ab
 202 */
 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));
 206
 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;
 214 }
 215
 216
 217 //
 218 // N.B.: l is the line, c is the circle
 219 //
 220 void Geometry::CheckLineToCircleIntersection(Object * l, Object * c)
 221 {
 222         // Set up global vars
 223         Global::numIntersectPoints = Global::numIntersectParams = 0;
 224
 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();
 230
 231         // Step 2: See if we have 0, 1, or 2 intersection points
 232
 233         // Case #1: No intersection points
 234         if (distance > c->radius[0])
 235                 return;
 236         // Case #2: One intersection point (possibly--tangent)
 237         else if (distance == c->radius[0])
 238         {
 239                 // Only intersects if the parameter is on the line segment!
 240                 if ((t >= 0.0) && (t <= 1.0))
 241                 {
 242                         Global::intersectPoint[0] = c->p[0] + radial;
 243                         Global::numIntersectPoints = 1;
 244                 }
 245
 246                 return;
 247         }
 248
 249         // Case #3: Two intersection points (possibly--secant)
 250
 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));
 257
 258         // Now, find the intersection points using the length...
 259         Vector lineUnit = Vector(l->p[0], l->p[1]).Unit();
 260         Point i1 = p + (lineUnit * perpendicularLength);
 261         Point i2 = p - (lineUnit * perpendicularLength);
 262
 263         // Next we need to see if they are on the line segment...
 264         double u = ParameterOfLineAndPoint(l->p[0], l->p[1], i1);
 265         double v = ParameterOfLineAndPoint(l->p[0], l->p[1], i2);
 266
 267         if ((u >= 0.0) && (u <= 1.0))
 268         {
 269                 Global::intersectPoint[Global::numIntersectPoints] = i1;
 270                 Global::numIntersectPoints++;
 271         }
 272
 273         if ((v >= 0.0) && (v <= 1.0))
 274         {
 275                 Global::intersectPoint[Global::numIntersectPoints] = i2;
 276                 Global::numIntersectPoints++;
 277         }
 278 }
 279
 280
 281 // should we just do common trig solves, like AAS, ASA, SAS, SSA?
 282 // Law of Cosines:
 283 // c² = a² + b² - 2ab * cos(C)
 284 // Solving for C:
 285 // cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab
 286 // Law of Sines:
 287 // a / sin A = b / sin B = c / sin C
 288
 289 // Solve the angles of the triangle given the sides. Angles returned are
 290 // opposite of the given sides (so a1 consists of sides s2 & s3, and so on).
 291 void Geometry::FindAnglesForSides(double s1, double s2, double s3, double * a1, double * a2, double * a3)
 292 {
 293         // Use law of cosines to find 1st angle
 294         double cosine1 = ((s2 * s2) + (s3 * s3) - (s1 * s1)) / (2.0 * s2 * s3);
 295
 296         // Check for a valid triangle
 297         if ((cosine1 < -1.0) || (cosine1 > 1.0))
 298                 return;
 299
 300         double angle1 = acos(cosine1);
 301
 302         // Use law of sines to find 2nd & 3rd angles
 303 // sin A / a = sin B / b
 304 // sin B = (sin A / a) * b
 305 // B = arcsin( sin A * (b / a))
 306 // ??? ==> B = A * arcsin(b / a)
 307 /*
 308 Well, look here:
 309 sin B = sin A * (b / a)
 310 sin B / sin A = b / a
 311 arcsin( sin B / sin A ) = arcsin( b / a )
 312
 313 hmm... dunno...
 314 */
 315
 316         double angle2 = asin(s2 * (sin(angle1) / s1));
 317         double angle3 = asin(s3 * (sin(angle1) / s1));
 318
 319         if (a1)
 320                 *a1 = angle1;
 321
 322         if (a2)
 323                 *a2 = angle2;
 324
 325         if (a3)
 326                 *a3 = angle3;
 327 }
 328
 329
 330 Point Geometry::GetPointForParameter(Object * obj, double t)
 331 {
 332         if (obj->type == OTLine)
 333         {
 334                 // Translate line vector to the origin, then add the scaled vector to
 335                 // initial point of the line.
 336                 Vector v = obj->p[1] - obj->p[0];
 337                 return obj->p[0] + (v * t);
 338         }
 339
 340         return Point(0, 0);
 341 }
 342
 343
 344 /*
 345 How to find the tangent of a point off a circle:
 346
 347  •  Calculate the midpoint on the point and the center of the circle
 348  •  Get the length of the line segment from the and the center divided by two
 349  •  Use that length to construct a circle with the point at the center and the
 350     radius equal to that length
 351  •  The intersection of the two circles are the tangent points
 352
 353 */
 354