% PCA example:  learning by trial and error
%       cao 11/3/06
%
% Data:  single slice (single surface coil) for a block-design experiment;
% flickering checkerboards off for 8 TRs, on for 8 TRs; 10 cycles
load slice
% First, traditional Fourier analysis
f = fft(slice,[],3); % FT 3rd dimension
% look at 10 cycles/scan component
figure(1)
subplot(1,2,1); imagesc(abs(f(:,:,11))); axis image
% and plot a voxel
subplot(1,2,2); plot(squeeze(slice(51,52,:))); axis tight

% Principle components analysis.  This relies on Matlab's SVD routine,
% which for a given matrix returns 3 matrices.  
%
% [U S V] = svd(X)
%
% When the input matrix is non-symmetric, there are 2 sets of eigenvectors 
% associated with each eigenvalue.  The left eigen vectors are returned in
% U, the right eigenvectors are returned in V, and the eigenvalues are
% returned in S (on the diagonal). X = U*S*V.  Just to play, check with a
% small matrix:

X1 = [1 2 3; 4 5 6]; % X is 2 x 3
[U1 S1 V1] = svd(X1); % U is 2 x 2 and V is 3 x 3

% We should be able to multiply these back together to get X:
U*S*V' % not sure why the primes ...
% X should be the eigenvalue-weighted sum of the products of the principle
% components
S(1,1)*U(:,1)*V(:,1)' + S(2,2)*U(:,2)*V(:,2)'
% This points toward why the rows of U and the columns of V are the paired
% eigenvectors of the matrix: the data can be recreated by summing
% component matrices, each of which is the product of a vector in U and a
% vector in V.

% Now with real data.
tic; [U S V] = svd(reshape(slice,[64*128 160])'); toc
% The reshaping puts voxels on the row dimension and time in the column
% dimension.  I thought SVD ran faster when the 1st dimension is the smallest,
% hence the transpose.  Tic and toc will tell us that this takes
% a bit over 2 minutes to complete; doing it the other way around (with
% voxels on the row dimension) takes about 2 minutes, as well.
%
% Now U has the temporal components and V has the spatial components.
% Following the example above, we'll take the rows of U as the principle
% temporal components and the columns of V as the principle spatial
% components.
figure(2)
% Look at the diagonal elements of S.  These usually fall off quite
% quickly.
imagesc(S(1:20,1:20));
%or 
plot(diag(S(2:end,2:end))); % skip the 1st since it's huge

% Then look at the first 5 principle components    
subplot(5,2,1); imagesc(reshape(V(:,1),[64 128])); axis image off; colorbar
subplot(5,2,2); plot(U(:,1)); title(S(1,1))
subplot(5,2,3); imagesc(reshape(V(:,2),[64 128])); axis image off; colorbar
subplot(5,2,4); plot(U(:,2)); title(S(2,2))
subplot(5,2,5); imagesc(reshape(V(:,3),[64 128])); axis image off; colorbar
subplot(5,2,6); plot(U(:,3)); title(S(3,3))
subplot(5,2,7); imagesc(reshape(V(:,4),[64 128])); axis image off; colorbar
subplot(5,2,8); plot(U(:,4)); title(S(4,4))
subplot(5,2,9); imagesc(reshape(V(:,5),[64 128])); axis image off; colorbar
subplot(5,2,10); plot(U(:,5)); title(S(5,5))

% It's not clear where the stimulus is.  But the first component is
% obvoiusly the magnitude map (coil sensitivity), and the 2nd component
% looks a whole lot like motion. 

% Repeat the exercise with normalized data (which is how we look at BOLD
% data, anyway).
normedData = reshape(slice,[64*128 160])';
normedData = normedData./repmat(mean(normedData,1),[160 1]) - 1;

tic; [U S V] = svd(normedData); toc
figure(3)
% Check out S
plot(diag(S))
% Then look at the first 5 principle components    
subplot(5,2,1); imagesc(reshape(V(:,1),[64 128])); axis image off; colorbar
subplot(5,2,2); plot(U(:,1)); title(S(1,1))
subplot(5,2,3); imagesc(reshape(V(:,2),[64 128])); axis image off; colorbar
subplot(5,2,4); plot(U(:,2)); title(S(2,2))
subplot(5,2,5); imagesc(reshape(V(:,3),[64 128])); axis image off; colorbar
subplot(5,2,6); plot(U(:,3)); title(S(3,3))
subplot(5,2,7); imagesc(reshape(V(:,4),[64 128])); axis image off; colorbar
subplot(5,2,8); plot(U(:,4)); title(S(4,4))
subplot(5,2,9); imagesc(reshape(V(:,5),[64 128])); axis image off; colorbar
subplot(5,2,10); plot(U(:,5)); title(S(5,5))

% I don't like how this turns out.  I checked and triple-checked that
% normedData didn't get messed up somehow, but PCA fails to find anything
% interesting here.

% One could reconstruct the original data without the first 2 components
% from the original analysis, and then go back to the Fourier analysis.
% By now, we've hogged my 2GB RAM, so need to clean up ... pack or
clear 
load slice
tic; [U S V] = svd(reshape(slice,[64*128 160])'); toc
cleanS = S; cleanS(1,1) = 0; cleanS(2,2) = 0;
cleanData = U*cleanS*V';
cleanData = reshape(cleanData',[64 128 160]);
f1 = fft(slice,[],3); % FT 3rd dimension
f2 = fft(cleanData,[],3); % FT 3rd dimension
% look at 10 cycles/scan component
figure(3)
subplot(2,2,1); imagesc(abs(f1(:,:,11))); axis image off; title('orig data')
subplot(2,2,3); imagesc(abs(f2(:,:,11))); axis image off; title('cleaned data')
% and plot a voxel
subplot(2,2,2); plot(squeeze(slice(51,52,:))); axis tight
subplot(2,2,4); plot(squeeze(cleanData(51,52,:))); axis tight

% This data, being high quality and already motion compensated, shows
% little change.  I wonder if other data would show a more clear difference
% ...



