X-Git-Url: http://shamusworld.gotdns.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fgeometry.cpp;h=b9158974ec97d902d4ff2e936a2cc4b43e469082;hb=0fcc2d879e1e0ca17eeaceae2159f5143a06586f;hp=6ed6fab1020c70449e64955f4bc9c6c0548f28b8;hpb=6a7baa2814a8b4d0b93df776a4c99689bcfb3ffa;p=architektonas diff --git a/src/geometry.cpp b/src/geometry.cpp index 6ed6fab..b915897 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -20,32 +20,6 @@ #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). @@ -105,7 +79,7 @@ double Geometry::Determinant(Point p1, Point p2) } -void 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) { Global::numIntersectPoints = Global::numIntersectParams = 0; @@ -113,6 +87,10 @@ void Geometry::Intersects(Object * obj1, Object * obj2)//, double * tp/*= 0*/, d 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); } @@ -136,17 +114,14 @@ So check if the above two numbers are both >=0 and <=1. // Finds the intersection between two lines (if any) -void 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; @@ -154,29 +129,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; -#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; @@ -188,47 +144,19 @@ printf(" -->(q - p) . s = %lf, (p - q) . r = %lf\n", v1.Dot(s), v2.Dot(r)); t = 1.0, u = 1.0; else return; -#endif } 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 0 - if (tp) - *tp = t; - - if (up) - *up = 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; -#else 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)) Global::numIntersectParams = 1; -#endif } @@ -251,15 +179,27 @@ void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2) // 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])) - || (d == fabs(c1->radius[0] - c2->radius[0]))) + if (d == (c1->radius[0] + c2->radius[0])) { 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; } + else if (d == fabs(c1->radius[0] - c2->radius[0])) + { + 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; + } +/* + 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)); @@ -274,188 +214,72 @@ void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2) } -#if 0 -// Finds the intersection between two lines (if any) -int Geometry::Intersects(Line * l1, Dimension * d1, double * tp/*= 0*/, double * up/*= 0*/) +// +// N.B.: l is the line, c is the circle +// +void Geometry::CheckLineToCircleIntersection(Object * l, Object * c) { - Line l2(d1->position, d1->endpoint); - return Intersects(l1, &l2, tp, up); -} + // Set up global vars + Global::numIntersectPoints = Global::numIntersectParams = 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(); -// 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*/) -{ -#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); + // Step 2: See if we have 0, 1, or 2 intersection points - double discriminant = ((c->radius * c->radius) * (dr * dr)) - (determinant * determinant); + // 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; + } - if (discriminant < 0) - return false; + return; + } - + // Case #3: Two intersection points (possibly--secant) - 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; + // 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 (up) - *up = Vector(c->position, p).Angle(); + // Now, find the points using the length, then check to see if they are on + // the line segment + Vector lineUnit = Vector(l->p[0], l->p[1]).Unit(); + Point i1 = p + (lineUnit * perpendicularLength); + Point i2 = p - (lineUnit * perpendicularLength); - return 1; - } + // Now we have our intersection points, 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); - // 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 ((u >= 0.0) && (u <= 1.0)) { - intersections++; + Global::intersectPoint[Global::numIntersectPoints] = i1; + Global::numIntersectPoints++; } - if ((t2 >= 0) && (t2 <= 1.0)) + if ((v >= 0.0) && (v <= 1.0)) { - intersections++; + Global::intersectPoint[Global::numIntersectPoints] = i2; + Global::numIntersectPoints++; } - - return intersections; -#endif } -// 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*/) -{ - // 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); - - // 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; - - // 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(); - - if (up) - *up = centerLine.Angle() + PI; - - return 1; - } - - // Find the distance from the center of c1 to the perpendicular chord - // (which contains the points of intersection) - // [N.B.: This is derived from Pythagorus by using the unknown distance - // from the center line to the point where the two radii coincide as - // a common unknown to two instances of the formula.] - 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; - - // 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); - - if (tp) - *tp = centerLine.Angle() - angleC1; - - if (up) - *up = (centerLine.Angle() + PI) - angleC2; - - 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)); - - return 2; -} -#endif - // should we just do common trig solves, like AAS, ASA, SAS, SSA? // Law of Cosines: // c² = a² + b² - 2ab * cos(C)