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 is a forward finite difference operator for a ith coordinate, and k is an iteration counter; calculating a first auxiliary variable s(k) from
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