% Define spatial component
Nx=32; Ny=Nx;Nvoxels=Nx*Ny;
thickness=18;
U=zeros(Nx,Ny);
x=[1:Nx]/Nx;
y=x.*(1-x);
xi=round(Nx*x);

for j=1:thickness,
    yi=round(Ny*y)+j+(ceil(Nx/8));
    for i=1:Nx,
        U(xi(i),yi(i))=1;
    end
end

% Vectorise spatial component
Uv=U(:)';

% Define temporal component
noise_dev=0.2;
s=boxcars(200,0.5,20);
Nscans=length(signal);
n=noise_dev*randn(Nscans,1);
cortex=s+n;

s2=boxcars(200,0.5,40);
n=noise_dev*randn(Nscans,1);
not_cortex=s2+n;

figure
subplot(2,1,1);
plot(cortex);
axis([0 Nscans -0.1 1.1]);
title('Cortical voxel');
subplot(2,1,2);
plot(not_cortex);
axis([0 Nscans -0.1 1.1]);
title('Non-cortical voxel');


M=s*Uv;
Uvcomp=ones(1,Nx*Ny)-Uv;
M=M+s2*Uvcomp;
M=M+noise_dev*randn(Nscans,Nvoxels);

figure
imagesc(U);
axis image
title('True spatial component');

figure
subplot(2,2,1);
imagesc(U);
axis image
title('True spatial component');

% Code from spm_regions.m
[m n]   = size(M);
[u s u] = svd(spm_atranspa(M'));
s       = diag(s);
u       = u(:,1);
v       = M'*u/sqrt(s(1));
d       = sign(sum(v));
u       = u*d;
v       = v*d;
Y       = u*sqrt(s(1)/n); % Scale by SD per voxel

% Get mean time series
mean_data=mean(M,2);

Cm=corrcoef(mean_data,cortex);
Ce=corrcoef(Y,cortex);

subplot(2,2,2);
plot(mean_data);
axis([0 Nscans -0.1 1.1]);
title(sprintf('Mean, r=%1.2f',Cm(1,2)));

subplot(2,2,4);
plot(Y);
axis([0 Nscans -0.1 1.1]);
title(sprintf('Eigenvariate, r=%1.2f',Ce(1,2)));

space_comp=v;
space_img=reshape(space_comp,Nx,Ny);
subplot(2,2,3);
imagesc(space_img);
axis image
title('Estimated spatial component');

ps=s./sum(s);
disp('First component explains:')
ps(1)

disp('Second component explains:')
ps(2)

disp('First two explain');
ps(1)+ps(2)