Principal Component Analisys (PCA)

Aplication of SVD: Computation of the Bounding Box of an object

With the developement of scanners, it is becoming easier to have a cloud of 3D points that correspond to a given object. In many cases it is important to know the space occupied by these points, so the envelope box or Bounding Box is calculated which can be of two types:

Aligned Bounding Box (ABB): The box is built with the planes parallel to the axes.
Oriented Bounding Box (OBB): The box is constructed appropriately at the given points.

Let's see how the calculations should be done in each case.

Load data points

We upload a file where the dots are recorded. Each row is a point and has three columns (each coordinate).

points = load('puntsObj.txt','-ascii');

numpoints=size(points,1); %the rows are the number of points

% Obs: in the plot window press the circular arrow and you can rotate the plot in 3D

plot3(points(:,1),points(:,2),points(:,3),'.'); % visualize the loaded points

Aligned Bounding Box:

we compute the maximum and minimum values associated with each coordinate, size of the box and the centroid point.

maxP=max(points); %it gives together: xmax,ymax,zmax

minP=min(points); %it gives together: xmin,ymin,zmin

% size of the box

boxsize=maxP-minP;

center=0.5*(maxP+minP); % centroid

%Now we visualize the box. For that we use a new function 'plotBB' that you can find at the end of this file.

boxAxis=eye(3); %for the ABB case, the axis of the box are the usual R^3 axis

plotBB(points, boxsize, boxAxis, center);

Oriented Bounding Box

To compute the OBB we need four steps:

To move the origin of coordinates to the centroid (express the points in local coordinates)
Compute the covariance matrix C, which is a 3x3 matrix defined below. N is the total number of points.
Compute the eigen vectors of the covariance matrix (Principal Components) that will be the axis of the OBB. In terms of Mechanics, these eigenvectors are called Principal Axes of Inertia. The eigenvector associated with the maximum eigenvalue is the one that gives less moment of inertia when rotating the body with respect to this axis. Because we are using SVD decomposition and the covariance matrix is symmetric, the eigenvectors are the columns of the V matrix and the eigen values are the singular values.
Compute the size of the box.

%
% local Coordinates (substract from the new centroid origin)
%
centroid=mean(points); 
localCoord=points-centroid; %substract to each point
% Compute the covariance matrix C (symmetric)
%
xx=localCoord(:,1);
yy=localCoord(:,2);
zz=localCoord(:,3);
 
C(1,1)=sum(xx.*xx); 
C(1,2)=sum(xx.*yy);
C(1,3)=sum(xx.*zz);
C(2,1)=C(1,2);
C(2,2)=sum(yy.*yy);
C(2,3)=sum(yy.*zz);
C(3,1)=C(1,3);
C(3,2)=C(2,3);
C(3,3)=sum(zz.*zz);
CC=C/(numpoints-1);

New reference vectors (PCA)

Compute new reference oriented according to the eigen vector basis (V matrix of the SVD decomposition).

Then, we compute the local coordinates of the points in the new reference base.

% Compute the principal components (PCA)
% (they are the axis of the OBB)
%
[U, D, V]=svd(CC);
%
% We express the 'local' points in the new base corresponding to the eigenvectors V making the product
% with the inverse. Because V is orthogonal, the inverse is V'.
%
inV = V';
% we transpose localC to perform matrix multiplication 
% and then we recover the matrix shape.
ptO=inV*localCoord'; 
ptO=ptO'; 

Compute the Oriented Bounding Box

% Now the computation of the OBB is equivalent to the ABB in the new reference

% (we repeat the computations)

maxP=max(ptO);

minP=min(ptO);

% size of the box

boxsize=maxP-minP;

% centroid

centLocal=0.5*(maxP+minP);

% Finally we need to express these points in the original reference:

% we change the local centroid (centL) to the original base and we add it

% the initial centroid.

centLWorld=V*centLocal';

centerReal=centroid+centLWorld';

%now we visualize the BB

boxAxis=V; %the axis are now the eigenvectors

plotBB(points, boxsize, boxAxis, centerReal);

view([-64.80 23.40])

Function to plot the BB

function plotBB(points, boxsize, eixos, centroid)
% Plot the BB of a set of points according to the directions and
% size of the Box.
%
% Input:
% points  -> nx3 matrix with the points coordinates
% boxsize -> 3 vector with the lengths of the 3 axis of the box 
% eixos   -> 3x3 matrix with the directions of the axis of the box
% centroid-> 1x3 vector with the centroid of the points (it's mean).
% 
% (c) Numerical Factory 2023
 
    figure() 
    plot3(points(:,1),points(:,2),points(:,3),'.'); 
    hold on;
    v1=eixos(:,1);
    v2=eixos(:,2);
    v3=eixos(:,3);
    cm=centroid';
    plot3(centroid(:,1),centroid(:,2),centroid(:,3),'ro'); 
 
    %-------------------------------------
    % Compute the 8 vertices of the box
    %-------------------------------------
    % lower face
    midaX=boxsize(1);
    midaY=boxsize(2);
    midaZ=boxsize(3);
    cv1=cm+0.5*(-midaX*v1-midaY*v2-midaZ*v3);
    cv2=cm+0.5*( midaX*v1-midaY*v2-midaZ*v3);
    cv3=cm+0.5*( midaX*v1+midaY*v2-midaZ*v3);
    cv4=cm+0.5*(-midaX*v1+midaY*v2-midaZ*v3);
    % upper face
    cv5=cm+0.5*(-midaX*v1-midaY*v2+midaZ*v3);
    cv6=cm+0.5*( midaX*v1-midaY*v2+midaZ*v3);
    cv7=cm+0.5*( midaX*v1+midaY*v2+midaZ*v3);
    cv8=cm+0.5*(-midaX*v1+midaY*v2+midaZ*v3);
 
    %-------------------------------------
    % draw the edges of the box
    %-------------------------------------
    % lower face
    aresta=[cv1';cv2'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv2';cv3'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv3';cv4'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv4';cv1'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
 
    % upper face
    aresta=[cv5';cv6'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv6';cv7'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv7';cv8'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv8';cv5'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    
    % vertical edges
    aresta=[cv1';cv5'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv2';cv6'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv3';cv7'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
    aresta=[cv4';cv8'];
    plot3(aresta(:,1),aresta(:,2),aresta(:,3),'r');
 
end