cad.net 合并样条曲线

合并样条曲线

cad的样条曲线也叫非均匀有理B样条,
数据结构:
利用控制点的(点,权重)集,阶数,非均匀节点向量,再通过构建多项式进行绘制.
也就是不会储存拟合点的.

连接多条样条曲线这个功能:
1,现有样条曲线通过采样获取拟合点,
采样策略
a:均匀采样 b:曲率采样 c:节点向量获取稠密性
2,按顺序链接多条曲线的拟合点点集,注意首尾点最近要把曲线点集逆序处理.
3,拟合点反算控制点.
4,动态阶数,二阶和二阶合并可能是二阶也可能三阶,三阶和三阶就必然是三阶了.
5,绘制完整的样条曲线.
6,如果两条曲线首尾点不吻合<=1e-6,那么这个情况下会实现插值法才对...

注明:
二阶三阶曲线转为直线不在本次范围,单独制作功能.

1,cad选择图元集合是无序的,
要通过判断两条曲线首尾点之间最近点作为链接点,
因此要点对最近算法.
https://www.cnblogs.com/JJBox/p/19003792
2,闭合曲线提示
3,单点曲线
曲线点集太少<2,会导致构造错误,提示并跳过.

// 在端点收集时跳过无效曲线
if (fitPoints.Count < 2) continue;

4,不连续曲线集

// 在距离检查时添加警告
if (next.Distance > tolerance * 10) {
    Console.WriteLine("警告：检测到不连续的曲线，距离=" + next.Distance);
    // 可以选择在此处创建多条独立曲线
}

5,曲线首尾两点,插入到另一曲线腰部.T字.
这种是几何错误,直接提示无法处理,然后退出.
注意要是无最近点的孤立首尾点,我们是允许8字形样条曲线连接的.

代码

下面的代码没有实现上面注意事项,因为懒...
并且还是AI生成的,逃...

using System;
using System.Collections.Generic;
using System.Linq;

namespace SplineConnector {
    /// <summary>专业多条样条曲线合并工具</summary>
    public static class ProfessionalSplineConnector {
        #region 公开主入口
        /// <summary>连接多条样条曲线，返回一条新的NURBS样条</summary>
        public static Spline ConnectSplines(IList<Spline> splinesToConnect, double tolerance = 1e-4) {
            try {
                if (splinesToConnect == null || splinesToConnect.Count == 0) throw new ArgumentException("必须输入至少一条样条");
                var fitPoints = ExtractAndCombineFitPoints(splinesToConnect, tolerance);
                int maxDegree = splinesToConnect.Max(s => s.Degree);
                int feasibleDegree = Math.Min(maxDegree, fitPoints.Count - 1);
                if (feasibleDegree < 1) feasibleDegree = 1;
                var controlPoints = CalculateControlPoints(fitPoints, feasibleDegree, tolerance);
                return CreateConnectedSpline(controlPoints, feasibleDegree, IsClosed(fitPoints, tolerance));
            } catch (Exception ex) {
                Console.WriteLine($"ConnectSplines错误: {ex.Message}");
                throw;
            }
        }
        #endregion

        #region 核心算法
        /// <summary>提取并合并多条样条的拟合点</summary>
        private static List<Point> ExtractAndCombineFitPoints(IList<Spline> splines, double tolerance) {
            try {
                var all = new List<Point>();
                foreach (var spline in splines) {
                    var pts = spline.IsReversed ? spline.GetFitPoints().Reverse() : spline.GetFitPoints();
                    foreach (var p in pts) {
                        if (all.Count == 0 || !ArePointsEqual(all.Last(), p, tolerance)) all.Add(p);
                    }
                }
                if (all.Count > 2 && ArePointsEqual(all.First(), all.Last(), tolerance)) all.RemoveAt(all.Count - 1);
                return all;
            } catch (Exception ex) {
                Console.WriteLine($"ExtractAndCombineFitPoints错误: {ex.Message}");
                throw;
            }
        }

        /// <summary>根据拟合点反算控制点</summary>
        private static List<Point> CalculateControlPoints(IList<Point> fitPoints, int degree, double tolerance) {
            try {
                int n = fitPoints.Count;
                if (n < 2) return new List<Point>(fitPoints);
                
                bool closed = IsClosed(fitPoints, tolerance);
                int ctrlCount = closed ? n + degree : n;
                
                var uk = ChordLengthParameterization(fitPoints);
                var U = GenerateKnotVector(uk, degree, ctrlCount, closed);
                var N = BuildObservationMatrix(uk, U, degree, n, ctrlCount);
                
                var Qx = fitPoints.Select(p => p.X).ToArray();
                var Qy = fitPoints.Select(p => p.Y).ToArray();
                
                var Px = SolveLeastSquares(N, Qx);
                var Py = SolveLeastSquares(N, Qy);
                
                var cps = new List<Point>();
                for (int i = 0; i < ctrlCount; i++) cps.Add(new Point(Px[i], Py[i]));
                
                if (closed) {
                    // For closed curves, ensure the first and last control points match
                    cps[ctrlCount - 1] = cps[0];
                }
                
                return cps;
            } catch (Exception ex) {
                Console.WriteLine($"CalculateControlPoints错误: {ex.Message}");
                throw;
            }
        }

        /// <summary>创建最终样条</summary>
        private static Spline CreateConnectedSpline(IList<Point> controlPoints, int degree, bool closed) {
            try {
                int n = controlPoints.Count;
                var knots = GenerateKnotVector(null, degree, n, closed);
                return new Spline {
                    Degree = degree,
                    ControlPoints = new List<Point>(controlPoints),
                    Knots = knots.ToList(),
                    IsClosed = closed
                };
            } catch (Exception ex) {
                Console.WriteLine($"CreateConnectedSpline错误: {ex.Message}");
                throw;
            }
        }
        #endregion

        #region 数学工具
        /// <summary>弦长参数化算法</summary>
        private static double[] ChordLengthParameterization(IList<Point> pts) {
            try {
                int n = pts.Count;
                var d = new double[n]; d[0] = 0.0;
                double total = 0;
                for (int i = 1; i < n; i++) total += Distance(pts[i - 1], pts[i]);
                if (total < 1e-12) {
                    for (int i = 0; i < n; i++) d[i] = (double)i / (n - 1);
                    return d;
                }
                double cum = 0;
                for (int i = 1; i < n; i++) {
                    cum += Distance(pts[i - 1], pts[i]);
                    d[i] = cum / total;
                }
                return d;
            } catch (Exception ex) {
                Console.WriteLine($"ChordLengthParameterization错误: {ex.Message}");
                throw;
            }
        }

        /// <summary>生成B样条节点向量</summary>
        private static double[] GenerateKnotVector(double[] uk, int p, int ctrlCount, bool closed = false) {
            try {
                int n = ctrlCount;
                int m = n + p + 1;
                var U = new double[m + 1];
                
                if (closed) {
                    // Uniform knot vector for closed curves
                    for (int i = 0; i <= m; i++) {
                        U[i] = (double)i / (m);
                    }
                } else {
                    // Open uniform knot vector
                    for (int i = 0; i <= p; i++) U[i] = 0.0;
                    for (int i = m - p; i <= m; i++) U[i] = 1.0;
                    
                    if (uk != null && uk.Length > 0) {
                        // Use parameter values for internal knots
                        for (int j = 1; j <= n - p; j++) {
                            double sum = 0;
                            for (int k = j; k < j + p; k++) {
                                if (k < uk.Length) sum += uk[k];
                            }
                            U[p + j] = sum / p;
                        }
                    } else {
                        // Default uniform distribution
                        for (int j = 1; j <= n - p; j++) {
                            U[p + j] = (double)j / (n - p + 1);
                        }
                    }
                }
                return U;
            } catch (Exception ex) {
                Console.WriteLine($"GenerateKnotVector错误: {ex.Message}");
                throw;
            }
        }

        /// <summary>构造观测矩阵N</summary>
        private static double[,] BuildObservationMatrix(double[] uk, double[] U, int p, int n, int ctrl) {
            try {
                var N = new double[n, ctrl];
                for (int i = 0; i < n; i++) {
                    for (int j = 0; j < ctrl; j++) {
                        N[i, j] = BasisFunction(p, U, j, uk[i]);
                    }
                }
                return N;
            } catch (Exception ex) {
                Console.WriteLine($"BuildObservationMatrix错误: {ex.Message}");
                throw;
            }
        }

        /// <summary>B样条基函数计算</summary>
        private static double BasisFunction(int p, double[] U, int i, double u) {
            try {
                // Special case for first and last control points
                if ((i == 0 && u == U[0]) || (i == U.Length - p - 2 && u == U[U.Length - 1])) {
                    return 1.0;
                }

                // Check if u is in the range of this basis function
                if (u < U[i] || u >= U[i + p + 1]) return 0.0;

                // Initialize zero-degree basis functions
                double[] N = new double[p + 1];
                for (int j = 0; j <= p; j++) {
                    N[j] = (u >= U[i + j] && u < U[i + j + 1]) ? 1.0 : 0.0;
                }

                // Compute triangular table
                for (int k = 1; k <= p; k++) {
                    double saved = (N[0] == 0) ? 0 : ((u - U[i]) * N[0]) / (U[i + k] - U[i]);
                    for (int j = 0; j < p - k + 1; j++) {
                        double Uleft = U[i + j + 1];
                        double Uright = U[i + j + k + 1];
                        
                        if (N[j + 1] == 0) {
                            N[j] = saved;
                            saved = 0;
                        } else {
                            double temp = N[j + 1] / (Uright - Uleft);
                            N[j] = saved + (Uright - u) * temp;
                            saved = (u - Uleft) * temp;
                        }
                    }
                }
                
                return N[0];
            } catch (Exception ex) {
                Console.WriteLine($"BasisFunction错误: i={i}, p={p}, u={u}, U=[{string.Join(",", U)}]\n{ex.Message}");
                throw;
            }
        }

        /// <summary>最小二乘求解</summary>
        private static double[] SolveLeastSquares(double[,] A, double[] b) {
            try {
                int rows = A.GetLength(0), cols = A.GetLength(1);
                
                // Check for underdetermined system
                if (rows < cols) {
                    throw new ArgumentException("系统欠定 - 方程数少于未知数");
                }

                var AtA = new double[cols, cols];
                for (int i = 0; i < cols; i++) {
                    for (int j = 0; j < cols; j++) {
                        double sum = 0;
                        for (int k = 0; k < rows; k++) sum += A[k, i] * A[k, j];
                        AtA[i, j] = sum;
                    }
                }
                
                var Atb = new double[cols];
                for (int i = 0; i < cols; i++) {
                    double sum = 0;
                    for (int k = 0; k < rows; k++) sum += A[k, i] * b[k];
                    Atb[i] = sum;
                }
                
                return SolveLinearSystem(AtA, Atb);
            } catch (Exception ex) {
                Console.WriteLine($"SolveLeastSquares错误: {ex.Message}");
                throw;
            }
        }

        /// <summary>线性方程组求解</summary>
        private static double[] SolveLinearSystem(double[,] A, double[] b) {
            try {
                int n = b.Length;
                double[] x = new double[n];
                
                // Gaussian elimination with partial pivoting
                for (int i = 0; i < n; i++) {
                    // Partial pivoting
                    int maxRow = i;
                    for (int k = i + 1; k < n; k++) {
                        if (Math.Abs(A[k, i]) > Math.Abs(A[maxRow, i])) {
                            maxRow = k;
                        }
                    }
                    
                    // Swap rows if needed
                    if (maxRow != i) {
                        for (int k = i; k < n; k++) {
                            (A[maxRow, k], A[i, k]) = (A[i, k], A[maxRow, k]);
                        }
                        (b[maxRow], b[i]) = (b[i], b[maxRow]);
                    }
                    
                    // Elimination
                    for (int k = i + 1; k < n; k++) {
                        double factor = A[k, i] / A[i, i];
                        for (int j = i; j < n; j++) {
                            A[k, j] -= factor * A[i, j];
                        }
                        b[k] -= factor * b[i];
                    }
                }
                
                // Back substitution
                for (int i = n - 1; i >= 0; i--) {
                    x[i] = b[i];
                    for (int j = i + 1; j < n; j++) {
                        x[i] -= A[i, j] * x[j];
                    }
                    x[i] /= A[i, i];
                }
                
                return x;
            } catch (Exception ex) {
                Console.WriteLine($"SolveLinearSystem错误: {ex.Message}");
                throw;
            }
        }
        #endregion

        #region 辅助工具
        private static double Distance(Point a, Point b) {
            try {
                return Math.Sqrt(Math.Pow(a.X - b.X, 2) + Math.Pow(a.Y - b.Y, 2));
            } catch (Exception ex) {
                Console.WriteLine($"Distance计算错误: {ex.Message}");
                throw;
            }
        }

        private static bool ArePointsEqual(Point a, Point b, double tol) {
            try {
                return Math.Abs(a.X - b.X) < tol && Math.Abs(a.Y - b.Y) < tol;
            } catch (Exception ex) {
                Console.WriteLine($"ArePointsEqual错误: {ex.Message}");
                throw;
            }
        }

        private static bool IsClosed(IList<Point> pts, double tol) {
            try {
                return pts.Count > 2 && ArePointsEqual(pts.First(), pts.Last(), tol);
            } catch (Exception ex) {
                Console.WriteLine($"IsClosed判断错误: {ex.Message}");
                throw;
            }
        }
        #endregion
    }

    #region 数据结构
    /// <summary>样条曲线定义</summary>
    public sealed class Spline {
        public int Degree { get; set; }
        public List<Point> ControlPoints { get; set; } = new List<Point>();
        public List<double> Knots { get; set; } = new List<double>();
        public bool IsReversed { get; set; } = false;
        public bool IsClosed { get; set; } = false;
        public IEnumerable<Point> GetFitPoints() {
            try {
                return ControlPoints;
            } catch (Exception ex) {
                Console.WriteLine($"GetFitPoints错误: {ex.Message}");
                throw;
            }
        }
    }

    /// <summary>二维点坐标（record不可变类型）</summary>
    public record Point(double X, double Y) {
        public override string ToString() => $"({X:F3}, {Y:F3})";
    }
    #endregion

    class Program {
        static void Main(string[] args) {
            try {
                var spline1 = new Spline {
                    Degree = 3,
                    ControlPoints = new List<Point> { new(0, 0), new(1, 2), new(3, 1), new(4, 3) }
                };
                var spline2 = new Spline {
                    Degree = 2,
                    ControlPoints = new List<Point> { new(4, 3), new(5, 5), new(7, 4) }
                };
                var connected = ProfessionalSplineConnector.ConnectSplines(new[] { spline1, spline2 });
                Console.WriteLine("合并后的样条控制点:");
                foreach (var pt in connected.ControlPoints) Console.WriteLine(pt);
            } catch (Exception ex) {
                Console.WriteLine($"程序运行错误: {ex.Message}");
            }
        }
    }
}