%% The Use of Random Projections for the Analysis of MALDI Mass Spectrometry Imaging Data
%% Andrew D. Palmer, Josephine Bunch, Iain B. Styles

%% Introduction to File Processing
% This is a script to perform Random Projection on
% hyperspectral imaging data
% The data in this example is publically available in full at
% imzmlconverter.co.uk (AM Race, IB Styles, J Bunch,Inclusive sharing of mass spectrometry imaging data requires a converter for all Journal of proteomics 75 (16), 2012, 5111-5112)
% 
% This code follows the algorithm steps presented in the paper and has been
% written with clarity in mind, rather than computational efficiency.
%% Load Data
% Sample data is provided in the supporting information in the form of a 2D
% matrix containing spectra at each pixel and acompanying metadata.
% Change this directory if this script and the data are in different
% directories
fileToLoad = ['sampleData_rebinned.mat'];
% Load sample data
load(fileToLoad)
% Hyperspectral data is loaded as a 2D matrix, X, with n rows and m
% columns (n = number of channels, m = number of pixels)
[nChannels, mPixels] = size(X);
% The spectral labellings, channelValues, is in the metadata
channelValues = fileInfo.mzAxis;
channelType = 'm/z';
% The image size is also in the metadata
nRows = fileInfo.nRows;
nColumns = fileInfo.nColumns;
%% Perform Compression
% choose k, the number of samplings 
k = 100; 
% generate a random matrix
randomMatrix = randn(k,nChannels);
% sample the data with the random matrix
S = randomMatrix*X;
%% Segmentation 
nClusters = 5;
[idx, C] = kmeans(S',nClusters);
figure, 
subplot(nClusters,2,[2 4])
imagesc(reshape(idx,nRows,nColumns)), axis image, axis off
colorbar
title('Segmentation map')
for c = 1 : nClusters
    subplot(nClusters,2,2*(c-1)+1)
    plot(channelValues,mean(X(:,idx==c)'))
    title(['Molecular Profile from Cluster ' num2str(c)])
    xlabel('m/z')
    ylabel('Mean Intensity')
end
set(gcf,'Color','w')