Referències ortonormals i matrius simètriques: Anàlisis de Components Principals (PCA)

El càcul de vectors propis en Matlab està simplificat per la funcio eig (eigen values son valors propis en anglès). Aquesta funció retorna els valors i també els vectors propis associats. En el cas que la matriu sigui simètrica, la funció eig retorna una b.o.n. formada per vectors propis.

Utilització de la funció eig de Matlab

A = [1,2,3; 2,5,7; 5,8,9]
A = 3×3
1 2 3 2 5 7 5 8 9
D = eig(A) % retorna en un vector els valors propis (l'ordre de magnitud pot variar)
D = 3×1
15.9853 0.2094 -1.1947
[V,D] = eig(A) % retorna en dues matrius els vectors i els valors propis. Les columnes de V
V = 3×3
-0.2339 -0.5871 -0.3235 -0.5524 0.7364 -0.6608 -0.8001 -0.3362 0.6772
D = 3×3
15.9853 0 0 0 0.2094 0 0 0 -1.1947
% son els vectors propis associats a cada valor propi.
% verifiquem que son vectors i valors propis dins la tolerància numèrica:
% apliquem la definició pel primer vector propi (1a columna de la matriu V)
error = A*V(:,1)-D(1,1)*V(:,1) %El que importa és l'exponent del resultat!!!
error = 3×1
10-13 ×
-0.1021 -0.1066 -0.1066
Per veure el canvi en l'ordre dels valors propis feu
A = magic(5) % matriu especial en la que totes les files i columnes sumen el mateix
A = 5×5
17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9
[V,D] = eig(A)
V = 5×5
-0.4472 0.0976 -0.6330 0.6780 -0.2619 -0.4472 0.3525 0.5895 0.3223 -0.1732 -0.4472 0.5501 -0.3915 -0.5501 0.3915 -0.4472 -0.3223 0.1732 -0.3525 -0.5895 -0.4472 -0.6780 0.2619 -0.0976 0.6330
D = 5×5
65.0000 0 0 0 0 0 -21.2768 0 0 0 0 0 -13.1263 0 0 0 0 0 21.2768 0 0 0 0 0 13.1263
%
% Per donar els vaps ordenats podem fer servir la comanda sort de Matlab
%
[variableInutil,permutacio] = sort(diag(D)); %determinem l'ordre creixent dins de permutacio
D = D(permutacio,permutacio) %reordenem les matrius corresponents
D = 5×5
-21.2768 0 0 0 0 0 -13.1263 0 0 0 0 0 13.1263 0 0 0 0 0 21.2768 0 0 0 0 0 65.0000
V = V(:,permutacio)
V = 5×5
0.0976 -0.6330 -0.2619 0.6780 -0.4472 0.3525 0.5895 -0.1732 0.3223 -0.4472 0.5501 -0.3915 0.3915 -0.5501 -0.4472 -0.3223 0.1732 -0.5895 -0.3525 -0.4472 -0.6780 0.2619 0.6330 -0.0976 -0.4472
A = [1,2,3; 2,5,7; 3,7,9] %matriu simètrica
A = 3×3
1 2 3 2 5 7 3 7 9
[V,D] = eig(A) %els vectors V formen b.o.n.
V = 3×3
0.3944 0.8855 0.2459 0.6712 -0.4603 0.5810 -0.6277 0.0641 0.7758
D = 3×3
-0.3707 0 0 0 0.1776 0 0 0 15.1931
V(:,1)'*V(:,2) %veiem per exemple que <V1,V2>=0
ans = 0

Referències ortonormals: Càlcul de la Bounding Box d'un objecte

Amb l'aparició dels escaners, cada vegada és més facil disposar d'un núvol de punts en 3D que es correspon amb un determinat objecte. En molts casos és important saber l'espai que ocupen aquests punts, per això es calcula la caixa envolvent o Bounding Box que pot ser de dos tipus:
Anem a veure com s'han de fer els càlculs en cada cas

Càrrega dels punts escanejats

Carreguem un fitxer on hi estan posats els punts. Cada fila és un punt i te tres columnes (cada coordenada).
punts = load('puntsObj.txt','-ascii');
numPunts=size(punts,1); %les files son el número de punts
plot3(punts(:,1),punts(:,2),punts(:,3),'.'); % visualitzem els punts carregats
% nota: en la finestra del dibuix premeu la fletxa circular i podreu rotar el dibuix en 3D

Aligned Bounding Box (ABB):

Calculem els valors màxim i mínim associats a cada coordenada i calculem el centre de la caixa.
maxX=max(punts(:,1));
minX=min(punts(:,1));
maxY=max(punts(:,2));
minY=min(punts(:,2));
maxZ=max(punts(:,3));
minZ=min(punts(:,3));
%mides de la caixa
midaX=maxX-minX;
midaY=maxY-minY;
midaZ=maxZ-minZ;
% el centre serà el punt mig de les mides de la caixa
centre=0.5*[maxX+minX,maxY+minY,maxZ+minZ];
%Ara visualitzem la caixa. Per això crearem una funció propia 'plotBB' que trobareu al
%final d'aquest guió
eixos=eye(3); %pel cas ABB els eixos de la caixa son els vectors de la base canònica
plotBB(punts, midaX, midaY, midaZ, eixos, centre);

Oriented Bounding Box (OBB):

Trobarem, en aquest ordre, els seus eixos, les seves mides i el seu centre.
Ho calcularem en sis passos:
  1. Trobar les coordenades dels punts en la referència R1 = {centroid; e1,e2,e3}, on e1,e2,e3 és la base canònica de
  2. Calcular la matriu de covariances C. És una matriu 3x3, simètrica i es defineix segons la fórmula de sota, usant les coordenades dels punts en referència R1, on N és el nombre total de punts.
  3. Trobar una b.o.n. {v1,v2,v3} formada per vectors propis de la matriu de covariances (Principal Components), que seran els eixos de la caixa OBB. En termes de Mecànica, aquests vectors propis s'anomenen Eixos Principals d'Inèrcia. El vector propi associat al valor propi màxim és el que dona menor moment d'inèrcia en rotar el cos respecte a aquest eix.
  4. Trobar les coordenades dels punts en la referència (ortonormal) R2 = {centroid; v1,v2,v3}. En aquesta referència, la caixa OBB es calcula com la caixa ABB (perquè la caixa OBB és la "caixa ABB respecte els eixos v1,v2,v3").
  5. Trobar les mides de la caixa OBB i les coordenades del seu centre en referència R2.
  6. Trobar les coordenades del centre de la caixa OBB en referència ordinària.
%
% Coordenades en referència R1 (restem respecte del nou origen que és el centroid)
%
centroid=mean(punts); %fa les mitjanes per columnes
coord1=[punts(:,1)-centroid(1), punts(:,2)-centroid(2), punts(:,3)-centroid(3)];
% Calculem la matriu de covariances (és simètrica)
%
xx=coord1(:,1);
yy=coord1(:,2);
zz=coord1(:,3);
C(1,1)=sum(xx.*xx); % aquí fem la suma i producte de les x per les x
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/numPunts;
% Calculem les components principals (PCA)
% (son els eixos de la OBB)
%
[V, D]=eig(CC);
%
% Expressem els punts en la referència R2 = {centroid; v1,v2,v3}. El canvi
% de referència R1 a referència R2 es fa, simplement, multiplicant per la
% inversa de V que coincideix amb la seva transposada: inv(V)=V' (perquè V
% és b.o.n.)
%
coord2=V'*coord1'; %transposem coord1 per poder multiplicar i després ho refem.
coord2=coord2'; %les mides de la matriu coord2 son com les de la matriu de punts original.
%
% Ara el càlcul de la OBB és equivalent al càlcul de la ABB en la nova
% referència (repetim les instruccions anteriors)
%
maxX=max(coord2(:,1));
minX=min(coord2(:,1));
maxY=max(coord2(:,2));
minY=min(coord2(:,2));
maxZ=max(coord2(:,3));
minZ=min(coord2(:,3));
%mides de la caixa
midaX=maxX-minX;
midaY=maxY-minY;
midaZ=maxZ-minZ;
% centre de la caixa en la referència R2
centre2=0.5*[maxX+minX,maxY+minY,maxZ+minZ];
%
% Finalmet cal expressar els punts en la referència ordinària.
% Trobem les coordenades del centre en referència ordinària.
%
centre = centroid + (V*centre2')';
%Ara visualitzem la caixa.
eixos=V; %els eixos son els de la nova base de Veps
plotBB(punts, midaX, midaY, midaZ, eixos, centre);
view([23.31 -2.66]) %punt de vista de la càmera 3D

Exercici 1:

Calculeu els volums de les respectives caixes ABB i OBB d'aquest exemple.
Sol: volum(ABB) = 1.9402 , volum(OBB) = 0.5939

Exercici 2:

Calculeu els centres de cadascuna de les 6 cares de la OBB d'aquest exemple i pinteu-los sobre la OBB.
Sol:
0.6682 0.6632 0.6714
0.3368 0.8288 0.8101
0.5135 0.9852 0.4815
0.4915 0.5069 1.0000
1.0825 1.3888 1.3584
-0.0776 0.1032 0.1231
puntsMitjos.png

Funció per pintar la BB

function plotBB(punts, midaX, midaY, midaZ, eixos, centre)
% Plot the BB of a set of points according to the directions and
% size of the Box
%
% (c) Numerical Factory 2020
%
figure() %obrim una nova finestra de dibuix
plot3(punts(:,1),punts(:,2),punts(:,3),'.'); % visualitzem els punts
hold on;
v1=eixos(:,1);
v2=eixos(:,2);
v3=eixos(:,3);
cm=centre';
plot3(centre(:,1),centre(:,2),centre(:,3),'ro'); % visualitzem els punts
%-------------------------------------
% Calculem els 8 vertex de la caixa
%-------------------------------------
% cara inferior
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);
% cara superior
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);
%-------------------------------------
%pintem les arestes de la caixa
%-------------------------------------
% cara inferior
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');
% cara superior
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');
% arestes verticals
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
(c) Numerical Factory 2020