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;

            }

      }

}