Алгоритм упаковки окружностей в окружность

using System.Collections.Generic;
 
using System;
using static System.Net.Mime.MediaTypeNames;
using Teigha.Geometry;
using Teigha.Runtime;
using HostMgd.ApplicationServices;
using Teigha.DatabaseServices;
using HostMgd.EditorInput;
 
namespace CirclePacking
 
{
 
    public class Commands
    {
 
        static double _minRad = 0.01;
 
        static short _color = 5;
 
        public class MyCircle
 
        {
 
            public Point2d Center;
 
            public double Radius;
 
            public short ColorIndex;
        }
 
        [CommandMethod("AGC")]
 
        public static void ApollonianGasket()
 
        {
 
        Document doc = HostMgd.ApplicationServices.Application.DocumentManager.MdiActiveDocument;
 
        Database db = doc.Database;
 
            Editor ed = doc.Editor;
 
            PromptEntityOptions peo = new PromptEntityOptions("\nSelect circle: ");
 
            peo.AllowNone = true;
 
            peo.SetRejectMessage("\nMust be a circle.");
 
            peo.AddAllowedClass(typeof(Circle), false);
 
            PromptEntityResult per = ed.GetEntity(peo);
 
            if (per.Status != PromptStatus.OK &&
 
                per.Status != PromptStatus.None
 
              )
 
                return;
 
            Transaction tr = doc.TransactionManager.StartTransaction();
 
            using (tr)
 
            {
 
                BlockTableRecord btr = (BlockTableRecord)tr.GetObject(
 
                    db.CurrentSpaceId, OpenMode.ForWrite
 
                  );
 
                // If the user selected a circle, use it
 
                Circle cir;
                if (per.Status != PromptStatus.None && per.ObjectId != null)
 
                {
 
                    cir =
 
                      tr.GetObject(per.ObjectId, OpenMode.ForRead) as Circle;
 
                }
                else
 
                {
 
                    // Otherwise create a new one at the default location
 
                    cir = new Circle(Point3d.Origin, Vector3d.ZAxis, 500);
 
                    btr.AppendEntity(cir);
 
                    tr.AddNewlyCreatedDBObject(cir, true);
 
                }
 
                List<MyCircle> circles =
 
                  FindApollonianPacking(cir.Radius * 2);
 
                // Offset each of the returned circles to fit our outer one
 
                Vector3d offset = cir.Center - Point3d.Origin - new Vector3d(cir.Radius, cir.Radius, 0);
 
                // Add each of the returned circles to the current space
 
                foreach (MyCircle myc in circles)
                {
                    Point3d cen = new Point3d(myc.Center.X, myc.Center.Y, 0) + offset;
 
                    Circle c = new Circle(cen, Vector3d.ZAxis, myc.Radius);
 
                    c.ColorIndex = myc.ColorIndex;
 
                    btr.AppendEntity(c);
 
                    tr.AddNewlyCreatedDBObject(c, true);
 
                }
 
                tr.Commit();
 
            }
 
        }
 
        // Find the Apollonian packing
 
        private static List<MyCircle> FindApollonianPacking(
 
          double width
 
        )
 
        {
 
            List<MyCircle> circles = new List<MyCircle>();
 
            // Create the three central tangent circles
 
            double radius = width * 0.232f, x = width / 2, gaskHeight = 2 * (radius + 2 * radius / Math.Sqrt(3)),
 
            y = (width - gaskHeight) / 2 + radius;
 
            MyCircle cir0 = CreateCircle(x, y, radius);
 
            circles.Add(cir0);
 
            x -= radius;
 
            y += radius * Math.Sqrt(3);
 
            MyCircle cir1 = CreateCircle(x, y, radius);
 
            circles.Add(cir1);
 
            x += 2 * radius;
 
            MyCircle cir2 = CreateCircle(x, y, radius);
 
            circles.Add(cir2);
 
            // Find the circle that contains them all
 
            MyCircle outer = FindApollonianCircle(cir0, cir1, cir2, -1, -1, -1);
 
            circles.Add(outer);
 
            // Set level to smaller values such as 3 to see partially
 
            // drawn gaskets
 
            int level = 10000;
 
            // Find the central circle
 
            FindCircleOutsideAll(circles, level, cir0, cir1, cir2);
 
            // Find circles tangent to the big circle
 
            FindCircleOutsideTwo(circles, level, cir0, cir1, outer);
 
            FindCircleOutsideTwo(circles, level, cir1, cir2, outer);
 
            FindCircleOutsideTwo(circles, level, cir2, cir0, outer);
 
            return circles;
 
        }
 
        // Draw a circle tangent to these three circles and that is
 
        // outside all three
 
        private static void FindCircleOutsideAll(
 
          List<MyCircle> circles,
 
          int level, MyCircle cir0, MyCircle cir1, MyCircle cir2
 
        )
 
        {
 
            MyCircle newCir = FindApollonianCircle(cir0, cir1, cir2, 1, 1, 1);
 
            if (newCir == null || newCir.Radius < _minRad)
 
                return;
 
            circles.Add(newCir);
 
            if (--level > 0)
 
            {
 
                FindCircleOutsideAll(circles, level, cir0, cir1, newCir);
 
                FindCircleOutsideAll(circles, level, cir0, cir2, newCir);
 
                FindCircleOutsideAll(circles, level, cir1, cir2, newCir);
 
            }
 
        }
 
        // Draw a circle tangent to these three circles and that is
 
        // outside two of them
 
        private static void FindCircleOutsideTwo(
 
          List<MyCircle> circles, int level,
 
          MyCircle cir0, MyCircle cir1, MyCircle contains
 
        )
 
        {
 
            MyCircle newCir =
 
              FindApollonianCircle(cir0, cir1, contains, 1, 1, -1);
 
            if (newCir == null || newCir.Radius < _minRad)
 
                return;
 
            circles.Add(newCir);
 
            if (--level > 0)
 
            {
 
                FindCircleOutsideTwo(circles, level, newCir, cir0, contains);
 
                FindCircleOutsideTwo(circles, level, newCir, cir1, contains);
 
                FindCircleOutsideAll(circles, level, cir0, cir1, newCir);
 
            }
 
        }
 
        private static MyCircle FindApollonianCircle(
 
          MyCircle c1, MyCircle c2, MyCircle c3, int s1, int s2, int s3
 
        )
 
        {
 
            // Make sure c2 doesn't have the same X or Y coordinate as
 
            // the others
 
            double tiny = Tolerance.Global.EqualPoint;
 
            if (
 
              Math.Abs(c2.Center.X - c1.Center.X) < tiny ||
 
              Math.Abs(c2.Center.Y - c1.Center.Y) < tiny
 
            )
 
            {
 
                MyCircle tempCir = c2;
 
                c2 = c3;
 
                c3 = tempCir;
 
                int temp_s = s2;
 
                s2 = s3;
 
                s3 = temp_s;
 
            }
 
            if (
 
              Math.Abs(c2.Center.X - c3.Center.X) < tiny ||
 
              Math.Abs(c2.Center.Y - c3.Center.Y) < tiny
 
            )
 
            {
 
                MyCircle tempCir = c2;
 
                c2 = c1;
 
                c1 = tempCir;
 
                int temp_s = s2;
 
                s2 = s1;
 
                s1 = temp_s;
 
            }
 
            double x1 = c1.Center.X,
 
                   y1 = c1.Center.Y,
 
                   r1 = c1.Radius,
 
                   x2 = c2.Center.X,
 
                   y2 = c2.Center.Y,
 
                   r2 = c2.Radius,
 
                   x3 = c3.Center.X,
 
                   y3 = c3.Center.Y,
 
                   r3 = c3.Radius;
 
            double v11 = 2 * x2 - 2 * x1,
 
                   v12 = 2 * y2 - 2 * y1,
 
                   v13 =
 
                     x1 * x1 - x2 * x2 + y1 * y1 -
 
                     y2 * y2 - r1 * r1 + r2 * r2,
 
                   v14 = 2 * s2 * r2 - 2 * s1 * r1;
 
            double v21 = 2 * x3 - 2 * x2,
 
                   v22 = 2 * y3 - 2 * y2,
 
                   v23 =
 
                     x2 * x2 - x3 * x3 + y2 * y2 - y3 * y3 -
 
                     r2 * r2 + r3 * r3,
 
                   v24 = 2 * s3 * r3 - 2 * s2 * r2;
 
            double w12 = v12 / v11,
 
                   w13 = v13 / v11,
 
                   w14 = v14 / v11;
 
            double w22 = v22 / v21 - w12,
 
                   w23 = v23 / v21 - w13,
 
                   w24 = v24 / v21 - w14;
 
            double P = -w23 / w22,
 
                   Q = w24 / w22,
 
                   M = -w12 * P - w13,
 
                   N = w14 - w12 * Q;
 
            double a = N * N + Q * Q - 1,
 
                   b =
 
                     2 * M * N - 2 * N * x1 + 2 * P * Q -
 
                     2 * Q * y1 + 2 * s1 * r1,
 
                   c =
 
                     x1 * x1 + M * M - 2 * M * x1 + P * P +
 
                     y1 * y1 - 2 * P * y1 - r1 * r1;
 
            // Find roots of a quadratic equation
 
            double[] solutions = QuadraticSolutions(a, b, c);
 
            if (solutions.Length < 1)
 
                return null;
 
            double rs = solutions[0],
 
                   xs = M + N * rs,
 
                   ys = P + Q * rs;
 
            if (
 
              Math.Abs(xs) < tiny ||
 
              Math.Abs(ys) < tiny ||
 
              Math.Abs(rs) < tiny
 
            )
 
                return null;
 
            return CreateCircle(xs, ys, Math.Abs(rs));
 
        }
 
        // Return solutions to a quadratic equation
 
        private static double[] QuadraticSolutions(
 
          double a, double b, double c
 
        )
 
        {
 
            double discriminant = b * b - 4 * a * c;
 
            // See if there are no real solutions
 
            if (discriminant < 0)
 
                return new double[] { };
 
            // See if there is one solution
 
            if (discriminant < Tolerance.Global.EqualPoint)
 
                return new double[] { -b / (2 * a) };
 
            // There are two solutions
 
            return new double[]
 
            {
 
        (-b + Math.Sqrt(discriminant)) / (2 * a),
 
        (-b - Math.Sqrt(discriminant)) / (2 * a),
 
            };
 
        }
 
        private static MyCircle CreateCircle(
        double x, double y, double radius
 
        )
 
        {
 
            MyCircle cir = new MyCircle();
 
            cir.Center = new Point2d(x, y);
 
            cir.Radius = radius;
 
            //cir.ColorIndex = _color;
 
            return cir;
 
        }
 
    }
 
}