This commit is contained in:
xyy
2026-05-25 13:56:11 +08:00
parent d4edefa050
commit 3c309d6470
3 changed files with 272 additions and 559 deletions

View File

@@ -1,263 +1,111 @@
using MathNet.Numerics.LinearAlgebra;
using MathNetMatrix = MathNet.Numerics.LinearAlgebra.Matrix<double>;
using MathNetVector = MathNet.Numerics.LinearAlgebra.Vector<double>;
using System.Drawing;
using System;
using System.Collections.Generic;
using System.Linq;
namespace .Services
{
class GetArea
public static class GetArea
{
// 设备固定参数
public static double R = 325; // 半球半径
// 定义参数(和你代码里一致)
public const int totalLights = 81;
public static double verticalAngleStep = 90.0 / (totalLights - 1); // 上下灯条用
public static double horizontalAngleStep = 180.0 / (totalLights - 1); // 左右灯条用
//public const double BlankTotalBaseArea = 4610;
//public const double lysmcdSrea = 780;
///// <summary>双目标准标定总面积无面罩空标准头模的双目总实测面积国标总视野保存率计算的基准值单位cm²</summary>
//public const double StandardTotal = 10360;
//空白视野面积计算
public static readonly double _standardTotalArea = (Math.PI * R * R) / 100;
public static double GetBlankViewArea(double binocularTotalArea)
{
// 公式:空白 = 标准总面积 - 双目总视野
return _standardTotalArea - binocularTotalArea;
}
//计算双目视野
public static double CalculateBinocularArea(
List<int> leftFinal,
List<int> rightFinal,
List<(int m, int n)> lightPositions)
{
// 取两者中较小的长度,防止一方数据被截断
int length = Math.Min(leftFinal.Count, rightFinal.Count);
int[] binocularData = new int[length];
for (int i = 0; i < length; i++)
{
binocularData[i] = leftFinal[i] & rightFinal[i];
}
//System.Diagnostics.Debug.WriteLine($"【双目亮灯】长度:{length}, 左眼亮灯:{leftFinal.Count}, 右眼亮灯:{rightFinal.Count}, 双目亮灯:{binocularData.Length}");
return CalculateEllipseArea(binocularData, lightPositions);
}
//下方视野 下方视野角度 = 人眼到「最低亮灯条」的夹角
// 最安全的写法:线程安全 + 不会空引用
private static readonly Random _random = new Random();
// 灯条参数
public const int totalLights = 81; // 每条灯条81个灯
public static double verticalAngleStep = 90.0 / (totalLights - 1); // 1.125°
public static double CalculateBottomViewAngle(int[] lightData, List<(int m, int n)> lightPositions)
{
List<double> bottomAngles = new List<double>();
for (int i = 0; i < lightData.Length; i++)
int bottomLampCount = 0;
for (int i = 0; i < lightData.Length && i < lightPositions.Count; i++)
{
if (lightData[i] == 1)
{
if (lightPositions.Count < lightData.Length)
{
return 0;
}
var (m, n) = lightPositions[i];
if (n == 1)
{
double angle = m * verticalAngleStep;
bottomAngles.Add(angle);
}
}
if (lightPositions[i].n == 1 && lightData[i] == 1)
bottomLampCount++;
}
if (bottomAngles.Count == 0)
return 0;
double baseAngle = bottomAngles.Max() - 13;
// ✅ 绝对不会报空引用
double fluctuation = (_random.NextDouble() * 4) - 2;
double finalAngle = baseAngle + fluctuation;
return finalAngle;
}
//视野保存率
public static double CalcVisionRate(double binocularRate)
{
// 1. 总视野保存率
double ratioTotal = binocularRate / _standardTotalArea;
double gammaTotal = GetVisionGamma(ratioTotal);
double totalRate = gammaTotal * ratioTotal * 100;
return (totalRate);
return (double)bottomLampCount / 81 * 90.0;
}
/// <summary>
/// GB2890-2022 自动获取 总视野/双目视野 校正系数γ
/// 根据径向灯条0/1数据计算边界角度极角
/// </summary>
/// <param name="ratio">实测面积/标准面积 比值(0~1)</param>
/// <returns>校正系数 γ</returns>
public static double GetVisionGamma(double ratio)
/// <param name="lightStates">从上到下的灯条状态索引0=顶部0°索引80=底部90°</param>
public static double ComputeBoundaryAngle(int[] lightStates)
{
// X视野残存率 Si/S0
double[] xData = { 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 };
// 总视野 γ 对应值
double[] gammaTotal = { 1.22, 1.18, 1.14, 1.10, 1.06, 1.03, 1.02, 1.01, 1.00 };
double[] yData = gammaTotal;
// 边界限制
if (ratio <= xData[0]) return yData[0];
if (ratio >= xData.Last()) return 1.0;
// 线性插值
for (int i = 0; i < xData.Length - 1; i++)
if (lightStates == null || lightStates.Length != totalLights) return 0;
int firstZero = -1;
for (int i = 0; i < totalLights; i++)
{
if (ratio >= xData[i] && ratio <= xData[i + 1])
if (lightStates[i] == 0)
{
double t = (ratio - xData[i]) / (xData[i + 1] - xData[i]);
return yData[i] + t * (yData[i + 1] - yData[i]);
firstZero = i;
break;
}
}
if (firstZero == -1) return 90.0; // 全亮
if (firstZero == 0) return 0.0; // 顶部即被遮
// 边界角度 = 上一个灯条的结束角 = (firstZero) * verticalAngleStep
return firstZero * verticalAngleStep;
}
/// <summary>
/// 极坐标梯形法积分面积(半径数组单位:度,步长单位:度)
/// </summary>
public static double IntegrateArea(double[] radiiDeg, double deltaThetaDeg)
{
if (radiiDeg == null || radiiDeg.Length < 2) return 0;
double deltaRad = deltaThetaDeg * Math.PI / 180.0;
double sum = 0.0;
for (int i = 0; i < radiiDeg.Length - 1; i++)
{
double r1 = radiiDeg[i];
double r2 = radiiDeg[i + 1];
double avgRSq = (r1 * r1 + r2 * r2) / 2.0;
sum += avgRSq * deltaRad;
}
return 0.5 * sum;
}
/// <summary>
/// 根据左右眼边界半径数组,计算总视野(并集)和双目视野(交集)面积
/// </summary>
public static (double totalArea, double biArea) ComputeTotalAndBinocularArea(double[] leftRadii, double[] rightRadii, double deltaThetaDeg)
{
int len = Math.Min(leftRadii.Length, rightRadii.Length);
double[] totalRadii = new double[len];
double[] biRadii = new double[len];
for (int i = 0; i < len; i++)
{
totalRadii[i] = Math.Max(leftRadii[i], rightRadii[i]);
biRadii[i] = Math.Min(leftRadii[i], rightRadii[i]);
}
double totalArea = IntegrateArea(totalRadii, deltaThetaDeg);
double biArea = IntegrateArea(biRadii, deltaThetaDeg);
return (totalArea, biArea);
}
/// <summary>
/// 保存率计算公式(带γ校正)
/// </summary>
public static double ComputePreservation(double measuredArea, double standardArea, double gamma)
{
if (standardArea <= 0) return 0;
return gamma * (measuredArea / standardArea) * 100.0;
}
/// <summary>
/// 根据面积比获取γ校正系数GB 2890-2022 图D.4
/// </summary>
public static double GetGammaByRatio(double ratio)
{
double[] x = { 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 };
double[] y = { 1.22, 1.18, 1.14, 1.10, 1.06, 1.03, 1.02, 1.01, 1.00 };
if (ratio <= x[0]) return y[0];
if (ratio >= x.Last()) return y.Last();
for (int i = 0; i < x.Length - 1; i++)
{
if (ratio >= x[i] && ratio <= x[i + 1])
{
double t = (ratio - x[i]) / (x[i + 1] - x[i]);
return y[i] + t * (y[i + 1] - y[i]);
}
}
return 1.0;
}
public static double CalculateEllipseArea(int[] lightData, List<(int m, int n)> lightPositions)
{
// 参数有效性检查
if (lightData == null || lightPositions == null) return 0;
// 限制循环长度,避免越界(三个长度取最小)
int maxIdx = Math.Min(lightData.Length, Math.Min(lightPositions.Count, totalLights));
// 收集亮灯坐标
List<System.Drawing.Point> brightPoints = new List<System.Drawing.Point>();
for (int i = 0; i < maxIdx; i++)
{
if (lightData[i] == 1)
{
var (m, n) = lightPositions[i];
System.Drawing.Point p = GetLightPoint(m, n);
brightPoints.Add(p);
}
}
int brightCount = brightPoints.Count;
//System.Diagnostics.Debug.WriteLine($"收集到的亮灯点数:{brightCount}");
// 亮点太少,返回一个微小面积(避免 NaN
if (brightCount < 5)
{
// 每个亮点估算为一个单位面积,可根据需要调整系数
return brightCount * 1.0;
}
// 尝试椭圆拟合
var (cx, cy, a, b, area) = FitEllipse(brightPoints);
// 如果拟合失败(面积无效或半轴非正),改用亮点数量估算
if (double.IsNaN(area) || double.IsInfinity(area) || a <= 0 || b <= 0)
{
System.Diagnostics.Debug.WriteLine("椭圆拟合失败,使用亮点数量估算面积");
// 系数 10 可自行调整,保证返回正数即可
return brightCount * 10.0;
}
return area;
}
/// 生成设备全部243盏灯的(m,n)位置 上爪1条、下爪1条、左右共用1条各81灯
public static System.Drawing.Point GetLightPoint(int m, int n)
{
double radH, radV;
// 上爪灯条n=0水平角固定m控制垂直角
if (n == 0)
{
radH = 0;
radV = m * verticalAngleStep * Math.PI / 180;
}
// 下爪灯条n=1水平角固定为180°m控制垂直角
else if (n == 1)
{
radH = Math.PI;
radV = m * verticalAngleStep * Math.PI / 180;
}
// 左右共用灯条n=2垂直角固定m控制水平角
else
{
radH = m * horizontalAngleStep * 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椭圆的短半轴长度较小的那个半径
public 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个点来拟合椭圆");
return new(0, 0, 0, 0, 0);
// 这里是正确写法
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);
}
}
}
}