% This is my Robust Super Resolution Method:
function [reconstPhoto,Ws] = robustSuperResolution(ims,Fs,photoWidth,photoHeight,X0,sigma,tol,maxIter);
    photoSize = photoWidth*photoHeight;
    
    X = X0(:);
    sigmaSqr = sigma^2;
    
    stillConverges =  1;
    lastResidual   = -1;
    lastX = X;
    
    i=1;
    while(i<=maxIter & stillConverges)
        A = sparse(photoSize,photoSize);
        b = sparse(photoSize,1);
        
        for j = 1:size(ims,2)
            xj = ims{j};
            Fj = Fs{j};
            
            imageSize = size(Fj,1);
            
            Wj = sigmaSqr ./ (sigmaSqr + (Fj*X-xj(:)).^2).^2;
            Ws{j} = reshape(Wj,size(xj,1),size(xj,2));
            
            Qj = (Fj)' * sparse(1:imageSize,1:imageSize,Wj);
            A = A + Qj * Fj;
            
            b = b + Qj * xj(:);
        end
        
        X = cgs(A,b,1.e-3);
        
        % Compute the residual:
        residual = 0;
        for j = 1:size(ims,2)
            xj = ims{j};
            Fj = Fs{j};
            residual = residual + norm(Fj*X-xj(:),2);
        end
        
        % Check if it still is converging:
        stillConverges = (lastResidual-residual > tol | lastResidual == -1);
        if(~stillConverges & residual >= lastResidual)
            X  = lastX;
            Ws = lastWs;
        else
            lastX  = X;
            lastWs = Ws;
            lastResidual = residual;
        end
        
        i = i+1;
    end
    
    reconstPhoto = reshape(X,photoHeight,photoWidth);
    