A method of compressed sensing imaging includes acquiring a sparse digital image b, said image comprising a plurality of intensities corresponding to an I-dimensional grid of points, initializing points (x(k), y(k)), wherein x(k) is an element of a first expanded image x defined by b=RΦ−1 x, wherein R is a Fourier transform matrix, Φ is a wavelet transform matrix, y(k) is a point in ∂ ( ∑ i = 1 l ⁢ ( ∇ i ⁢ Φ - 1 ⁢ x ( k ) ) 2 ) 1 / 2 , ∇i is a forward finite difference operator for a ith coordinate, and k is an iteration counter; calculating a first auxiliary variable s(k) from x ( k ) - τ 1 ( αΦ ⁢ ∑ n ⁢ L n * ⁢ y n ( k ) + Φ ⁢ ⁢ R * ( R ⁢ ⁢ Φ - 1 ⁢ x ( k ) - b ) ) , wherein τ1,α are predetermined positive scalar constants, the sum is over all points n in x, and L* is an adjoint of operator L=(∇1, . . . , ∇I); calculating a second auxiliary variable tn(k) from yn(k)+τ2LnΦ−1x(k), wherein τ2 is a predetermined positive scalar constant; updating x(k+1) from sign(s(k))max{0,|s(k)|−τ1β}, wherein β is a predetermined positive scalar constant; and updating yn(k+1) from min ⁢ { 1 τ 2 ,  t n ( k )  2 } ⁢ t n ( k )  t n ( k )  2 .