function patches = sampleIMAGES(IMAGES,patchsize,norm_ena)
% sampleIMAGES
% inputs
% numpatches: for example 1E4
% patchsize% for example 8*8
% IMAGES: images in a 3D matrix
% output
% patches % a vector of randomly chosen patches
if nargin==2
    norm_ena=1;
end
visibleSize = patchsize*patchsize;   % number of input units 

% get size and number of images
[xn, yn, zn]=size(IMAGES);

scale=patchsize/xn;
imgs = imresize(IMAGES, scale);
patches=(reshape(imgs,visibleSize,zn));

%% ---------------------------------------------------------------
% For the autoencoder to work well we need to normalize the data
% Specifically, since the output of the network is bounded between [0,1]
% (due to the sigmoid activation function), we have to make sure 
% the range of pixel values is also bounded between [0,1]
if norm_ena==1
patches = normalizeData(patches);
end
end

%% ---------------------------------------------------------------
function patches = normalizeData(patches)

% Squash data to [0.1, 0.9] since we use sigmoid as the activation
% function in the output layer

% Remove DC (mean of images). 
patches = bsxfun(@minus, patches, mean(patches));

% Truncate to +/-3 standard deviations and scale to -1 to 1
pstd = 3 * std(patches(:));
patches = max(min(patches, pstd), -pstd) / pstd;

% Rescale from [-1,1] to [0.1,0.9]
patches = (patches + 1) * 0.4 + 0.1;

end