/* Computational Geometry for Circles */
#include "geom-circle.h"
#include <math.h>
/****************************************************************
* Calculates the intersection points of two overlapping circles.
* Returns true if the circles intersect.
* Returns false if the circles do not intersect.
****************************************************************/
int getCircleCircleIntersections(/* Circle 0 Radius */
double r0,
/* Circle 0 Center */
double px0, double py0,
/* Circle 1 Radius */
double r1,
/* Circle 1 Center */
double px1, double py1,
/* Intersection Point */
double* px3, double* py3,
/* Intersection Point */
double* px4, double* py4)
{
double a, h, px2, py2, mx, my;
double dx = px1-px0;
double dy = py1-py0;
double d = sqrt(dx*dx + dy*dy); /* Distance between centers */
/*Circles share centers. This results in division by zero,
infinite solutions or one circle being contained within the other. */
if(d == 0.0)
return 0;
/* Circles do not touch each other */
else if(d > (r0 + r1))
return 0;
/* One circle is contained within the other */
else if(d < (r0 - r1))
return 0;
/*
* Considering the two right triangles p0p2p3 and p1p2p3 we can write:
* a^2 + h^2 = r0^2 and b^2 + h^2 = r1^2
*
* BEGIN PROOF
*
* Remove h^2 from the equation by setting them equal to themselves:
* r0^2 - a^2 = r1^2 - b^2
* Substitute b with (d - a) since it is proven that d = a + b:
* r0^2 - a^2 = r1^2 - (d - a)^2
* Complete the square:
* r0^2 - a^2 = r1^2 - (d^2 -2da + a^2)
* Subtract r1^2 from both sides:
* r0^2 - r1^2 - a^2 = -(d^2 -2da + a^2)
* Invert the signs:
* r0^2 - r1^2 - a^2 = -d^2 + 2da - a^2
* Adding a^2 to each side cancels them out:
* r0^2 - r1^2 = -d^2 + 2da
* Add d^2 to both sides to shift it to the other side:
* r0^2 - r1^2 + d^2 = 2da
* Divide by 2d to finally solve for a:
* a = (r0^2 - r1^2 + d^2)/ (2d)
*
* END PROOF
*/
a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d);
/* Solve for h by substituting a into a^2 + h^2 = r0^2 */
h = sqrt((r0*r0) - (a*a));
/*Find point p2 by adding the a offset in relation to line d to point p0 */
px2 = px0 + (dx * a/d);
py2 = py0 + (dy * a/d);
/* Tangent circles have only one intersection */
if(d == (r0 + r1))
{
*px3 = *py4 = px2;
*py3 = *py4 = py2;
return 1;
}
/* Get the perpendicular slope by multiplying by the negative reciprocal
* Then multiply by the h offset in relation to d to get the actual offsets */
mx = -(dy * h/d);
my = (dx * h/d);
/* Add the offsets to point p2 to obtain the intersection points */
*px3 = px2 + mx;
*py3 = py2 + my;
*px4 = px2 - mx;
*py4 = py2 - my;
return 1;
}
/****************************************************************
* Calculates the tangent points on circle from a given point.
* Returns true if the given point lies outside the circle.
* Returns false if the given point is inside the circle.
****************************************************************/
int getCircleTangentPoints(/* Circle Radius and Center */
double cr, double cx, double cy,
/* Point to determine tangency */
double px, double py,
/* Tangent Point 0 */
double* tx0, double* ty0,
/* Tangent Point 1 */
double* tx1, double* ty1)
{
double pr;
double dx = px-cx;
double dy = py-cy;
double hyp = sqrt(dx*dx + dy*dy); /* Distance to center of circle */
/* Point is inside the circle */
if(hyp < cr)
return 0;
/* Point is lies on the circle, so there is only one tangent point */
else if(hyp == cr)
{
*tx0 = *tx1 = px;
*ty0 = *ty1 = py;
return 1;
}
/* Since the tangent lines are always perpendicular to the radius, so
* we can use the Pythagorean theorem to solve for the missing side */
pr = sqrt((hyp*hyp) - (cr*cr));
return getCircleCircleIntersections(cr,
cx, cy,
pr,
px, py,
tx0, ty0,
tx1, ty1);
}
/* NOTE: Uncomment the #define below to compile test
* gcc -std=c89 -o geom-circle geom-circle.c -lm
*/
/* #define TEST_GEOM_CIRCLE */
#ifdef TEST_GEOM_CIRCLE
#include <stdio.h>
void testTangentPoints(double cr, double cx, double cy,
double px, double py)
{
double tx0, ty0, tx1, ty1;
if(!getCircleTangentPoints(cr, cx, cy,
px, py,
&tx0, &ty0,
&tx1, &ty1))
printf("Error calculating tangent points.\n");
else
printf("Circle : cr=%f, cx=%f, cy=%f\n"
"Point : px=%f, py=%f\n"
"Tangent: tx0=%f, ty0=%f\n"
"Tangent: tx1=%f, ty1=%f\n\n",
cr, cx, cy,
px, py,
tx0, ty0,
tx1, ty1);
}
int main(void)
{
testTangentPoints(3.0, 0.0, 0.0, 4.0, 0.0);
/* Solution: tx0 = 2.2500, ty0 = 1.9843 */
/* Solution: tx1 = 2.2500, ty1 = -1.9843 */
testTangentPoints(6.8221, 20.1762, 10.7170, 24.3411, 18.2980);
/* Solution: tx0 = 19.0911, ty0 = 17.4522 */
/* Solution: tx1 = 26.4428, ty1 = 13.4133 */
return 0;
}
#endif