using System;
using System.Drawing;
namespace Geometry
{
/// <summary>
///
Geometric solutions for circles
/// </summary>
public
class CircleGeometry
{
private
CircleGeometry()
{
}
/// <summary>
/// Find circles tangent to a line and two circles
/// </summary>
/// <param
name="pt1">First point on line</param>
/// <param
name="pt2">Second point on line</param>
/// <param
name="ptCentre1">Centre of first
circle</param>
/// <param
name="nRad2">Radius of first
circle</param>
/// <param
name="ptCentre2">Centre of second
circle</param>
/// <param
name="nRad3">Radius of second
circle</param>
/// <param
name="Centres">Array to receive list of centre points</param>
/// <param
name="Radii">Array to receive
list of radii</param>
/// <returns>Number of solutions</returns>
public
static int CircleL2C(Point pt1, Point pt2, Point
ptCentre1, int nRad2,
Point ptCentre2, int nRad3, ref Point[] Centres, ref int[] Radii)
//***************************************************************************
// Find tangents to a
line and two circles
{
// Initialise
int nSolutions
= 0;
Centres = new Point[0];
Radii = new int[0];
// Translate
so first circle is on origin
PointF pt1N = new PointF(pt1.X - ptCentre1.X, pt1.Y - ptCentre1.Y);
PointF pt2N = new PointF(pt2.X - ptCentre1.X, pt2.Y - ptCentre1.Y);
PointF ptCentre2N = new PointF(ptCentre2.X - ptCentre1.X, ptCentre2.Y - ptCentre1.Y);
// Find angle to rotate
second circle to x-axis
double
ang = Math.Atan2(-ptCentre2N.Y, ptCentre2N.X);
double
angCos = Math.Cos(ang);
double
angSin = Math.Sin(ang);
// Rotate
circle centre
double
x3 = (ptCentre2N.X * angCos) - (ptCentre2N.Y * angSin);
// Rotate first point on
line
double
iTemp = (pt1N.X * angCos) -
(pt1N.Y * angSin);
pt1N.Y = (float)((pt1N.X * angSin) + (pt1N.Y * angCos));
pt1N.X = (float)iTemp;
// Rotate second point on
line
iTemp
= (pt2N.X * angCos) - (pt2N.Y * angSin);
pt2N.Y = (float)((pt2N.X * angSin) + (pt2N.Y * angCos));
pt2N.X = (float)iTemp;
// Derive
solutions
PointF[] fCentres = new PointF[0];
double[]
fRadii = new double[0];
int nfSolutions
= NormCircleL2C(pt1N, pt2N, nRad2, x3, nRad3, ref
fCentres, ref fRadii);
// Transform solutions
back to original coordinate system
double
xc, yc;
Point ptCentre = new Point();
double
angNCos = Math.Cos(-ang);
double
angNSin = Math.Sin(-ang);
for(int iSolution=0;
iSolution<nfSolutions;
++iSolution)
{
xc
= fCentres[iSolution].X;
yc
= fCentres[iSolution].Y;
ptCentre.X = (int)Math.Round(((xc * angNCos - yc * angNSin) + ptCentre1.X));
ptCentre.Y = (int)Math.Round(((xc * angNSin + yc * angNCos) + ptCentre1.Y));
nSolutions
= AddSolution(ptCentre, (int)Math.Round(fRadii[iSolution]), ref Centres, ref Radii);
}
// Return number of
solutions
return
nSolutions;
}
private
static int NormCircleL2C(PointF
pt1, PointF pt2, int nRad2, double
x3, int nRad3,
ref PointF[]
fCentres, ref double[] fRadii)
//************************************************************************
// Find circles tangent
to a line and 2 circles. The first circle is centred
on the origin,
// the second at (x3, 0)
{
double
a1, b1, c1, t, r2, r3, a, b, c, u, s;
double
A, B, C, b24ac, b24acRoot, fRadius, xc, yc, Ac, Bc,
Cc;
// Initialise
int nSolutions
= 0;
fCentres
= new PointF[0];
fRadii
= new double[0];
// Get line equation
a1 = 0.0; b1
= 0.0; c1 = 0.0;
if(!LineGeometry.Equation(pt1, pt2, ref
a1, ref b1, ref
c1)) return 0;
// Take all combinations
of t to get circles on other side of line,
// and radii to get
circles on other side of given circles
for(int iCase=0;
iCase<8; ++iCase)
{
t = ((iCase &
1) == 0) ? 1 : -1;
r2 = ((iCase &
2) == 0) ? nRad2 : -nRad2;
r3 = ((iCase &
4) == 0) ? nRad3 : -nRad3;
// Get constants
a = 2*(a1*(r2 - r3) - x3*t);
b = 2*b1*(r2 - r3);
c = 2*c1*(r2 - r3) + t*(r2*r2 - r3*r3 +
x3*x3);
if(b
!= 0.0)
{
u = b1*c - b*c1;
s = a1*b - a*b1;
A = t*t*b*b*(a*a + b*b) - b*b*s*s;
B = 2*(u*t*b*(a*a + b*b) + a*c*s*t*b -
b*b*s*s*r2);
C = u*u*(a*a + b*b) + 2*a*c*s*u + c*c*s*s
- b*b*s*s*r2*r2;
}
else
{
u = a1*c - a*c1;
s = a*b1;
A = a*a*(t*t*a*a - s*s);
B = 2*a*a*(u*t*a - s*s*r2);
C = u*u*a*a + c*c*s*s - a*a*s*s*r2*r2;
}
// Calculate radius
b24ac = (B * B) - (4 * A * C);
if((b24ac < 0.0) || (A ==
0.0)) continue;
b24acRoot = Math.Sqrt(b24ac);
for(int iRoot=0;
iRoot<2; ++iRoot)
{
if((b24ac == 0.0) && (iRoot > 0)) continue;
if(iRoot == 0) fRadius = (-B + b24acRoot) / (A + A);
else fRadius = (-B -
b24acRoot) / (A + A);
if((fRadius <=
0.0) || (fRadius > (double)Int32.MaxValue))
continue;
// Derive x-coordinate of centre
double[] xSols = new double[2];
int nxSols =
0;
if(x3 != 0.0)
{
xc = ((r2+fRadius)*(r2+fRadius) -
(r3+fRadius)*(r3+fRadius) + x3*x3) / (2 * x3);
if(Math.Abs(xc) > (double)Int32.MaxValue)
continue;
xSols[nxSols++]
= xc;
}
else // If circles
are concentric there are 2 solutions for x
{
Ac = (a1*a1 + b1*b1);
Bc = -2*a1*(fRadius*t -
c1);
Cc = (fRadius*t
- c1)*(fRadius*t - c1) - b1*b1*(r2 + fRadius)*(r2 + fRadius);
b24ac = Bc*Bc
- 4*Ac*Cc;
if((b24ac < 0.0) || (Ac == 0.0)) continue;
b24acRoot = Math.Sqrt(b24ac);
for(int xCase=0; xCase<2;
++xCase)
{
//
Calculate x-coordinate
if(xCase == 0) xc = (-Bc + b24acRoot) / (Ac +
Ac);
else xc = (-Bc - b24acRoot) / (Ac + Ac);
if(Math.Abs(xc) > (double)Int32.MaxValue)
continue;
xSols[nxSols++]
= xc;
}
}
// Now derive
y-coordinate(s)
for(int ixSol=0; ixSol<nxSols; ++ixSol)
{
xc = xSols[ixSol];
if(b1 != 0.0)
{
yc = (-a1*xc - c1 + fRadius*t) / b1;
}
else
{
double ycSquare
= (r2+fRadius)*(r2+fRadius)- (xc*xc);
if(ycSquare <
0.0) continue;
yc = Math.Sqrt(ycSquare);
}
if(Math.Abs(yc) > (double)Int32.MaxValue)
continue;
// Add solution
PointF ptCentre = new PointF((float)xc, (float)yc);
nSolutions = AddSolution(ptCentre, fRadius, ref fCentres, ref fRadii);
if(b1 == 0.0)
{
ptCentre = new PointF((float)xc, (float)-yc);
nSolutions = AddSolution(ptCentre, fRadius, ref fCentres, ref fRadii);
}
}
}
}
// Return number of
solutions found
return
nSolutions;
}
private
static int AddSolution(Point
ptCentre, int
nRadius,
ref
Point[] Centres, ref int[] Radii)
//*****************************************************
// Add
a solution to a list if not already present
// In
this overload, all values are integral
{
Point[] CentreCopy = new Point[0];
int[] RadiiCopy = new int[0];
// Don't
accept zero or negative radius solutions
int nSolutions
= Centres.Length;
if(nRadius <= 0) return nSolutions;
// Check if present
bool bPresent
= false;
for(int iSolution=0;
iSolution<nSolutions; ++iSolution)
{
if((ptCentre == Centres[iSolution]) &&
(nRadius == Radii[iSolution]))
{
bPresent
= true;
break;
}
}
// If found, keep current
solutions
if(bPresent) return nSolutions;
// Otherwise,
copy existing solutions to temporary arrays
if(nSolutions > 0)
{
CentreCopy = new Point[nSolutions];
RadiiCopy = new int[nSolutions];
Array.Copy(Centres, 0, CentreCopy,
0, nSolutions);
Array.Copy(Radii, 0, RadiiCopy, 0, nSolutions);
}
// Set up new solution
arrays and load
Centres = new Point[nSolutions + 1];
Radii = new int[nSolutions + 1];
if(nSolutions > 0)
{
Array.Copy(CentreCopy, 0, Centres, 0, nSolutions);
Array.Copy(RadiiCopy, 0, Radii, 0, nSolutions);
}
Centres[nSolutions] = ptCentre;
Radii[nSolutions]
= nRadius;
// Return updated number of solutions
return
(nSolutions + 1);
}
private
static int AddSolution(PointF ptCentre, double fRadius,
ref
PointF[] Centres, ref double[] Radii)
//*****************************************************
// Add
a solution to a list if not already present
// In
this overload, all values are floating-point
{
PointF[] CentreCopy = new PointF[0];
double[]
RadiiCopy = new double[0];
// Don't
accept zero or negative radius solutions
int nSolutions
= Centres.Length;
if(fRadius <= 0.0) return nSolutions;
// Check if present
bool bPresent
= false;
for(int iSolution=0;
iSolution<nSolutions; ++iSolution)
{
if((ptCentre == Centres[iSolution]) &&
(fRadius
== Radii[iSolution]))
{
bPresent = true;
break;
}
}
// If found, keep current
solutions
if(bPresent) return nSolutions;
// Otherwise,
copy existing solutions to temporary arrays
if(nSolutions > 0)
{
CentreCopy = new PointF[nSolutions];
RadiiCopy = new double[nSolutions];
Array.Copy(Centres, 0, CentreCopy,
0, nSolutions);
Array.Copy(Radii, 0, RadiiCopy, 0, nSolutions);
}
// Set up new solution
arrays and load
Centres = new PointF[nSolutions + 1];
Radii = new double[nSolutions
+ 1];
if(nSolutions > 0)
{
Array.Copy(CentreCopy, 0, Centres, 0, nSolutions);
Array.Copy(RadiiCopy, 0, Radii, 0, nSolutions);
}
Centres[nSolutions] = ptCentre;
Radii[nSolutions]
= fRadius;
// Return updated number
of solutions
return
(nSolutions + 1);
}
}
}
using System;
using System.Drawing;
namespace Geometry
{
/// <summary>
///
Geometric solutions for straight lines
/// </summary>
public
class LineGeometry
{
private
LineGeometry()
{
}
/// <summary>
/// Get the normalised equation of a
line between two points in the form: Ax + By + C = 0
/// </summary>
/// <param
name="pt1">First point</param>
/// <param
name="pt2">Second point</param>
/// <param
name="A">X factor</param>
/// <param
name="B">Y factor</param>
/// <param
name="C">Constant</param>
/// <returns>true on success</returns>
public
static bool Equation(PointF
pt1, PointF pt2, ref double A, ref double B, ref double C)
//****************************************************************************************
// Derive
the equation of the line between 2 points in the form Ax + By + C = 0
// If the points are coincident, return
false
// In
this overload, coordinates are floating-point values
{
// Validation
if((pt1
== pt2)) return false;
// Get constants
double
xDiff = pt2.X - pt1.X;
double
yDiff = pt2.Y - pt1.Y;
double
rSquare = (xDiff * xDiff) + (yDiff * yDiff);
double
rInv = 1.0 / rSquare;
// Derive
parameters
A = -yDiff * rInv;
B = xDiff * rInv;
C = ((double)pt1.X * pt2.Y - (double)pt2.X
* pt1.Y) * rInv;
// Normalize
the equation for convenience
double
sMult = 1.0 / Math.Sqrt( A
* A + B * B);
A *= sMult;
B *= sMult;
C *= sMult;
// Return success
return
true;
}
}
}