function encoder = trainEncoder(images, varargin)
% TRAINENCODER   Train image encoder: LASC
%   ENCODER = TRAINENCOER(IMAGES, MAXIMSIZE) trains a LASC encoder from the
%   specified list of images IMAGES.
%
%   TRAINENCODER(..., 'OPT', VAL, ...) accepts the following options:
%
%   Type:: 'lasc'
%     Locality-constrained Affine Subspace Coding.
%
%   numPcaDimension:: +inf
%     The subspace dimensionality. Use +inf to preserve all PCs.
% 
%   Subdivisions:: []
%     A list of spatial subdivisions. Each column is a rectangle
%     [XMIN YMIN XMAX YMAX WEIGHT]. The spatial subdivisions are
%
%   Layouts:: {'1x1'}
%     A list of strings representing regular spatial subdivisions
%     in the format MxN, where M is the number of vertical
%     subdivisions and N the number of horizontal ones. For
%     example {'1x1', 2x2'} uses 5 partitions: the whole image and
%     four quadrants. The subdivisions are appended to the ones
%     specified by the SUBDIVISIONS option.
%
%   ReadImageFn:: @readImage
%     The function used to load an image.
%
%   ExtractorFn:: @getDenseSIFT
%     The function used to extract the feature frames and
%     descriptors from an image.
% 
% Reference: VLFeat open source library (url: http://www.vlfeat.org/).

opts.type = 'lasc';
opts.seed = 1 ;
opts.numPcaDimensions = +inf ;
opts.numWords = 256 ;
opts.numSamplesPerWord = 4000 ;
opts.layouts = {'1x1'} ;
opts.subdivisions = zeros(5,0) ;
opts.readImageFn = @readImage ;
opts.extractorFn = @getDenseSIFT ;
opts = vl_argparse(opts, varargin) ;

pgrid = zeros(1,numel(opts.layouts));
for i = 1:numel(opts.layouts)
  t = sscanf(opts.layouts{i},'%dx%d') ;
  pgrid(i) = t(1) * t(2);
end
weights = (1./pgrid);
weights = weights/sum(weights);

for i = 1:numel(opts.layouts)
  t = sscanf(opts.layouts{i},'%dx%d') ;
  m = t(1) ;
  n = t(2) ;
  [x,y] = meshgrid( linspace(0,1,n+1), linspace(0,1,m+1) ) ;
  x1 = x(1:end-1,1:end-1) ;
  y1 = y(1:end-1,1:end-1) ;
  x2 = x(2:end,2:end) ;
  y2 = y(2:end,2:end) ;
  w = repmat(weights(i), 1, m*n);
  opts.subdivisions = cat(2, opts.subdivisions, [x1(:)'; y1(:)'; x2(:)'; y2(:)'; w(:)']) ;
end

disp(opts) ;
encoder.type = opts.type ;
encoder.numWords = opts.numWords ;
encoder.readImageFn = opts.readImageFn ;
encoder.extractorFn = opts.extractorFn ;
encoder.subdivisions = opts.subdivisions ;

%% Step 0 obtain sample image descriptors
numImages = numel(images) ;
numDescrsPerImage = ceil(opts.numWords * opts.numSamplesPerWord / numImages) ;
fprintf('%s: obtaining sample image descrptors...\n', mfilename) ;
parfor i = 1:numImages
  im = encoder.readImageFn(images{i}) ;
  features = encoder.extractorFn(im) ;
  randn('state',0) ;
  rand('state',0) ;
  sel = vl_colsubset(1:size(features.descr,2), single(numDescrsPerImage)) ;
  descrs{i} = features.descr(:,sel) ;
end
descrs = cat(2, descrs{:}) ;

%% Step 1: learn a k-means vocabulary
fprintf('%s: running k-means...\n', mfilename) ; 
vl_twister('state', opts.seed) ;
nzidx = sum(abs(descrs),1) > 0;
descrs = descrs(:,nzidx);
% running k-means by sp_kmeans %
options = foptions;
options(1) = 1;     % 1: display 0: no display
options(2) = 1;
options(3) = 1e-5;  % precision
options(5) = 1;     % initialization
options(14) = 10;   % maximum iterations
centers = zeros(opts.numWords, size(descrs,1));
[encoder.words,index] = sp_kmeans(centers, descrs', options);
encoder.words = encoder.words';

%% Step 2: build affine subspaces of the k-means dictionary
encoder.pca = SubspacePCA(encoder.words, min(size(descrs,1),opts.numPcaDimensions), descrs, index);

