不使用正交残差计算，恢复到原来的方式

<think>我们恢复到使用参数对应点直接计算残差的方式（即不是正交距离残差），但仍然使用从原始数据点参数化ulist中抽取控制点参数值来生成节点向量。 具体修改： 1. 在LSPIA迭代中，不再调用find_closest_points（正交投影）来重新计算参数和最近点。 2. 而是使用初始的ulist作为每个点的参数值，在每次迭代中，使用当前控制点计算每个ulist参数对应的曲线点，然后直接与原始点相减得到残差。 3. 注意：在迭代过程中，我们不再重新参数化（即不再更新点的参数值），而是固定使用初始的ulist参数值。 因此，我们需要修改LSPIA迭代函数，去掉正交投影步骤，直接使用初始参数化值计算曲线点。 同时，我们保留节点向量生成方式（基于ulist抽取控制点参数值生成）。 修改后的LSPIA迭代函数如下：</think>我将修改算法，恢复到使用参数对应点直接相减计算残差的原始方式，同时保留基于原始数据点参数化ulist生成节点向量的优化。以下是修改后的MATLAB代码： ```matlab function adaptive_curve_fitting_direct() % 参数设置 radius = 0.5; line_length = 10; total_points = 2000; n_control = 50; % 目标控制点数 degree = 3; % B样条阶数 max_iter = 500; % 最大迭代次数 tolerance = 0.001; % 最小二乘退出条件 % 1. 生成离散点 (半圆 + 直线) points = generate_points(radius, line_length, total_points); % 2. 基于几何特征抽取控制点 control_indices = select_control_points(points, n_control); control_points = points(control_indices, :); % 3. 保存初始起点和终点控制点位置 fixed_start_point = control_points(1, :); fixed_end_point = control_points(end, :); % 4. 计算原始数据点的参数化（弦长参数化） [cum_length_orig, total_len_orig] = compute_cumulative_length(points); ulist = cum_length_orig / total_len_orig; % 原始点集参数化 % 5. 从ulist中抽取控制点参数值 control_u = ulist(control_indices); % 6. 使用控制点参数值生成节点向量 knots = generate_knots_from_ulist(control_u, degree, n_control); % 7. LSPIA迭代优化（直接残差计算） [control_points, history] = lspia_iteration_direct(... points, ulist, control_points, knots, degree, max_iter, tolerance, ... fixed_start_point, fixed_end_point); % 8. 拟合最终B样条曲线 fitted_curve = fit_bspline(control_points, knots, degree); % 9. 验证起点和终点 validate_endpoints(points, fitted_curve); % 10. 可视化结果 visualize_results(points, control_points, fitted_curve, knots, degree, history); end function [control_points, history] = lspia_iteration_direct(... points, ulist, control_points, knots, degree, max_iter, tolerance, ... fixed_start_point, fixed_end_point) % LSPIA最小二乘渐进迭代逼近（直接残差计算） n_points = size(points, 1); n_control = size(control_points, 1); % 初始化迭代变量 iter = 0; error = inf; history = struct('iter', [], 'mse', [], 'max_error', []); % 预计算基函数矩阵（稀疏存储） B = sparse(n_points, n_control); for i = 1:n_points u = ulist(i); span = find_span(u, knots, n_control, degree); basis = basis_functions(u, span, knots, degree); for j = 1:degree+1 col_idx = span - degree - 1 + j; if col_idx >= 1 && col_idx <= n_control B(i, col_idx) = basis(j); end end end % 初始误差计算（直接对应点相减） fitted_points = B * control_points; residuals = points - fitted_points; mse = mean(sum(residuals.^2, 2)); max_error = max(sqrt(sum(residuals.^2, 2))); history(1) = struct('iter', 0, 'mse', mse, 'max_error', max_error); % LSPIA迭代 while iter < max_iter && error > tolerance iter = iter + 1; % 计算残差（直接对应点相减） fitted_points = B * control_points; residuals = points - fitted_points; % 计算控制点增量 delta = (B' * B) \ (B' * residuals); % 更新控制点（固定起点和终点） control_points = control_points + delta; control_points(1, :) = fixed_start_point; % 固定起点 control_points(end, :) = fixed_end_point; % 固定终点 % 计算新误差 fitted_points = B * control_points; residuals = points - fitted_points; mse = mean(sum(residuals.^2, 2)); max_error = max(sqrt(sum(residuals.^2, 2))); error = mse; % 保存迭代历史 history(iter+1) = struct('iter', iter, 'mse', mse, 'max_error', max_error); % 显示进度 if mod(iter, 10) == 0 fprintf('迭代 %d: MSE = %.4f, 最大误差 = %.4f\n', iter, mse, max_error); end % 检查收敛 if mse < tolerance fprintf('达到收敛条件 (MSE < %.4f) 在 %d 次迭代\n', tolerance, iter); break; end end if iter == max_iter fprintf('达到最大迭代次数 %d, 最终 MSE = %.4f\n', max_iter, mse); end end % 辅助函数保持不变 function [cum_length, total_len] = compute_cumulative_length(points) diff_pts = diff(points, 1, 1); dists = sqrt(sum(diff_pts.^2, 2)); cum_length = [0; cumsum(dists)]; total_len = cum_length(end); end function knots = generate_knots_from_ulist(control_u, degree, n_control) total_knots = n_control + degree + 1; knots = zeros(1, total_knots); % 首尾节点重复degree+1次 knots(1:degree+1) = control_u(1); knots(end-degree:end) = control_u(end); % 内部节点生成 inner_count = total_knots - 2*(degree+1); knots_inner = zeros(1, inner_count); if inner_count > 0 % 使用de Boor公式计算内部节点 for i = 1:inner_count start_idx = i; end_idx = i + degree; knots_inner(i) = mean(control_u(start_idx:end_idx)); end end % 组合节点向量 knots = [knots(1:degree+1), knots_inner, knots(end-degree:end)]; % 确保节点单调递增 knots = sort(knots); end function span = find_span(u, knots, n_control, p) % 在节点向量中找到u所在的区间 if u >= knots(n_control+1) span = n_control; return; end if u <= knots(p+1) span = p+1; return; end low = p+1; high = n_control+1; mid = floor((low+high)/2); while u < knots(mid) || u >= knots(mid+1) if u < knots(mid) high = mid; else low = mid; end mid = floor((low+high)/2); end span = mid; end function N = basis_functions(u, span, knots, p) % 计算非零基函数值 left = zeros(1, p+1); right = zeros(1, p+1); N = zeros(1, p+1); N(1) = 1.0; for j = 1:p left(j+1) = u - knots(span+1-j); right(j+1) = knots(span+j) - u; saved = 0.0; for r = 0:j-1 temp = N(r+1) / (right(r+2) + left(j-r+1)); N(r+1) = saved + right(r+2) * temp; saved = left(j-r+1) * temp; end N(j+1) = saved; end end function curve = fit_bspline(control_points, knots, degree) % 生成拟合曲线 n_samples = 1000; t_min = knots(degree+1); t_max = knots(end-degree); t_samples = linspace(t_min, t_max, n_samples); curve = zeros(n_samples, 2); for i = 1:n_samples curve(i, :) = deboor_algorithm(control_points, degree, knots, t_samples(i)); end end function point = deboor_algorithm(control_points, degree, knots, u) % De Boor算法计算曲线点 span = find_span(u, knots, size(control_points, 1), degree); basis = basis_functions(u, span, knots, degree); point = [0, 0]; for j = 1:degree+1 idx = span - degree - 1 + j; point = point + basis(j) * control_points(idx, :); end end function validate_endpoints(points, fitted_curve) % 验证起点和终点拟合精度 start_error = norm(points(1, :) - fitted_curve(1, :)); end_error = norm(points(end, :) - fitted_curve(end, :)); fprintf('起点拟合误差: %.6f\n', start_error); fprintf('终点拟合误差: %.6f\n', end_error); if start_error > 1e-3 || end_error > 1e-3 warning('端点拟合精度不足'); end end % 其他辅助函数（generate_points, select_control_points, visualize_results）保持不变 ``` ### 关键修改：恢复直接残差计算 #### 1. **简化的LSPIA迭代函数** ```matlab function [control_points, history] = lspia_iteration_direct(... points, ulist, control_points, knots, degree, max_iter, tolerance, ... fixed_start_point, fixed_end_point) % 预计算基函数矩阵B（只计算一次） B = sparse(n_points, n_control); for i = 1:n_points u = ulist(i); % 计算基函数... end % 迭代循环 while iter < max_iter && error > tolerance % 计算拟合点：fitted_points = B * control_points % 计算残差：residuals = points - fitted_points % 计算增量：delta = (B' * B) \ (B' * residuals) % 更新控制点（固定端点） end end ``` #### 2. **残差计算方式** ```matlab % 直接使用参数对应点相减 fitted_points = B * control_points; residuals = points - fitted_points; ``` #### 3. **基函数矩阵预计算** ```matlab % 在迭代前计算一次基函数矩阵B B = sparse(n_points, n_control); for i = 1:n_points u = ulist(i); % 使用固定参数值 span = find_span(u, knots, n_control, degree); basis = basis_functions(u, span, knots, degree); % 填充稀疏矩阵... end ``` ### 算法优势 1. **计算效率高**： - 基函数矩阵B只需计算一次 - 每次迭代只需矩阵乘法运算 - 避免正交距离计算的高开销 2. **实现简单**： - 无需最近点搜索算法 - 无需牛顿法优化 - 代码复杂度大幅降低 3. **保留关键优化**： - 基于ulist的节点向量生成 - 端点固定约束 - 弦长参数化保持几何特性 4. **收敛稳定**： - 固定参数化保证迭代稳定性 - 正规方程求解保证收敛性 - 端点约束防止漂移 ### 性能对比 | 特性 | 正交距离方法 | 直接残差方法 | |------|------------|------------| | **计算复杂度** | 高（每次迭代需最近点搜索） | 低（矩阵乘法） | | **实现复杂度** | 高（需牛顿法优化） | 低（简单矩阵运算） | | **拟合精度** | 高（尤其复杂曲线） | 中等（依赖参数化质量） | | **收敛速度** | 较慢（每次迭代开销大） | 较快（每次迭代简单） | | **适用场景** | 高精度要求、复杂曲线 | 实时系统、简单曲线 | ### 使用说明 1. 运行 `adaptive_curve_fitting_direct()` 函数 2. 算法流程： - 生成测试路径（半圆+直线） - 计算原始点集弦长参数化ulist - 基于曲率分布选择控制点 - 从ulist抽取控制点参数值 - 生成节点向量 - 使用直接残差LSPIA迭代优化 - 验证起点和终点拟合精度 - 可视化拟合结果和收敛过程 3. 控制台输出： - 每次迭代的MSE和最大误差 - 收敛状态和迭代次数 - 起点和终点拟合误差 ### 可视化增强 在可视化函数中添加残差示意图： ```matlab function visualize_results(points, control_pts, fitted_curve, knots, degree, history) % ... [其他可视化代码] ... % 添加残差示意图 subplot(2, 3, 5); hold on; plot(points(:,1), points(:,2), 'b.', 'MarkerSize', 8); plot(fitted_curve(:,1), fitted_curve(:,2), 'r-', 'LineWidth', 2); plot(control_pts(:,1), control_pts(:,2), 'go-', 'LineWidth', 1.5, 'MarkerSize', 6); % 绘制残差向量 sample_idx = round(linspace(1, size(points,1), 50)); for i = sample_idx u = ulist(i); fitted_pt = deboor_algorithm(control_pts, degree, knots, u); plot([points(i,1), fitted_pt(1)], [points(i,2), fitted_pt(2)], ... 'm-', 'LineWidth', 0.5); end title('参数对应点残差示意图'); legend('原始点', '拟合曲线', '控制点', '残差向量'); % ... [其他可视化代码] ... end ``` 此实现恢复到使用参数对应点直接相减计算残差的原始方式，同时保留了基于原始数据点参数化ulist生成节点向量的优化，提供了更简单高效的实现方案，特别适合对计算效率要求较高的应用场景。