GaborScaleListLastIteration = para.GaborScaleListLastIteration;
DoGScaleListLastIteration = para.DoGScaleListLastIteration;
isLocalNormalizeLastIteration = para.isLocalNormalizeLastIteration;
numWaveletLastIteration = para.numWaveletLastIteration;

if isLocalNormalizeLastIteration
    DoGScaleListLastIteration=[];
end

nScaleGabor=length(GaborScaleListLastIteration);
nScaleDoG=length(DoGScaleListLastIteration);


%%% create filter bank
filters=[];
for iScale=1:nScaleGabor
    f = MakeFilter(GaborScaleListLastIteration(iScale),nOrient);
    for i=1:nOrient
        f_r{i} =single(real(f{i}));
        f_i{i} =single(imag(f{i}));
    end
    
    filters = [filters f_r f_i];
end

for iScale=1:nScaleDoG
    
    f0 = single(dog(DoGScaleListLastIteration(iScale),0));
    filters = [filters f0];
    
end

numFilter = length(filters);
halfFilterSizes = zeros(size(filters));
for iF = 1:numFilter
    halfFilterSizes(iF)=(size(filters{iF},1)-1)/2;
end
%overAllArea = sum((sx-2*halfFilterSizes).*(sy-2*halfFilterSizes));
Corr = CorrFilterFrame(filters); %calculate correlation among filters

%%%%%%%% optional: for drawing sketch in the learning process
[allFilterR, ~, filterSymbol] = MakeFilterBank(GaborScaleListLastIteration, DoGScaleListLastIteration, nOrient); % filter bank for function "drawSketch"
half = zeros(1, nOrient*nScaleGabor+nScaleDoG);
for i=1:(nOrient*nScaleGabor+nScaleDoG)
    half(i) = (size(allFilterR{1, i}, 1)-1)/2;
end
%%%%%%%%%%

for iFilter = 1:numFilter
    model.rHat{iFilter}=zeros(sx, sy,'single');
end

for iImg = 1:numImage
    
    tic
    % load cropped images for training
    imageLoaded = load(fullfile(detectedCroppedSavingFolder,['detected-img-' num2str(iImg,'%04d')]), 'cropedImage');
    
    cropedImage = imageLoaded.cropedImage;
    
    cropedImage = cropedImage - mean(cropedImage(:));
    cropedImage = cropedImage/std(cropedImage(:))*sqrt(sigsq);
    
    
    % compute feature map to learn
    [~,MAX1] = applyfilterBank_MultiResolution_sparseV4({cropedImage}, filters, halfFilterSizes, nOrient, locationShiftLimit,orientShiftLimit,...
        isLocalNormalizeLastIteration,isSeparate,localNormScaleFactor,thresholdFactor,nScaleGabor,nScaleDoG, sqrt(sigsq));  % if using local normalization, the last parameter is important.
        
    % sum over the observed statistics (within cluster)
    for iFilter = 1:numFilter
        model.rHat{iFilter}=model.rHat{iFilter}+MAX1{iFilter}/numImage;
    end
    
    disp(['cropping time for image ' num2str(iImg) ': ' num2str(toc) ' seconds']);
end% iImage

savingFolder=fullfile(resultPath,['iteration' num2str(it+1) '_refine/']);
if ~exist(savingFolder)
    mkdir(savingFolder);
end

clusterSavingFolder=fullfile(savingFolder,'model/');
if ~exist(clusterSavingFolder)
    mkdir(clusterSavingFolder);
end

[template, currSample, logZ]=sparseFRAMElearnGibbs_multipleSelection(filters, nScaleGabor, nScaleDoG, nOrient, filterSymbol, half, model.rHat, ...
    sx, sy, halfFilterSizes, locationShiftLimit, nTileRow, nTileCol, lambdaLearningRate, numFilter, numWaveletLastIteration, interval, Corr, threshold_corrBB, c_val_list, ...
    lower_bound_rand, upper_bound_rand, nPartCol, nPartRow, part_sx, part_sy, gradient_threshold_scale, clusterSavingFolder);

