    

    %% solution of least squares min |I_syn - B_selected * C|^2
    %% C = inv(B'B)B'I;  I_reconstruct = B*C = B* inv(B'B)B' I_syn, we compute and store R=B inv(B'B)B' as reconstruction matrix
    tic;
    B=[];
    for i=1:template.numSelected    
       Bi = filterTiling(filters{template.selectedFilter(i)}, sx, sy, template.selectedRow(i), template.selectedCol(i));
       B = [B,Bi(:)];
    end
    R=B*(pinv(B'*B)*(B'));

    disp(['preparing the reconstruction matrix for least square: ' num2str(toc) ' seconds']);
    %%%% save each chain of synthesized images respectively
    EachChainSynFolder=[Synfolder2 'EachChainSyn_numWavelet_'  num2str(template.numSelected) '_numSketch_' num2str(nStroke) '/'];
    if exist(EachChainSynFolder)
        rmdir(EachChainSynFolder,'s'); 
        mkdir(EachChainSynFolder);
    else
        mkdir(EachChainSynFolder);
    end
    %%%% reconstruct each chain of synthesized images respectively
    EachChainReconstructionFolder=[Synfolder2 'EachChainReconSyn_numWavelet_' num2str(template.numSelected) '_numSketch_' num2str(nStroke) '/'];
    if exist(EachChainReconstructionFolder)
        rmdir(EachChainReconstructionFolder,'s'); 
        mkdir(EachChainReconstructionFolder);
    else
        mkdir(EachChainReconstructionFolder);
    end


    I_multi_chain_projected = zeros(nTileRow*sx, nTileCol*sy);
    I_multi_chain = zeros(nTileRow*sx, nTileCol*sy);

    index=0;
    for j=1:nTileCol
       for i=1:nTileRow
         index=index+1;
         crop=I_syn(1+(i-1)*sx:sx+(i-1)*sx, 1+(j-1)*sy:sy+(j-1)*sy);
         gLow = min(crop(:));
         gHigh = max(crop(:));
         img = (crop-gLow)/(gHigh-gLow);       
         imwrite(img,[EachChainSynFolder 'original_chain-' num2str(index,'%04d') '.png']);
         I_multi_chain(1+(i-1)*sx: sx+(i-1)*sx, 1+(j-1)*sy: sy+(j-1)*sy) = img; 
         %% project the synthesized images onto the domain specified by the selected basis
         img_projected = reshape(R*crop(:),sx,sy);  % R = B*(inv(B'*B)*(B'));
         gLow = min(img_projected(:));
         gHigh = max(img_projected(:));
         img_projected = (img_projected-gLow)/(gHigh-gLow);       
         imwrite(img_projected,[EachChainReconstructionFolder 'reconstruct_chain-' num2str(index,'%04d') '.png']);
         I_multi_chain_projected(1+(i-1)*sx: sx+(i-1)*sx, 1+(j-1)*sy: sy+(j-1)*sy) = img_projected; 
       end
    end
    imwrite(I_multi_chain,[EachChainSynFolder 'raw.png']);
    imwrite(I_multi_chain_projected,[EachChainReconstructionFolder 'reconstruct.png']);

            
            