using System;

using System.Drawing;

namespace Geometry

{

/// <summary>

/// Geometric solutions for circles

/// </summary>

public class CircleGeometry

{

private CircleGeometry()

{

}

/// <summary>

/// Find circles tangent to three other circles

/// </summary>

/// <param name="ptCentre1">Centre of first circle</param>

/// <param name="nRadius1">Radius of first circle</param>

/// <param name="ptCentre2">Centre of second circle</param>

/// <param name="nRadius2">Radius of second circle</param>

/// <param name="ptCentre3">Centre of third circle</param>

/// <param name="nRadius3">Radius of third 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 Circle3C(Point ptCentre1, int nRadius1, Point ptCentre2, int nRadius2,

Point ptCentre3, int nRadius3, ref Point[] Centres, ref int[] Radii)

//*************************************************************************************

// Find tangents to three circles

{

// Initialise

int nSolutions = 0;

Centres = new Point[0];

Radii = new int[0];

// If all circles concentric, there are no solutions

if((ptCentre1 == ptCentre2) && (ptCentre2 == ptCentre3)) return 0;

// Ensure first 2 circles are not concentric

if(ptCentre1 == ptCentre2)

{

Point ptTemp = ptCentre2;

ptCentre2 = ptCentre3;

ptCentre3 = ptTemp;

int nTemp = nRadius2;

nRadius2 = nRadius3;

nRadius3 = nTemp;

}

// Translate so first circle is on origin

Point ptCentre2N = new Point();

ptCentre2N.X = ptCentre2.X - ptCentre1.X;

ptCentre2N.Y = ptCentre2.Y - ptCentre1.Y;

Point ptCentre3N = new Point();

ptCentre3N.X = ptCentre3.X - ptCentre1.X;

ptCentre3N.Y = ptCentre3.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 second and third circles

double x2 = ptCentre2N.X * angCos - ptCentre2N.Y * angSin;

double x3 = ptCentre3N.X * angCos - ptCentre3N.Y * angSin;

double y3 = ptCentre3N.X * angSin + ptCentre3N.Y * angCos;

// Derive solutions

PointF[] fCentres = new PointF[0];

double[] fRadii = new double[0];

int nfSolutions = NormCircle3C(nRadius1, x2, nRadius2, x3, y3, nRadius3, 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 NormCircle3C(int nRadius1, double x2, int nRadius2, double x3, double y3,

int nRadius3, ref PointF[] fCentres, ref double[] fRadii)

//********************************************************************************

// Find circles tangent to 3 circles. The first circle is centered at the origin, the

// second at (x2, 0) and the third at (x3, y3)

// NB x2 must not be zero. y3 must not be zero.

{

double r1, r2, r3, a, b, c, t, A, B, C, b24ac, b24acRoot;

double fRadius, xc, yc;

// Initialise

int nSolutions = 0;

fCentres = new PointF[0];

fRadii = new double[0];

// Negate the radii to get all combinations

for(int iCase=0; iCase<8; ++iCase)

{

r1 = ((iCase & 1) == 0) ? nRadius1 : -nRadius1;

r2 = ((iCase & 2) == 0) ? nRadius2 : -nRadius2;

r3 = ((iCase & 4) == 0) ? nRadius3 : -nRadius3;

// Special case where radii of first 2 circles are equal

if(r1 == r2)

{

CoRadialSol(r1, x2, r2, x3, y3, r3, ref fCentres, ref fRadii);

continue;

}

// Get constants

a = 2*(x2*(r3 - r1) - x3*(r2 - r1));

b = 2*y3*(r1 - r2);

c = (r2 - r1)*(x3*x3 + y3*y3 -

(r3 - r1)*(r3 - r1)) - (r3 - r1)*(x2*x2 - (r2 - r1)*(r2 - r1));

t = (x2*x2 + r1*r1 - r2*r2)/2.0;

A = (r1 - r2)*(r1 - r2)*(a*a + b*b) - (x2*x2*b*b);

B = 2*(t*(r1 - r2)*(a*a + b*b) + a*c*x2*(r1 - r2) - (r1*x2*x2*b*b));

C = t*t*(a*a + b*b) + (2*a*c*x2*t) + (c*c*x2*x2) - (r1*r1*x2*x2*b*b);

// Calculate radius

b24ac = B*B - 4*A*C;

if((b24ac > 0.0) && (A != 0.0))

{

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 (x2 may not be zero)

xc = (fRadius*(r1 - r2) + t) / x2;

if(Math.Abs(xc) > (double)Int32.MaxValue) continue;

// Derive y-coordinate of centre. b should never be 0, as

// r1=r2 is special case and y3 may not be zero

yc = (-a*xc - c) / b;

if(Math.Abs(yc) > (double)Int32.MaxValue) continue;

// Load results

PointF ptCentre = new PointF((float)xc, (float)yc);

nSolutions = AddSolution(ptCentre, fRadius, ref fCentres, ref fRadii);

}

}

}

// Return number of valid results

return fCentres.Length;

}

private static void CoRadialSol(double r1, double x2, double r2, double x3, double y3,

double r3, ref PointF[] fCentres, ref double[] fRadii)

//*******************************************************************************

// Find circle tangent to 3 circles, where the first 2 circles have the same radius

// The first circle is at (0,0), radius r1, the second at (x2, 0)

// and the third at (x3,y3), radius r3

{

double b24ac, b24acRoot, yc, fRadius, A, B, C;

PointF ptCentre;

// Calculate x-cordinate of centre

double xc = x2 / 2.0;

// If all radii are equal, there will be only one solution

if(r1 == r3)

{

if(y3 == 0.0) return;

// Calculate y-coordinate of centre

yc = (x3*x3 - 2*xc*x3 + y3*y3) / (y3 + y3);

if(Math.Abs(yc) > (double)Int32.MaxValue) return;

// Derive radius

A = 1;

B = 2 * r1;

C = r1*r1 - xc*xc - yc*yc;

b24ac = B*B - 4*A*C;

if(b24ac < 0.0) return;

b24acRoot = Math.Sqrt(b24ac);

fRadius = (-B + b24acRoot) / (A + A);

if((fRadius <= 0.0) || (Math.Abs(fRadius) > (double)Int32.MaxValue))

{

fRadius = (-B - b24acRoot) / (A + A);

if((fRadius <= 0.0) || (Math.Abs(fRadius) > (double)Int32.MaxValue)) return;

}

// Add solution

ptCentre = new PointF((float)xc, (float)yc);

AddSolution(ptCentre, fRadius, ref fCentres, ref fRadii);

}

else

{

// Evaluate constants

double k = r1*r1 - r3*r3 + x3*x3 + y3*y3 - 2*xc*x3;

A = 4*((r1 - r3)*(r1 - r3) - y3*y3);

B = 4*(k*(r1 - r3) - 2*y3*y3*r1);

C = 4*xc*xc*y3*y3 + k*k - 4*y3*y3*r1*r1;

b24ac = B*B - 4*A*C;

if((b24ac >= 0.0) && (A != 0.0))

{

// Calculate radii

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;

// Evaluate y-coordinate

yc = (2 * fRadius * (r1 - r3) + k)/(2 * y3);

// Add solution

ptCentre = new PointF((float)xc, (float)yc);

AddSolution(ptCentre, fRadius, ref fCentres, ref fRadii);

}

}

}

}

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

}

}

}