Here is a quick demo:
function pca_point_cloud_demo
% Get random rotation matrix
r=randn(3,1);
r=r/norm(r); % direction of rotation vector
t=(50*rand(1)+25)*(pi/180); % rotation amount; between 25 and 75 degrees
r=t*r; % random rotation vector
K=zeros(3);
K(1,2)=-r(3);
K(1,3)= r(2);
K(2,3)=-r(1);
K=K-K'; % skew symmetric matrix
R=expm(K); % rotation matrix corresponding to r
% Simulate uniformly distributed point cloud
N=1E3;
X=2*rand(N,3)-1;
idx=sum(X.^2,2)>1;
X(idx,:)=[];
% Non-uniformly scale, rotate, and translate X
S=diag(sort(5*rand(1,3),'descend'));
t=10*randn(1,3);
X=bsxfun(@plus,(R*S*X')',t);
% Do PCA on X
X_ave=mean(X,1); % centroid
dX=bsxfun(@minus,X,X_ave);
C=dX'*dX; % covariance matrix
[U,~]=svd(C);
U(:,3)=cross(U(:,1),U(:,2)); % make sure there is no reflection
% Transformation that aligns centroid of X with the origin and its
% prinicial axes with Cartesian basis vectors
T1=eye(4); T1(1:3,4)=-X_ave(:);
T2=eye(4); T2(1:3,1:3)=U';
T=T2*T1;
% Apply T to X to get Y
Y=X;
Y(:,4)=1;
Y=(T*Y')';
Y(:,4)=[];
% PCA-based bounding box (for better visualization)
BBo=unit_cube_mesh;
L=max(Y)-min(Y);
V=bsxfun(@times,L,BBo.vertices);
V=bsxfun(@plus,V,min(Y));
BBo.vertices=V; % BB around Y
V(:,4)=1;
V=(T\V')';
V(:,4)=[];
BB=BBo;
BB.vertices=V; % BB around X; same as BBo but rotated and traslated
% Visualize X and Y
figure('color','w')
subplot(1,2,1)
plot3(X(:,1),X(:,2),X(:,3),'.k','MarkerSize',20)
axis equal
set(get(gca,'Title'),'String','Original Point Cloud','FontSize',20)
view([20 20])
hold on
h=patch(BB);
set(h,'FaceColor','b','FaceAlpha',0.25,'EdgeColor','r')
subplot(1,2,2)
plot3(Y(:,1),Y(:,2),Y(:,3),'.k','MarkerSize',20)
axis equal
set(get(gca,'Title'),'String','After PCA Normalization','FontSize',20)
view([20 20])
hold on
h=patch(BBo);
set(h,'FaceColor','b','FaceAlpha',0.25,'EdgeColor','r')
function fv=unit_cube_mesh
% Construct quadrilateral mesh of a unit cube with edges along x-, y-, and
% z-axes, and one corner at the origin.
X=[0 0 0; ...
1 0 0; ...
1 1 0; ...
0 1 0; ...
0 0 1; ...
1 0 1; ...
1 1 1; ...
0 1 1];
F=[1 4 3 2;
5 6 7 8;
2 3 7 6;
3 4 8 7;
1 5 8 4;
1 2 6 5];
fv.faces=F;
fv.vertices=X;
