function xest = cspocs_K(Phi, aPhi, K, y, niter)
% Syntax : 
%  xest = cspocs_K(Phi, aPhi, K, y, niter)
%
% Description : Variation on the Candes and Romberg POCS (alternate
% Projection Onto Convex Sets) [1] for CS recovery, aka recovering of x from
% y = Phi*x when size(Phi,1) < size(Phi,2), using a **sparsity** criterion
% as the one of Lu Gan [2] (Stage 1 only of the proposed method without
% Wiener filtering). The soft-tresholding used for projecting signal onto
% the l1 ball of radius ||x||_1 in [1], is replaced in [2] by hard
% thresholding keeping the K largest elements. The signal is there projected
% onto the non-convex (!) set of sparse signals.
% 
% In :
% * Phi, aPhi : Measurement matrix and its reconstruction
% * K : the assumed sparsity of x 
% * niter : maximun number of alternate projections.
%
% Out :
% * xest : the recovered (or estimated) signal x
%
% Author of this mfile : L. Jacques, LTS2/EPFL, 2008.
%    
% References : 
% [1] Lu Gan, Block compressed sensing of natural
% images. (Proc. Int. Conf. on Digital Signal Processing (DSP), Cardiff, UK,
% July 2007) http://www.dsp.ece.rice.edu/cs/block_CS.pdf
% [2] E. J. Candès and J. Romberg. Practical signal recovery from random
% projections. Wavelet Applications in Signal and Image Processing XI,
% Proc. SPIE
% Conf. 5914. http://www.acm.caltech.edu/~emmanuel/papers/PracticalRecovery.pdf
%
% Example :
% >> N=128; K=20;
% >> x=[rand(1,K) zeros(1,N-K)]'; x=x(randperm(128));    
% >> m=floor(3.5*K);Phi=randn(m,N)/sqrt(m);aPhi=Phi';          
% >> y=Phi*x;
% >> figure; plot(x);
% >> nx=cspocs_K(Phi,aPhi,K,y,10000);
% Tolerance reached. Stop at n=90
% Final score: ||y - Phi*xest|| = 8.321048e-11
% >> hold on; plot(nx,'ro');
% 
% Remark: Even if it is less obvious why this algo works (since the K-sparse
% signals set is not convex, excepted perhaps locally), (1) this algo stops
% after much less iterations than cspocs_l1, (2) this algo is also much less
% sensible to K than cspocs_l1 is sensitive to the a priori l1 norm of
% x. Recovery/NoRecovery transition points seems located around 3.2K
%
% This script/program is released under the Commons Creative Licence
% with Attribution Non-commercial Share Alike (by-nc-sa)
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% Short Disclaimer: this script is for educational purpose only.
% Longer Disclaimer see  http://igorcarron.googlepages.com/disclaimer
            
    [m,N] = size(Phi);
    
    pPhi = pinv(Phi);
    pPhiPhi = pPhi*Phi;
    
    %% First projector (on the hyperplane Phi*x=y)    
    P = @(u) u - (pPhiPhi)*u + pPhi*y;
    
    %% Second projector on the l1 ball of radius l1val
    H = @(u) u .* (abs(u) >= abs(K2th(u,K)));
    
    %% Initializations
    xest = pPhi*y;
    Tol = 1e-10;
    score = 0;
    
    for n = 1:niter;
        xest = P(xest);
        xest = H(xest);
        
        oldscore = score;
        score = norm(Phi*xest - y) / norm(xest);
        
        if (score < Tol)
            fprintf('Tolerance reached. Stop at n=%i\n',n);
            break
        end
    end

    fprintf('Final score: ||y - Phi*xest|| = %e\n', score);
    
function out = K2th(x,K)
% Pick the magnitude the Kth largest element of x.
% Used by Hard Thresolding.
[ans, ind] = sort(abs(x), 'descend');
out = x(ind(K));