disp(['runing on testing images']);
working_add='test_working';    % change working space
mkdir(working_add);
imageFolder = 'testImage'; % folder of training images  
imageName = dir([imageFolder '/*.jpg']);
numImage = size(imageName, 1); % number of training images 
Ioriginal = cell(1, numImage);
for img = 1 : numImage         % read images
    tmpIm = imread([imageFolder '/' imageName(img).name]); 
    if size(tmpIm,3) == 3
        tmpIm = rgb2gray(tmpIm);
    end

    sx = size(tmpIm,1); sy = size(tmpIm,2);
    if resizeTrainingImages
        if (resizeFactor < 0.)
            resizeFactor = sqrt(constantImageArea/(sx*sy))
        end
    	tmpIm = imresize( tmpIm, resizeFactor, 'bilinear' );
    end
    Ioriginal{img} = single(tmpIm);
    J0 = Ioriginal{img};
    J = cell(1, numResolution);
    for r=1:numResolution
       J{r} = imresize(J0, allResolution(r), 'nearest');  % scaled images
    end
    multipleResolutionImageName = [working_add '/multipleResolutionImage' num2str(img)]; 
    save(multipleResolutionImageName, 'J');
end
% calculate Sum1Max1
SUM1MAX1;
%
test_activations = []; % 3 by N matrix, where N is an unknown large number
for i = 1:numImage
    % compute SUM2
    SUM1MAX1mapName = [working_add '/SUM1MAX1map' 'image' num2str(i) 'scale' num2str(1)];
    load(SUM1MAX1mapName, 'SUM1map', 'MAX1map', 'M1RowShift', 'M1ColShift',...
                'M1OriShifted', 'J');
    for iRes = 1:numResolution
        SUM2map(:,iRes) = mexc_ComputeSUM2( numOrient, MAX1map(iRes,:), TransformedTemplate, subsampleS2, DELTA );
    end
    % random perturbation (to break ties arbitrarily for MAX2)
    for ii = 1:numel(SUM2map)
        SUM2map{ii}(:) = SUM2map{ii}(:) + 1e-3 * ( rand(numel(SUM2map{ii}),1) - .5 );
    end
    % exclude the near-boundary region
    for ii = 1:numel(SUM2map)
        SUM2map{ii}(:,[1:floor(templateSize(2)/4) end-floor(templateSize(2)/4):end]) = min(-1001,S2Thres-1);
        SUM2map{ii}([1:floor(templateSize(1)/4) end-floor(templateSize(1)/4):end],:) = min(-1001,S2Thres-1);
    end

    % compute MAX2, perform surround supression and get activations
    if iter == 1
        if strcmp('MatchingPursuit',supressionModeInEStep) == 1
            tmpActivations = mexc_ComputeMAX2MP( SUM2map, int32(locationPerturbationFraction*partSize/subsampleS2), -1000 );
        elseif strcmp('LocalSurroundSurpression',supressionModeInEStep) == 1
            subsampleM2 = 1;
            [MAX2map M2LocationTrace M2TemplateTrace M2RowColShift tmpActivations] = ...
                mexc_ComputeMAX2( templateAffinityMatrix, SUM2map, locationPerturbationFraction, ...
                int32(partSize/subsampleS2*ones(numCluster*nTransform,1)), subsampleM2 );
            tmpActivations = tmpActivations( :,tmpActivations(4,:) > -1000 );
        end
    else
        % discard the activated instances that have a low S2 score
        if iter > floor(numIter/2) % for the later iterations, increase the sparsity
            locationPerturbationFraction = locationPerturbationFraction_final;
        end
        if strcmp('MatchingPursuit',supressionModeInEStep) == 1
            tmpActivations = mexc_ComputeMAX2MP( SUM2map, int32(locationPerturbationFraction*partSize/subsampleS2), 0. );
        elseif strcmp('LocalSurroundSurpression',supressionModeInEStep) == 1
            subsampleM2 = 1;
            [MAX2map M2LocationTrace M2TemplateTrace M2RowColShift tmpActivations] = ...
                mexc_ComputeMAX2( templateAffinityMatrix, SUM2map, locationPerturbationFraction, ...
                int32(partSize/subsampleS2*ones(numCluster*nTransform,1)), subsampleM2 );
            tmpActivations = tmpActivations( :,tmpActivations(4,:) > S2Thres );
        end
    end
    test_activations = [test_activations,[single(i*ones(1,size(tmpActivations,2)));tmpActivations]];
end
test_activations(2:3,:) = test_activations(2:3,:) * subsampleS2;