RectangleCenterCalculator.cs 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430
  1. using MathNet.Numerics.LinearAlgebra;
  2. using System;
  3. using System.Collections.Generic;
  4. using System.Linq;
  5. using System.Text;
  6. using System.Threading.Tasks;
  7. namespace TeamAAS_VP.Core
  8. {
  9. /// <summary>
  10. /// 提供多种方法用于计算由四个角点定义的矩形中心点、长轴角度以及相关矩形参数的辅助类。
  11. /// 角点顺序约定:左上(0) → 右上(1) → 右下(2) → 左下(3)
  12. /// 长轴方向:从左上角(0)指向右上角(1)
  13. /// 角度范围:[-180, 180) 度,0°指向X轴正方向,逆时针为正
  14. /// </summary>
  15. public class RectangleCenterCalculator
  16. {
  17. /// <summary>
  18. /// 矩形中心及方向信息的综合结果
  19. /// </summary>
  20. public class RectangleResult
  21. {
  22. /// <summary>
  23. /// 矩形中心点坐标
  24. /// </summary>
  25. public Vector<double> Center { get; set; }
  26. /// <summary>
  27. /// 矩形长轴方向向量(从左上指向右上的单位向量)
  28. /// </summary>
  29. public Vector<double> MajorAxis { get; set; }
  30. /// <summary>
  31. /// 矩形短轴方向向量(单位向量)
  32. /// </summary>
  33. public Vector<double> MinorAxis { get; set; }
  34. /// <summary>
  35. /// 矩形长度(沿长轴方向,左上到右上的距离)
  36. /// </summary>
  37. public double Length { get; set; }
  38. /// <summary>
  39. /// 矩形宽度(沿短轴方向)
  40. /// </summary>
  41. public double Width { get; set; }
  42. /// <summary>
  43. /// 长轴方向角度[-180, 180)度
  44. /// 从左上指向右上的向量在XY平面的投影角度
  45. /// 0° = 指向X轴正方向,正角度逆时针,负角度顺时针
  46. /// </summary>
  47. public double MajorAxisAngle { get; set; }
  48. /// <summary>
  49. /// 四个角点的坐标(保持输入顺序)
  50. /// </summary>
  51. public List<Vector<double>> Corners { get; set; }
  52. /// <summary>
  53. /// 左上角点坐标
  54. /// </summary>
  55. public Vector<double> TopLeft => Corners != null && Corners.Count >= 1 ? Corners[0] : null;
  56. /// <summary>
  57. /// 右上角点坐标
  58. /// </summary>
  59. public Vector<double> TopRight => Corners != null && Corners.Count >= 2 ? Corners[1] : null;
  60. /// <summary>
  61. /// 右下角点坐标
  62. /// </summary>
  63. public Vector<double> BottomRight => Corners != null && Corners.Count >= 3 ? Corners[2] : null;
  64. /// <summary>
  65. /// 左下角点坐标
  66. /// </summary>
  67. public Vector<double> BottomLeft => Corners != null && Corners.Count >= 4 ? Corners[3] : null;
  68. /// <summary>
  69. /// 矩形的面积
  70. /// </summary>
  71. public double Area { get; set; }
  72. /// <summary>
  73. /// 拟合误差(均方根误差)
  74. /// </summary>
  75. public double FitError { get; set; }
  76. }
  77. /// <summary>
  78. /// 验证角点顺序并进行必要的调整
  79. /// 确保顺序为:左上 → 右上 → 右下 → 左下
  80. /// </summary>
  81. private static List<Vector<double>> ValidateAndOrderCorners(List<Vector<double>> corners)
  82. {
  83. if (corners == null || corners.Count != 4)
  84. throw new ArgumentException("需要4个角点");
  85. // 这里可以添加自动排序逻辑,但既然您已知顺序,直接返回
  86. // 如果未来需要自动排序,可以使用以下逻辑:
  87. // 1. 找到最左边的两个点作为左边界
  88. // 2. 根据Y坐标区分左上和左下
  89. // 3. 计算中心点,确定其他点位置
  90. return corners; // 假设输入已经按正确顺序
  91. }
  92. /// <summary>
  93. /// 计算从左上指向右上的长轴方向角度
  94. /// </summary>
  95. private static double CalculateMajorAxisAngle(Vector<double> topLeft, Vector<double> topRight)
  96. {
  97. // 计算方向向量(从左上指向右上)
  98. Vector<double> direction = topRight - topLeft;
  99. // 只考虑XY平面上的投影
  100. double x = direction[0];
  101. double y = direction[1];
  102. // 使用Math.Atan2计算角度(弧度)
  103. double angleRad = Math.Atan2(y, x);
  104. // 转换为度,得到有符号角度 [-180, 180)
  105. double angleDeg = angleRad * 180.0 / Math.PI;
  106. // Math.Atan2的结果已经是[-180, 180],但确保一下边界
  107. if (angleDeg >= 180.0)
  108. angleDeg -= 360.0;
  109. else if (angleDeg < -180.0)
  110. angleDeg += 360.0;
  111. return angleDeg;
  112. }
  113. /// <summary>
  114. /// 使用最小二乘法思想计算最优矩形中心点及相关参数(对外统一入口)
  115. /// 推荐使用改进的PCA方法,因为它能利用所有点信息进行优化
  116. /// </summary>
  117. public static RectangleResult CalculateOptimalCenter(List<Vector<double>> corners)
  118. {
  119. // 验证和调整角点顺序
  120. var orderedCorners = ValidateAndOrderCorners(corners);
  121. // 使用改进的PCA方法(考虑已知顺序)
  122. return CalculateByPCAWithKnownOrder(orderedCorners);
  123. }
  124. /// <summary>
  125. /// 方法1:基于已知角点顺序的直接计算(最准确)
  126. /// 直接使用左上和右上点计算长轴,然后计算其他参数
  127. /// </summary>
  128. public static RectangleResult CalculateByDirectMethod(List<Vector<double>> corners)
  129. {
  130. var orderedCorners = ValidateAndOrderCorners(corners);
  131. var topLeft = orderedCorners[0];
  132. var topRight = orderedCorners[1];
  133. var bottomRight = orderedCorners[2];
  134. var bottomLeft = orderedCorners[3];
  135. // 计算长轴(从左上指向右上)
  136. Vector<double> majorAxis = (topRight - topLeft).Normalize(2);
  137. // 计算短轴(从左上指向左下)
  138. Vector<double> minorAxis = (bottomLeft - topLeft).Normalize(2);
  139. // 确保短轴与长轴垂直(处理非正交情况)
  140. minorAxis = (minorAxis - majorAxis * minorAxis.DotProduct(majorAxis)).Normalize(2);
  141. int dim = orderedCorners[0].Count;
  142. // 计算中心点(四个角点的平均)
  143. var center = Vector<double>.Build.DenseOfEnumerable(
  144. Enumerable.Range(0, dim)
  145. .Select(i => orderedCorners.Select(c => c[i]).Average()));
  146. // 计算长度和宽度
  147. double length = Distance(topLeft, topRight);
  148. double width = Distance(topLeft, bottomLeft);
  149. // 计算长轴角度
  150. double angle = CalculateMajorAxisAngle(topLeft, topRight);
  151. // 计算拟合误差
  152. double fitError = CalculateFitErrorDirect(center, majorAxis, minorAxis, length, width, orderedCorners);
  153. return new RectangleResult
  154. {
  155. Center = center,
  156. MajorAxis = majorAxis,
  157. MinorAxis = minorAxis,
  158. Length = length,
  159. Width = width,
  160. MajorAxisAngle = angle,
  161. Corners = orderedCorners,
  162. Area = length * width,
  163. FitError = fitError
  164. };
  165. }
  166. /// <summary>
  167. /// 方法2:使用PCA但考虑已知角点顺序(推荐)
  168. /// 结合PCA的鲁棒性和已知顺序的准确性
  169. /// </summary>
  170. public static RectangleResult CalculateByPCAWithKnownOrder(List<Vector<double>> corners)
  171. {
  172. var orderedCorners = ValidateAndOrderCorners(corners);
  173. int dim = orderedCorners[0].Count;
  174. // 步骤1:使用所有点进行PCA得到初步估计
  175. var matrix = Matrix<double>.Build.DenseOfRows(orderedCorners);
  176. var mean = Vector<double>.Build.DenseOfEnumerable(
  177. Enumerable.Range(0, dim).Select(i => matrix.Column(i).Average()));
  178. var centered = matrix.Clone();
  179. for (int i = 0; i < 4; i++)
  180. {
  181. centered.SetRow(i, centered.Row(i) - mean);
  182. }
  183. var covariance = centered.Transpose() * centered / 3.0;
  184. var evd = covariance.Evd();
  185. // 获取特征向量
  186. var eigenvectors = evd.EigenVectors;
  187. // 步骤2:使用已知顺序确定长轴方向
  188. var topLeft = orderedCorners[0];
  189. var topRight = orderedCorners[1];
  190. // 计算实际的长轴方向(从左上指向右上)
  191. Vector<double> actualMajorDirection = (topRight - topLeft).Normalize(2);
  192. // 从PCA的特征向量中找到最接近实际方向的那个
  193. Vector<double> bestAxis = eigenvectors.Column(0);
  194. double bestSimilarity = Math.Abs(actualMajorDirection.DotProduct(bestAxis));
  195. for (int i = 1; i < eigenvectors.ColumnCount; i++)
  196. {
  197. var axis = eigenvectors.Column(i);
  198. double similarity = Math.Abs(actualMajorDirection.DotProduct(axis));
  199. if (similarity > bestSimilarity)
  200. {
  201. bestSimilarity = similarity;
  202. bestAxis = axis;
  203. }
  204. }
  205. // 确保PCA轴的方向与实际方向一致(点积为正)
  206. Vector<double> majorAxis = bestAxis;
  207. if (actualMajorDirection.DotProduct(majorAxis) < 0)
  208. {
  209. majorAxis = -majorAxis;
  210. }
  211. // 计算短轴(与长轴垂直)
  212. Vector<double> minorAxis = CalculatePerpendicularAxis(majorAxis);
  213. // 将点投影到轴上计算中心、长度和宽度
  214. var projMajor = orderedCorners.Select(p => (p - mean).DotProduct(majorAxis)).ToList();
  215. var projMinor = orderedCorners.Select(p => (p - mean).DotProduct(minorAxis)).ToList();
  216. double centerMajor = (projMajor.Min() + projMajor.Max()) / 2.0;
  217. double centerMinor = (projMinor.Min() + projMinor.Max()) / 2.0;
  218. var center = mean + majorAxis * centerMajor + minorAxis * centerMinor;
  219. double length = Math.Abs(projMajor.Max() - projMajor.Min());
  220. double width = Math.Abs(projMinor.Max() - projMinor.Min());
  221. // 计算长轴角度(使用实际的左上-右上方向)
  222. double angle = CalculateMajorAxisAngle(topLeft, topRight);
  223. // 计算拟合误差
  224. double fitError = CalculateFitError(center, majorAxis, minorAxis, length, width, orderedCorners);
  225. return new RectangleResult
  226. {
  227. Center = center,
  228. MajorAxis = majorAxis,
  229. MinorAxis = minorAxis,
  230. Length = length,
  231. Width = width,
  232. MajorAxisAngle = angle,
  233. Corners = orderedCorners,
  234. Area = length * width,
  235. FitError = fitError
  236. };
  237. }
  238. /// <summary>
  239. /// 方法3:对角线优化方法(考虑已知顺序)
  240. /// </summary>
  241. public static RectangleResult CalculateByDiagonalOptimization(List<Vector<double>> corners)
  242. {
  243. var orderedCorners = ValidateAndOrderCorners(corners);
  244. var topLeft = orderedCorners[0];
  245. var topRight = orderedCorners[1];
  246. // 直接使用已知顺序计算长轴角度
  247. double angle = CalculateMajorAxisAngle(topLeft, topRight);
  248. // 计算长轴方向向量
  249. Vector<double> majorAxis = (topRight - topLeft).Normalize(2);
  250. Vector<double> minorAxis = CalculatePerpendicularAxis(majorAxis);
  251. // 使用所有对角线组合计算最佳中心
  252. var diagonalCombinations = new[]
  253. {
  254. (orderedCorners[0], orderedCorners[2]), // 左上-右下
  255. (orderedCorners[1], orderedCorners[3]) // 右上-左下
  256. };
  257. var centers = diagonalCombinations.Select(p => (p.Item1 + p.Item2) / 2.0).ToList();
  258. var center = (centers[0] + centers[1]) / 2.0;
  259. // 计算长度和宽度
  260. double length = Distance(topLeft, topRight);
  261. double width = Distance(topLeft, orderedCorners[3]); // 左上-左下
  262. // 计算拟合误差
  263. double fitError = CalculateFitError(center, majorAxis, minorAxis, length, width, orderedCorners);
  264. return new RectangleResult
  265. {
  266. Center = center,
  267. MajorAxis = majorAxis,
  268. MinorAxis = minorAxis,
  269. Length = length,
  270. Width = width,
  271. MajorAxisAngle = angle,
  272. Corners = orderedCorners,
  273. Area = length * width,
  274. FitError = fitError
  275. };
  276. }
  277. /// <summary>
  278. /// 辅助方法:计算与给定向量垂直的单位向量
  279. /// </summary>
  280. private static Vector<double> CalculatePerpendicularAxis(Vector<double> axis)
  281. {
  282. if (axis == null) throw new ArgumentNullException(nameof(axis));
  283. int n = axis.Count;
  284. if (n == 0) throw new ArgumentException("axis must have positive dimension", nameof(axis));
  285. // 专门处理2D:(-y, x) 是垂直向量
  286. if (n == 2)
  287. {
  288. var perp2 = Vector<double>.Build.DenseOfArray(new[] { -axis[1], axis[0] });
  289. if (perp2.L2Norm() < 1e-12) throw new ArgumentException("axis is zero vector", nameof(axis));
  290. return perp2.Normalize(2);
  291. }
  292. // 一般n维:选一个与axis不共线的标准基向量
  293. int idx = 0;
  294. for (int i = 0; i < n; i++)
  295. {
  296. if (Math.Abs(axis[i]) < 0.9) { idx = i; break; }
  297. }
  298. var temp = Vector<double>.Build.Dense(n, i => i == idx ? 1.0 : 0.0);
  299. var perpendicular = temp - axis * temp.DotProduct(axis);
  300. if (perpendicular.L2Norm() < 1e-12) throw new InvalidOperationException("failed to compute perpendicular vector");
  301. return perpendicular.Normalize(2);
  302. }
  303. /// <summary>
  304. /// 辅助方法:计算两点距离
  305. /// </summary>
  306. private static double Distance(Vector<double> a, Vector<double> b)
  307. {
  308. return (a - b).L2Norm();
  309. }
  310. /// <summary>
  311. /// 辅助方法:计算拟合误差
  312. /// </summary>
  313. private static double CalculateFitError(Vector<double> center, Vector<double> majorAxis, Vector<double> minorAxis, double length, double width, List<Vector<double>> corners)
  314. {
  315. // 生成理想矩形的四个角点
  316. var idealCorners = new List<Vector<double>>
  317. {
  318. center + majorAxis * (length / 2) + minorAxis * (width / 2), // 右上
  319. center + majorAxis * (length / 2) - minorAxis * (width / 2), // 右下
  320. center - majorAxis * (length / 2) - minorAxis * (width / 2), // 左下
  321. center - majorAxis * (length / 2) + minorAxis * (width / 2) // 左上
  322. };
  323. // 注意:理想角点顺序可能与输入不同,需要重新排序匹配
  324. // 这里简单计算所有对应点的最小距离和
  325. double sumSquared = 0;
  326. for (int i = 0; i < 4; i++)
  327. {
  328. // 找到最近的理想角点
  329. double minDist = double.MaxValue;
  330. foreach (var ideal in idealCorners)
  331. {
  332. double dist = Distance(corners[i], ideal);
  333. if (dist < minDist) minDist = dist;
  334. }
  335. sumSquared += minDist * minDist;
  336. }
  337. return Math.Sqrt(sumSquared / 4);
  338. }
  339. /// <summary>
  340. /// 直接计算方法的拟合误差
  341. /// </summary>
  342. private static double CalculateFitErrorDirect(Vector<double> center, Vector<double> majorAxis, Vector<double> minorAxis, double length, double width, List<Vector<double>> corners)
  343. {
  344. // 对于直接方法,我们知道每个角点的理想位置
  345. var idealCorners = new List<Vector<double>>
  346. {
  347. center - majorAxis * (length / 2) + minorAxis * (width / 2), // 左上
  348. center + majorAxis * (length / 2) + minorAxis * (width / 2), // 右上
  349. center + majorAxis * (length / 2) - minorAxis * (width / 2), // 右下
  350. center - majorAxis * (length / 2) - minorAxis * (width / 2) // 左下
  351. };
  352. double sumSquared = 0;
  353. for (int i = 0; i < 4; i++)
  354. {
  355. sumSquared += Math.Pow(Distance(corners[i], idealCorners[i]), 2);
  356. }
  357. return Math.Sqrt(sumSquared / 4);
  358. }
  359. }
  360. }