znlgis 博客

GIS开发与技术分享

第四章:数学库详解

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的图元系统和渲染管线至关重要。


上一章第三章:核心架构与分层设计

下一章第五章:核心数据模型详解