X-Git-Url: http://shamusworld.gotdns.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fgeometry.cpp;h=95d12c65068b16715086191e1d70c7e01483dcbd;hb=bf5a50feb0f84a4627a65c5b82c3ca2d2eefe54b;hp=b67c64e30e6211d8cd9dafefbd1204671101f02b;hpb=ea7712f342020baf61cf33ba98b12140da6aecf7;p=architektonas diff --git a/src/geometry.cpp b/src/geometry.cpp index b67c64e..95d12c6 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -16,35 +16,10 @@ #include "geometry.h" #include #include +#include "global.h" #include "mathconstants.h" -// This is unused -#if 0 -Point Geometry::IntersectionOfLineAndLine(Point p1, Point p2, Point p3, Point p4) -{ - // Find the intersection of the lines by formula: - // px = (x1y2 - y1x2)(x3 - x4) - (x1 - x2)(x3y4 - y3x4) - // py = (x1y2 - y1x2)(y3 - y4) - (y1 - y2)(x3y4 - y3x4) - // d = (x1 - x2)(y3 - y4) - (y1 - y2)(x3 - x4) = 0 if lines are parallel - // Intersection is (px / d, py / d) - - double d = ((p1.x - p2.x) * (p3.y - p4.y)) - ((p1.y - p2.y) * (p3.x - p4.x)); - - // Check for parallel lines, and return sentinel if so - if (d == 0) - return Point(0, 0, -1); - - double px = (((p1.x * p2.y) - (p1.y * p2.x)) * (p3.x - p4.x)) - - ((p1.x - p2.x) * ((p3.x * p4.y) - (p3.y * p4.x))); - double py = (((p1.x * p2.y) - (p1.y * p2.x)) * (p3.y - p4.y)) - - ((p1.y - p2.y) * ((p3.x * p4.y) - (p3.y * p4.x))); - - return Point(px / d, py / d, 0); -} -#endif - - // Returns the parameter of a point in space to this vector. If the parameter // is between 0 and 1, the normal of the vector to the point is on the vector. // Note: lp1 is the tail, lp2 is the head of the line (vector). @@ -64,6 +39,29 @@ double Geometry::ParameterOfLineAndPoint(Point tail, Point head, Point point) } +double Geometry::DistanceToLineFromPoint(Point tail, Point head, Point point) +{ + // Interpretation: given a line in the form x = a + tu, where u is the + // unit vector of the line, a is the tail and t is a parameter which + // describes the line, the distance of a point p to the line is given by: + // || (a - p) - ((a - p) . u) u || + // We go an extra step: we set the sign to reflect which side of the line + // it's on (+ == to the left if head points away from you, - == to the + // right) + Vector line(tail, head); + Vector u = line.Unit(); + Vector a_p = tail - point; + Vector dist = a_p - (u * (a_p).Dot(u)); + + double angle = Vector::Angle(tail, point) - line.Angle(); + + if (angle < 0) + angle += TAU; + + return dist.Magnitude() * (angle < HALF_TAU ? +1.0 : -1.0); +} + + Point Geometry::MirrorPointAroundLine(Point point, Point tail, Point head) { // Get the vector of the intersection of the line and the normal on the @@ -104,12 +102,18 @@ double Geometry::Determinant(Point p1, Point p2) } -int Geometry::Intersects(Object * obj1, Object * obj2, double * tp/*= 0*/, double * up/*= 0*/, double * vp/*= 0*/, double * wp/*= 0*/) +void Geometry::Intersects(Object * obj1, Object * obj2) { - if ((obj1->type == OTLine) && (obj2->type == OTLine)) - return CheckLineToLineIntersection(obj1, obj2, tp, up); + Global::numIntersectPoints = Global::numIntersectParams = 0; - return 0; + if ((obj1->type == OTLine) && (obj2->type == OTLine)) + CheckLineToLineIntersection(obj1, obj2); + else if ((obj1->type == OTCircle) && (obj2->type == OTCircle)) + CheckCircleToCircleIntersection(obj1, obj2); + else if ((obj1->type == OTLine) && (obj2->type == OTCircle)) + CheckLineToCircleIntersection(obj1, obj2); + else if ((obj1->type == OTCircle) && (obj2->type == OTLine)) + CheckLineToCircleIntersection(obj2, obj1); } @@ -133,15 +137,14 @@ So check if the above two numbers are both >=0 and <=1. // Finds the intersection between two lines (if any) -int Geometry::CheckLineToLineIntersection(Object * l1, Object * l2, double * tp, double * up) +void Geometry::CheckLineToLineIntersection(Object * l1, Object * l2) { + Global::numIntersectPoints = Global::numIntersectParams = 0; + Vector r(l1->p[0], l1->p[1]); Vector s(l2->p[0], l2->p[1]); Vector v1 = l2->p[0] - l1->p[0]; // q - p -#if 0 - Vector v2 = l1->p[0] - l2->p[0]; // p - q -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); -#endif + double rxs = (r.x * s.y) - (s.x * r.y); double t, u; @@ -149,29 +152,10 @@ printf("l1: (%lf, %lf) (%lf, %lf), l2: (%lf, %lf) (%lf, %lf)\n", l1->p[0].x, l1- { double qpxr = (v1.x * r.y) - (r.x * v1.y); -#if 0 -printf(" --> R x S = 0! (q - p) x r = %lf\n", qpxr); -printf(" -->(q - p) . r = %lf, r . r = %lf\n", v1.Dot(r), r.Dot(r)); -printf(" -->(p - q) . s = %lf, s . s = %lf\n", v2.Dot(s), s.Dot(s)); -printf(" -->(q - p) . s = %lf, (p - q) . r = %lf\n", v1.Dot(s), v2.Dot(r)); -#endif - // Lines are parallel, so no intersection... if (qpxr != 0) - return 0; + return; -#if 0 -//this works IFF the vectors are pointing in the same direction. everything else -//is fucked! - // If (q - p) . r == r . r, t = 1, u = 0 - if (v1.Dot(r) == r.Dot(r)) - t = 1.0, u = 0; - // If (p - q) . s == s . s, t = 0, u = 1 - else if (v2.Dot(s) == s.Dot(s)) - t = 0, u = 1.0; - else - return 0; -#else // Check to see which endpoints are connected... Four possibilities: if (l1->p[0] == l2->p[0]) t = 0, u = 0; @@ -182,226 +166,146 @@ printf(" -->(q - p) . s = %lf, (p - q) . r = %lf\n", v1.Dot(s), v2.Dot(r)); else if (l1->p[1] == l2->p[1]) t = 1.0, u = 1.0; else - return 0; -#endif + return; } else { t = ((v1.x * s.y) - (s.x * v1.y)) / rxs; u = ((v1.x * r.y) - (r.x * v1.y)) / rxs; } -/* -Now there are five cases (NOTE: only valid if vectors face the same way!): - -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. - -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. - -3. If r × s = 0 and (q − p) × r ≠ 0, then the two lines are parallel and non-intersecting. -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. - -5. Otherwise, the two line segments are not parallel but do not intersect. -*/ - // Return parameter values, if user passed in valid pointers - if (tp) - *tp = t; - - if (up) - *up = u; + Global::intersectParam[0] = t; + Global::intersectParam[1] = u; // If the parameters are in range, we have overlap! if ((t >= 0) && (t <= 1.0) && (u >= 0) && (u <= 1.0)) - return 1; - - return 0; + Global::numIntersectParams = 1; } -#if 0 -// Finds the intersection between two lines (if any) -int Geometry::Intersects(Line * l1, Dimension * d1, double * tp/*= 0*/, double * up/*= 0*/) -{ - Line l2(d1->position, d1->endpoint); - return Intersects(l1, &l2, tp, up); -} - - -// Finds the intersection(s) between a line and a circle (if any) -int Geometry::Intersects(Line * l, Circle * c, double * tp/*= 0*/, double * up/*= 0*/, double * vp/*= 0*/, double * wp/*= 0*/) +void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2) { -#if 0 - Vector center = c->position; - Vector v1 = l->position - center; - Vector v2 = l->endpoint - center; - Vector d = v2 - v1; - double dr = d.Magnitude(); - double determinant = (v1.x * v2.y) - (v1.y * v2.x); + // Set up global vars + Global::numIntersectPoints = Global::numIntersectParams = 0; - double discriminant = ((c->radius * c->radius) * (dr * dr)) - (determinant * determinant); - - if (discriminant < 0) - return false; - - - - return true; -#else -/* -I'm thinking a better approach to this might be as follows: - --- Get the distance of the circle's center from the line segment. If it's - > the radius, it doesn't intersect. --- If the parameter is off the line segment, check distance to endpoints. (Not sure - how to proceed from here, it's different than the following.) - [Actually, you can use the following for all of it. You only know if you have - an intersection at the last step, which is OK.] --- If the radius == distance, we have a tangent line. --- If radius > distance, use Pythagorus to find the length on either side of the - normal to the spots where the hypotenuse (== radius' length) contact the line. --- Use those points to find the parameter on the line segment; if they're not on - the line segment, no intersection. -*/ - double t = ParameterOfLineAndPoint(l->position, l->endpoint, c->position); -//small problem here: it clamps the result to the line segment. NOT what we want -//here! !!! FIX !!! [DONE] - Vector p = l->GetPointAtParameter(t); - double distance = Vector::Magnitude(c->position, p); - - // If the center of the circle is farther from the line than the radius, fail. - if (distance > c->radius) - return 0; - - // Now we have to check for intersection points. - // Tangent case: (needs to return something) - if ((distance == c->radius) && (t >= 0.0) && (t <= 1.0)) - { - // Need to set tp & up to something... !!! FIX !!! - if (tp) - *tp = t; - - if (up) - *up = Vector(c->position, p).Angle(); + // Get the distance between the centers of the circles + Vector centerLine(c1->p[0], c2->p[0]); + double d = centerLine.Magnitude(); + double clAngle = centerLine.Angle(); - return 1; - } + // If the distance between centers is greater than the sum of the radii or + // less than the difference between the radii, there is NO intersection + if ((d > (c1->radius[0] + c2->radius[0])) + || (d < fabs(c1->radius[0] - c2->radius[0]))) + return; - // The line intersects the circle in two points (possibly). Use Pythagorus - // to find them for testing. - double offset = sqrt((c->radius * c->radius) - (distance * distance)); -//need to convert distance to paramter value... :-/ -//t = position on line / length of line segment, so if we divide the offset by length, -//that should give us what we want. - double length = Vector::Magnitude(l->position, l->endpoint); - double t1 = t + (offset / length); - double t2 = t - (offset / length); - -//need to find angles for the circle... - Vector cp1 = l->position + (Vector(l->position, l->endpoint) * (length * t1)); - Vector cp2 = l->position + (Vector(l->position, l->endpoint) * (length * t2)); - double a1 = Vector(c->position, cp1).Angle(); - double a2 = Vector(c->position, cp2).Angle(); - -//instead of this, return a # which is the # of intersections. [DONE] - int intersections = 0; - - // Now check for if the parameters are in range - if ((t1 >= 0) && (t1 <= 1.0)) + // If the distance between centers is equal to the sum of the radii or + // equal to the difference between the radii, the intersection is tangent + // to both circles. + if (d == (c1->radius[0] + c2->radius[0])) { - intersections++; + Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0]); + Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0]); + Global::numIntersectPoints = 1; + return; } - - if ((t2 >= 0) && (t2 <= 1.0)) + else if (d == fabs(c1->radius[0] - c2->radius[0])) { - intersections++; + double sign = (c1->radius[0] > c2->radius[0] ? +1 : -1); + Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle) * c1->radius[0] * sign); + Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle) * c1->radius[0] * sign); + Global::numIntersectPoints = 1; + return; } - return intersections; -#endif +/* + c² = a² + b² - 2ab·cos µ +2ab·cos µ = a² + b² - c² + cos µ = (a² + b² - c²) / 2ab +*/ + // Use the Law of Cosines to find the angle between the centerline and the + // radial line on Circle #1 + double a = acos(((c1->radius[0] * c1->radius[0]) + (d * d) - (c2->radius[0] * c2->radius[0])) / (2.0 * c1->radius[0] * d)); + + // Finally, find the points of intersection by using +/- the angle found + // from the centerline's angle + Global::intersectPoint[0].x = c1->p[0].x + (cos(clAngle + a) * c1->radius[0]); + Global::intersectPoint[0].y = c1->p[0].y + (sin(clAngle + a) * c1->radius[0]); + Global::intersectPoint[1].x = c1->p[0].x + (cos(clAngle - a) * c1->radius[0]); + Global::intersectPoint[1].y = c1->p[0].y + (sin(clAngle - a) * c1->radius[0]); + Global::numIntersectPoints = 2; } -// Finds the intersection(s) between a circle and a circle (if any) -// There can be 0, 1, or 2 intersections. -// Returns the angles of the points of intersection in tp thru wp, with the -// angles returned as c1, c2, c1, c2 (if applicable--in the 1 intersection case, -// only the first two angles are returned: c1, c2). -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*/) +// +// N.B.: l is the line, c is the circle +// +void Geometry::CheckLineToCircleIntersection(Object * l, Object * c) { - // Get the distance between centers. If the distance plus the radius of the - // smaller circle is less than the radius of the larger circle, there is no - // intersection. If the distance is greater than the sum of the radii, - // there is no intersection. If the distance is equal to the sum of the - // radii, they are tangent and intersect at one point. Otherwise, they - // intersect at two points. - Vector centerLine(c1->position, c2->position); - double d = centerLine.Magnitude(); -//printf("Circle #1: pos=<%lf, %lf>, r=%lf\n", c1->position.x, c1->position.y, c1->radius); -//printf("Circle #2: pos=<%lf, %lf>, r=%lf\n", c2->position.x, c2->position.y, c2->radius); -//printf("Distance between #1 & #2: %lf\n", d); + // Set up global vars + Global::numIntersectPoints = Global::numIntersectParams = 0; - // Check to see if we actually have an intersection, and return failure if not - if ((fabs(c1->radius - c2->radius) > d) || ((c1->radius + c2->radius) < d)) - return 0; + // Step 1: Find shortest distance from center of circle to the infinite line + double t = ParameterOfLineAndPoint(l->p[0], l->p[1], c->p[0]); + Point p = l->p[0] + (Vector(l->p[0], l->p[1]) * t); + Vector radial = Vector(c->p[0], p); + double distance = radial.Magnitude(); - // There are *two* tangent cases! - if (((c1->radius + c2->radius) == d) || (fabs(c1->radius - c2->radius) == d)) - { - // Need to return something in tp & up!! !!! FIX !!! [DONE] - if (tp) - *tp = centerLine.Angle(); + // Step 2: See if we have 0, 1, or 2 intersection points - if (up) - *up = centerLine.Angle() + PI; + // Case #1: No intersection points + if (distance > c->radius[0]) + return; + // Case #2: One intersection point (possibly--tangent) + else if (distance == c->radius[0]) + { + // Only intersects if the parameter is on the line segment! + if ((t >= 0.0) && (t <= 1.0)) + { + Global::intersectPoint[0] = c->p[0] + radial; + Global::numIntersectPoints = 1; + } - return 1; + return; } - // Find the distance from the center of c1 to the perpendicular chord - // (which contains the points of intersection) - double x = ((d * d) - (c2->radius * c2->radius) + (c1->radius * c1->radius)) - / (2.0 * d); - // Find the the length of the perpendicular chord -// Not needed...! - double a = sqrt((-d + c2->radius - c1->radius) * (-d - c2->radius + c1->radius) * (-d + c2->radius + c1->radius) * (d + c2->radius + c1->radius)) / d; + // Case #3: Two intersection points (possibly--secant) - // Now, you can use pythagorus to find the length of the hypotenuse, but we - // already know that length, it's the radius! :-P - // What's needed is the angle of the center line and the radial line. Since - // there's two intersection points, there's also four angles (two for each - // circle)! - // We can use the arccos to find the angle using just the radius and the - // distance to the perpendicular chord...! - double angleC1 = acos(x / c1->radius); - double angleC2 = acos((d - x) / c2->radius); + // So, we have the line, and the perpendicular from the center of the + // circle to the line. Now figure out where the intersection points are. + // This is a right triangle, though do we really know all the sides? + // Don't need to, 2 is enough for Pythagoras :-) + // Radius is the hypotenuse, so we have to use c² = a² + b² => a² = c² - b² + double perpendicularLength = sqrt((c->radius[0] * c->radius[0]) - (distance * distance)); - if (tp) - *tp = centerLine.Angle() - angleC1; + // Now, find the intersection points using the length... + Vector lineUnit = Vector(l->p[0], l->p[1]).Unit(); + Point i1 = p + (lineUnit * perpendicularLength); + Point i2 = p - (lineUnit * perpendicularLength); - if (up) - *up = (centerLine.Angle() + PI) - angleC2; + // Next we need to see if they are on the line segment... + double u = ParameterOfLineAndPoint(l->p[0], l->p[1], i1); + double v = ParameterOfLineAndPoint(l->p[0], l->p[1], i2); - if (vp) - *vp = centerLine.Angle() + angleC1; - - if (wp) - *wp = (centerLine.Angle() + PI) + angleC2; - - if (p1) - *p1 = c1->position + (centerLine.Unit() * x) + (Vector::Normal(Vector(), centerLine) * (a / 2.0)); - - if (p2) - *p2 = c1->position + (centerLine.Unit() * x) - (Vector::Normal(Vector(), centerLine) * (a / 2.0)); + if ((u >= 0.0) && (u <= 1.0)) + { + Global::intersectPoint[Global::numIntersectPoints] = i1; + Global::numIntersectPoints++; + } - return 2; + if ((v >= 0.0) && (v <= 1.0)) + { + Global::intersectPoint[Global::numIntersectPoints] = i2; + Global::numIntersectPoints++; + } } -#endif + // should we just do common trig solves, like AAS, ASA, SAS, SSA? // Law of Cosines: -// c^2 = a^2 + b^2 -2ab*cos(C) +// c² = a² + b² - 2ab * cos(C) // Solving for C: -// cos(C) = (c^2 - a^2 - b^2) / -2ab = (a^2 + b^2 - c^2) / 2ab +// cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab // Law of Sines: // a / sin A = b / sin B = c / sin C @@ -459,3 +363,22 @@ Point Geometry::GetPointForParameter(Object * obj, double t) return Point(0, 0); } + +Point Geometry::Midpoint(Line * line) +{ + return Point((line->p[0].x + line->p[1].x) / 2.0, + (line->p[0].y + line->p[1].y) / 2.0); +} + + +/* +How to find the tangent of a point off a circle: + + • Calculate the midpoint on the point and the center of the circle + • Get the length of the line segment from the and the center divided by two + • Use that length to construct a circle with the point at the center and the + radius equal to that length + • The intersection of the two circles are the tangent points + +*/ +