[Matlab]椭球拟合+正圆修正+Qt转换

% 已知三组数据(三个坐标平面上的截线) 椭球拟合
x1 = [-267, -215, -140, -99, -65, -90, -110, -128, -164, -222, -258, -390, -456, -506, -542, -559, -549, -509, -462, -395, -335];
y1 = [-276, -251, -213, -147, -74, 11, 43, 64, 106, 137, 154, 141, 114, 67, 3, -34, -110, -177, -225, -259, -273];

y2 = [-265, -227, -150, -74, -11, 43, 110, 140, 152, 138, 115, 78, 56, -8, -52, -131, -182, -238, -263, -281, -274];
z2 = [197, 286, 346, 365, 358, 340, 260, 214, 165, 36, -16, -58, -84, -123, -116, -98, -68, -13, 47, 109, 158];

x3 = [-514, -528, -506, -484, -472, -425, -388, -365, -335, -296, -241, -151, -126, -114, -66, -72, -99, -135, -214, -325, -383, -469, -518];
z3 = [113, 170, 251, 272, 300, 339, 362, 372, 363, 377, 354, 344, 318, 288, 208, 71, 17, -43, -73, -92, -60, -9, 65];

fprintf(‘开始椭球拟合…\n’);

% 方法:使用数值优化直接拟合椭球参数
% 椭球方程: (x-x0)²/a² + (y-y0)²/b² + (z-z0)²/c² = 1

% 初始参数估计
x0_est = mean([x1, x3]);
y0_est = mean([y1, y2]);
z0_est = mean([z2, z3]);
a_est = std([x1, x3]);
b_est = std([y1, y2]);
c_est = std([z2, z3]);

initial_params = [x0_est, y0_est, z0_est, a_est, b_est, c_est];
fprintf(‘初始参数: 中心(%.1f,%.1f,%.1f), 半径(%.1f,%.1f,%.1f)\n’, initial_params);

% 使用优化方法
options = optimset(‘Display’, ‘iter’, ‘MaxFunEvals’, 10000, ‘TolX’, 1e-6);
optimal_params = fminsearch(@(p) ellipsoid_error(p, x1, y1, y2, z2, x3, z3), initial_params, options);

x0 = optimal_params(1);
y0 = optimal_params(2);
z0 = optimal_params(3);
a = optimal_params(4);
b = optimal_params(5);
c = optimal_params(6);

fprintf(‘\n=== 椭球拟合结果 ===\n’);
fprintf(‘椭球方程: (x-%.2f)²/%.2f² + (y-%.2f)²/%.2f² + (z-%.2f)²/%.2f² = 1\n’, x0, a, y0, b, z0, c);
fprintf(‘椭球中心: (%.2f, %.2f, %.2f)\n’, x0, y0, z0);
fprintf(‘半轴长度: a=%.2f, b=%.2f, c=%.2f\n’, a, b, c);

% 生成椭球面(确保实数)
[x_ell, y_ell, z_ell] = generate_simple_ellipsoid(x0, y0, z0, a, b, c, 30);

% 计算三个平面的截线椭圆
theta = linspace(0, 2*pi, 100);

% XY平面截线 (z=0)
xy_x = x0 + a * cos(theta);
xy_y = y0 + b * sin(theta);

% YZ平面截线 (x=0)
yz_y = y0 + b * cos(theta);
yz_z = z0 + c * sin(theta);

% XZ平面截线 (y=0)
xz_x = x0 + a * cos(theta);
xz_z = z0 + c * sin(theta);

% 可视化结果
figure(‘Position’, [100, 100, 1400, 1000]);

% 子图1: 三维椭球与数据点
subplot(2,3,1);
surf(x_ell, y_ell, z_ell, ‘FaceAlpha’, 0.2, ‘EdgeColor’, ‘none’, ‘FaceColor’, ‘cyan’);
hold on;

% 绘制数据点
scatter3(x1, y1, zeros(size(x1)), 60, ‘r’, ‘filled’, ‘MarkerEdgeColor’, ‘k’);
scatter3(zeros(size(y2)), y2, z2, 60, ‘g’, ‘filled’, ‘MarkerEdgeColor’, ‘k’);
scatter3(x3, zeros(size(x3)), z3, 60, ‘b’, ‘filled’, ‘MarkerEdgeColor’, ‘k’);

scatter3(x0, y0, z0, 120, ‘k’, ‘filled’, ‘Marker’, ‘x’, ‘LineWidth’, 3);
title(‘三维椭球与截线数据点’);
xlabel(‘X’); ylabel(‘Y’); zlabel(‘Z’);
legend(‘椭球’, ‘XY平面数据’, ‘YZ平面数据’, ‘XZ平面数据’, ‘椭球中心’);
grid on; axis equal;

% 子图2: XY平面截线
subplot(2,3,2);
plot(xy_x, xy_y, ‘r-‘, ‘LineWidth’, 3);
hold on;
scatter(x1, y1, 50, ‘r’, ‘filled’);
scatter(x0, y0, 80, ‘k’, ‘x’, ‘LineWidth’, 2);
title(‘XY平面截线 (z=0)’);
xlabel(‘X’); ylabel(‘Y’);
legend(‘拟合椭圆’, ‘数据点’, ‘中心’);
grid on; axis equal;

% 子图3: YZ平面截线
subplot(2,3,3);
plot(yz_y, yz_z, ‘g-‘, ‘LineWidth’, 3);
hold on;
scatter(y2, z2, 50, ‘g’, ‘filled’);
scatter(y0, z0, 80, ‘k’, ‘x’, ‘LineWidth’, 2);
title(‘YZ平面截线 (x=0)’);
xlabel(‘Y’); ylabel(‘Z’);
legend(‘拟合椭圆’, ‘数据点’, ‘中心’);
grid on; axis equal;

% 子图4: XZ平面截线
subplot(2,3,4);
plot(xz_x, xz_z, ‘b-‘, ‘LineWidth’, 3);
hold on;
scatter(x3, z3, 50, ‘b’, ‘filled’);
scatter(x0, z0, 80, ‘k’, ‘x’, ‘LineWidth’, 2);
title(‘XZ平面截线 (y=0)’);
xlabel(‘X’); ylabel(‘Z’);
legend(‘拟合椭圆’, ‘数据点’, ‘中心’);
grid on; axis equal;

% 计算拟合误差
fprintf(‘\n=== 拟合质量分析 ===\n’);

% XY平面误差
errors_xy = zeros(size(x1));
for i = 1:length(x1)
errors_xy(i) = ((x1(i)-x0)^2/a^2 + (y1(i)-y0)^2/b^2) – 1;
end

% YZ平面误差
errors_yz = zeros(size(y2));
for i = 1:length(y2)
errors_yz(i) = ((y2(i)-y0)^2/b^2 + (z2(i)-z0)^2/c^2) – 1;
end

% XZ平面误差
errors_xz = zeros(size(x3));
for i = 1:length(x3)
errors_xz(i) = ((x3(i)-x0)^2/a^2 + (z3(i)-z0)^2/c^2) – 1;
end

fprintf(‘XY平面 – 平均绝对误差: %.4f, 最大误差: %.4f\n’, mean(abs(errors_xy)), max(abs(errors_xy)));
fprintf(‘YZ平面 – 平均绝对误差: %.4f, 最大误差: %.4f\n’, mean(abs(errors_yz)), max(abs(errors_yz)));
fprintf(‘XZ平面 – 平均绝对误差: %.4f, 最大误差: %.4f\n’, mean(abs(errors_xz)), max(abs(errors_xz)));

% 子图5: 拟合误差分布
subplot(2,3,5);
all_errors = [errors_xy, errors_yz, errors_xz];
histogram(all_errors, 20, ‘FaceColor’, ‘m’, ‘FaceAlpha’, 0.7);
xlabel(‘拟合误差’); ylabel(‘频数’);
title(‘拟合误差分布’);
grid on;

% 子图6: 数据点与拟合椭圆距离
subplot(2,3,6);
plot(1:length(errors_xy), abs(errors_xy), ‘ro-‘, ‘LineWidth’, 1.5, ‘MarkerSize’, 4);
hold on;
plot(length(errors_xy)+1:length(errors_xy)+length(errors_yz), abs(errors_yz), ‘go-‘, ‘LineWidth’, 1.5, ‘MarkerSize’, 4);
plot(length(errors_xy)+length(errors_yz)+1:length(all_errors), abs(errors_xz), ‘bo-‘, ‘LineWidth’, 1.5, ‘MarkerSize’, 4);
xlabel(‘数据点索引’); ylabel(‘绝对误差’);
title(‘各数据点拟合误差’);
legend(‘XY平面’, ‘YZ平面’, ‘XZ平面’);
grid on;

% 显示拟合优度
fprintf(‘总体平均绝对误差: %.4f\n’, mean(abs(all_errors)));

% 辅助函数:椭球误差计算
function error = ellipsoid_error(params, x1, y1, y2, z2, x3, z3)
x0 = params(1); y0 = params(2); z0 = params(3);
a = params(4); b = params(5); c = params(6);

total_error = 0;

% XY平面数据误差 (z=0)
for i = 1:length(x1)
error_val = ((x1(i)-x0)^2/a^2 + (y1(i)-y0)^2/b^2) – 1;
total_error = total_error + error_val^2;
end

% YZ平面数据误差 (x=0)
for i = 1:length(y2)
error_val = ((y2(i)-y0)^2/b^2 + (z2(i)-z0)^2/c^2) – 1;
total_error = total_error + error_val^2;
end

% XZ平面数据误差 (y=0)
for i = 1:length(x3)
error_val = ((x3(i)-x0)^2/a^2 + (z3(i)-z0)^2/c^2) – 1;
total_error = total_error + error_val^2;
end

error = total_error;
end

% 辅助函数:生成简单椭球面
function [x, y, z] = generate_simple_ellipsoid(x0, y0, z0, a, b, c, n)
theta = linspace(0, 2*pi, n);
phi = linspace(0, pi, n);
[Theta, Phi] = meshgrid(theta, phi);

% 生成椭球坐标(确保实数)
x = x0 + a * sin(Phi) .* cos(Theta);
y = y0 + b * sin(Phi) .* sin(Theta);
z = z0 + c * cos(Phi);

% 确保所有坐标都是实数
x = real(x);
y = real(y);
z = real(z);
end

//正圆修正

% 原始椭球参数
x0_orig = -303.51;
y0_orig = -60.65;
z0_orig = 132.76;
a_orig = 237.59;
b_orig = 215.39;
c_orig = 241.03;

fprintf(‘=== 原始椭球参数 ===\n’);
fprintf(‘中心: (%.2f, %.2f, %.2f)\n’, x0_orig, y0_orig, z0_orig);
fprintf(‘半轴长度: a=%.2f, b=%.2f, c=%.2f\n’, a_orig, b_orig, c_orig);

% 步骤1: 将椭球平移到原点
% 新坐标 = 原始坐标 – 中心坐标
fprintf(‘\n=== 坐标平移 ===\n’);
fprintf(‘平移向量: (%.2f, %.2f, %.2f)\n’, -x0_orig, -y0_orig, -z0_orig);

% 步骤2: 计算平均半径并生成正圆
r_avg = (a_orig + b_orig + c_orig) / 3;
fprintf(‘\n=== 半径归一化 ===\n’);
fprintf(‘原始平均半径: %.2f\n’, r_avg);
fprintf(‘目标半径: %.2f (正圆)\n’, r_avg);

% 生成原始椭球(用于对比)
[x_orig, y_orig, z_orig] = generate_ellipsoid(x0_orig, y0_orig, z0_orig, a_orig, b_orig, c_orig, 40);

% 生成平移后的椭球(中心在原点)
[x_centered, y_centered, z_centered] = generate_ellipsoid(0, 0, 0, a_orig, b_orig, c_orig, 40);

% 生成正圆球体
[x_sphere, y_sphere, z_sphere] = generate_sphere(0, 0, 0, r_avg, 40);

% 可视化结果
figure(‘Position’, [100, 100, 1500, 500]);

% 子图1: 原始椭球
subplot(1,3,1);
surf(x_orig, y_orig, z_orig, ‘FaceAlpha’, 0.3, ‘EdgeColor’, ‘none’, ‘FaceColor’, ‘blue’);
hold on;
scatter3(x0_orig, y0_orig, z0_orig, 100, ‘r’, ‘filled’, ‘Marker’, ‘x’, ‘LineWidth’, 3);
scatter3(0, 0, 0, 100, ‘k’, ‘filled’, ‘Marker’, ‘o’);
title(‘原始椭球’);
xlabel(‘X’); ylabel(‘Y’); zlabel(‘Z’);
legend(‘原始椭球’, ‘椭球中心’, ‘坐标原点’);
grid on; axis equal;

% 子图2: 平移后的椭球(中心在原点)
subplot(1,3,2);
surf(x_centered, y_centered, z_centered, ‘FaceAlpha’, 0.3, ‘EdgeColor’, ‘none’, ‘FaceColor’, ‘green’);
hold on;
scatter3(0, 0, 0, 100, ‘r’, ‘filled’, ‘Marker’, ‘x’, ‘LineWidth’, 3);
title(‘平移后椭球(中心在原点)’);
xlabel(‘X’); ylabel(‘Y’); zlabel(‘Z’);
legend(‘平移椭球’, ‘中心(0,0,0)’);
grid on; axis equal;

% 子图3: 正圆球体
subplot(1,3,3);
surf(x_sphere, y_sphere, z_sphere, ‘FaceAlpha’, 0.3, ‘EdgeColor’, ‘none’, ‘FaceColor’, ‘red’);
hold on;
scatter3(0, 0, 0, 100, ‘r’, ‘filled’, ‘Marker’, ‘x’, ‘LineWidth’, 3);
title(sprintf(‘正圆球体 (半径=%.2f)’, r_avg));
xlabel(‘X’); ylabel(‘Y’); zlabel(‘Z’);
legend(‘正圆球体’, ‘中心(0,0,0)’);
grid on; axis equal;

% 坐标变换函数
fprintf(‘\n=== 坐标变换公式 ===\n’);
fprintf(‘1. 坐标平移:\n’);
fprintf(‘ x_new = x_original – (%.2f)\n’, x0_orig);
fprintf(‘ y_new = y_original – (%.2f)\n’, y0_orig);
fprintf(‘ z_new = z_original – (%.2f)\n’, z0_orig);

fprintf(‘\n2. 半径归一化(生成正圆):\n’);
fprintf(‘ 目标半径 R = %.2f\n’, r_avg);
fprintf(‘ 对于任意点 (x, y, z),正圆坐标为:\n’);
fprintf(‘ magnitude = sqrt(x² + y² + z²)\n’);
fprintf(‘ scale = %.2f / magnitude\n’, r_avg);
fprintf(‘ x_sphere = x * scale\n’);
fprintf(‘ y_sphere = y * scale\n’);
fprintf(‘ z_sphere = z * scale\n’);

% 验证变换
fprintf(‘\n=== 变换验证 ===\n’);
test_points = [x0_orig, y0_orig, z0_orig; % 原始中心点
x0_orig + a_orig, y0_orig, z0_orig; % X轴端点
x0_orig, y0_orig + b_orig, z0_orig; % Y轴端点
x0_orig, y0_orig, z0_orig + c_orig]; % Z轴端点

for i = 1:size(test_points, 1)
x_orig = test_points(i, 1);
y_orig = test_points(i, 2);
z_orig = test_points(i, 3);

% 坐标平移
x_centered = x_orig – x0_orig;
y_centered = y_orig – y0_orig;
z_centered = z_orig – z0_orig;

% 半径归一化
magnitude = sqrt(x_centered^2 + y_centered^2 + z_centered^2);
scale = r_avg / magnitude;
x_sphere = x_centered * scale;
y_sphere = y_centered * scale;
z_sphere = z_centered * scale;
sphere_radius = sqrt(x_sphere^2 + y_sphere^2 + z_sphere^2);

fprintf(‘点%d:\n’, i);
fprintf(‘ 原始: (%.1f, %.1f, %.1f)\n’, x_orig, y_orig, z_orig);
fprintf(‘ 平移: (%.1f, %.1f, %.1f), 距离=%.1f\n’, x_centered, y_centered, z_centered, magnitude);
fprintf(‘ 正圆: (%.1f, %.1f, %.1f), 半径=%.1f\n’, x_sphere, y_sphere, z_sphere, sphere_radius);
end

% 完整的坐标变换函数
fprintf(‘\n=== 完整变换函数 ===\n’);
fprintf(‘function [x_sphere, y_sphere, z_sphere] = transformToSphere(x, y, z)\n’);
fprintf(‘ %% 输入: 原始坐标\n’);
fprintf(‘ %% 输出: 正圆球体上的坐标\n’);
fprintf(‘ \n’);
fprintf(‘ %% 椭球参数\n’);
fprintf(‘ x0 = %.4f;\n’, x0_orig);
fprintf(‘ y0 = %.4f;\n’, y0_orig);
fprintf(‘ z0 = %.4f;\n’, z0_orig);
fprintf(‘ R = %.4f;\n’, r_avg);
fprintf(‘ \n’);
fprintf(‘ %% 坐标平移\n’);
fprintf(‘ x_centered = x – x0;\n’);
fprintf(‘ y_centered = y – y0;\n’);
fprintf(‘ z_centered = z – z0;\n’);
fprintf(‘ \n’);
fprintf(‘ %% 半径归一化\n’);
fprintf(‘ magnitude = sqrt(x_centered^2 + y_centered^2 + z_centered^2);\n’);
fprintf(‘ if magnitude > 0\n’);
fprintf(‘ scale = R / magnitude;\n’);
fprintf(‘ x_sphere = x_centered * scale;\n’);
fprintf(‘ y_sphere = y_centered * scale;\n’);
fprintf(‘ z_sphere = z_centered * scale;\n’);
fprintf(‘ else\n’);
fprintf(‘ x_sphere = R; y_sphere = 0; z_sphere = 0;\n’);
fprintf(‘ end\n’);
fprintf(‘end\n’);

% 辅助函数:生成椭球面
function [x, y, z] = generate_ellipsoid(x0, y0, z0, a, b, c, n)
theta = linspace(0, 2*pi, n);
phi = linspace(0, pi, n);
[Theta, Phi] = meshgrid(theta, phi);

x = x0 + a * sin(Phi) .* cos(Theta);
y = y0 + b * sin(Phi) .* sin(Theta);
z = z0 + c * cos(Phi);
end

% 辅助函数:生成球体
function [x, y, z] = generate_sphere(x0, y0, z0, r, n)
theta = linspace(0, 2*pi, n);
phi = linspace(0, pi, n);
[Theta, Phi] = meshgrid(theta, phi);

x = x0 + r * sin(Phi) .* cos(Theta);
y = y0 + r * sin(Phi) .* sin(Theta);
z = z0 + r * cos(Phi);
end

//Qt转换

QVector3D MainWindow::transformToSphere(double x, double y, double z)
{
// 椭球参数
const double centerX = -303.51;
const double centerY = -60.65;
const double centerZ = 132.76;
const double radius = 231.34; // 平均半径

// 坐标平移
double x_centered = x – centerX;
double y_centered = y – centerY;
double z_centered = z – centerZ;

// 半径归一化
double magnitude = std::sqrt(x_centered * x_centered +
y_centered * y_centered +
z_centered * z_centered);

if (magnitude > 1e-10) {
double scale = radius / magnitude;
return QVector3D(x_centered * scale,
y_centered * scale,
z_centered * scale);
} else {
return QVector3D(radius, 0, 0);
}
}

发表评论

邮箱地址不会被公开。 必填项已用*标注