%求三角形法向量
function [x,y,z]=fa_vector(nodes,x1,x2,x3)
%三个点的坐标
p1=[nodes(x1,2),nodes(x1,3),nodes(x1,4)];
p2=[nodes(x2,2),nodes(x2,3),nodes(x1,4)];
p3=[nodes(x3,2),nodes(x3,3),nodes(x3,4)];
%两个边向量
a=p2-p1;
b=p3-p1;
%求法向量
c=cross(b,a);
%归一化
 norm = sqrt(c(1,1)^2+c(1,2)^2+c(1,3)^2);
 x=c(1,1)/norm;
 y=c(1,2)/norm;
 z=c(1,3)/norm;
end

思路还是很容易想到的：

1.首先使用KD树寻找当前点邻域的N个点，这里取了10个，直接调用了vlfeat。

2.用最小二乘估计当前邻域点组成的平面，得到法向量。

3.根据当前邻域点平均值确定邻域质心，通常质心会在弯曲表面的内部，反方向即为法线方向。

vlfeat在这里下载，配置参考这里，rabbit.pcd下载地址

处理效果如下：

原始点云：



点云表面法向量，做了降采样处理：



兔子果断变刺猬。

matlab代码如下：



clear all;
close all;
clc;
warning off;

pc = pcread('rabbit.pcd');
pc=pcdownsample(pc,'random',0.3);       %0.3倍降采样
pcshow(pc);

pc_point = pc.Location';                %得到点云数据
kdtree = vl_kdtreebuild(pc_point);      %使用vlfeat建立kdtree

normE=[];
for i=1:length(pc_point)
    
    p_cur = pc_point(:,i);
    [index, distance] = vl_kdtreequery(kdtree, pc_point, p_cur, 'NumNeighbors', 10);    %寻找当前点最近的10个点
    p_neighbour = pc_point(:,index)';
    p_cent = mean(p_neighbour);     %得到局部点云平均值，便于计算法向量长度和方向
    
    %最小二乘估计平面
    X=p_neighbour(:,1);
    Y=p_neighbour(:,2);
    Z=p_neighbour(:,3);
    XX=[X Y ones(length(index),1)];
    YY=Z;
    %得到平面法向量
    C=(XX'*XX)\XX'*YY;
    
    %局部平面指向局部质心的向量
    dir1 = p_cent-p_cur';
    %局部平面法向量
    dir2=[C(1) C(2) -1];
    
    %计算两个向量的夹角
    ang = sum(dir1.*dir2) / (sqrt(dir1(1)^2 +dir2(1)^2) + sqrt(dir1(2)^2 +dir2(2)^2)+sqrt(dir1(3)^2 +dir2(3)^2) );
    
    %根据夹角判断法向量正确的指向
    flag = acos(ang);
    dis = norm(dir1);
    if flag<0
        dis = -dis;
    end
    
    %画出当前点的表面法向量
    t=(0:dis/100:dis)';
    x = p_cur(1) + C(1)*t;
    y = p_cur(2) + C(2)*t;
    z = p_cur(3) + (-1)*t;
    
    normE =[normE;x y z];
    i
end
pcshowpair(pc,pointCloud(normE));

原理

PCA原理

主元分析法PCA学习笔记

点云法向量与点云平面拟合的关系（PCA）

Estimating Surface Normals in a PointCloud

3D【24】PCA点云法向量估计

利用PCA计算点云的法线

3D点云法向量估计(最小二乘拟合平面)

为什么用PCA做点云法线估计？

利用PCA求点云的法向量

pca_demo.m

clc
clear
close all
n=50;
z=peaks(n);
x=1:n;y=1:n;
[x,y]=meshgrid(x,y);
P=[x(:),y(:),z(:)];

%读入屋顶点云
% P=load('building_Example2.txt');
% P=P(:,1:1:3);

%读入Building1.txt
% P=load('point_w.txt');
% P=P(:,1:1:3);


P=P';%P=3x2500
k=8;
[pn,pw] = pca(P, k);
pp = P+3.*pn;%pn为法向量，pw为评价曲率的一个参数

%求取曲率最大的500个点
P_=P';%P_ 是2500x3
[pw_,indice]=sort(pw);
indice=indice';

%找出曲率最大（或最小）的几个点，验证
% pw_first=indice(size(indice,1)-100:1:size(indice,1),:);%曲率最大值排前1000的点
pw_first=indice(1:1:700,:);%曲率最小值（最平坦）排前1000的点

for i=1:size(pw_first)
    Point_pw(i,1)=P_(pw_first(i),1);
    Point_pw(i,2)=P_(pw_first(i),2);
    Point_pw(i,3)=P_(pw_first(i),3);%Point_pw为曲率的前500个点
end


figure;
% scatter3( P(1,:)',P(2,:)',P(3,:)','.');
plot3(P(1,:)',P(2,:)',P(3,:)','g.');
hold on
plot3(Point_pw(:,1),Point_pw(:,2),Point_pw(:,3),'r*');
title("最平坦的前780个点");




% 使用matlab工具箱计算的法向量
figure(2);
P=P';
pt=pointCloud(P);
pcshow(pt);
hold on;

normals=pcnormals(pt,8);
u = normals(1:5:end,1);
v = normals(1:5:end,2);
w = normals(1:5:end,3);

x=P(1:5:end,1);
y=P(1:5:end,2);
z=P(1:5:end,3);
title('Matlab点云工具箱计算的 Normals of Point Cloud')
hold on
quiver3(x, y, z, u, v, w);
view(-37.5,45);
hold off

%使用pca计算的法向量
figure(3)
pcshow(pt);
hold on;
pn=pn';%对pn进行转置

u_p = pn(1:5:end,1);
v_p = pn(1:5:end,2);
w_p = pn(1:5:end,3);


title('PCA计算的 Normals of Point Cloud')
hold on
quiver3(x, y, z, u_p, v_p, w_p);
view(-37.5,45);
hold off


pca.m

% PCA主元分析法求法向量
% 输入：
% p:3*n的数值矩阵
% k:k近邻参数
% neighbors = transpose(knnsearch(transpose(p), transpose(p), 'k', k+1));
% neighbors一般可缺省。若之前做过k邻域求取操作也可直接调用，提高运算效率
% 输出
% n:法矢，已规定方向由邻域拟合出的平面指向查询点
% w:用于评估曲率的参数，详见：Mark P,et al. Multi-scale Feature Extraction on Point-Sampled Surfaces[J]. Computer Graphics Forum, 2010, 22(3): 281-289.

function [n,w] = pca(p, k, neighbors)%p为输入的点云
if nargin < 2
    error('no bandwidth specified')
end
if nargin < 3
    neighbors = transpose(knnsearch(transpose(p), transpose(p), 'k', k+1));%neighbor为一个索引矩阵，第一行代表第几个点，后8行代表K近邻的点。记录每个点及其周围的8个点
end
m = size(p,2);%返回第2维的维度
n = zeros(3,m); %存放法线的矩阵
w = zeros(1,m);
for i = 1:m
    x = p(:,neighbors(2:end, i));%x为8个邻域点    ，3x8的矩阵
    p_bar=mean(x,2);%每一行求均值(一共三行)
    
%     P =  (x - repmat(p_bar,1,k)) * transpose(x - repmat(p_bar,1,k));%中心化样本矩阵，再计算协方差矩阵
%     P = 1/(k) * (x - repmat(p_bar,1,k)) * transpose(x - repmat(p_bar,1,k)); %邻域协方差矩阵P
    P=(x - repmat(p_bar,1,k)) * transpose(x - repmat(p_bar,1,k))./(size(x,2)-1);
    
    [V,D] = eig(P);%求P的特征值、特征向量。 D是对应的特征值对角矩阵，V是特征向量(因为协方差矩阵为实对称矩阵，故特征向量为单位正交向量)
    
    [d0, idx] = min(diag(D)); %d0为最小特征值  idx为特征值的列数索引。diag():创建对角矩阵或获取矩阵的对角元素

    
    n(:,i) = V(:,idx);   % 最小特征值对应的特征向量为法矢，即法向量    
    
    %规定法矢方向指向
    flag = p(:,i) - p_bar;%由近邻点的平均点指向对应点的向量
    if dot(n(:,i),flag)<0%如果这个向量与法向量的数量积为负数（反向）
        n(:,i)=-n(:,i);%法向量取反向
    end
    if nargout > 1 
        w(1,i)=abs(d0)./sum(abs(diag(D)));%最小特征值的绝对值在协方差矩阵特征值绝对值的总和中占的比重
    end
end


效果 

在pca.m中是用到了matlab内置的K-NN查找K近邻点

近邻点查找效果如下，这里K取8，



使用手写的PCA计算出的法向量和matlab点云处理工具箱计算的法向量的效果对比如下，



 

可以看到，基本上是一样的。

 

 

参考周志华老师的西瓜书，对PCA的求解构成进行细节的推导，希望能有所帮助。
主成分分析简称PCA，是一种在尽可能减少信息损失的情况下找到某种方式降低数据的维度的方法。对于正交属性空间中的样本点，如何用一个超平面对所有样本进行恰当的表达？这里这个超平面应该有两个性质：
I.最近重构性：重构后的样本映射回原空间，与原样本的距离都足够的近；
II.最大可分性：样本在这个超平面上的投影尽可能分开。
本文主要基于这个两个性质，从PCA的公式推导的角度，来看优化在其中的作用。
一、最近重构性
设在n维空间中有m个样本点矩阵X：  
现在对这些点进行压缩，使其投影到k维空间中，其中k<n，使其损失的信息最小。假设数据样本进行了去中心化，即 
 ， 再假定矩阵 
 ，是投影变换后的得到的新坐标系， 是投影矩阵，其中 
 是标准的正交基向量，即满足：
若丢弃新坐标系中的部分坐标，即将维度降低到k<n。
由矩阵的乘法可知，将样本点 
 低维空间的坐标为 
 ，其中 
 ， j<=k。
所以，利用该坐标系重构数据，即把数据集合Z从k维空间重新映射回n维空间，得到新的坐标点 
这里假设新得到的数据与原始的数据点之间的距离最小，即PCA可转化为求解约束最优化问题：
            1.1
根据F范数与矩阵迹的关系:
           1.2
利用1.2对1.1进行简化:
  1.3
由于已知项，不会影响到最优化的结果；有负号转为最大值，所以1.4等式简化为：
       1.4
1.4式变换为：
     1.5
最终PCA的最优化问题就简化为：
     1.6
利用拉个朗日乘子法来求解1.6的最优化问题，引入乘子矩阵 
 ，若此时仅考虑约束
  ，则拉格朗日乘子矩阵 
 此时为对角矩阵 ， 另新的拉格朗日乘子矩阵为 
，1.6式转化为求解拉格朗日函数的极值：
     1.7
对1.7式求W的偏导数，在导数为0处取极值：
由矩阵微分公式 
 ：
         1.8
从1.8式根据特征向量的定义可知， 
 分别表示由协方差矩阵 
 的特征值和特征向量，特征值越大，数据在其对应额特征向量的方向上所包含的信息越丰富。
显然，此时只需要令
 和 
 分别为 协方差
 的前 
 个最大的特征值和单位特征向量就能使得目标函数达到最优值。
二、最大可分性
从最大可分性出发，能得到PCA的另一种解释。样本点 
 在新空间中超平面上的投影是 
 ，若所有样本点的投影能尽可能分开，则应该使投影后样本点的方差最大化：
投影后样本点的方差是 
 ，优化目标可以转化为矩阵的迹：
后面的求解构成与第一种的最近重构性一致。
PCA的算法步骤：
输入：n维空间的样本集合其中；映射到k维空间
1、归一化，将X中样本变换为标准正态分布 
：

     1.1         

     1.2        

2、计算协方差矩阵 

3、对协方差矩阵 
 进行特征分解

4、求取最大的k个特征值以及对应的特征向量，依次记录为 

输出，其中 
三、PCA的应用
（一）、数据降维。数据降维是处理高维度数据的基础。使用PCA降维有什么意义？
数据在低维下更容易进行处理与使用，算法的开销也将大大减少，比如在研究高维度数据分群中，使用无监督聚类的距离公式计算相似度不准确且开销大；
相关特征，重要特征通过降维能在数据中显现出来；同时，降至2维或3维也能进行可视化；
去除数据噪声，当数据受到噪声影响时，最小的特征值所对应的特征向量往往与噪声有关，将它们舍弃能再一定程度上起到去噪的作用。
（二）、人脸识别与手写数字识别。PCA在这方面的应用这几年随着深度学习技术的发展，而不断弱化。这里可以看下相关的blog文章，以人脸识别为例，计算人脸图像库中的“平均脸”并提取前K为特征向量，用于做映射矩阵。最终就是计算残差的大小，来判断人脸。残差公式就是上面的1.1式
PCA检测人脸的简单示例_matlab实现 - CSDN博客​blog.csdn.net
手写数字识别的例子没有找到，可以参考下PRML的截图：
南瓜书PumpkinBook_机器学习公式推导​datawhalechina.github.io阿泽：【机器学习】降维——PCA（非常详细）​zhuanlan.zhihu.com
矩阵乘法的本质是什么？​www.zhihu.com