%% demo_overdet_hyper_by2.m
%% 
%% Script to show the effect of using 21 wavelengths rather than 2
%% in solving single-layer hyperspectral unmixing problem from the lab
%% problem from HW 2, ECEU692, SP04.
%% 
%% Script first solves an overdetermined 21x2 system and plots the
%% results for the rectangular patch whose color varied vertically.
%% It then adds a 3rd column of 'grey' (all ones in its spectrum)' to
%% make a 21x3 system, solves it, and plots those results. Next it
%% chooses pairs of consecutive wavelengths to make a bunch of 2x2
%% systems and plots those results---user needs to hit return to step
%% through these plots. Finally it chooses pairs of wavelengths spaced
%% half-way through the 21, solves those 2x2 systems, and plots those
%% results, again requiring the user to hit return to step through the
%% results. Spectral coefficients for system matrix are obtained by
%% averaging blocks of pixels in the single-color regions.
%% 
%% Script assumes that
%%
%% data is loaded, R is reflection data cube
%% u = 1./ R -1 is computed
%% 
%% one pure color is in: u(121:250 171;220,:)
%% next pure color is in: u(300:430 171;220,:)
%% 
%% rectangular patch: 
%% one color: u(121:230,401:450,:)
%% mixed color: u(280:300,401:450,:)
%% one color: u(330:430,401:450,:)
%%
%% D.H. Brooks
%% 22 Feb 04
%% 

%% create 21x2 matrix from absorption data and solve that system for
%% relative concentration of chromophores.

A1(:,1) = squeeze(mean(mean(u(160:190,180:210,:))));
A1(:,2) = squeeze(mean(mean(u(340:370,180:210,:))));
x1 = A1\(squeeze(u(121:430,425,:)).');

%% Now add a 3rd column of all ones and solve for relative
%% concentrations of 3 chromophores from the same data

A2 =[A1 ones(21,1)];
x2 = A2\(squeeze(u(121:430,425,:)).');

%% plot the results so far

figure, plot(x1.')
title('LS, 21 wavelengths by 2 chromophores')
figure,plot(x2.')
title('LS, 21 wavelengths by 2 chromophores + grey')

%% now choose consecutive pairs of wavelengths, solve for chromophore
%% concentrations, and plot the results.

figure
for indx = 1:2:19
  A3 = A1(indx:indx+1,:);
  x_3consec = A3\(squeeze(u(121:430,425,indx:indx+1)).');
  plot(x_3consec.')
  title(['Using 2 consecutive wavelengths starting from number ' num2str(indx)])
  pause
end

%% and now repeat choosing wavelengths spaced 10 wavelengths apart.

figure
for indx = 1:10
  A3 = A1([indx indx+10],:);
  x_3distrib = A3\(squeeze(u(121:430,425,[indx indx+10])).');
  plot(x_3distrib.')
  title(['Using 2 distributed wavelengths starting from number ' num2str(indx)])
  pause
end