GUI functionality fixes.

[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 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

*/