function [ X ] = B_HTP( Phi,y,s,dx1,dx2,N0 )
%B Summary of this function goes here
%   Detailed explanation goes here

    X = zeros(dx1,dx2);
    k = 1;
%     Phi=conj(Phi);
    while k <= N0
        
        %Residual
        Xtemp = X + reshape(Phi'*(y-Phi*X(:)),dx1,dx2);
        
        %find indices of the s largest rows in norm of Xtemp
        nX = sqrt(sum(Xtemp.^2,2));
%         for l = 1:dx1
%             
%             nX(l) = norm(Xtemp(l,:),2);
%             
%         end
        [~,I] = sort(nX,'descend');
        J = sort(I(1:s));
        if k > 1
            if I == I_old
                k = N0;
            end
%             disp(norm(I-I_old));
        end
%         disp(norm(y-Phi*X(:)))
        %solve LeastSquares on support J
        X = LeastSquares(Phi,y,J,dx1,dx2);
%         disp(norm(y-Phi*X(:)))
        I_old = I;
        k = k+1;
    end

end

