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

            }

      }

}