Projeccions sobre un subespai o el seu ortogonal.

Estudiarem com fer projeccions d'un vector sobre un subespai (podeu consultar el document d'Atenea: problemesOctaveMatlabTema2.pdf)

(recordeu a fer un clear all si heu estat treballant previament amb Matlab)

Projecció d'un vector sobre un subespai
Projecció ortogonal d'un punt sobre una varietat lineal.
Projecció NO ortogonal d'un punt sobre un pla.
Aplicació: Ombres d'un objecte
Exercici:

Projecció d'un vector sobre un subespai

Segons les dimensions del subespai es pot fer el càlcul de diferents formes.

En primer lloc farem la projecció d'un vector v d'R^4 sobre un subespai F 3-dim. v=(2,1,0,3) i F={(x,y,z,t) d'R^4 tal que 2x+y+3z+t=0}

%-----------------------------------------------
% Primer mètode
%-----------------------------------------------
v=[2;1;0;3]; %normalment sempre es consideren vectors columna
% busquem la base del subespai
B=[2,1,3,1]; %escribim l'equació de F
A=null(B) %la base de vectors de F està formada per les columnes d'A (nucli de B)
% comprovem que efectivament és un subespai 3-dim
rank(A)

A =

   -0.2582   -0.7746   -0.2582
    0.9560   -0.1319   -0.0440
   -0.1319    0.6043   -0.1319
   -0.0440   -0.1319    0.9560


ans =

     3

% fem la projecció resolent el sistema A'A p = A'v
%(sempre que es tingui la base en comptes de l'equació es fa igual)
%
p=A'*A\A'*v;
proj1=A*p

%-----------------------------------------------
% segon mètode:
%-----------------------------------------------
% Aprofitant que dim(F)+dim(F_ort) = 4, vol dir que dim(F_ort)=1 i
% llavors és més fàcil projectar sobre ortF i calcular proj(v_sobre_F)=v-proj(v_sobre_F_ort)
%
F_ort=[2;1;3;1]; %coeficients de l'equació
% projectem sobre l'ortogonal
projF_ort=(dot(F_ort,v)/dot(F_ort,F_ort))*F_ort;
proj2=v-projF_ort
compError =(proj1-proj2)'

proj2 =

    0.9333
    0.4667
   -1.6000
    2.4667


compError =

   1.0e-15 *

    0.5551    0.0555         0         0

Projecció ortogonal d'un punt sobre una varietat lineal.

En el cas de considerar un punt de R^3, en comptes d'un vector, llavors podem passar el problema al cas anterior a partir d'algun punt de la varietat lineal que farà d'origen. Exemple: Calcular la projecció ortogonal de P=(2,1,3) sobre el pla definit de la forma x=2+a+b; y=3+2a+2b; z=4+3a+2b;

P=[2;1;3];
puntPla=[2;3;4]; %punt del pla
v=P-puntPla;
u1=[1;2;3]; %vectors directors del pla
u2=[1;2;2];
A=[u1,u2]; %els posem com a columnes de la matriu
pp=A'*A\A'*v; %
projv=A*pp;
projP=puntPla+projv

projP =

    1.2000
    1.4000
    3.0000

Projecció NO ortogonal d'un punt sobre un pla.

En el cas de fer la projecció en una direcció NO ortogonal d'un determinat punt sobre un pla, el problema es redueix, donada la direcció de projecció, a calcular la recta que passa pel punt, te aquesta direcció i fer la intersecció amb el pla.

Veiem un exemple: Projecteu el punt Q=(2,1,3) sobre el pla de l'exemple anterior segons la direcció v=(1,1,1) El millor és treballar en forma paràmètrica de la recta r = Q+mu*v i la forma general del pla: ax+by+cz+d=0. Fent la substituicó de r en el pla obtenim el valor del paràmetre mu

Q=[1;2;3];
% calculem l'equació general del pla vNormal=[a,b,c]
normalPla=cross(u1,u2); %el vector normal és el producte vectorial dels dos vectors directors
d = -normalPla'*puntPla; %calculem el terme independent a partir de la normal i el punt del pla
% Calculem el valor del paràmetre a partir de la substitució de r en
% l'equació del pla i aillant mu
mu=(-normalPla'*Q-d)/(normalPla'*v);
projQ=Q+mu*v %verifiqueu que és un punt que pertany al pla !!

projQ =

    1.0000
    1.0000
    2.5000

Aplicació: Ombres d'un objecte

L'ombra que projecta un objecte sobre un pla és un cas de projecció. Els punts d'ombra es troben a partir de la projecció dels punts de l'objecte sobre el pla. Segons la direcció de la llum, hi ha dos possibles casos: * Una llum puntual: Seria el cas d'il·luminar un objecte amb una llanterna. La direcció de projecció és la que uneix el focus de llum amb cada punt de l'objecte. * Una llum direccional: Si el focus de llum està molt lluny (el Sol) tots els rajos de llum tindran sempre la mateixa direcció (depenent de l'hora solar)

% Llum puntual: apliquem-ho als punts d'un objecte quan el focus de llum es troba
% en el punt Q=(0,0,10) i projectem sobre el pla que està a la base de
% l'objecte.

punts=load('Bunny.txt','-ascii'); %carreguem un conjunt de punts 3D
% es tracta d'una taula de files: numPunts i colum: 3 coordenades.
figure('units','normalized','outerposition',[0 0 0.3 1]) %fem la figura full screen
plot3(punts(:,1),punts(:,2),punts(:,3),'.');
hold on;
axis equal;
%
% Q és el punt on està la llum
%
Q = [0,0,10]';
plot3(Q(1),Q(2),Q(3),'or'); %pintem el punt de llum
plot3(Q(1),Q(2),Q(3),'*y');
%
% Definim el pla d'ombra de la forma ax+by+cz+d=0
%
base=min(punts(:,3)); %calculem la z de la base del model
pla = [0,0,1,-base];  %pla horitzontal situat a la base
normalPla=[pla(1),pla(2),pla(3)]';
d=pla(4);
%
% és com el cas anterior: projectarem Q sobre el pla, però per cada punt P de
% l'objecte la direcció de projecció serà v=P-Q
normPlaQ = normalPla'*Q; %aquest terme és constant en totes les direccions
numPunts = size(punts,1);
ombra=[];
for i=1:numPunts
    puntObj=punts(i,:)';
    v=puntObj-Q;
%    plaP = pla(1)*punts(i,1)+pla(2)*punts(i,2)+pla(3)*punts(i,3);
    mu=(-normPlaQ-d)/(normalPla'*v);
%    projQ=Q+mu*v %verifiqueu que és un punt que pertany al pla !!

%     lambda=(-plaQ-pla(4))/(plaP-plaQ);
%    if (mu > 0)
        ptOmbra=Q+mu*v;
        ombra=[ombra; ptOmbra'];
%    end
end
plot3(ombra(:,1),ombra(:,2),ombra(:,3),'k.');

Exercici:

Implementeu el cas d'una ombra generada per una llum direccional. Per exemple, per la direcció v=(-1,-1,-1), la projecció del punt 425 del model sobre el pla de la base és el punt:

-4.3196e+00  -4.9877e+00  -1.4484e+00