平面度误差计算（matlab）：最小二乘法、对角线法、最小区域法

function [flatness_error, best_planes] = calculate_flatness_error(points) % 基于最小包容区域法计算平面度误差 flatness_error = inf; best_planes = []; num_points = size(points, 1); % 第一步：使用三角形准则 for i = 1:num_points-2 for j = i+1:num_points-1 for k = j+1:num_points [normal, D, current_error] = triangle_criterion(points, i, j, k); if current_error < flatness_error projection_inside = check_projection_inside_triangle(points, normal, D, i, j, k); if projection_inside flatness_error = current_error; best_planes = [{normal, D};{normal, D + current_error}]; end end end end end % 第二步：如果极点不在任何三角形内，则使用线交叉准则 if isempty(best_planes) for i = 1:num_points-3 for j = i+1:num_points-2 for k = j+1:num_points-1 for l = k+1:num_points [current_error, plane1, plane2] = line_cross_criterion(points, i, j, k, l); if current_error < flatness_error flatness_error = current_error; best_planes = [plane1; plane2]; end end end end end end % 绘制最佳包容平面和点 plot_results(points, best_planes); % 输出平面度误差 fprintf('平面度误差: %f\n', flatness_error); end function [normal, D, current_error] = triangle_criterion(points, i, j, k) % 三角形准则计算平面度误差 p1 = points(i, :); p2 = points(j, :); p3 = points(k, :); normal = cross(p2 - p1, p3 - p1); if norm(normal) < eps current_error = inf; D = []; normal = []; return; end normal = normal / norm(normal); D = -dot(normal, p1); distances = points * normal' + D; current_error = max(abs(distances)); end function projection_inside = check_projection_inside_triangle(points, normal, D, i, j, k) % 检查极点投影是否在三角形内部 distances = points * normal' + D; [~, idx] = max(abs(distances)); extreme_point = points(idx, :); distances(idx) = 0; % Remove the extreme point from the distances % Determine if the extreme point projects inside the triangle triangle_points = [points(i, :); points(j, :); points(k, :)]; extreme_proj = extreme_point - dot(extreme_point, normal) * normal; projection_inside = point_in_triangle(extreme_proj, triangle_points); end function inside = point_in_triangle(pt, triangle_pts) % 使用重心坐标来检查点是否在三角形内部 v0 = triangle_pts(3, :) - triangle_pts(1, :); v1 = triangle_pts(2, :) - triangle_pts(1, :); v2 = pt - triangle_pts(1, :); dot00 = dot(v0, v0); dot01 = dot(v0, v1); dot02 = dot(v0, v2); dot11 = dot(v1, v1); dot12 = dot(v1, v2); invDenom = 1 / (dot00 * dot11 - dot01 * dot01); u = (dot11 * dot02 - dot01 * dot12) * invDenom; v = (dot00 * dot12 - dot01 * dot02) * invDenom; inside = (u >= 0) && (v >= 0) && (u + v <= 1); end function [current_error, plane1, plane2] = line_cross_criterion(points, i, j, k, l) % 线交叉准则计算平面度误差 p1 = points(i, :); p2 = points(j, :); p3 = points(k, :); p4 = points(l, :); line1 = p2 - p1; line2 = p4 - p3; normal = cross(line1, line2); if norm(normal) < eps current_error = inf; plane1 = []; plane2 = []; return; end normal = normal / norm(normal); D = -dot(normal, p1); distances = points * normal' + D; current_error = max(abs(distances)); plane1 = {normal, D}; plane2 = {normal, D + current_error}; end function plot_results(points, planes) % 绘制点和最佳包容平面 scatter3(points(:,1), points(:,2), points(:,3), 'filled'); hold on; for i = 1:length(planes) plot_plane(planes{i,1}, planes{i,2}, points); end hold off; xlabel('X'); ylabel('Y'); zlabel('Z'); title('包容平面与平面度误差'); grid on; end function plot_plane(normal, D, points) % 绘制一个平面 [x, y] = meshgrid(linspace(min(points(:,1)), max(points(:,1)), 10), ... linspace(min(points(:,3)), max(points(:,3)), 10)); z = (-D - normal(1) * x - normal(2) * y) / normal(3); surf(x, y, z, 'FaceAlpha', 0.5, 'EdgeColor', 'none'); end