%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.
imageInfo=[ 42 0268
 42 2291
 42 2319
 84 2326
 42 2328
 42 2342
 84 2344
 42 2347
 42 2351
 84 2358
 42 2360
 80 2374
 84 2434
 42 2435
 84 2464
 82 2465
 84 2475
 84 2492
 42 2513
 84 2521
 42 2531
 84 2550
 84 2569
 84 2590
 82 2591
 84 2592
 83 2593
 42 2595
 42 2596
 83 2597
126 2598
 84 2600
 42 2609
 42 2614
 84 2617
 84 2618
 42 2620
 42 2621
 42 2622
 84 2623
 42 2624
 42 2637
 42 2638
 42 2640
 42 2641
 84 2642
 42 2643
 42 2644
 42 2646
 42 2647
 42 2648
 84 2649
 42 2650
 84 2651
 84 2652
 42 2653
 42 2654
 81 2655
 84 2656
 42 2657
 42 2658
 84 2661
		42 2662];
imageInfo = imageInfo';
faceNum = size(imageInfo,2);
imagesNum = imageInfo(1,:);
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/equinox/longwave/' num2str(imageInfo(2,jj),'%04d') ...
						   '/L2-2-' 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 %d from image %d (%s)...\n', getsample, 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











