第四章:数学库详解
4.1 数学库概述
4.1.1 LightCAD.MathLib的定位
LightCAD.MathLib是整个LightCAD平台的数学基础层,为上层所有模块提供几何计算和数学运算能力。作为最底层的模块,它具有以下特点:
- 零CAD业务依赖:不依赖任何CAD业务逻辑,纯粹的数学计算
- 高精度:全部使用双精度浮点数(double),满足工程计算精度需求
- 外部依赖精简:仅依赖MathNet.Numerics 5.0.0用于高级数学运算
- 性能优化:关键路径上的计算进行了内联和缓存优化
4.1.2 模块结构
LightCAD.MathLib/
├── LightCAD.MathLib.csproj # 项目文件
├── Basic/ # 基本数学类型
│ ├── Point2d.cs # 二维点
│ ├── Point3d.cs # 三维点
│ ├── Vector2d.cs # 二维向量
│ ├── Vector3d.cs # 三维向量
│ ├── Matrix3d.cs # 3×3矩阵
│ ├── Matrix4d.cs # 4×4变换矩阵
│ ├── BoundingBox2d.cs # 二维包围盒
│ ├── BoundingBox3d.cs # 三维包围盒
│ ├── Plane.cs # 平面
│ └── Transform.cs # 变换组合
├── Curve/ # 曲线数学
│ ├── BezierCurve.cs # 贝塞尔曲线
│ ├── NurbsCurve.cs # NURBS曲线
│ ├── ArcMath.cs # 圆弧数学
│ ├── EllipseMath.cs # 椭圆数学
│ └── SplineMath.cs # 样条数学
├── Intersection/ # 交集计算
│ ├── LineLineIntersect.cs # 直线-直线
│ ├── LineArcIntersect.cs # 直线-圆弧
│ ├── LineCircleIntersect.cs # 直线-圆
│ ├── ArcArcIntersect.cs # 圆弧-圆弧
│ └── CurveCurveIntersect.cs # 曲线-曲线
└── Utilities/ # 工具函数
├── Tolerance.cs # 精度容差
├── AngleUtils.cs # 角度工具
├── GeometryUtils.cs # 几何工具
└── MathConstants.cs # 数学常量
4.2 基本数学类型
4.2.1 二维点(Point2d)
Point2d是二维空间中点的表示,是所有二维几何运算的基础:
public struct Point2d
{
public double X { get; set; }
public double Y { get; set; }
// 构造函数
public Point2d(double x, double y)
{
X = x;
Y = y;
}
// 原点
public static readonly Point2d Origin = new Point2d(0, 0);
// 距离计算
public double DistanceTo(Point2d other)
{
var dx = X - other.X;
var dy = Y - other.Y;
return Math.Sqrt(dx * dx + dy * dy);
}
// 中点
public Point2d MidPoint(Point2d other)
{
return new Point2d(
(X + other.X) / 2.0,
(Y + other.Y) / 2.0
);
}
// 点到向量的转换
public Vector2d ToVector()
{
return new Vector2d(X, Y);
}
// 运算符重载
public static Point2d operator +(Point2d p, Vector2d v)
=> new Point2d(p.X + v.X, p.Y + v.Y);
public static Vector2d operator -(Point2d a, Point2d b)
=> new Vector2d(a.X - b.X, a.Y - b.Y);
// 容差比较
public bool IsEqualTo(Point2d other, double tolerance = 1e-10)
{
return Math.Abs(X - other.X) < tolerance
&& Math.Abs(Y - other.Y) < tolerance;
}
}
使用示例:
var p1 = new Point2d(0, 0);
var p2 = new Point2d(3, 4);
// 计算距离
double dist = p1.DistanceTo(p2); // 5.0
// 计算中点
var mid = p1.MidPoint(p2); // (1.5, 2.0)
// 点的偏移
var offset = p1 + new Vector2d(10, 20); // (10, 20)
// 容差比较
bool equal = p1.IsEqualTo(new Point2d(1e-11, 1e-11)); // true
4.2.2 三维点(Point3d)
Point3d扩展了二维点到三维空间:
public struct Point3d
{
public double X { get; set; }
public double Y { get; set; }
public double Z { get; set; }
public Point3d(double x, double y, double z)
{
X = x; Y = y; Z = z;
}
public static readonly Point3d Origin = new Point3d(0, 0, 0);
// 距离计算
public double DistanceTo(Point3d other)
{
var dx = X - other.X;
var dy = Y - other.Y;
var dz = Z - other.Z;
return Math.Sqrt(dx * dx + dy * dy + dz * dz);
}
// 投影到XY平面
public Point2d ToPoint2d()
{
return new Point2d(X, Y);
}
// 矩阵变换
public Point3d TransformBy(Matrix4d matrix)
{
return matrix.Transform(this);
}
}
4.2.3 二维向量(Vector2d)
public struct Vector2d
{
public double X { get; set; }
public double Y { get; set; }
// 预定义单位向量
public static readonly Vector2d XAxis = new Vector2d(1, 0);
public static readonly Vector2d YAxis = new Vector2d(0, 1);
public static readonly Vector2d Zero = new Vector2d(0, 0);
// 长度
public double Length => Math.Sqrt(X * X + Y * Y);
// 单位化
public Vector2d Normalize()
{
var len = Length;
if (len < 1e-15) return Zero;
return new Vector2d(X / len, Y / len);
}
// 点积
public double Dot(Vector2d other)
{
return X * other.X + Y * other.Y;
}
// 叉积(返回Z分量)
public double Cross(Vector2d other)
{
return X * other.Y - Y * other.X;
}
// 旋转
public Vector2d Rotate(double angle)
{
var cos = Math.Cos(angle);
var sin = Math.Sin(angle);
return new Vector2d(
X * cos - Y * sin,
X * sin + Y * cos
);
}
// 垂直向量
public Vector2d Perpendicular()
{
return new Vector2d(-Y, X);
}
// 两向量之间的角度
public double AngleTo(Vector2d other)
{
return Math.Atan2(Cross(other), Dot(other));
}
}
4.2.4 三维向量(Vector3d)
public struct Vector3d
{
public double X { get; set; }
public double Y { get; set; }
public double Z { get; set; }
public static readonly Vector3d XAxis = new Vector3d(1, 0, 0);
public static readonly Vector3d YAxis = new Vector3d(0, 1, 0);
public static readonly Vector3d ZAxis = new Vector3d(0, 0, 1);
public static readonly Vector3d Zero = new Vector3d(0, 0, 0);
public double Length => Math.Sqrt(X * X + Y * Y + Z * Z);
// 单位化
public Vector3d Normalize()
{
var len = Length;
if (len < 1e-15) return Zero;
return new Vector3d(X / len, Y / len, Z / len);
}
// 点积
public double Dot(Vector3d other)
{
return X * other.X + Y * other.Y + Z * other.Z;
}
// 叉积
public Vector3d Cross(Vector3d other)
{
return new Vector3d(
Y * other.Z - Z * other.Y,
Z * other.X - X * other.Z,
X * other.Y - Y * other.X
);
}
// 绕轴旋转
public Vector3d RotateAround(Vector3d axis, double angle)
{
var normalizedAxis = axis.Normalize();
var cos = Math.Cos(angle);
var sin = Math.Sin(angle);
var dot = Dot(normalizedAxis);
return this * cos
+ normalizedAxis.Cross(this) * sin
+ normalizedAxis * dot * (1 - cos);
}
}
4.2.5 变换矩阵(Matrix4d)
Matrix4d是LightCAD中最核心的数学类型之一,用于表示三维空间中的仿射变换:
public struct Matrix4d
{
// 4x4矩阵数据(行优先存储)
private double[] data; // 16个元素
// 单位矩阵
public static readonly Matrix4d Identity = new Matrix4d(
1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1
);
// 创建平移矩阵
public static Matrix4d CreateTranslation(double x, double y, double z)
{
var m = Identity;
m[0, 3] = x;
m[1, 3] = y;
m[2, 3] = z;
return m;
}
// 创建缩放矩阵
public static Matrix4d CreateScale(double sx, double sy, double sz)
{
var m = Identity;
m[0, 0] = sx;
m[1, 1] = sy;
m[2, 2] = sz;
return m;
}
// 创建绕X轴旋转矩阵
public static Matrix4d CreateRotationX(double angle)
{
var cos = Math.Cos(angle);
var sin = Math.Sin(angle);
var m = Identity;
m[1, 1] = cos; m[1, 2] = -sin;
m[2, 1] = sin; m[2, 2] = cos;
return m;
}
// 创建绕Y轴旋转矩阵
public static Matrix4d CreateRotationY(double angle)
{
var cos = Math.Cos(angle);
var sin = Math.Sin(angle);
var m = Identity;
m[0, 0] = cos; m[0, 2] = sin;
m[2, 0] = -sin; m[2, 2] = cos;
return m;
}
// 创建绕Z轴旋转矩阵
public static Matrix4d CreateRotationZ(double angle)
{
var cos = Math.Cos(angle);
var sin = Math.Sin(angle);
var m = Identity;
m[0, 0] = cos; m[0, 1] = -sin;
m[1, 0] = sin; m[1, 1] = cos;
return m;
}
// 矩阵乘法
public static Matrix4d operator *(Matrix4d a, Matrix4d b)
{
// 标准4x4矩阵乘法
var result = new Matrix4d();
for (int i = 0; i < 4; i++)
for (int j = 0; j < 4; j++)
{
double sum = 0;
for (int k = 0; k < 4; k++)
sum += a[i, k] * b[k, j];
result[i, j] = sum;
}
return result;
}
// 变换点
public Point3d Transform(Point3d point)
{
var x = data[0] * point.X + data[1] * point.Y + data[2] * point.Z + data[3];
var y = data[4] * point.X + data[5] * point.Y + data[6] * point.Z + data[7];
var z = data[8] * point.X + data[9] * point.Y + data[10] * point.Z + data[11];
return new Point3d(x, y, z);
}
// 求逆矩阵
public Matrix4d Inverse()
{
// 使用高斯消元或伴随矩阵法求逆
// ...
}
}
变换组合示例:
// 创建一个复合变换:先缩放,再旋转,最后平移
var scale = Matrix4d.CreateScale(2, 2, 2);
var rotate = Matrix4d.CreateRotationZ(Math.PI / 4); // 旋转45度
var translate = Matrix4d.CreateTranslation(10, 20, 0);
// 变换的应用顺序是从右到左(矩阵乘法的特性)
var transform = translate * rotate * scale;
// 应用到点
var original = new Point3d(1, 0, 0);
var result = transform.Transform(original);
4.2.6 包围盒(BoundingBox)
包围盒用于快速的空间查询和碰撞检测:
// 二维包围盒
public struct BoundingBox2d
{
public Point2d Min { get; set; }
public Point2d Max { get; set; }
// 从点集创建
public static BoundingBox2d FromPoints(IEnumerable<Point2d> points)
{
var minX = double.MaxValue;
var minY = double.MaxValue;
var maxX = double.MinValue;
var maxY = double.MinValue;
foreach (var p in points)
{
minX = Math.Min(minX, p.X);
minY = Math.Min(minY, p.Y);
maxX = Math.Max(maxX, p.X);
maxY = Math.Max(maxY, p.Y);
}
return new BoundingBox2d
{
Min = new Point2d(minX, minY),
Max = new Point2d(maxX, maxY)
};
}
// 包含测试
public bool Contains(Point2d point)
{
return point.X >= Min.X && point.X <= Max.X
&& point.Y >= Min.Y && point.Y <= Max.Y;
}
// 相交测试
public bool Intersects(BoundingBox2d other)
{
return Min.X <= other.Max.X && Max.X >= other.Min.X
&& Min.Y <= other.Max.Y && Max.Y >= other.Min.Y;
}
// 合并
public BoundingBox2d Union(BoundingBox2d other)
{
return new BoundingBox2d
{
Min = new Point2d(
Math.Min(Min.X, other.Min.X),
Math.Min(Min.Y, other.Min.Y)),
Max = new Point2d(
Math.Max(Max.X, other.Max.X),
Math.Max(Max.Y, other.Max.Y))
};
}
// 宽度和高度
public double Width => Max.X - Min.X;
public double Height => Max.Y - Min.Y;
public Point2d Center => Min.MidPoint(Max);
}
4.3 曲线数学
4.3.1 圆弧数学(ArcMath)
public static class ArcMath
{
/// <summary>
/// 通过三点计算圆弧(起点、经过点、终点)
/// </summary>
public static ArcData FromThreePoints(Point2d start, Point2d through, Point2d end)
{
// 通过三点求圆心
// 使用两条弦的垂直平分线的交点
var mid1 = start.MidPoint(through);
var mid2 = through.MidPoint(end);
var dir1 = (through - start).Perpendicular();
var dir2 = (end - through).Perpendicular();
var center = LineLineIntersect.Compute(mid1, dir1, mid2, dir2);
var radius = center.DistanceTo(start);
// 计算起始角和终止角
var startAngle = Math.Atan2(start.Y - center.Y, start.X - center.X);
var endAngle = Math.Atan2(end.Y - center.Y, end.X - center.X);
return new ArcData(center, radius, startAngle, endAngle);
}
/// <summary>
/// 计算圆弧上的点
/// </summary>
public static Point2d PointOnArc(Point2d center, double radius,
double angle)
{
return new Point2d(
center.X + radius * Math.Cos(angle),
center.Y + radius * Math.Sin(angle)
);
}
/// <summary>
/// 计算圆弧长度
/// </summary>
public static double ArcLength(double radius, double startAngle,
double endAngle)
{
var sweepAngle = endAngle - startAngle;
if (sweepAngle < 0) sweepAngle += 2 * Math.PI;
return radius * sweepAngle;
}
/// <summary>
/// 将圆弧离散为点序列
/// </summary>
public static List<Point2d> Tessellate(Point2d center, double radius,
double startAngle, double endAngle,
int segments)
{
var points = new List<Point2d>(segments + 1);
var sweep = endAngle - startAngle;
if (sweep < 0) sweep += 2 * Math.PI;
for (int i = 0; i <= segments; i++)
{
var t = (double)i / segments;
var angle = startAngle + sweep * t;
points.Add(PointOnArc(center, radius, angle));
}
return points;
}
}
4.3.2 贝塞尔曲线(BezierCurve)
public class BezierCurve
{
private Point2d[] controlPoints;
public BezierCurve(Point2d[] controlPoints)
{
this.controlPoints = controlPoints;
}
/// <summary>
/// 使用De Casteljau算法计算曲线上的点
/// </summary>
public Point2d Evaluate(double t)
{
var points = (Point2d[])controlPoints.Clone();
int n = points.Length;
for (int r = 1; r < n; r++)
{
for (int i = 0; i < n - r; i++)
{
points[i] = new Point2d(
(1 - t) * points[i].X + t * points[i + 1].X,
(1 - t) * points[i].Y + t * points[i + 1].Y
);
}
}
return points[0];
}
/// <summary>
/// 计算曲线在参数t处的切线方向
/// </summary>
public Vector2d Tangent(double t)
{
int n = controlPoints.Length - 1;
var tangentPoints = new Point2d[n];
for (int i = 0; i < n; i++)
{
tangentPoints[i] = new Point2d(
n * (controlPoints[i + 1].X - controlPoints[i].X),
n * (controlPoints[i + 1].Y - controlPoints[i].Y)
);
}
var bezier = new BezierCurve(tangentPoints);
var point = bezier.Evaluate(t);
return new Vector2d(point.X, point.Y).Normalize();
}
/// <summary>
/// 将贝塞尔曲线离散为折线段
/// </summary>
public List<Point2d> Tessellate(int segments = 32)
{
var points = new List<Point2d>(segments + 1);
for (int i = 0; i <= segments; i++)
{
var t = (double)i / segments;
points.Add(Evaluate(t));
}
return points;
}
}
4.3.3 NURBS曲线(NurbsCurve)
NURBS(Non-Uniform Rational B-Spline)是CAD中最通用的曲线表示方法:
public class NurbsCurve
{
public Point3d[] ControlPoints { get; }
public double[] Weights { get; }
public double[] KnotVector { get; }
public int Degree { get; }
public NurbsCurve(Point3d[] controlPoints, double[] weights,
double[] knotVector, int degree)
{
ControlPoints = controlPoints;
Weights = weights;
KnotVector = knotVector;
Degree = degree;
}
/// <summary>
/// 计算NURBS曲线上参数t处的点
/// </summary>
public Point3d Evaluate(double t)
{
int n = ControlPoints.Length - 1;
double x = 0, y = 0, z = 0, w = 0;
for (int i = 0; i <= n; i++)
{
double basis = BSplineBasis(i, Degree, t) * Weights[i];
x += basis * ControlPoints[i].X;
y += basis * ControlPoints[i].Y;
z += basis * ControlPoints[i].Z;
w += basis;
}
return new Point3d(x / w, y / w, z / w);
}
/// <summary>
/// 计算B样条基函数(Cox-de Boor递推公式)
/// </summary>
private double BSplineBasis(int i, int p, double t)
{
if (p == 0)
{
return (t >= KnotVector[i] && t < KnotVector[i + 1]) ? 1.0 : 0.0;
}
double left = 0, right = 0;
var denom1 = KnotVector[i + p] - KnotVector[i];
if (Math.Abs(denom1) > 1e-15)
left = (t - KnotVector[i]) / denom1 * BSplineBasis(i, p - 1, t);
var denom2 = KnotVector[i + p + 1] - KnotVector[i + 1];
if (Math.Abs(denom2) > 1e-15)
right = (KnotVector[i + p + 1] - t) / denom2 * BSplineBasis(i + 1, p - 1, t);
return left + right;
}
}
4.4 交集计算
4.4.1 直线-直线求交
public static class LineLineIntersect
{
/// <summary>
/// 计算两条直线的交点
/// </summary>
/// <param name="p1">直线1上的点</param>
/// <param name="d1">直线1的方向</param>
/// <param name="p2">直线2上的点</param>
/// <param name="d2">直线2的方向</param>
/// <returns>交点,如果平行则返回null</returns>
public static Point2d? Compute(Point2d p1, Vector2d d1,
Point2d p2, Vector2d d2)
{
var cross = d1.Cross(d2);
if (Math.Abs(cross) < 1e-10)
return null; // 平行或重合
var diff = p2 - p1;
var t = diff.Cross(d2) / cross;
return p1 + d1 * t;
}
/// <summary>
/// 计算两条线段的交点(考虑端点范围)
/// </summary>
public static Point2d? ComputeSegments(Point2d a1, Point2d a2,
Point2d b1, Point2d b2)
{
var d1 = a2 - a1;
var d2 = b2 - b1;
var cross = d1.Cross(d2);
if (Math.Abs(cross) < 1e-10)
return null;
var diff = b1 - a1;
var t1 = diff.Cross(d2) / cross;
var t2 = diff.Cross(d1) / cross;
if (t1 >= 0 && t1 <= 1 && t2 >= 0 && t2 <= 1)
return a1 + d1 * t1;
return null;
}
}
4.4.2 直线-圆求交
public static class LineCircleIntersect
{
/// <summary>
/// 计算直线与圆的交点
/// </summary>
public static List<Point2d> Compute(Point2d linePoint, Vector2d lineDir,
Point2d center, double radius)
{
var result = new List<Point2d>();
var dir = lineDir.Normalize();
// 直线参数方程:P = linePoint + t * dir
// 圆方程:|P - center|² = radius²
var oc = linePoint - center;
var a = dir.Dot(dir);
var b = 2.0 * oc.Dot(dir);
var c = oc.Dot(oc) - radius * radius;
var discriminant = b * b - 4 * a * c;
if (discriminant < -1e-10)
return result; // 无交点
if (Math.Abs(discriminant) < 1e-10)
{
// 相切,一个交点
var t = -b / (2 * a);
result.Add(linePoint + dir * t);
}
else
{
// 两个交点
var sqrtD = Math.Sqrt(discriminant);
var t1 = (-b + sqrtD) / (2 * a);
var t2 = (-b - sqrtD) / (2 * a);
result.Add(linePoint + dir * t1);
result.Add(linePoint + dir * t2);
}
return result;
}
}
4.4.3 圆弧-圆弧求交
public static class ArcArcIntersect
{
/// <summary>
/// 计算两个圆的交点
/// </summary>
public static List<Point2d> ComputeCircles(Point2d c1, double r1,
Point2d c2, double r2)
{
var result = new List<Point2d>();
var d = c1.DistanceTo(c2);
// 判断是否有交点
if (d > r1 + r2 + 1e-10) return result; // 相离
if (d < Math.Abs(r1 - r2) - 1e-10) return result; // 内含
if (d < 1e-10 && Math.Abs(r1 - r2) < 1e-10) return result; // 重合
// 计算交点
var a = (r1 * r1 - r2 * r2 + d * d) / (2 * d);
var h = Math.Sqrt(Math.Max(0, r1 * r1 - a * a));
var px = c1.X + a * (c2.X - c1.X) / d;
var py = c1.Y + a * (c2.Y - c1.Y) / d;
if (Math.Abs(h) < 1e-10)
{
// 相切
result.Add(new Point2d(px, py));
}
else
{
// 两个交点
var ox = h * (c2.Y - c1.Y) / d;
var oy = h * (c2.X - c1.X) / d;
result.Add(new Point2d(px + ox, py - oy));
result.Add(new Point2d(px - ox, py + oy));
}
return result;
}
}
4.5 工具函数
4.5.1 精度容差(Tolerance)
在CAD计算中,浮点数的精度问题是一个永恒的挑战。LightCAD通过容差系统来处理这个问题:
public static class Tolerance
{
// 全局容差值
public const double Global = 1e-10;
// 距离容差
public const double Distance = 1e-6;
// 角度容差(弧度)
public const double Angle = 1e-10;
// 参数容差(曲线参数)
public const double Parameter = 1e-10;
// 容差比较方法
public static bool IsZero(double value, double tolerance = Global)
{
return Math.Abs(value) < tolerance;
}
public static bool IsEqual(double a, double b, double tolerance = Global)
{
return Math.Abs(a - b) < tolerance;
}
public static int Compare(double a, double b, double tolerance = Global)
{
if (a < b - tolerance) return -1;
if (a > b + tolerance) return 1;
return 0;
}
}
4.5.2 角度工具(AngleUtils)
public static class AngleUtils
{
// 角度转弧度
public static double DegToRad(double degrees)
{
return degrees * Math.PI / 180.0;
}
// 弧度转角度
public static double RadToDeg(double radians)
{
return radians * 180.0 / Math.PI;
}
// 将角度规范化到[0, 2π)范围
public static double Normalize(double angle)
{
angle = angle % (2 * Math.PI);
if (angle < 0) angle += 2 * Math.PI;
return angle;
}
// 将角度规范化到[-π, π)范围
public static double NormalizeSigned(double angle)
{
angle = Normalize(angle);
if (angle > Math.PI) angle -= 2 * Math.PI;
return angle;
}
// 计算从角度a到角度b的最短旋转角度
public static double ShortestRotation(double from, double to)
{
var diff = NormalizeSigned(to - from);
return diff;
}
}
4.5.3 几何工具(GeometryUtils)
public static class GeometryUtils
{
/// <summary>
/// 计算点到直线的距离
/// </summary>
public static double PointToLineDistance(Point2d point,
Point2d lineStart,
Vector2d lineDir)
{
var v = point - lineStart;
var projected = v.Dot(lineDir.Normalize());
var closest = lineStart + lineDir.Normalize() * projected;
return point.DistanceTo(closest);
}
/// <summary>
/// 计算点到线段的最近点
/// </summary>
public static Point2d ClosestPointOnSegment(Point2d point,
Point2d segStart,
Point2d segEnd)
{
var dir = segEnd - segStart;
var len = dir.Length;
if (len < 1e-15) return segStart;
var normalized = dir.Normalize();
var t = (point - segStart).Dot(normalized);
if (t <= 0) return segStart;
if (t >= len) return segEnd;
return segStart + normalized * t;
}
/// <summary>
/// 判断点是否在多边形内部(射线法)
/// </summary>
public static bool PointInPolygon(Point2d point, List<Point2d> polygon)
{
int n = polygon.Count;
bool inside = false;
for (int i = 0, j = n - 1; i < n; j = i++)
{
if (((polygon[i].Y > point.Y) != (polygon[j].Y > point.Y)) &&
(point.X < (polygon[j].X - polygon[i].X) *
(point.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) +
polygon[i].X))
{
inside = !inside;
}
}
return inside;
}
/// <summary>
/// 计算多边形面积(Shoelace公式)
/// </summary>
public static double PolygonArea(List<Point2d> polygon)
{
double area = 0;
int n = polygon.Count;
for (int i = 0; i < n; i++)
{
int j = (i + 1) % n;
area += polygon[i].X * polygon[j].Y;
area -= polygon[j].X * polygon[i].Y;
}
return Math.Abs(area) / 2.0;
}
/// <summary>
/// 计算凸包(Graham扫描法)
/// </summary>
public static List<Point2d> ConvexHull(List<Point2d> points)
{
// Graham扫描算法实现
if (points.Count < 3) return new List<Point2d>(points);
// 找到最下方的点作为起点
var start = points.OrderBy(p => p.Y).ThenBy(p => p.X).First();
// 按极角排序
var sorted = points
.Where(p => !p.IsEqualTo(start))
.OrderBy(p => Math.Atan2(p.Y - start.Y, p.X - start.X))
.ToList();
var hull = new List<Point2d> { start };
foreach (var point in sorted)
{
while (hull.Count > 1)
{
var a = hull[hull.Count - 2];
var b = hull[hull.Count - 1];
var cross = (b - a).Cross(point - b);
if (cross <= 0)
hull.RemoveAt(hull.Count - 1);
else
break;
}
hull.Add(point);
}
return hull;
}
}
4.6 MathNet.Numerics集成
4.6.1 高级数学运算
LightCAD通过MathNet.Numerics库提供高级数学能力:
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.Interpolation;
// 矩阵运算
var matrix = Matrix<double>.Build.DenseOfArray(new double[,]
{
{ 1, 2, 3 },
{ 4, 5, 6 },
{ 7, 8, 10 }
});
var inverse = matrix.Inverse();
// 最小二乘拟合
var x = new double[] { 0, 1, 2, 3, 4 };
var y = new double[] { 0.1, 1.0, 3.9, 9.1, 16.0 };
var polyFit = MathNet.Numerics.Fit.Polynomial(x, y, 2);
// 样条插值
var spline = CubicSpline.InterpolateAkima(x, y);
double interpolated = spline.Interpolate(2.5);
4.6.2 方程求解
// 线性方程组求解
var A = Matrix<double>.Build.DenseOfArray(new double[,]
{
{ 3, 2, -1 },
{ 2, -2, 4 },
{ -1, 0.5, -1 }
});
var b = Vector<double>.Build.Dense(new double[] { 1, -2, 0 });
var solution = A.Solve(b);
4.7 数学库的测试
LightCAD为MathLib提供了独立的单元测试项目LightCAD.MathLib.Tests,确保数学计算的正确性:
[TestClass]
public class Point2dTests
{
[TestMethod]
public void Distance_ShouldCalculateCorrectly()
{
var p1 = new Point2d(0, 0);
var p2 = new Point2d(3, 4);
Assert.AreEqual(5.0, p1.DistanceTo(p2), 1e-10);
}
[TestMethod]
public void MidPoint_ShouldCalculateCorrectly()
{
var p1 = new Point2d(0, 0);
var p2 = new Point2d(10, 20);
var mid = p1.MidPoint(p2);
Assert.AreEqual(5.0, mid.X, 1e-10);
Assert.AreEqual(10.0, mid.Y, 1e-10);
}
}
[TestClass]
public class IntersectionTests
{
[TestMethod]
public void LineLineIntersect_PerpendicularLines()
{
var result = LineLineIntersect.Compute(
new Point2d(0, 0), new Vector2d(1, 0), // X轴
new Point2d(5, -5), new Vector2d(0, 1) // x=5的竖直线
);
Assert.IsNotNull(result);
Assert.AreEqual(5.0, result.Value.X, 1e-10);
Assert.AreEqual(0.0, result.Value.Y, 1e-10);
}
[TestMethod]
public void LineLineIntersect_ParallelLines_ReturnsNull()
{
var result = LineLineIntersect.Compute(
new Point2d(0, 0), new Vector2d(1, 0),
new Point2d(0, 1), new Vector2d(1, 0)
);
Assert.IsNull(result);
}
}
4.8 性能优化策略
4.8.1 结构体vs类
MathLib中的基本类型(Point2d、Vector2d等)使用struct而非class:
// 使用struct避免堆分配
public struct Point2d // 值类型
{
public double X;
public double Y;
}
优势:
- 避免GC压力
- 在栈上分配,访问更快
- 更好的内存局部性
4.8.2 内联优化
对于频繁调用的小方法,使用MethodImplOptions.AggressiveInlining提示JIT编译器进行内联:
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public double DistanceTo(Point2d other)
{
var dx = X - other.X;
var dy = Y - other.Y;
return Math.Sqrt(dx * dx + dy * dy);
}
4.8.3 避免不必要的分配
// 避免:创建中间数组
// var points = curve.Tessellate().ToArray();
// 推荐:使用Span或预分配缓冲区
Span<Point2d> buffer = stackalloc Point2d[segmentCount];
curve.TessellateInto(buffer);
4.9 本章小结
本章全面介绍了LightCAD.MathLib数学库的各个组成部分。从基本的点、向量、矩阵类型,到曲线数学和交集计算,再到容差处理和性能优化。MathLib作为LightCAD的数学基础,为上层所有几何计算提供了可靠、高效的支持。理解这些数学基础对于后续学习LightCAD的图元系统和渲染管线至关重要。
上一章:第三章:核心架构与分层设计
下一章:第五章:核心数据模型详解