MATLAB中，曲面点云法向量可视化，各点法向量方向不一致，如何修改代码使得各点法向量朝向统一？

Question

问题如下图：

希望效果图如下：

或者这样的，只要朝向统一就好：

点云坐标p.txt存在百度云中：链接：https://pan.baidu.com/s/13DFDJi7184CjJvigpj1BNQ
提取码：nhu4
代码如下：

load p.txt

%MYNORMALVECTOR 采用PCA算法求取点云矩阵的每个点的法向量。
%Input: 
%      p是Nx3的矩阵，包含点云的三维坐标信息[X,Y,Z]。
%Output:
%      n是Nx3的矩阵，对应p中每一点的单位法向量。

N = size(p,1);
n = zeros(N,3); % 存放法线的矩阵
% 寻找每个点K近邻的索引
k = 5; % k近邻
[NeiborIndex,~] = knnsearch(p,p,'Distance','euclidean','NSMethod','kdtree','K',k);
for i=1:N
    % 当前点邻域点集的索引
    IndexOfCurrentNeibors = NeiborIndex(i,:);
    % 当前点邻域点集
    CurrentNeighbors = p(IndexOfCurrentNeibors,:);
    PtCent = mean(CurrentNeighbors); %得到局部点云平均值，便于计算法向量长度和方向    
    % PCA法计算法向量
    B = CurrentNeighbors - repmat(PtCent,k,1); % 去中心化点集矩阵
    PtCov = transpose(B)*B; %协方差矩阵
    % 求特征值、特征向量。 
    % D是对应的特征值对角矩阵，V是特征向量
    % (因为协方差矩阵为实对称矩阵，故特征向量为单位正交向量)
    [V,D] = eig(PtCov); 
    % d0为最小特征值  idx为特征值的列数索引。
    [~, idx] = min(diag(D)); % diag():创建对角矩阵或获取矩阵的对角元素
    n(i,:) = V(:,idx)';   % 最小特征值对应的特征向量为法矢，即单位法向量  
    % 朝向统一
    n1 = n(i,:);
    n2 = p(i,:);
    if dot(n1,n2)<0 % 如果这个向量与法向量的数量积为负数（反向）
        n(i,:) = -n(i,:);  % 法向量取反向
    end
end
%法向量可视化
ptCloud = pointCloud(p);
figure()
pcshow(ptCloud);
hold on 
quiver3(p(:,1),p(:,2),p(:,3),n(:,1),n(:,2),n(:,3));
hold off 
请赐教

Answer 1

你把n1的方向换成y的正方向[0,1,0]，然后n2换成n(i,:)，就可以啦，最后的view(2)是为了显示xOy面设置的，可以看见箭头朝向确实一致了

load p.txt
%MYNORMALVECTOR 采用PCA算法求取点云矩阵的每个点的法向量。
%Input: 
%      p是Nx3的矩阵，包含点云的三维坐标信息[X,Y,Z]。
%Output:
%      n是Nx3的矩阵，对应p中每一点的单位法向量。
N = size(p,1);
n = zeros(N,3); % 存放法线的矩阵
% 寻找每个点K近邻的索引
k = 5; % k近邻
[NeiborIndex,~] = knnsearch(p,p,'Distance','euclidean','NSMethod','kdtree','K',k);
for i=1:N
    % 当前点邻域点集的索引
    IndexOfCurrentNeibors = NeiborIndex(i,:);
    % 当前点邻域点集
    CurrentNeighbors = p(IndexOfCurrentNeibors,:);
    PtCent = mean(CurrentNeighbors); %得到局部点云平均值，便于计算法向量长度和方向    
    % PCA法计算法向量
    B = CurrentNeighbors - repmat(PtCent,k,1); % 去中心化点集矩阵
    PtCov = transpose(B)*B; %协方差矩阵
    % 求特征值、特征向量。 
    % D是对应的特征值对角矩阵，V是特征向量
    % (因为协方差矩阵为实对称矩阵，故特征向量为单位正交向量)
    [V,D] = eig(PtCov); 
    % d0为最小特征值  idx为特征值的列数索引。
    [~, idx] = min(diag(D)); % diag():创建对角矩阵或获取矩阵的对角元素
    n(i,:) = V(:,idx)';   % 最小特征值对应的特征向量为法矢，即单位法向量  
    % 朝向统一
    n1 = [0,1,0]; %% n1用的是y轴正方向
    n2 = n(i,:); %% n2用的是n的i个方向
    if dot(n1,n2)<0 % 如果这个向量与法向量的数量积为负数（反向）
        n(i,:) = -n(i,:);  % 法向量取反向
    end
end
%法向量可视化
ptCloud = pointCloud(p);
figure()
pcshow(ptCloud);
hold on 
quiver3(p(:,1),p(:,2),p(:,3),n(:,1),n(:,2),n(:,3));
hold off 
xlabel('x'); ylabel('y');zlabel('z')
view(2)