function c1Database = build_c1_layer(s1Database,param)
% This function is to extract mid-level features with max pooling followed 
% by spatial grouping.
% Input:   s1Database -- the database of low-level sparse feature maps
%          param      -- see dedault settings below
% Output:
% Written by Jimei Yang
% Oct. 28, 2009
if ~isfield(param,'poolingS1')
    param.poolingS1 = 4;
end
if ~isfield(param,'biasC1')
    param.biasC1 = 1;
end
if ~isfield(param,'cellSizeC1')
    param.cellSizeC1 = 1;
end

c1Database = s1Database;
nSmp = length(s1Database.path);

for i = 1:nSmp
    fprintf('--extract the c1 feature...');
    tic;
    fpath = s1Database.path{i};
    load(fpath);
    bSize = size(s1Fea.feaData,1);
    
    % extract c1 features

    % find to which block (bin) each feature belongs
    xBin = ceil(s1Fea.x / param.poolingS1);
    yBin = ceil(s1Fea.y / param.poolingS1);
    idxBin = (yBin - 1)*max(xBin) + xBin;
    tBins = max(idxBin);
    feaData1 = zeros(bSize,tBins);

    for k = 1:tBins
        sidxBin = find(idxBin == k);
        if isempty(sidxBin)
            continue;
        end
        % max pooling of absolute value
%         feaData(:,k) = max(abs(s1Fea.feaData(:, sidxBin)), [], 2);
        tmp = max(abs(s1Fea.feaData(:, sidxBin)), [], 2);
%         tmp = tanh(tmp);
        fnorm = norm(tmp);
        if fnorm > 0
            feaData1(:,k) = tmp/fnorm;
        else
            feaData1(:,k) = tmp;
        end
    end
    [gridX,gridY] = meshgrid(unique(xBin),unique(yBin));
    gridX = gridX(:);
    gridY = gridY(:);
    width = max(gridX);
    height = max(gridY);
    rangeX = 1:param.biasC1:width-param.cellSizeC1+1;
    rangeY = 1:param.biasC1:height-param.cellSizeC1+1;
    feaData2 = zeros(bSize*param.cellSizeC1*param.cellSizeC1,length(rangeX)*length(rangeY));
    fCount = 1;
    for ii = rangeX
        for jj = rangeY
            idx = find(gridX>=ii&gridX<ii+param.cellSizeC1&gridY>=jj&gridY<jj+param.cellSizeC1);
            tmp = feaData1(:,idx);
            tmp = tmp(:);
            fnorm = norm(tmp);
            if fnorm > 0
                feaData2(:,fCount) = tmp/fnorm;
                fCount = fCount+1;
            else
                feaData2(:,fCount) = tmp;
                fCount = fCount+1;
            end
        end
    end
    c1Fea.feaData = feaData2; 
    [gridX,gridY] = meshgrid(rangeX,rangeY);
    c1Fea.x = gridX(:);
    c1Fea.y = gridY(:);
    c1Fea.width = length(rangeX);
    c1Fea.height = length(rangeY);

    
    % save c1Database
    [pdir,fname] = fileparts(s1Database.path{i});
    seps = findstr(pdir,'\');
    fdir = ['data\c1_',pdir(seps(1)+1:end)];
    if ~isdir(fdir)
        mkdir(fdir);
    end
    fpath = fullfile(fdir,[fname,'.mat']);
    save(fpath,'c1Fea');
    c1Database.path{i} = fpath;
    disp([num2str(i),' in ',num2str(toc),' seconds;']);
end
