for (c = 1 : numCluster)
    t = 0; 
    for (img = 1 : numImage)
       if (MAX2scoreAll(img, c) == max(MAX2scoreAll(img, :)))
           SUM1mapName = [outputFolder 'working/SUM1map' 'image' num2str(img)];   
           load(SUM1mapName); 
           t = t + 1;
           for (orient = 1 : numOrient)
               SUM1mapLearn{t, orient} = SUM1mapFind{originalResolution, orient} + single(0.); 
           end
       end
    end
    disp(['**************** start learning template ' num2str(c)]); tic; 
    CsharedSketch(numOrient, locationShiftLimit, orientShiftLimit, ... % about active basis  
           numElement, numImage, t, sizeTemplatex, sizeTemplatey, SUM1mapLearn, ... % about training images 
           halfFilterSize, Correlation, allSymbol(1, :), ... % about filters
           numStoredPoint, storedlambda, storedExpectation, storedLogZ, ... % about exponential model 
           selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ... % learned parameters
           commonTemplate, deformedTemplate); % learned templates 
    disp(['mex-C learning time: ' num2str(toc) ' seconds']);
    Outepsgif(-double(commonTemplate), ['nonalign_sym0' num2str(c)]); 
    RotateTemplate;   
    templateName = [outputFolder 'working/template' 'cluster' num2str(c)];   
    save(templateName, 'commonTemplate', 'allSelectedOrient', 'allSelectedx', 'allSelectedy', 'selectedlambda', 'selectedLogZ', 'centerx', 'centery');                    
end
%% Iteration part 1: detect the object in each image
for(it = 1:numIteration)
    close all; 
    disp(['**************** detection for iteration ' num2str(it)]);
    for (img=1:numImage)
      disp(['======> start detecting in image ' num2str(img)]); 
      MAX1mapName = [outputFolder 'working/MAX1map' 'image' num2str(img)]; 
      load(MAX1mapName); 
      for (c = 1 : numCluster) 
        %disp(['matching template ' num2str(c)]); tic; 
        templateName = [outputFolder 'working/template' 'cluster' num2str(c)];   
        load(templateName);                    
        MMAX2 = -1e10; 
        for (flip = 0 : flipOrNot)
          for (rot = -rotateShiftLimit : rotateShiftLimit)
            r = rot+rotateShiftLimit+1 + (rotateShiftLimit*2+1)*flip;  
            CgetSUM2(numResolution, allSizex, allSizey, numOrient, ...
                     numElement, allSelectedOrient(r, :), allSelectedx(r, :), allSelectedy(r, :), selectedlambda, selectedLogZ, ...
                     MAX1map, SUM2map, MAX2score, allFx, allFy, rangePercent);    
            [maxOverResolution ind] = max(MAX2score);   % most likely resolution
            if (MMAX2 < maxOverResolution) 
               MMAX2 = maxOverResolution; Mrot = rot; Mflip = flip; 
               Mind = ind; MFx = allFx(ind); MFy = allFy(ind); 
            end
          end
        end
        MAX2scoreAll(img, c) = MMAX2; 
        MrotAll(img, c) = Mrot; MflipAll(img, c) = Mflip; 
        MindAll(img, c) = Mind; MFxAll(img, c) = MFx; MFyAll(img, c) = MFy; 
        %disp(['mex-C finding time: ' num2str(toc) ' seconds']);
      end
    end
%% Iteration part 2: re-learn the template for each cluster
    for (c = 1 : numCluster)
      disp(['**************** starting cropping images for cluster ' num2str(c) ' in iteration ' num2str(it)]); tic; 
      templateName = [outputFolder 'working/template' 'cluster' num2str(c)];   
      load(templateName);       
      t = 0; 
      for (img = 1 : numImage)
        if (MAX2scoreAll(img, c) == max(MAX2scoreAll(img, :)))
           SUM1mapName = [outputFolder 'working/SUM1map' 'image' num2str(img)];   
           load(SUM1mapName); 
           Mrot = MrotAll(img, c); Mflip = MflipAll(img, c); 
           Mind = MindAll(img, c); MFx = MFxAll(img, c); MFy = MFyAll(img, c); 
           t = t + 1; 
           for (orient = 1 : numOrient)
              orient1 = orient - Mrot; 
              if (orient1 > numOrient)
                 orient1 = orient1 - numOrient; 
              end
              if (orient1 <= 0)
                 orient1 = orient1 + numOrient; 
              end                    
             Ccopy(SUM1mapLearn{t, orient}, SUM1mapFind{Mind, orient1}, MFx, MFy, centerx, centery, sizeTemplatex, sizeTemplatey, allSizex(Mind), allSizey(Mind), Mrot*pi/numOrient);  
           end
           if (Mflip > 0)
             for (orient = 1 : numOrient)
                 SUM1mapLearnTmp{1, orient} = SUM1mapLearn{t, orient}+0.; 
             end
             for (orient = 1 : numOrient)
                 if (orient-1>0)
                     orient2 = numOrient - (orient-1) + 1; 
                 else
                     orient2 = orient;   
                 end
                 Cflip(SUM1mapLearn{t, orient2}, SUM1mapLearnTmp{1, orient}, sizeTemplatex, sizeTemplatey); 
             end
           end           
           if (it == numIteration)
              r = Mrot+rotateShiftLimit+1 + (rotateShiftLimit*2+1)*Mflip;  
              CdrawTemplate(numResolution, allSizex, allSizey, numOrient, ...
                locationShiftLimit, orientShiftLimit, ...
                halfFilterSize, SUM1mapFind, allSymbol(1, :), ...
                numElement, allSelectedOrient(r, :), allSelectedx(r, :), allSelectedy(r, :), ...
                translatedTemplate, Mind, MFx, MFy);
              Jshow{img} = ImageMultiResolution{Mind}+translatedTemplate{Mind}*100.;  % overlaid sketch 
           end
        end
      end
      disp(['cropping time: ' num2str(toc) ' seconds']);
      disp(['**************** re-learning template ' num2str(c) ' in iteration ' num2str(it)]); tic; 
      CsharedSketch(numOrient, locationShiftLimit, orientShiftLimit, ... % about active basis  
           numElement, numImage, t, sizeTemplatex, sizeTemplatey, SUM1mapLearn, ... % about training images 
           halfFilterSize, Correlation, allSymbol(1, :), ... % about filters
           numStoredPoint, storedlambda, storedExpectation, storedLogZ, ... % about exponential model 
           selectedOrient, selectedx, selectedy, selectedlambda, selectedLogZ, ... % learned parameters
           commonTemplate, deformedTemplate); % learned templates 
    Outepsgif(-double(commonTemplate), ['nonalign_sym' 'iteration' num2str(it) 'cluster' num2str(c)]); 
    RotateTemplate;   
    templateName = [outputFolder 'working/template' 'cluster' num2str(c)];   
    save(templateName, 'commonTemplate', 'allSelectedOrient', 'allSelectedx', 'allSelectedy', 'selectedlambda', 'selectedLogZ', 'centerx', 'centery');                    
    disp(['mex-C learning time: ' num2str(toc) ' seconds']);
   end
end