Decomposition Method

Singular Value Decomposition

A singular value decomposition (SVD) is a factorization of a matrix into several constitutive components all of which have a specific meaning in applications. Also, the SVD is a transformation that stretches/compresses and rotates a vector. SVD decomposition is given by Equation (1).

(1) $\begin{equation*} A=U \Sigma V^*\end{equation*}$

Where $U\in \mathbb{C}^{m\times m}$ and $V\in \mathbb{C}^{n\times n}$ are unitary matrices. $\Sigma \in \mathbb{R}^{m\times n}$ is diagonal, which components are non-negative and ordered from the largest to smallest.

In this page, we will use it to combine with the Principal Component Analysis.

Principal Component Analysis

When we have a matrix $X\in \mathbb{R}^{m\times n}$ , by implement SVD of the matrix $X$ , the data is projected into $Y$ in terms of the new orthonormal basis $U$ . So, in different measurements, the measurements have zeros covariance thus, the principal components were picking out. As a result, the principle information was extracted, and the noise was eliminated.

The calculating procedure is below, it can show how to erase noise and get the new diagonal matrix.

First, we project matrix $X$ in terms of basis $U$ :

As a result, get Equation (2).

(2) $\begin{equation*} U^{-1}X=Y\end{equation*}$

Since $U$ is a unitary matrix, it should have $U^{-1}=U^{*}$ .

(3) $\begin{equation*} U^{-1}X=U^{\ast }X=Y\end{equation*}$

Then, follows that to get the covariance of $Y$ .

(4) $\begin{equation*} \begin{align*}C_{Y}&=\frac{1}{n-1} YY^{T} \\&=\frac{1}{n-1} (U^{\ast }X){(U^{\ast }X)}^{T} \\&=\frac{1}{n-1} U^{\ast }(XX^{T})U \\&=\frac{1}{n-1} U^{\ast }U\Sigma^{2} UU^{\ast } \\&=\frac{1}{n-1} \Sigma^{2} \end{align*}\end{equation*}$

(5) $\begin{equation*} \Sigma =\left[ \begin{matrix}\sigma_{1} && && \\&& \ddots && \\&& && \sigma_{n}\end{matrix} \right]\end{equation*}$

As a result, get equation (6).

(6) $\begin{equation*}C_Y = \frac{1}{n-1} \Sigma^{2}\end{equation*}$

clear all; clc; close all;

load('cam1_1.mat')
load('cam1_2.mat')
load('cam1_3.mat')
load('cam1_4.mat')
load('cam2_1.mat')
load('cam2_2.mat')
load('cam2_3.mat')
load('cam2_4.mat')  
load('cam3_1.mat')
load('cam3_2.mat')
load('cam3_3.mat')
load('cam3_4.mat')

%%
% Tracking for camera 1 test 1
[x11, y11] = trackPG(vidFrames1_1,[321,228],[20,20,20,20]);
% Tracking for camera 2 test 1
[x21, y21] = trackPG(vidFrames2_1,[274,274],[20,20,20,20]);
% Tracking for camera 3 test 1
[x31, y31] = trackPG(vidFrames3_1,[318,271],[15,15,15,15]);

% SVD test 1
A1(1,:) = x11 - mean(x11);
A1(2,:) = y11 - mean(y11);
A1(3,:) = x21(10:235) - mean(x21(10:235));
A1(4,:) = y21(10:235) - mean(y21(10:235));
A1(5,:) = x31(1:226) - mean(x31(1:226));
A1(6,:) = y31(1:226) - mean(y31(1:226));
[U1,s1,V1] = svd(A1.',0);

sigment1 = diag(s1);
en11 = sigment1(1)/sum(sigment1);
en21 = sum(sigment1(1:2))/sum(sigment1);
en31 = sum(sigment1(1:3))/sum(sigment1);

figure(1)
subplot(2,1,1)
plot(1:length(sigment1),sigment1,'ro')
title('Singular Values')

subplot(2,1,2)
plot(1:6,V1(:,1),'.-',1:6,V1(:,2),'.-',1:6,V1(:,3),'.-')
legend('Mode 1','Mode 2','Mode 3','Location','southeast')
xticks([1 2 3 4 5 6])
title('Linear POD Modes')


%%
% Tracking for camera 1 test 2
[x12, y12] = trackP(vidFrames1_2,[325,308],[17,17*0.9,15,15*0.9]);
% Tracking for camera 2 test 2
[x22, y22] = trackP(vidFrames2_2,[314,357],[34,34,34,35]);
% Tracking Camera 2 test 2
% Tracking for camera 3 test 2
[x32, y32] = trackPG(vidFrames3_2,[349,245],[44,45,44,44]);

% SVD test 2
A2 = zeros(6,314);
A2(1,:) = x12 - mean(x12);
A2(2,:) = y12 - mean(y12);
A2(3,:) = x22(15:328) - mean(x22(15:328));
A2(4,:) = y22(15:328) - mean(y22(15:328));
A2(5,:) = x32(1:314) - mean(x32(1:314));
A2(6,:) = y32(1:314) - mean(y32(1:314));
[U2,s2,V2] = svd(A2.',0);

sigment2 = diag(s2);
en12 = sigment2(1)/sum(sigment2);
en22 = sum(sigment2(1:2))/sum(sigment2);
en32 = sum(sigment2(1:3))/sum(sigment2);

figure(2)
subplot(2,1,1)
plot(1:length(sigment2),sigment2,'ro')
title('Singular Values')

subplot(2,1,2)
plot(1:6,V2(:,1),'.-',1:6,V2(:,2),'.-',1:6,V2(:,3),'.-')
legend('Mode 1','Mode 2','Mode 3','Location','southeast')
xticks([1 2 3 4 5 6])
title('Linear POD Modes')


% Tracking for camera 1 test 3, tracking pink

[x13, y13] = trackPP(vidFrames1_3,[330,290],[20,20,20,20]);
% Tracking for camera 2 test 3, tracking pink
[x23, y23] = trackPP(vidFrames2_3,[252,294],[20,20,20,20]);
% Tracking for camera 3 test 3
[x33, y33] = trackPG(vidFrames3_3,[353,229],[17,17,17,17]);

% SVD test 3
A3 = zeros(6,206);
A3(1,:) = x13(20:225) - mean(x13(20:225));
A3(2,:) = y13(20:225) - mean(y13(20:225));
A3(3,:) = x23(5:210) - mean(x23(5:210));
A3(4,:) = y23(5:210) - mean(y23(5:210));
A3(5,:) = x33(10:215) - mean(x33(10:215));
A3(6,:) = y33(10:215) - mean(y33(10:215));
[U3,s3,V3] = svd(A3.',0);
sigment3 = diag(s3);

en13 = sigment3(1)/sum(sigment3);
en23 = sum(sigment3(1:2))/sum(sigment3);
en33 = sum(sigment3(1:3))/sum(sigment3);

figure(3)
subplot(2,1,1)
plot(1:length(sigment3),sigment3,'ro')
title('Singular Values')

subplot(2,1,2)
plot(1:6,V3(:,1),'.-',1:6,V3(:,2),'.-',1:6,V3(:,3),'.-')
legend('Mode 1','Mode 2','Mode 3','Location','southeast')
xticks([1 2 3 4 5 6])
title('Linear POD Modes')

% Tracking for camera 1 test 4
[x14, y14] = trackPP(vidFrames1_4,[402,263],[44,45,44,44]);
% Tracking for camera 2 test 4
[x24, y24] = trackPP(vidFrames2_4,[244,243],[44,45,44,44]);
% Tracking for camera 3 test 4
[x34, y34] = trackP(vidFrames3_4,[363,244],[40,40,30,30]);

% SVD test 4
A4 = zeros(6,392);
A4(1,:) = x14 - mean(x14);
A4(2,:) = y14 - mean(y14);
A4(3,:) = x24(1:392) - mean(x24(1:392));
A4(4,:) = y24(1:392) - mean(y24(1:392));
A4(5,:) = x34(1:392) - mean(x34(1:392));
A4(6,:) = y34(1:392) - mean(y34(1:392));
[U4,s4,V4] = svd(A4.',0);

sigment4 = diag(s4);
en14 = sigment4(1)/sum(sigment4);
en24 = sum(sigment4(1:2))/sum(sigment4);
en34 = sum(sigment4(1:3))/sum(sigment4);

figure(4)
subplot(2,1,1)
plot(1:length(sigment4),sigment4,'ro')
title('Singular Values')

subplot(2,1,2)
plot(1:6,V4(:,1),'.-',1:6,V4(:,2),'.-',1:6,V4(:,3),'.-')
legend('Mode 1','Mode 2','Mode 3','Location','southeast')
title('Linear POD Modes')

% Plot all position vectors

figure(5)
plot(1:226,y11,'r')
hold on
plot(1:226,y21(10:235),'k')
plot(1:226,x31(1:226),'b')
hold off
title('Test 1: Ideal Case')
lgd = legend('Camera 1', 'Camera 2', 'Camera 3','Location','southeast');
lgd.FontSize =7;
ylabel('Displacement')
xlabel('Frame')

figure(6)
plot(1:314,y12,'r')
hold on
plot(1:314  ,y22(15:328),'k')
plot(1:314,x32(1:314),'b')
lgd = legend('Camera 1', 'Camera 2', 'Camera 3','Location','southeast');
lgd.FontSize =7;
hold off
title('Test 2: Noisy Case')
ylabel('Displacement')
xlabel('Frame')


figure(7)
plot(1:206,y13(20:225),'r')
hold on
plot(1:206  ,y23(5:210),'k')
plot(1:206,x33(10:215),'b')
hold off
title('Test 3: Horizontal Displacement')
ylabel('Displacement')
xlabel('Frame')
lgd = legend('Camera 1', 'Camera 2', 'Camera 3','Location','southeast');
lgd.FontSize =7;

figure(8)
plot(1:392,y14,'r')
hold on
plot(1:392  ,y24(1:392),'k')
plot(1:392,x34(1:392),'b')
hold off
title('Test 4: Horizontal Displacement and Rotation')
ylabel('Displacement')
xlabel('Frame')
lgd = legend('Camera 1', 'Camera 2', 'Camera 3','Location','southeast');
lgd.FontSize =7;

%------------------------------------------------------------------------------

function [xp,yp] = trackPG(vidFrames,guess,width)
dim = size(vidFrames);
L = dim(4); %length of video in frames
xp = zeros(L,1);
yp = zeros(L,1);
maxloc = [0,0];
maxval = 0;

init_guess = guess;
xp(1) = init_guess(1);
yp(1) = init_guess(2);

%center the first search window to the initial point
center = init_guess;

for f = 2:L
    frame = double(rgb2gray(vidFrames(:,:,:,f)));
    for x = center(1)-width(1):center(1)+width(2)
        for y = center(2)-width(3):center(2)+width(4)
            point = frame(y,x,:);
            if point > maxval
                maxloc = [x,y];
                maxval = point;
            end   
        end
    end
    xp(f) = maxloc(1);
    yp(f) = maxloc(2);
    maxloc = [0,0];
    maxval = 0;
    
    center(1) = xp(f);
    center(2) = yp(f);
end
end
%-------------------------------------------------------------
function [xp,yp] = trackP(vidFrames,guess,width)
dim = size(vidFrames);
L = dim(4);     %length of video in frames
xp = zeros(L,1);
yp = zeros(L,1);
maxloc = [0,0];
maxval = 0;

init_guess = guess;
xp(1) = init_guess(1);
yp(1) = init_guess(2);

%center the first search window to the initial point
center = init_guess;

for f = 2:L
    frame = double(vidFrames(:,:,:,f));
     for x = center(1)-width(1):center(1)+width(2)
        for y = center(2)-width(3):center(2)+width(4)
            point = frame(y,x,:);
            colorsum = point(1,1,1)+point(1,1,2)+point(1,1,3);
            if colorsum > maxval
                maxloc = [x,y];
                maxval = colorsum;
            end   
        end
    end
    xp(f) = maxloc(1);
    yp(f) = maxloc(2);
    maxloc = [0,0];
    maxval = 0;
    
    center(1) = xp(f);
    center(2) = yp(f);
end
end
%-------------------------------------------------------------
function [xp,yp] = trackPP(vidFrames,guess,width)
dim = size(vidFrames);
L = dim(4); %length of video in frames
xp = zeros(L,1);
yp = zeros(L,1);
maxloc = [0,0];
maxred = 0;

init_guess = guess;
xp(1) = init_guess(1);
yp(1) = init_guess(2);

%center the first search window to the initial point
center = init_guess;

for f = 2:L
    frame = double(vidFrames(:,:,:,f));
     for x = center(1)-width(1):center(1)+width(2)
        for y = center(2)-width(3):center(2)+width(4)
            point = frame(y,x,:);
            if ((point(1,1,1) > maxred)&(point(1,1,2)<150)&(point(1,1,3)<150))
                maxloc = [x,y];
                maxred = point;
            end   
        end
    end
    xp(f) = maxloc(1);
    yp(f) = maxloc(2);
    maxloc = [0,0];
    maxred = 0;
    
    center(1) = xp(f);
    center(2) = yp(f);
end
end

Welcome to Intelligent Additive Manufacturing Lab