% *****************************************************
% display activated templates on selected images
% *****************************************************

storeExponentialModelName = ['storedExponentialModel' num2str(1)];   
load(storeExponentialModelName);


colors = colormap(hsv(numCluster));

showPartBoundingBox = true;
partSizeX = templateSize(1);
partSizeY = templateSize(2);

%% preparation
% transform the templates
S2Templates = bestS2Templates;
TransformedTemplate = cell(nTransform,numCluster);
selectedScale = zeros(1,numElement,'single');
for cc = 1:numCluster
	for iT = 1:nTransform
		templateScaleInd = templateTransform{iT}(1);
		rowScale = templateTransform{iT}(2);
		colScale = templateTransform{iT}(3);
		rotation = templateTransform{iT}(4);
		[tmpSelectedRow tmpSelectedCol tmpSelectedOri tmpSelectedScale] = ...
			mexc_TemplateAffineTransform( templateScaleInd, rowScale,...
			colScale, rotation, S2Templates{cc}.selectedRow, S2Templates{cc}.selectedCol,...
			S2Templates{cc}.selectedOri, selectedScale, numOrient );
		TransformedTemplate{iT,cc}.selectedRow = tmpSelectedRow;
		TransformedTemplate{iT,cc}.selectedCol = tmpSelectedCol;
		TransformedTemplate{iT,cc}.selectedOri = tmpSelectedOri;
		TransformedTemplate{iT,cc}.selectedScale = tmpSelectedScale;
		TransformedTemplate{iT,cc}.selectedLambda = S2Templates{cc}.selectedLambda;
		TransformedTemplate{iT,cc}.selectedLogZ = S2Templates{cc}.selectedLogZ;
	end
end



%% begin displaying
% pruning invalid assignments due to empty clusters
iCluster = ceil( ( activations(5,:) + 1 ) / nTransform );
invClu=ismember(iCluster,non_empty_cluster);  
inds=find(invClu==1);
activations=activations(:,inds);

buf_length = 0;
for img = selected_img

	for k = 1:buf_length
		fprintf(1,'\b');
	end
	str = sprintf('%d',img);
	fprintf(1,str);
	buf_length = length(str);

	% load SUM1/MAX1 map
	SUM1MAX1mapName = [working_add '/SUM1MAX1map' 'image' num2str(img) 'scale' num2str(1)];
	load(SUM1MAX1mapName, 'SUM1map', 'MAX1map', 'ARGMAX1map', 'M1RowShift', 'M1ColShift',...
		'M1OriShifted', 'J');
	% find activations
	ind = find( activations(1,:) == img );
	
	% initialize the mask of bounding boxes and sketched image
	matchedBoundingBox = zeros([size(J{end}),3]); % the bounding box map is of the highest resolution
	matchedSym = zeros(size(J{end}),'uint8'); % the sketched image is of the highest resolution
    colorSketch = 254*ones([size(J{end}),3],'uint8');
	matchedArea = zeros(size(J{end}),'uint8'); % areas to containing templates
    
    % Gabor basis elements locations
    for t = ind(end:-1:1) % for each activated template (or codeword, or part)
	
		% track Gabor elements of activated templates for this template only
		gaborXX = [];
		gaborYY = [];
		gaborOO = [];
		gaborMM = [];
	
    	iCluster = ceil( ( activations(5,t) + 1 ) / nTransform );
    	iTransform = activations(5,t) + 1 - (iCluster-1) * nTransform;
		iRes = activations(4,t) + 1;
    	for j = 1:numElementAll(iCluster)
		    gaborX = floor(activations(2,t) + TransformedTemplate{iTransform,iCluster}.selectedRow(j));
		    gaborY = floor(activations(3,t) + TransformedTemplate{iTransform,iCluster}.selectedCol(j));
		    gaborO = TransformedTemplate{iTransform,iCluster}.selectedOri(j);
		    if gaborX > 0 && gaborX <= size(MAX1map{iRes,1},1) && gaborY > 0 && gaborY <= size(MAX1map{iRes,1},2)
		        trace = ARGMAX1map{iRes,gaborO+1}(gaborX,gaborY) + 1;
		        dx = M1RowShift{gaborO+1}(trace);
		        dy = M1ColShift{gaborO+1}(trace);
		        shiftedo = M1OriShifted{gaborO+1}(trace);
		        gaborX = floor(.5 + gaborX + single(dx));
		        gaborY = floor(.5 + gaborY + single(dy));
		        gaborO = single(shiftedo);
		    end
		    gaborXX = [gaborXX;gaborX];
		    gaborYY = [gaborYY;gaborY];
		    gaborOO = [gaborOO;gaborO];
		    if gaborX > 0 && gaborX <= size(MAX1map{iRes,1},1) && gaborY > 0 && gaborY <= size(MAX1map{iRes,1},2)
		        val = SUM1map{iRes,gaborO+1}(gaborX,gaborY);
		    else
		        val = 0;
		    end
		    gaborMM = [gaborMM; max(0,sqrt(val)-.2)];
			gaborMM = [gaborMM; 1];
        end
		
		% render the sketches for this activated template
		tmpMatchedSym = displayMatchedTemplate([size(J{iRes},1) size(J{iRes},2)],gaborXX,...
			gaborYY,gaborOO,zeros(length(gaborXX),1,'single'),gaborMM,allSymbol,numOrient);
        
		scaling = double(size(J{end},1)) / double(size(J{iRes},1));
		tmpMatchedSym = imresize(tmpMatchedSym,size(matchedSym),'nearest');
		matchedSym = max( matchedSym, tmpMatchedSym );
        tmpColorSketch = repmat(double(tmpMatchedSym),[1 1 3]);
        ind = find(tmpColorSketch);
        temid=find(non_empty_cluster(tempOrder)==iCluster);
        for cc = 1:3
            tmpColorSketch(:,:,cc) = tmpColorSketch(:,:,cc) * colors(temid,cc);
        end
        colorSketch(ind) = uint8(tmpColorSketch(ind));
		
        if showPartBoundingBox
		    margin = 2;
			largerPartSizeX = floor(partSizeX * scaling+.5);
			largerPartSizeY = floor(partSizeY * scaling+.5);
		    xx = repmat((1:largerPartSizeX),1,margin*2);
		    yy = [];
		    for y = [1:margin largerPartSizeY-margin+1:largerPartSizeY]
		        yy = [yy,ones(1,largerPartSizeX)*y];
		    end
		    yy = [yy,repmat((1:largerPartSizeY),1,margin*2)];
		    for x = [1:margin largerPartSizeX-margin+1:largerPartSizeX]
		        xx = [xx,ones(1,largerPartSizeY)*x];
		    end
		    inRow = single(xx-floor(largerPartSizeX/2)); inCol = single(yy-floor(largerPartSizeY/2));
		    tScale = 0; rScale = 1; cScale = 1; inO = zeros(numel(inRow),1,'single'); inS = zeros(numel(inRow),1,'single');
		   % actualPartRotation = iTransform-1;
           actualPartRotation = rotationRange(iTransform);
		    [outRow, outCol] = ...
		        mexc_TemplateAffineTransform(tScale,rScale,cScale,...
		            actualPartRotation,inRow,inCol,inO,inS,numOrient);

		    % directly overwrite the corresponding pixels
            pos=[];
		    for p = 1:length(outRow)
		        x = floor(.5 + outRow(p) + activations(2,t)*scaling); y = floor(.5 + outCol(p) + activations(3,t)*scaling);
		        if x > 0 && x <= size(matchedBoundingBox,1) && y > 0 && y <= size(matchedBoundingBox,2)
		            matchedBoundingBox(x,y,:) = [.5 .9 .6];
                    pos=[pos; x y];
                end
                if x<=0  && y > 0 && y <= size(matchedBoundingBox,2)
                    pos=[pos; 1 y];
                end;
                 if x > size(matchedBoundingBox,1)  && y > 0 && y <= size(matchedBoundingBox,2)
                    pos=[pos; size(matchedBoundingBox,1) y];
                end;
                if x > 0 && x <= size(matchedBoundingBox,1) && y <= 0
                    pos=[pos; x 1];
                end;
                if x > 0 && x <= size(matchedBoundingBox,1) && y > size(matchedBoundingBox,2)
                    pos=[pos; x size(matchedBoundingBox,2)];
                end;
            end
            minx=min(pos(:,1));
            maxx=max(pos(:,1));
            for minii=minx:maxx
                inds=find(pos(:,1)==minii);
                matchedArea(minii,min(pos(inds,2)): max(pos(inds,2)))=1;
                if (strcmp(working_add, 'working'))
                    totalSumR(img)=totalSumR(img)+max(pos(inds,2))-min(pos(inds,2))+1;
                end;
            end;
        end;
    end;
	if (strcmp(working_add, 'working'))
        totalR(img)=sum(matchedArea(:));
    end;
    
    % overlay
    matchedSym = repmat(-single(matchedSym),[1 1 3]);
    matchedSym = 1 * (matchedSym-min(matchedSym(:)))/(max(matchedSym(:))-min(matchedSym(:)));
    if showPartBoundingBox
		for y = 1:size(matchedSym,2)
		    for x = 1:size(matchedSym,1)
		        if sum(abs(matchedBoundingBox(x,y,:))) > 0
		            matchedSym(x,y,:) = matchedBoundingBox(x,y,:);
		        end
		    end
		end
    end
	
	tmp = single(repmat(J{end},[1 1 3])) / 255;
	alpha = .6;
    tmp = tmp * alpha + matchedSym * (1-alpha);
    imwrite( tmp, sprintf('%s/%s_overlayed_image%d.png',destFolder,tr_or_tt,img) );
    imwrite( matchedSym, sprintf('%s/%s_sketch_image%d.png',destFolder,tr_or_tt,img) );
    imwrite( colorSketch, sprintf('%s/%s_colorsketch_image%d.png',destFolder,tr_or_tt,img) );
end
fprintf(1,'\n');

if (strcmp(working_add, 'working'))
    ratio=sum(totalR)/sum(totalSumR);
end;
