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Φ−1x, wherein R is a Fourier transform matrix, Φ is a wavelet transform matrix, y(k) is a point in∂(∑i=1l(∇iΦ-1x(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) fromx(k)-τ1(αΦ∑nLn*yn(k)+ΦR*(RΦ-1x(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) frommin{1τ2,tn(k)2}tn(k)tn(k)2.