今天遇到了球面拟合的问题,正好在上coursera上的机器学习的课程,干脆就用梯度下降拟合一个球面试试吧
%% clear
clear;
clc;
close all;
%% generate data
theta = 0:2*pi/100:2*pi;
beta = pi/6:2*pi/100:5*pi/6;
r = 1;
x = r*cos(beta)'*cos(theta) + 1;
y = r*cos(beta)'*sin(theta) + 2;
z = r*sin(beta)'*ones(1,length(theta)) + 3;
%% add noise
n = randn(length(beta),length(theta));
n = n/max(max(n)) * 0.05;
z = z+n;
figure
mesh(x,y,z)
%% construct cost function
% J(a,b,c,R) = 1/(2*length(theta)*length(beta)) * Sigma((x-a)^2+(y-b)^2+(z-c)^2-R^2)^2
a = x(1);
b = y(1);
c = z(1);
R = 1.1;
alpha = 0.35;
N = 100;
J = zeros(1,N);
x = reshape(x, length(beta)*length(theta), 1);
y = reshape(y, length(beta)*length(theta), 1);
z = reshape(z, length(beta)*length(theta), 1);
tic
% iterate
for i=1:N
deltaA = alpha * 1/(length(beta)*length(theta)) * sum(((x-a).^2+(y-b).^2+(z-c).^2-R^2).*(x-a)*-2);
deltaB = alpha * 1/(length(beta)*length(theta)) * sum(((x-a).^2+(y-b).^2+(z-c).^2-R^2).*(y-b)*-2);
deltaC = alpha * 1/(length(beta)*length(theta)) * sum(((x-a).^2+(y-b).^2+(z-c).^2-R^2).*(z-c)*-2);
deltaR = alpha * 1/(length(beta)*length(theta)) * sum(((x-a).^2+(y-b).^2+(z-c).^2-R^2)*R*-2);
a = a-deltaA;
b = b-deltaB;
c = c-deltaC;
R = R - deltaR;
J(i) = 1/(2*length(beta)*length(theta)) * sum(((x-a).^2 + (y-b).^2 + (z-c).^2-R^2).^2);
end
toc
figure
plot(J)
效果也还可以,100次迭代耗时0.044s,圆心x,y轴坐标误差很小,轴坐标差0.04, 半径差0.04。