using MathNet.Numerics.LinearAlgebra;
using MathNetMatrix = MathNet.Numerics.LinearAlgebra.Matrix<double>;
using MathNetVector = MathNet.Numerics.LinearAlgebra.Vector<double>;
using System.Drawing;
namespace 头罩视野.Services
{
    class GetArea
    {

        // 设备固定参数
        private static double R = 330;         // 半球半径
        private static double angleStep = 10;  // 每格角度

        // 定义参数（和你代码里一致）
        private const int totalLights = 81;

        // 传入：72个灯的亮灭数据（0=灭，1=亮）
        // 返回：椭圆面积
        public static double CalculateEllipseArea(int[] lightData, List<(int m, int n)> lightPositions)
        {
            //if (lightData.Length != totalLights || lightPositions.Count != totalLights)
            //    throw new Exception("必须是81个灯的数据");

            // 第一步：收集所有亮灯坐标
            List<System.Drawing.Point> brightPoints = new List<System.Drawing.Point>();
            for (int i = 0; i < totalLights; i++)
            {
                if (lightData[i] == 1)
                {
                    var (m, n) = lightPositions[i];
                    System.Drawing.Point p = GetLightPoint(m, n);
                    brightPoints.Add(p);
                }
            }

            // 第二步：用亮点拟合椭圆
            var (cx, cy, a, b, area) = FitEllipse(brightPoints);

            // 返回面积
            return area;
        } 
        /// 生成设备全部243盏灯的(m,n)位置   上爪1条、下爪1条、左右共用1条，各81灯
        private static System.Drawing.Point GetLightPoint(int m, int n)
        {
            double radH, radV;

            // 上爪灯条（n=0）：水平角固定，m控制垂直角
            if (n == 0)
            {
                radH = 0;
                radV = m * angleStep * Math.PI / 180;
            }
            // 下爪灯条（n=1）：水平角固定为180°，m控制垂直角
            else if (n == 1)
            {
                radH = Math.PI;
                radV = m * angleStep * Math.PI / 180;
            }
            // 左右共用灯条（n=2）：垂直角固定，m控制水平角
            else
            {
                radH = m * angleStep * Math.PI / 180;
                radV = 0;
            }

            // 保留你原来的投影公式
            double x = R * Math.Tan(radH);
            double y = R * Math.Tan(radV);

            return new System.Drawing.Point((int)Math.Round(x), (int)Math.Round(y));
        }

        // 最小二乘法拟合椭圆（核心算法）cx：椭圆中心点的 X 坐标 cy：椭圆中心点的 Y 坐标  a：椭圆的长半轴长度（较大的那个半径）b：椭圆的短半轴长度（较小的那个半径）


        private static (double cx, double cy, double a, double b, double area) FitEllipse(List<Point> points)
            {
                int n = points.Count;
                if (n < 5)
                    throw new Exception("至少需要5个点来拟合椭圆");

            // 这里是正确写法
            var M = MathNetMatrix.Build.Dense(n, 5);
            var Y = MathNetVector.Build.Dense(n, i => -1.0);

            for (int i = 0; i < n; i++)
                {
                    double x = points[i].X;
                    double y = points[i].Y;
                    M[i, 0] = x * x;
                    M[i, 1] = x * y;
                    M[i, 2] = y * y;
                    M[i, 3] = x;
                    M[i, 4] = y;
                }

                // 求解
                Vector<double> sol = M.QR().Solve(Y);
                double A = sol[0], B = sol[1], C = sol[2], D = sol[3], E = sol[4], F = 1;

                // 椭圆中心
                double cx = (2 * C * D - B * E) / (B * B - 4 * A * C);
                double cy = (2 * A * E - B * D) / (B * B - 4 * A * C);

                // 半轴
                double term1 = 2 * (A * E * E + C * D * D - B * D * E + (B * B - 4 * A * C) * F);
                double term2 = (A + C) + Math.Sqrt((A - C) * (A - C) + B * B);
                double term3 = (A + C) - Math.Sqrt((A - C) * (A - C) + B * B);

                double a = Math.Sqrt(Math.Abs(term1 / ((B * B - 4 * A * C) * term3)));
                double b = Math.Sqrt(Math.Abs(term1 / ((B * B - 4 * A * C) * term2)));

                if (a < b) (a, b) = (b, a);
                double area = Math.PI * a * b;

            return (cx, cy, a, b, area);
        }



  




    }


}