function find_S3(imageIDList,partOnOff) % if S3T is fixed, then the input should be the image, and just the index to S3T
% find_S3 - Find the best localization and deformation of S3 template [S3T] (shape script) on image [im].

load('working/S2Dict.mat','S2Dict','S2TNeighbors','canonicalS2IndByLoc', 'S2DictPaths', 'numCandPart', 'PartList');
load Config
save_harddisk_mode = true;
numOrient = 16;
load working/LearnedPart_1_1 TransformedTemplate
nTransform = length(TransformedTemplate);

%% compose the S3 template based on partOnOff
S3T = [];
selectedLogZ_S3 = [];
for iPartLoc = 1:length(S2DictPaths) % each .mat file contains transformed S2 templates 
    % write it into S3T
    if partOnOff(iPartLoc) > 0
        startx = PartList(1,iPartLoc);
        starty = PartList(2,iPartLoc);
        
        oneS3Element = single( [floor(startx+partSize/2 - outerBBSize/2-1);
            floor(starty+partSize/2-outerBBSize/2-1);
            canonicalS2IndByLoc(iPartLoc)-1] ); % to make it starting from 0
        S3T = [S3T, oneS3Element];
    end
end

% convert S3T to a detailed model
S3Template = cell(1,1);
nFS3 = size(S3T,2);
subsampleS2 = 1;
S3Template{1} = struct('selectedRow',single(floor(S3T(1,:)/subsampleS2)),'selectedCol',single(floor(S3T(2,:)/subsampleS2)),...
    'selectedInd',single(S3T(3,:)),'selectedTransform',single(zeros(nFS3,1)),'selectedLambda',single(ones(nFS3,1))...
    ,'selectedLogZ',single(zeros(nFS3,1)));


%% compute SUM MAX maps

for imageID = imageIDList
    
    disp( sprintf('processing image %d',imageID) );
    
    %% check if SUM1 MAX1 maps are computed
    pathSUM1MAX1map = sprintf('working/SUM1MAX1mapimage%d.mat',imageID);
    if ~exist(pathSUM1MAX1map,'file')
        disp('error: no SUM1 MAX1 map');
        return;
    end
    load(pathSUM1MAX1map,'MAX1map','M1RowShift','M1ColShift','M1OriShifted','ARGMAX1map','SUM1map');
    
    %% SUM2 MAX2
    S2M2mapfile = sprintf('working/SUM2MAX2mapimage%d.mat',imageID);
    if ~exist(S2M2mapfile,'file')

        SUM2map = cell( numCandPart * nTransform, 1 );
        for iPartLoc = 1:numCandPart
            startx = PartList(1,iPartLoc);
            starty = PartList(2,iPartLoc);
            load(sprintf( 'working/LearnedPart_%d_%d.mat', startx, starty ), 'TransformedTemplate');
            subsampleS2 = 1;
            SUM2map((iPartLoc-1)*nTransform+(1:nTransform), 1) = ...
                mexc_ComputeSUM2( numOrient, MAX1map, TransformedTemplate, subsampleS2 );
        end

        locationPerturbationFraction = double(partLocRange)/double(partSize);
        templateAffinityMatrix = cell( numCandPart * nTransform, 1 );
        for iPartLoc = 1:numCandPart
            from = (iPartLoc-1)*nTransform + 1;
            to = iPartLoc * nTransform;
            for jj = 1:nTransform
                templateAffinityMatrix{jj+(iPartLoc-1)*nTransform} = int32((from:to)-1);
            end
        end
        subsampleM2 = 1;
        [MAX2map M2LocationTrace M2TemplateTrace M2RowColShift] = ...
             mexc_ComputeMAX2( templateAffinityMatrix, SUM2map, locationPerturbationFraction, ...
                 int32(partSize*ones(numCandPart*nTransform,1)), subsampleM2 );

         % save S2Map M2Map
        save(S2M2mapfile, 'SUM2map', 'MAX2map', 'M2LocationTrace', 'M2TemplateTrace',...
            'M2RowColShift', 'templateAffinityMatrix');

    end

    load(S2M2mapfile, 'SUM2map', 'MAX2map', 'M2LocationTrace', 'M2TemplateTrace',...
            'M2RowColShift', 'templateAffinityMatrix');
        
    %% compute S3 map
    subsampleS3 = 1;
    SUM3map = mexc_ComputeSUM3( MAX2map, S3Template, subsampleS3, 1 );

    %% combine SUM/MAX maps for mexc_findS3.cpp
    feamapfile = sprintf('working/FeaMapImage%d.mat',imageID);
    % create ARGMAX2map
    ARGMAX2map = cell( size( MAX2map ) );
    for j = 1:numel(ARGMAX2map)
        S2CanonicalTemplateInd = floor(double(M2TemplateTrace{j}) / double(nTransform));
        S2TransformInd = double(M2TemplateTrace{j}) - S2CanonicalTemplateInd * double(nTransform);
        ARGMAX2map{j} = int32( double(M2LocationTrace{j}) + S2TransformInd * size(M2RowColShift,1) );
    end
    save(feamapfile,'SUM1map','MAX1map','SUM2map','MAX2map','ARGMAX1map','ARGMAX2map','SUM3map');

    %% deformation neighborhood for S2 templates
    if ~exist('working/S2TNeighbors.mat','file')
        S2TNeighbors = {};
        for fullIndS2 = 1:length(templateAffinityMatrix)
            SparseT = [];
            neighbors = templateAffinityMatrix{ fullIndS2 };
            % outer product of (rowshift,colshift) and neighbors
            for t = 1:length(neighbors)
                for jj = 1:size(M2RowColShift,1)
                    oneElement = [ single(M2RowColShift(jj,1) * partSize);...
                        single(M2RowColShift(jj,2) * partSize);...
                        single(neighbors(t)) ];
                    SparseT = [SparseT, oneElement]; % concatenate column vector
                end
            end
            S2TNeighbors{fullIndS2} = SparseT;
        end
        save('working/S2TNeighbors.mat','S2TNeighbors');
    else
        load('working/S2TNeighbors.mat','S2TNeighbors');
    end
    
    %% prepare "dictionaries" and "neighborhoods"
    dictionaries = cell(3,1);	% 3 is the number of levels
    dictionaries{1} = cell(16,1);		% first level are Gabor filters
    for j = 1:numel(S2Dict)
        S2Dict{j} = single(S2Dict{j});
    end
    dictionaries{2} = S2Dict;
    dictionaries{3} = {S3T};

    neighborhoods = cell(3,1);
    % create S1 neighborhood
    for iS1Filter = 1:numel(M1RowShift)
        nshift = numel( M1RowShift{iS1Filter} );
        NBS1{iS1Filter} = ...
            [single(reshape(M1RowShift{iS1Filter},1,nshift));...
            single(reshape(M1ColShift{iS1Filter},1,nshift));...
            single(reshape(M1OriShifted{iS1Filter},1,nshift))];
    end
    neighborhoods{1} = NBS1;
    neighborhoods{2} = S2TNeighbors;
    neighborhoods{3} = {};

    featureFolder = {'working'};
    paths = { {'positiveImage', sprintf('FeaMapImage%d.mat',imageID)} }; % the first item is not used 
    % for debug:
    % [DT1 DT2 DT3 dict2 nb2] = mexc_find_S3(paths,featureFolder,dictionaries,neighborhoods);
    % save(sprintf('working/DTimage%d.mat',imageID),'DT1','DT2','DT3','NBS1','dict2','nb2');
    [DT1 DT2 DT3] = mexc_FindS3(paths,featureFolder,dictionaries,neighborhoods);
    save(sprintf('working/DTimage%d.mat',imageID),'DT1','DT2','DT3','NBS1');

    %% render the symbolic template
    load Config outerBBSize partSize 
    nGaborOri = length(SUM1map);
    % shape motif level
    bbIm = zeros( partSize, partSize );
    margin = 5;
    bbIm([1:margin, partSize-margin+1:partSize], :) = 1;
    bbIm(:, [1:margin, partSize-margin+1:partSize]) = 1;
    bbIm(1:margin, (-margin:margin)+floor(partSize/2)) = 0;
    j = 1;
    rotationRange = -3:3;
    for angle = rotationRange
        transformedMotifBB{j} = imrotate( bbIm, angle * 180 / nGaborOri );
        j = j + 1;
    end
    MatchedMotifTemplate = zeros(outerBBSize,outerBBSize);
    for k = 1:length(DT2{1}(1,:))
        indNonRotated = floor((double(DT2{1}(3,k)) / double(nTransform) ) );
        indRotation = DT2{1}(3,k) - indNonRotated * double(nTransform) + 1;

        col = int32(DT2{1}(2,k) + 1);
        row = int32(DT2{1}(1,k) + 1);
        if row < 1 || row > size(SUM2map{DT2{1}(3,k)+1},1) || col < 1 || col > size(SUM2map{DT2{1}(3,k)+1},2)
        	continue;
        end
        val = SUM2map{DT2{1}(3,k)+1}( row, col );
        MatchedMotifTemplate = drawGaborSymbol( MatchedMotifTemplate,...
            transformedMotifBB, row, col, indRotation,...
            0, 0, max(min(20,val),0.1) );
    end

    % leaf level
    selectedRow = int32(DT1{1}(1,:));
    selectedCol = int32(DT1{1}(2,:));
    selectedOri = int32(DT1{1}(3,:));
    selectedS = zeros(size(DT1{1},2),1);
    selectedMean = ones(size(DT1{1},2),1);
    for j = 1:length(selectedRow)
        if selectedRow(j) + 1 < 1 || selectedRow(j) + 1 > size(SUM1map{1},1) ...
                || selectedCol(j) + 1 < 1 || selectedCol(j) + 1 > size(SUM1map{1},2)
            continue;
        end
        val = sqrt( SUM1map{selectedOri(j)+1}(selectedRow(j)+1,selectedCol(j)+1) );
        val = max(0, val-0.5);
        selectedMean(j) = min( 3, val );
    end
    load storedExponentialModel.mat allSymbol
    nGaborOri = numel(allSymbol);

    latticeSize = [outerBBSize outerBBSize];
    load(sprintf('working/multipleResolutionImage%d.mat',imageID),'J0');
    % inscribe
    srcImage = inscribe([outerBBSize outerBBSize], J0, 255);
    matchedSym = displayMatchedTemplate(latticeSize,selectedRow,...
        selectedCol,selectedOri, selectedS,selectedMean,allSymbol,nGaborOri);
    % overlay
    srcImage = uint8(255 * (srcImage-min(srcImage(:)))/(max(srcImage(:))-min(srcImage(:))));
    imwrite( srcImage, sprintf('output/matchedS3image%d.png',imageID)  );
    matchedSym = -double(matchedSym) + double(MatchedMotifTemplate) * 0;
    matchedSym = uint8(255 * (matchedSym-min(matchedSym(:)))/(max(matchedSym(:))-min(matchedSym(:))));
    % add a bounding frame
    matchedSym([1:5,end-4:end],:) = 0;
    matchedSym(:,[1:5,end-4:end]) = 0;
    imwrite( matchedSym, sprintf('output/matchedS3T%d.png',imageID)  );
    
    if save_harddisk_mode
    	delete(feamapfile);
        delete(S2M2mapfile);
    end
    
end



