308 lines
12 KiB
C#
308 lines
12 KiB
C#
|
||
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
|
||
{
|
||
// 设备固定参数
|
||
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 static readonly double _standardTotalArea = (Math.PI * R * R) / 100;
|
||
|
||
//计算双目视野
|
||
|
||
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);
|
||
}
|
||
|
||
////下方视野 下方视野角度 = 人眼到「最低亮灯条」的夹角
|
||
|
||
//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++)
|
||
// {
|
||
// // 只处理亮灯的情况
|
||
// if (lightData[i] == 1)
|
||
// {
|
||
// if (lightPositions.Count < lightData.Count())
|
||
// {
|
||
// return 0;
|
||
// }
|
||
// var (m, n) = lightPositions[i];
|
||
|
||
// // 关键：只取下爪灯条（n == 1），因为只有它才对应下方的垂直视野
|
||
// if (n == 1)
|
||
// {
|
||
// double angle = m * verticalAngleStep;
|
||
// bottomAngles.Add(angle);
|
||
// }
|
||
// }
|
||
// }
|
||
|
||
// // 没有亮灯，返回0
|
||
// if (bottomAngles.Count == 0)
|
||
// return 0;
|
||
|
||
// // 最大角度 = 最下方的亮灯条，也就是下方视野的边界
|
||
// return bottomAngles.Max();
|
||
//}
|
||
/// <summary>
|
||
/// 下方视野角度 = 从正下方（90°）向顶部扫描，第一个被遮挡的边界角度（度）
|
||
/// 使用下爪灯条（n==1）的径向数据
|
||
/// </summary>
|
||
public static double CalculateBottomViewAngle(int[] lightData, List<(int m, int n)> lightPositions)
|
||
{
|
||
// 提取下爪灯条（n==1）的数据，共81个灯
|
||
int radialCount = 81;
|
||
int[] radialData = new int[radialCount];
|
||
int idx = 0;
|
||
for (int i = 0; i < lightData.Length && idx < radialCount; i++)
|
||
{
|
||
var (m, n) = lightPositions[i];
|
||
if (n == 1)
|
||
{
|
||
radialData[idx++] = lightData[i];
|
||
}
|
||
}
|
||
if (idx != radialCount) return 0; // 数据不完整
|
||
|
||
// 反转数组（使索引0对应物理底部90°，索引80对应顶部0°）
|
||
Array.Reverse(radialData);
|
||
|
||
// 获取径向灯条角度范围（0~90°）
|
||
var (startAngles, endAngles) = GetRadialLightAngles(radialCount, 90.0);
|
||
double boundary = ComputeBoundaryAngle(radialData, startAngles, endAngles);
|
||
// 转换为从顶部开始的极角（即标准下方视野角度）
|
||
return 90 - boundary;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 根据灯条0/1状态计算径向边界角度（度）
|
||
/// </summary>
|
||
public static double ComputeBoundaryAngle(int[] lightStates, double[] startAngles, double[] endAngles)
|
||
{
|
||
if (lightStates == null || startAngles == null || endAngles == null) return 0;
|
||
int n = lightStates.Length;
|
||
if (n == 0 || startAngles.Length != n || endAngles.Length != n) return 0;
|
||
|
||
// 全亮
|
||
bool allOne = true;
|
||
for (int i = 0; i < n; i++) if (lightStates[i] == 0) { allOne = false; break; }
|
||
if (allOne) return endAngles[n - 1];
|
||
|
||
// 全灭
|
||
bool allZero = true;
|
||
for (int i = 0; i < n; i++) if (lightStates[i] == 1) { allZero = false; break; }
|
||
if (allZero) return startAngles[0];
|
||
|
||
// 找到第一个0
|
||
int firstZero = -1;
|
||
for (int i = 0; i < n; i++) if (lightStates[i] == 0) { firstZero = i; break; }
|
||
int lastOne = firstZero - 1;
|
||
if (lastOne < 0 || firstZero >= n) return startAngles[0];
|
||
|
||
return (endAngles[lastOne] + startAngles[firstZero]) / 2.0;
|
||
}
|
||
|
||
/// <summary>
|
||
/// 获取径向灯条的角度范围（等分0~maxAngle，灯条数为lightCount）
|
||
/// </summary>
|
||
public static (double[] start, double[] end) GetRadialLightAngles(int lightCount, double maxAngle = 90.0)
|
||
{
|
||
double step = maxAngle / (lightCount - 1);
|
||
double[] start = new double[lightCount];
|
||
double[] end = new double[lightCount];
|
||
for (int i = 0; i < lightCount; i++)
|
||
{
|
||
start[i] = i * step;
|
||
end[i] = (i + 1) * step;
|
||
}
|
||
end[lightCount - 1] = maxAngle;
|
||
return (start, end);
|
||
}
|
||
|
||
/// <summary>
|
||
/// GB2890-2022 自动获取 总视野/双目视野 校正系数γ
|
||
/// </summary>
|
||
/// <param name="ratio">实测面积/标准面积 比值(0~1)</param>
|
||
/// <returns>校正系数 γ</returns>
|
||
public static double GetVisionGamma(double ratio)
|
||
{
|
||
// 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 (ratio >= xData[i] && ratio <= xData[i + 1])
|
||
{
|
||
double t = (ratio - xData[i]) / (xData[i + 1] - xData[i]);
|
||
return yData[i] + t * (yData[i + 1] - yData[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);
|
||
}
|
||
|
||
}
|
||
|
||
|
||
}
|