圆柱体

1.绘制圆柱体，计算圆柱体表面某一点的法向量。

r=5;
h=10;
n=100;%分点数
radius = 1.0;%搜索半径
min_neighbors = 8;
%-------------------------------
n_h= floor(h/(r*2*pi/n));
X=zeros(n,n_h);
Y=zeros(n,n_h);
Z=zeros(n,n_h);
for j=1:n_h
for i=1:n
    X(i,j)=r*cos(i*2*pi/n);
    Y(i,j)=r*sin(i*2*pi/n);
    Z(i,j)=j*r*2*pi/n;
end
end

data=zeros(n*n_h,3);
data(:,1)=reshape(X,[n*n_h,1]);
data(:,2)=reshape(Y,[n*n_h,1]);
data(:,3)=reshape(Z,[n*n_h,1]);
tree=KDTreeSearcher(data);
query=[X(50,20),Y(50,20),Z(50,20)];
hold on
plot3(X,Y,Z,'r.','MarkerSize',1)
plot3(query(:,1),query(:,2),query(:,3),'b.','MarkerSize',4)
idx = rangesearch(tree, query, radius);                                                  
idxs = idx{1};
neighbors = [data(idxs(:),1) data(idxs(:),2) data(idxs(:),3)];
if size(neighbors) < min_neighbors
    normal = {1};
    return;
end
C = cov(neighbors);% 计算协方差矩阵
[v, lambda] = eig(C);
[~, i_min] = min(diag(lambda));
[~, i_max] = max(diag(lambda));
normal = v(:,i_min) ./ norm(v(:,i_min));
normal_max = v(:,i_max) ./ norm(v(:,i_max));
query1=[normal(1,1)+ query(1,1),normal(2,1)+ query(1,2),normal(3,1)+ query(1,3)];
query2=[normal_max(1,1)+ query(1,1),normal_max(2,1)+ query(1,2),normal_max(3,1)+ query(1,3)];
plot3([query1(1,1) query(1,1)],[query1(1,2) query(1,2)],[query1(1,3) query(1,3)],'LineWidth',1);
plot3([query2(1,1) query(1,1)],[query2(1,2) query(1,2)],[query2(1,3) query(1,3)],'LineWidth',1);


hold off
axis equal