function featureArray = BatchExtractFeatures(imgCell,codeBook, para)

% if exist('./working','dir')%working files
%     delete('./working/*.mat'); 
% else
%     mkdir('working');
% end
% if exist('./sketch','dir')%image sketches
%     delete('./sketch/*.*'); 
% else
%     mkdir('sketch');
% end

featureArray=[];
%% set parameters

%ParameterCodeImage;
numResolution=para.numResolution;
allResolution=para.allResolution;
working_add=para.working_add;
numFilter=para.numFilter;
resizeTrainingImages=para.resizeTrainingImages;
minHeightOrWidth=para.minHeightOrWidth;
isGlobalNormalization = para.isGlobalNormalization;
sigsq=para.sigsq;
filters=para.filters;
halfFilterSizes=para.halfFilterSizes;
numOrient=para.numOrient;
locationShiftLimit=para.locationShiftLimit;
orientShiftLimit=para.orientShiftLimit;
isLocalNormalize=para.isLocalNormalize;
isSeparate=para.isSeparate;
localNormScaleFactor=para.localNormScaleFactor;
thresholdFactor=para.thresholdFactor;
nScaleGabor=para.nScaleGabor;
nScaleDoG=para.nScaleDoG;

numImage = numel(imgCell); % number of training images 

for img = 1 : numImage
    tmpIm = imgCell{img};
%      tmpIm = imread([imageFolder '/' imageName(img).name]); 
    if size(tmpIm,3) == 3
        tmpIm = im2single(rgb2gray(tmpIm));
    else
        tmpIm = im2single(tmpIm);
    end

    sx = size(tmpIm,1); 
    sy = size(tmpIm,2);
    if resizeTrainingImages
       resizeFactor = minHeightOrWidth/min(sx,sy);
       tmpIm = imresize( tmpIm, resizeFactor, 'bilinear' );
    end
    
%     if isGlobalNormalization
%        tmpIm=tmpIm-mean(tmpIm(:));
%        tmpIm = tmpIm/std(tmpIm(:))*sqrt(sigsq);
%     end
    
    Ioriginal{img} = tmpIm;
    J0 = Ioriginal{img};
    J = cell(1, numResolution);
    
    J_raw = cell(1, numResolution); %% new !!
    
    allSizex = zeros(1, numResolution); 
    allSizey = zeros(1, numResolution);
    for r=1:numResolution
       img_scaled = imresize(J0, allResolution(r), 'nearest');  % scaled images
       
       J_raw{r} = img_scaled; %% new !!
       
       if isGlobalNormalization
           img_scaled = img_scaled-mean(img_scaled(:));
           img_scaled = img_scaled/std(img_scaled(:))*sqrt(sigsq);
       end
       
       J{r} = img_scaled;
       [sizex, sizey] = size(J{r}); 
       allSizex(r) = sizex; 
       allSizey(r) = sizey; 
    end
    multipleResolutionImageName = [working_add '/multipleResolutionImage' num2str(img)]; 
    save(multipleResolutionImageName, 'J','J_raw','allSizex','allSizey');
end

CreatingFeaturesMaps;

[numTransform, numCluster]=size(codeBook);

 for i=1:numel(imgCell)
    disp(['starting to compute SUM2map for Images' num2str(i)]);
          
    ticID=tic;
          
    % compute SUM2
    SUM1MAX1mapName = [working_add '/SUM1MAX1map' 'image' num2str(i)];
    load(SUM1MAX1mapName, 'MAX1map', 'allSizex','allSizey');
          
    SUM2map = cell(numTransform*numCluster,numResolution);
          
    for iCluster=1:numCluster           
       for r=1:numTransform
          [SUM2] = sparseFRAME_SUM2_logZ( single(numResolution), single(allSizex), single(allSizey), single(numFilter), codeBook{r,iCluster}, MAX1map);  
          SUM2map(r+(iCluster-1)*numTransform,:)=SUM2;
       end
    end
   
            
    % exclude the near-boundary region
%     for ii = 1:numel(SUM2map)
%           SUM2map{ii}(:,[1:floor(templateSize(2)/4) end-floor(templateSize(2)/4):end]) = -1001; %min(-1,S2Thres-1);
%           SUM2map{ii}([1:floor(templateSize(1)/4) end-floor(templateSize(1)/4):end],:) = -1001; %min(-1,S2Thres-1);
%     end
                    
    disp(['SUM2map computational time: ' num2str(toc(ticID)) ' seconds']);           
    
    
    %% collapse SUM3map along resolution dimension
    
    num_total_features = size(SUM2map,1);
    denseFeatureMaps=cell(num_total_features,1);
    
    for i_feature = 1:size(SUM2map,1)
        localMax = SUM2map{i_feature,numResolution};
        for i_res = 1:numResolution
            tmpMap = imresize(SUM2map{i_feature,i_res},size(localMax),'nearest');
            localMax = max(localMax,tmpMap);
        end
        denseFeatureMaps{i_feature}=localMax;
    end
        
    %% SPM
    
    numLayers = 3;
    threshold = 0;
    SPM_feature = featureSPM(denseFeatureMaps,numLayers,threshold);
    featureArray = [featureArray; SPM_feature];    
         
end
        
 


end % end of function


