Commit ada351be authored by Felix Kong's avatar Felix Kong

dmd_heart.m: rm old "script" code

parent 37be30ea
%% Try DMD on it
%% DMD for Hamlyn Center heart video dataset
% I could NOT get this to work well.
%
clear variables
load data/HamlynHeartImages.mat;
warning('This is not supposed to work well. Try Kernel DMD (coming soon!) for something that does work.');
disp('Data preparation...')
num_images = 100;
......@@ -12,9 +17,6 @@ for k = 1:num_images % first 100 images
Xall(:,k) = double(image_k(:));
end
X = Xall(:,1:end-1);
X2 = Xall(:,2:end);
% Step 1: initialize DMD object
verbose = true;
heart_dmd = dmd(Xall,verbose);
......@@ -28,59 +30,33 @@ initial_condition = Xall(:,1);
num_timesteps = size(Xall,2) + 100; % predict data and 100 frames into the future
xhat = heart_dmd.predict(initial_condition,num_timesteps);
disp('DMD Regression...')
[U0,S0,V0] = svd(X, 'econ');
figure(2);
clf
semilogy(diag(S0),'o');
%% Compute DMD (Phi are eigenvectors)
% r = size(X,2); % truncate at 5 modes
r = 50;
U = U0(:,1:r);
S = S0(1:r,1:r);
V = V0(:,1:r);
Atilde = U'*X2*(V/S);
[W,Lambda] = eig(Atilde);
Phi = X2*V*inv(S)*W;
eigs = diag(Lambda);
%% No calculations after this line.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% everything after this line is plotting
% there are no calculations after this line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b0 = Phi\X(:,1);
%% Calculate DMD output
disp('Reconstructing...')
% steps = size(X,2);
steps = 200;
xkvec = zeros(steps,size(X,1));
% calculate the output states
for j = 1:steps
xkvec(j,:) = real(Phi*Lambda^(j-1)*b0)';
if(norm(imag(xkvec(j,:))) > 1e-3)
error('xkvec(j,:) isnt real (has non-negligble imaginary component)')
end
reconstruction(:,:,j) = uint8((reshape(xkvec(j,:),size(grayImages,1),size(grayImages,2)))); % im2uint8 doesnt work
%% convert back into image format
for j = 1:size(xhat,2)
reconstruction(:,:,j) = uint8((reshape(xhat(:,j),size(grayImages,1),size(grayImages,2)))); % im2uint8 doesnt work
end
disp(norm(xkvec' - xhat))
%% replay
reconsructionHdl = implay(reconstruction);
set(reconsructionHdl.Parent,'Name',sprintf('Reconstruction with %d modes.',r));
set(reconsructionHdl.Parent,'Name',sprintf('Reconstruction with %d modes.',heart_dmd.r));
originalHdl = implay(grayImages(:,:,1:num_images));
set(originalHdl.Parent,'Name','Original data');
%%
%% show 1st mode
figure(3);
clf
mode_number = 1;
% dont use uint8 and imshow, U0 is unitary so columns have very small
% values
singleMode = surf(flipud(reshape(U0(:,mode_number),size(grayImages,1),size(grayImages,2))),'LineStyle','none');
singleMode = surf(flipud(reshape(heart_dmd.U0(:,mode_number),size(grayImages,1),size(grayImages,2))),'LineStyle','none');
colormap(flipud(gray))
modeEig = Lambda(mode_number,mode_number);
modeEig = heart_dmd.Lambda(mode_number,mode_number);
title(sprintf('Mode %d, eig %.3f + %.3fi', mode_number, real(modeEig), imag(modeEig)));
% view from above
axis([0 360 -0.001 288 -1 0])
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment