%function [X, whiteningMatrix, dewhiteningMatrix] = ...
%				data( samples, winsize, rdim );
samples = 40000;
winsize = 21;
% data - gathers patches from the gray-scale images
%
% INPUT variables:
% samples            total number of patches to take
% winsize            patch width in pixels
% rdim               reduced dimensionality
%
% OUTPUT variables:
% X                  the whitened data as column vectors
% whiteningMatrix    transformation of patch-space to X-space
% dewhiteningMatrix  inverse transformation
%

rand('seed', 0);
  
%----------------------------------------------------------------------
% Gather rectangular image windows
%----------------------------------------------------------------------

% We have a total of 13 images.
faceNum = 9;
imagesNum = [24 31 23 23 25 22 23 23 22];
imageInc = 2;
dataNum = 0;
for i=1:length(imagesNum),
 dataNum = dataNum+length(1:imageInc:imagesNum(i));
end
disp([sprintf('There are %d images in total to be sampled from.', dataNum)]);
% This is how many patches to take per image
getsample = floor(samples/dataNum);

% Initialize the matrix to hold the patches
X = zeros(winsize^2,getsample*dataNum);

sampleNum = 1;  
i = 1;
for jj=(1:faceNum)
     for k=(1:imageInc:imagesNum(jj)),
  % Even things out (take enough from last image)
  if i==dataNum, getsample = samples-sampleNum+1; end
  
  % Load the image
 % I = imread(['/brain2/liux/data/icabasis/textures/pami-21-4-11-h/Nat-10_' num2str(i) ...
  %  '.tif']);
 %I = imread(['/usr/home/liux/basis/textures/pami-21-4-11-i/Nat-10_' num2str(i) ...
  %  '.pgm']);

  imName =['/brain2/liux/data/ired/' char('A'+jj) ...
						   '/PP0000' num2str(k) '.PGM'];
  % I = imread(imName);
  I = readp5(imName);
  % Normalize to zero mean and unit variance
  I = double(I);
  I = I-mean(mean(I));
  I = I/sqrt(mean(mean(I.^2)));
  
  % Sample 
  fprintf('Sampling image %d (%s)...\n',i, imName);
  sizex = size(I,2); sizey = size(I,1);
  posx = floor(rand(1,getsample)*(sizex-winsize-2))+1;
  posy = floor(rand(1,getsample)*(sizey-winsize-1))+1;
  
  for j=1:getsample
    X(:,sampleNum) = reshape( I(posy(1,j):posy(1,j)+winsize-1, ...
			posx(1,j):posx(1,j)+winsize-1),[winsize^2 1]);
    sampleNum=sampleNum+1;
  end 
    i=i+1;
end
end











