function [tmpImage, errV] = reconstructImage(resF, sortRes, gaborBank, powerMapF, indFiltersToUse, im, imageSize, sizeRes, param, plotta)

%init matrix indicating if for each position a wavelet
%with orient and scale already used
tmpRes = initTmpRes(gaborBank, powerMapF);

%init error vector
errV = zeros(param.numWavelets, 1);

%init reconstructed image
tmpImage = zeros(imageSize);

%create figure
if plotta
    fsfigure
end %if plotta

%init counters
countW = 1; %wavelet counter
ind = 1; %sorted response counter
%loop on wavelet responses
while (countW <= param.numWavelets && ind <= sizeRes)
    
    %get information of corresponding wavelet
    currScale = sortRes(ind, 2);
    currTheta = sortRes(ind, 3);
    currX = sortRes(ind, 4);
    currY = sortRes(ind, 5);
    currO = sortRes(ind, 6);
    
    %check if current wavelet at current position already
    %used
    if (tmpRes(currO).value(currY, currX) == 1)
        %increment counter
        ind = ind + 1;
        continue
    end %if (tmpRes(currY, currX, currO)
    
    %check if current wavelet among the best
    if ~ismember(currO,  indFiltersToUse)
        %increment counter and skip
        ind = ind + 1;
        continue
    end %if ~ismember(currO,  ind_o_counter_sort)
    
    %get filters with current parameters
    currGaborEven = gaborBank(currO).even;
    currGaborOdd = gaborBank(currO).odd;
    
    %get weights for gabor filters
    weightEven = resF(currO).even;
    weightOdd = resF(currO).odd;
    
    %compute filtersize
    filterSize = computeFilterSizeFromScale(currScale);
    %compute step
    step = computeStepFromScale(currScale);
    
    %get filter position
    %check x,y correct
    [posX, posY, posX_filter, posY_filter] = findFilterPosition(currX, currY, step, imageSize(1), filterSize);
    
    %shift filter to center filter at x,y
    %currGaborEven = adjustGaborPos(currGaborEven, currX, currY, imageSize);
    %currGaborOdd = adjustGaborPos(currGaborOdd, currX, currY, imageSize);
    
    %add to image
    %tmpImage = tmpImage + (weightEven(currY, currX) * currGaborEven);
    %tmpImage = tmpImage + (weightOdd(currY, currX) * currGaborOdd);
    
    %add to image
    tmpImage(posY, posX) = tmpImage(posY, posX) + (weightEven(currY, currX) * currGaborEven(posY_filter, posX_filter));
    tmpImage(posY, posX) = tmpImage(posY, posX) + (weightOdd(currY, currX) * currGaborOdd(posY_filter, posX_filter));
    
    %Mean Squared Error between original image and
    %reconstructed
    imOrigNorm = normalizzaImg(im);
    err = immse(imOrigNorm, normalizzaImg(tmpImage));
    errV(countW) = err;
        
    if plotta && mod(countW, 1000) == 0
        
        %get image limits
        climG = getImageLimits(currGaborEven);
        climI = getImageLimits(tmpImage);
        
        
        subplot(2,2,1);
        imshow(im,[])
        title('Original Image')
        
        subplot(2,2,2);
        imshow(currGaborEven, climG);
        title('Current wavelet')
        
        subplot(2,2,3);
        imshow(tmpImage, climI)
        title({['Reconstructed image with ' num2str(numel(indFiltersToUse)) ' filters, after ' num2str(countW) ' wavelets'], ['MSE: ' num2str(err)]});
        
        subplot(2,2,4);
        plot(errV(1:10:countW), 'LineWidth', 2);
        xlabel('Number of wavelets added')
        ylabel('MSE')
        title('MSE of reconstruction')
        
        pause(0.01)
        
    end %if mod(countW, 100) == 0
    
    %increment counters
    countW = countW + 1;
    ind = ind + 1;
    
    %threshold
    %if (err < param.minMSE) && (countW > param.minCountW)
    %break
    %end %if err < param.minMSE
    
end %while (countW <= param.numWavelets && ind <= sizeRes)