function featureArray = BatchExtractFeatures(imgCell,codeBook)

% 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;
MakeFiltersBank;
CreatingFeaturesMaps;

(imgCell

 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(nTransform*numCluster,numResolution);
          
    for iCluster=1:numCluster           
       for r=1:numTransform
          [SUM2] = sparseFRAME_SUM2_logZ( single(numResolution), single(allSizex), single(allSizey), single(numFilter), AllTemplate{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:numImgResolution
            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


