[architektonas] / src / geometry.cpp

//
// geometry.cpp: Algebraic geometry helper functions
//
// Part of the Architektonas Project
// (C) 2011 Underground Software
// See the README and GPLv3 files for licensing and warranty information
//
// JLH = James Hammons <jlhamm@acm.org>
//
// WHO  WHEN        WHAT
// ---  ----------  ------------------------------------------------------------
// JLH  08/31/2011  Created this file
//
// NOTE: All methods in this class are static.
//

#include "geometry.h"
#include <math.h>
#include <stdio.h>
#include "global.h"
#include "mathconstants.h"

// 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).
double Geometry::ParameterOfLineAndPoint(Point tail, Point head, Point point)
{
        // Geometric interpretation:
        // The parameterized point on the vector lineSegment is where the normal of
        // the lineSegment to the point intersects lineSegment. If the pp < 0, then
        // the perpendicular lies beyond the 1st endpoint. If pp > 1, then the
        // perpendicular lies beyond the 2nd endpoint.

        Vector lineSegment = head - tail;
        double magnitude = lineSegment.Magnitude();
        Vector pointSegment = point - tail;
        double t = lineSegment.Dot(pointSegment) / (magnitude * magnitude);

        return t;
}

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
        // line to the point in question.
        double t = ParameterOfLineAndPoint(tail, head, point);
        Vector v = Vector(tail, head) * t;

        // Get the point normal to point to the line passed in
        Point normalOnLine = tail + v;

        // Make our mirrored vector (head - tail)
        Vector mirror = -(point - normalOnLine);

        // Find the mirrored point
        Point mirroredPoint = normalOnLine + mirror;

        return mirroredPoint;
}

//
// point: The point we're rotating
// rotationPoint: The point we're rotating around
//
Point Geometry::RotatePointAroundPoint(Point point, Point rotationPoint, double angle)
{
        Vector v = Vector(rotationPoint, point);
        double px = (v.x * cos(angle)) - (v.y * sin(angle));
        double py = (v.x * sin(angle)) + (v.y * cos(angle));

        return Vector(rotationPoint.x + px, rotationPoint.y + py, 0);
}

double Geometry::Determinant(Point p1, Point p2)
{
        return (p1.x * p2.y) - (p2.x * p1.y);
}

void Geometry::Intersects(Object * obj1, Object * obj2)
{
        Global::numIntersectPoints = Global::numIntersectParams = 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);
}

/*
Intersecting line segments:
An easier way:
Segment L1 has edges A=(a1,a2), A'=(a1',a2').
Segment L2 has edges B=(b1,b2), B'=(b1',b2').
Segment L1 is the set of points tA'+(1-t)A, where 0<=t<=1.
Segment L2 is the set of points sB'+(1-s)B, where 0<=s<=1.
Segment L1 meet segment L2 if and only if for some t and s we have
tA'+(1-t)A=sB'+(1-s)B
The solution of this with respect to t and s is

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?))

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?))

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)
{
        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

        double rxs = (r.x * s.y) - (s.x * r.y);
        double t, u;

        if (rxs == 0)
        {
                double qpxr = (v1.x * r.y) - (r.x * v1.y);

                // Lines are parallel, so no intersection...
                if (qpxr != 0)
                        return;

                // Check to see which endpoints are connected... Four possibilities:
                if (l1->p[0] == l2->p[0])
                        t = 0, u = 0;
                else if (l1->p[0] == l2->p[1])
                        t = 0, u = 1.0;
                else if (l1->p[1] == l2->p[0])
                        t = 1.0, u = 0;
                else if (l1->p[1] == l2->p[1])
                        t = 1.0, u = 1.0;
                else
                        return;
        }
        else
        {
                t = ((v1.x * s.y) - (s.x * v1.y)) / rxs;
                u = ((v1.x * r.y) - (r.x * v1.y)) / rxs;
        }

        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;
}

void Geometry::CheckCircleToCircleIntersection(Object * c1, Object * c2)
{
        // Set up global vars
        Global::numIntersectPoints = Global::numIntersectParams = 0;

        // 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();

        // 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;

        // 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]))
        {
                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));

        // 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;
}

//
// N.B.: l is the line, c is the circle
//
void Geometry::CheckLineToCircleIntersection(Object * l, Object * c)
{
        // 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();

        // Step 2: See if we have 0, 1, or 2 intersection points

        // 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;
        }

        // Case #3: Two intersection points (possibly--secant)

        // 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));

        // 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);

        // 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 ((u >= 0.0) && (u <= 1.0))
        {
                Global::intersectPoint[Global::numIntersectPoints] = i1;
                Global::numIntersectPoints++;
        }

        if ((v >= 0.0) && (v <= 1.0))
        {
                Global::intersectPoint[Global::numIntersectPoints] = i2;
                Global::numIntersectPoints++;
        }
}

// should we just do common trig solves, like AAS, ASA, SAS, SSA?
// Law of Cosines:
// c² = a² + b² - 2ab * cos(C)
// Solving for C:
// cos(C) = (c² - a² - b²) / -2ab = (a² + b² - c²) / 2ab
// Law of Sines:
// a / sin A = b / sin B = c / sin C

// Solve the angles of the triangle given the sides. Angles returned are
// opposite of the given sides (so a1 consists of sides s2 & s3, and so on).
void Geometry::FindAnglesForSides(double s1, double s2, double s3, double * a1, double * a2, double * a3)
{
        // Use law of cosines to find 1st angle
        double cosine1 = ((s2 * s2) + (s3 * s3) - (s1 * s1)) / (2.0 * s2 * s3);

        // Check for a valid triangle
        if ((cosine1 < -1.0) || (cosine1 > 1.0))
                return;

        double angle1 = acos(cosine1);

        // Use law of sines to find 2nd & 3rd angles
// sin A / a = sin B / b
// sin B = (sin A / a) * b
// B = arcsin( sin A * (b / a))
// ??? ==> B = A * arcsin(b / a)
/*
Well, look here:
sin B = sin A * (b / a)
sin B / sin A = b / a
arcsin( sin B / sin A ) = arcsin( b / a )

hmm... dunno...
*/

        double angle2 = asin(s2 * (sin(angle1) / s1));
        double angle3 = asin(s3 * (sin(angle1) / s1));

        if (a1)
                *a1 = angle1;

        if (a2)
                *a2 = angle2;

        if (a3)
                *a3 = angle3;
}

Point Geometry::GetPointForParameter(Object * obj, double t)
{
        if (obj->type == OTLine)
        {
                // Translate line vector to the origin, then add the scaled vector to
                // initial point of the line.
                Vector v = obj->p[1] - obj->p[0];
                return obj->p[0] + (v * 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 of the point and the center of the circle
 •  Get the length of the line segment from point 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

Another way:

 •  Find angle between radius and line between point and center
 •  Angle +/- from line (point to center) are the tangent points

*/
void Geometry::FindTangents(Object * c, Point p)
{
        // Set up global vars
        Global::numIntersectPoints = Global::numIntersectParams = 0;

        // Get the distance between the point and the center of the circle, the
        // length, and the angle.
        Vector centerLine(c->p[0], p);
        double d = centerLine.Magnitude();
        double clAngle = centerLine.Angle();

        // If the point is on or inside the circle, there are no tangents.
        if (d <= c->radius[0])
                return;

        // Use 'cos(a) = adjacent / hypotenuse' to get the tangent angle
        double a = acos(c->radius[0] / d);

        // Finally, find the points of intersection by using +/- the angle found
        // from the centerline's angle
        Global::intersectPoint[0].x = c->p[0].x + (cos(clAngle + a) * c->radius[0]);
        Global::intersectPoint[0].y = c->p[0].y + (sin(clAngle + a) * c->radius[0]);
        Global::intersectPoint[1].x = c->p[0].x + (cos(clAngle - a) * c->radius[0]);
        Global::intersectPoint[1].y = c->p[0].y + (sin(clAngle - a) * c->radius[0]);
        Global::numIntersectPoints = 2;
}

/*
The extangents can be found by collapsing the circle of smaller radius to a point, and diminishing the larger circle by the radius of the smaller, then treating it like the point-circle case (you have to translate the tangent back out by the length of the radius of the smaller circle once found).

Not sure what the analogous method for finding the intangents are...  :-/
Looks like the intangent can be found by augmenting the larger circle by the smaller radius, and doing the center of the smaller as the point-circle case again, only translating the line back by the radius of the smaller.
*/
void Geometry::FindTangents(Object * c1, Object * c2)
{
        // Set up global vars
        Global::numIntersectPoints = Global::numIntersectParams = 0;

        // Find the larger and smaller of the two:
        Object * cLg = c1, * cSm = c2;

        if (c2->radius[0] > c1->radius[0])
                cLg = c2, cSm = c1;

        // Get the distance between the point and the center of the circle, the
        // length, and the angle.
        Vector centerLine(cLg->p[0], cSm->p[0]);
        double d = centerLine.Magnitude();
        double clAngle = centerLine.Angle();

        // If one circle is completely inside the other, there are no tangents.
        if ((d + cSm->radius[0]) <= cLg->radius[0])
                return;

        // Subtract the radius of the smaller from the larger, and use the point-circle method to find the extangent:
        double a = acos((cLg->radius[0] - cSm->radius[0]) / d);

        // Finally, find the points of intersection by using +/- the angle found
        // from the centerline's angle
        Global::intersectPoint[0].x = cLg->p[0].x + (cos(clAngle + a) * cLg->radius[0]);
        Global::intersectPoint[0].y = cLg->p[0].y + (sin(clAngle + a) * cLg->radius[0]);
        Global::intersectPoint[1].x = cSm->p[0].x + (cos(clAngle + a) * cSm->radius[0]);
        Global::intersectPoint[1].y = cSm->p[0].y + (sin(clAngle + a) * cSm->radius[0]);

        Global::intersectPoint[2].x = cLg->p[0].x + (cos(clAngle - a) * cLg->radius[0]);
        Global::intersectPoint[2].y = cLg->p[0].y + (sin(clAngle - a) * cLg->radius[0]);
        Global::intersectPoint[3].x = cSm->p[0].x + (cos(clAngle - a) * cSm->radius[0]);
        Global::intersectPoint[3].y = cSm->p[0].y + (sin(clAngle - a) * cSm->radius[0]);
        Global::numIntersectPoints = 4;

        // If the circles overlap, there are no intangents.
        if (d <= (cLg->radius[0] + cSm->radius[0]))
                return;

        // Add the radius of the smaller from the larger, and use the point-circle method to find the intangent:
        a = acos((cLg->radius[0] + cSm->radius[0]) / d);

        // Finally, find the points of intersection by using +/- the angle found
        // from the centerline's angle
        Global::intersectPoint[4].x = cLg->p[0].x + (cos(clAngle + a) * cLg->radius[0]);
        Global::intersectPoint[4].y = cLg->p[0].y + (sin(clAngle + a) * cLg->radius[0]);
        Global::intersectPoint[5].x = cSm->p[0].x + (cos(clAngle - a) * cSm->radius[0]);
        Global::intersectPoint[5].y = cSm->p[0].y + (sin(clAngle - a) * cSm->radius[0]);

        Global::intersectPoint[6].x = cLg->p[0].x + (cos(clAngle - a) * cLg->radius[0]);
        Global::intersectPoint[6].y = cLg->p[0].y + (sin(clAngle - a) * cLg->radius[0]);
        Global::intersectPoint[7].x = cSm->p[0].x + (cos(clAngle + a) * cSm->radius[0]);
        Global::intersectPoint[7].y = cSm->p[0].y + (sin(clAngle + a) * cSm->radius[0]);
        Global::numIntersectPoints = 8;
}

//
// Parameter 1: point in question
// Parameter 2, 3: points we are comparing to
//
Point Geometry::NearestTo(Point point, Point p1, Point p2)
{
        double l1 = Vector::Magnitude(point, p1);
        double l2 = Vector::Magnitude(point, p2);

        return (l1 < l2 ? p1 : p2);
}