function [V, D]=pca_image(input_file, exclude, p, ...
    mask_file, mask_thresh, X_remove, X_interest, output_file_base)

%PCA_IMAGE PCA of a multivariate linear model for image data 
% as described in 
% Worsley, K.J., Poline, J.B., Friston, K.J. and Evans, A.C. (1998). 
% Characterizing the response of PET and fMRI data using 
% Multivariate Linear Models (MLM). NeuroImage, 6:305-319,
% but ignoring the temporal correlation. 
%
% The temporal or file components are normalized so that their sum of
% squares is 1, and the spatial components are divided by their sd
% over voxels. Produces nice figures of time and space components if an 
% output file base is not given. Slices are numbered starting at 0.
% 
% [V, D]=PCA_IMAGE(INPUT_FILE [, EXCLUDE [,P [, MASK_FILE [, MASK_THRESH 
%         [, X_REMOVE [, X_INTEREST [, OUTPUT_FILE_BASE ]]]]]]])
%
% INPUT_FILE is the name of a single image file with multiple frames, 
% e.g. fMRI data, or a matrix of image file names, each with a single,
% frame, e.g. effect files from multiple runs, sessions or subjects.
%
% EXCLUDE is a row vector of frames or files to be excluded, empty by default.
%
% P is number of desired components, not more than numX, the number
% of frames or files. Default is 1.
%
% MASK_FILE is a mask file. If empty, it is ignored. Default is [].
%
% MASK_THRESH defines the search volume as the first frame of MASK_FILE 
% > MASK_THRESH. If MASK_THRESH is a vector [a b], a<=b, then mask 
% is a < MASK_FILE <= b. Default is 0.
%
% X_REMOVE is a design matrix with a row for each frame or file of 
% the covariates to be removed from the data before doing the PCA. 
% Default is ones(numX,1), i.e. removing the mean image over time or 
% files; use [ones(numX,1) (1:numX)'] to remove a linear drift as well.
%
% X_INTEREST is a design matrix with a row for each frame or file of 
% the covariates of interest. The PCA is done on the effects for these 
% covariates. Default is eye(numX), i.e. a PCA on all the frames or  
% files; use repmat(eye(12),numX/12,1) for a PCA on the average time 
% course of each block of 12 consecutive frames. 
%
% OUTPUT_FILE_BASE_pca.mnc or .img is output file of the first P 
% spatial components as frames, none if empty (default), in which 
% case a figure of the temporal and spatial components is produced.
%
% V is numframes x P matrix of the first P temporal or file components.
%
% D is a diagonal matrix of all the eigen values of the sum of squares matrix. 

%############################################################################
% COPYRIGHT:   Copyright 2002 K.J. Worsley, 
%              Department of Mathematics and Statistics,
%              McConnell Brain Imaging Center, 
%              Montreal Neurological Institute,
%              McGill University, Montreal, Quebec, Canada. 
%              worsley@math.mcgill.ca
%
%              Permission to use, copy, modify, and distribute this
%              software and its documentation for any purpose and without
%              fee is hereby granted, provided that this copyright
%              notice appears in all copies. The author and McGill University
%              make no representations about the suitability of this
%              software for any purpose.  It is provided "as is" without
%              express or implied warranty.
%############################################################################

% Defaults:

if nargin < 2
   exclude=[]
end
if nargin < 3
   p=1
end
if nargin < 4
   mask_file=[]
end
if nargin < 5
   mask_thresh=0
end
if nargin < 8
   output_file_base=[]
end

numfiles=size(input_file,1)
d=fmris_read_image(input_file(1,:),0,0);
numframes=d.dim(4)
numslices=d.dim(3)
numys=d.dim(2)
numxs=d.dim(1)
numpix=numxs*numys;
file_path=d.file_path;
ext=input_file(1,(size(input_file,2)-2):size(input_file,2));
isanalyze= all(ext=='img')

if ~isempty(mask_file)
   mask_thresh1=mask_thresh(1);
   if length(mask_thresh)>=2
      mask_thresh2=mask_thresh(2);
   else
      mask_thresh2=Inf;
   end
   d=fmris_read_image(mask_file,1:numslices,1);
   mask= (d.data>mask_thresh1 & d.data<=mask_thresh2);
   xf=find(max(max(mask,[],2),[],3));
   xr=min(xf):max(xf);
   yf=find(max(max(mask,[],1),[],3));
   yr=min(yf):max(yf);
   mask=reshape(mask,numpix,numslices);
   N=sum(sum(mask));
else
   xr=1:numxs;
   yr=1:numys;
   N=prod(d.dim(1:3));
end

numX=numfiles*(numfiles>1)+numframes*(numfiles==1);
allpts = 1:numX;
allpts(exclude) = zeros(1,length(exclude));
keep = allpts( find( allpts ) );
n=length(keep)

if ~isempty(output_file_base)
   if strcmp(computer,'PCWIN')
      if strcmp(output_file_base(2:3),':\');file_path='';end
   else
      if strcmp(output_file_base(1),'/');file_path='';end
   end
   out_svd.file_name=[deblank(output_file_base) '_pca.' ext];
   out_svd.file_path=file_path;
   out_svd.dim=[numxs numys numslices p];
   out_svd.parent_file=deblank(input_file(1,:));
   out_svd.precision='float'; 
   out_svd.descrip='';
end

if nargin < 6
   X_remove=ones(numX,1);
end
if nargin < 7
   X_interest=eye(numX);
end

X=X_interest(keep,:);
if ~isempty(X_remove)
   Z=X_remove(keep,:);
   XZ=X-Z*pinv(Z)*X;
else
   XZ=X;
end
[UX,SX,VX]=svd(XZ);
n1=rank(XZ);
UX=UX(:,1:n1);

A=zeros(n1);
for slice=1:numslices
   slice
   if numfiles==1
      d=fmris_read_image(input_file,slice,1:numframes);
      Y=reshape(d.data(:,:,keep),numpix,n);
   else
      Y=zeros(numpix,n);
      for i=1:n
         d=fmris_read_image(input_file(keep(i),:),slice,1);
         Y(:,i)=d.data(:);
      end
   end
   if ~isempty(mask_file)
      Y=Y(mask(:,slice),:);
   end
   YX=Y*UX;
   A=A+YX'*YX;
end
[Vs, D]=eig(A);
[ds,is]=sort(-diag(D));
ds=-ds;
sd=sqrt(ds(1:p)/N)
pcntvar=ds(1:p)/sum(ds)*100

VX=UX*Vs(:,is(1:p));
U=zeros(numxs,numys,numslices,p);
for slice=1:numslices
   slice
   if numfiles==1
      d=fmris_read_image(input_file,slice,1:numframes);
      Y=reshape(d.data(:,:,keep),numpix,n);
   else
      Y=zeros(numpix,n);
      for i=1:n
         d=fmris_read_image(input_file(keep(i),:),slice,1);
         Y(:,i)=d.data(:);
      end
   end
   if ~isempty(mask_file)
      Y=Y.*repmat(mask(:,slice),1,n);
   end
   u=Y*VX;
   out_svd.data(:,:,slice,:)=reshape(u,[numxs numys p]);
end

% change the signs so that max=max(abs) of each image:
umins=squeeze(min(min(min(out_svd.data,[],1),[],2),[],3));
umaxs=squeeze(max(max(max(out_svd.data,[],1),[],2),[],3));
ok=(umaxs>-umins);
for i=1:p
   if ~ok(i)
      out_svd.data(:,:,:,i)=-out_svd.data(:,:,:,i)/sd(i);
      VX(:,i)=-VX(:,i);
   else
      out_svd.data(:,:,:,i)=out_svd.data(:,:,:,i)/sd(i);
   end
end

if ~isempty(mask_file)
   mask=reshape(mask,[numxs,numys,numslices]);
   mask(mask==0)=NaN;
   for i=1:p
      out_svd.data(:,:,:,i)=out_svd.data(:,:,:,i).*mask;
   end
end

end
V=zeros(numX,p);
V(keep,1:p)=VX;

if isempty(output_file_base)
   subplot(2,1,1);
   a=max(abs(VX),[],1);
   plot([0 numX],repmat(-(1:p),2,1),'k'); 
   hold on; plot(keep,VX./repmat(a,n,1)*0.5+repmat(-(1:p),n,1)); hold off;
   xlim([0 numX]*1.17);
   ylabel('-Component'); 
   if numfiles==1
      xlabel('Frame');
      title('Temporal components (sd, % variance explained)');
   else
      xlabel('File');
      title('File components (sd, % variance explained)');
   end
   for i=1:p
      text(numX,-i,[num2str(round(sd(i)*10)/10), ', ', ...
            num2str(round(pcntvar(i)*10)/10) '%']);
   end
   
   if d.vox(1)<0 & isanalyze
      xr=flipdim(xr,2);
   end
   if d.vox(2)<0 & isanalyze
      yr=flipdim(yr,2);
   end
   
   subplot(2,1,2);
   imagesc(-0.5:numslices-0.5,0.5:p+0.5, ...
      flipud(reshape(permute(out_svd.data(xr,yr,:,p:-1:1),[1 3 2 4]), ...
      length(xr)*numslices,length(yr)*p)'));
   xlabel('Slice'); ylabel('Component');
   title('Spatial components (in sd units)'); colorbar;
else
   fmris_write_image(out_svd,1:numslices,1:p); 
end

return