Depth Sensing Using Optical Pulses and Fixed Coded Aperature

[0001] This application is a related to U.S. patent application Ser. No. 13/405,504, “Depth Sensing Using Active Coherent Signals,” filed by Boufounos on Feb. 27, 2012, and incorporated herein by reference. That application provides a method for sensing depths in a scene using active coherent sensing of reflected radar pulses using a single transmitter and multiple receivers arranged in a plane, using back projection.

[0002] This invention relates generally to reconstructing a scene using optical pulses, and more particularly to reconstructing the scene as a three dimensional (3D) depth map.

[0003] Compressive Sensing (CS) uses fewer linear measurements than implied by the dimensionality of an acquired signal. To reduce the acquisition rate, CS reconstruction methods exploit the structure of signals. To capture the structure, the most commonly used signal model is sparsity.

[0004] As used herein, sparse is a well known term of art in signal processing, and not a relative indefinite term. A “sparse” approximation estimates a sparse vector, i.e., a vector with most coefficients equal or approximately equal to zero, satisfying a linear system of equations given high-dimensional measured data and an acquisition matrix.

[0005] A CS-based acquisition system can be model as

r=A(s)  (1)

where A(·) is a linear function, s is in some appropriate signal space, and r is in a measurement space. The measurement space has a much lower dimension than the signal space. A number of possible properties of A(·), such as low coherence, the restricted isometry property, depending on the model, guarantee that the reconstruction is possible using an appropriate method.

[0006] CS is successfully applied in many imaging systems, in which, typically, the signal s to be acquired is a 2D image in N_x×N_y.

[0007] Using compressive sensing approaches, it has been shown that images can be acquired with measurements as few as 10% of the number of pixels N_xN_y. These gains are not as relevant in conventional visible-light imaging, where charge-coupled devices and metal-oxide-semiconductor sensor technology make measurements extremely inexpensive, which has had a significant impact in other modalities, such as medical imaging, low-light imaging, hyperspectral imaging, and depth sensing. However, the gains are significant for imaging systems that depend on complex costly sensors, such as “time-of-flight” (TOF) light or radar sensors.

[0008] Model-Based Compressive Sensing

[0009] A model-based CS framework provides a general approach to developing a large number of signal models and characterizing their suitability for CS acquisition. Models under this framework are created by imposing restrictions on the signal support. A fundamental operation is the projection of a general signal to the set of signals that satisfy support constraints on the model. As long as such a projection can be computed, common greedy CS reconstruction methods, such as Compressive SAmpling Matching Pursuit (CoSaMP) and Iterative Hard

[0010] Thresholding (IHT), can be modified to reconstruct signals in the model. Furthermore, it has been shown that a pair of approximate projections with different approximation properties is sufficient to guarantee accurate reconstruction, instead of an exact projection, see Hegde et al., “Approximation-tolerant model-based compressive sensing,” Proc. ACM Symposium on Discrete Algorithms (SODA). SIAM, pp. 1544-1561, January 2014.

[0011] One signal model, approximation tolerant model-based compressive sensing, is motivated by signals measured using 2D seismic imaging. A signal Sε^N×T according to that model is a matrix with N rows and T columns. Each row of the signal only has S non-zero entries, which are spatially close to the S non-zero entries of the row above or below. This is enforced by restricting an earth-mover distance (EMD) between the support of subsequent rows of the signal. The projection under the approximation tolerant model is performed by solving a sequence of network flow problems.

[0012] Optical Sensing

[0013] Optical sensing measures distances to a scene by illuminating the scene with optical signal and analyzing the reflected light. One optical sensing technique, e.g., light radar (Lidar), can be used with applications such as geoscience, remote sensing, airborne laser swath mapping (ALSM), laser altimetry, contour mapping, and vehicle navigation.

[0014] Conventional high-resolution, high frame-rate optical based systems typically use an expensive array of precision TOF sensors and illuminate the scene with singular pulses transmitted by a stationary laser. Alternatively, at the expense of reduced frame-rate, the laser scans the scene. A smaller sensor array, which can also scan the scene, acquires the reflection. The resulting system significantly lowers the cost, but requires the use of mechanical components, which can be prone to failure and cumbersome in mobile applications.

[0015] Compressive sensing can exploit significant gains in computational power due to the reduced sensing cost, and allow elaborate signal models and reconstruction methods, which, in turn, enable reduced sensor complexity. For example, some compressive depth sensing systems use a single sensor combined with a spatial light modulator and multiple pulses illuminating the scene.

[0016] A spatial light modulator is used to implement a variable coded aperture that changes with every pulse. However, that restricts those system to static scenes, see Howland et al., “Compressive Sensing Lidar for 3D imaging,” CLEO: Science and Innovations. Optical Society of America, 2011, CMG3, Kirmani et al., “Exploiting sparsity in time-of-flight range acquisition using a single time-resolved sensor,” Opt. Express, vol. 19, no. 22, pp. 21485-21507, October 2011, and Kirmani et al., “CoDAC: a compressive depth acquisition camera framework,” IEEE Int. Conf. Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2012, pp. 5425-5428, and U.S. 20130088726. In addition, spatial light modulators (SLM) are relative expensive, e.g., $1000 or more depending on the desired resolution, see Holoeye® or Texas Instruments DLP® products.

[0017] With S=1, a restricted version of the above approximation tolerant model is a good model for 2D imaging (one spatial and one depth dimension) with light pulses as the model minimizes a 1-D total variation of the depth map. However, to extend the model to 3-D volumes, i.e., signals sεN_x××N_y×T (two spatial and one depth dimension) is unfortunately not obvious.

[0018] Continuous improvement in optical sensing components, as well as recent advances in signal acquisition theory, provide an opportunity to reduce the cost and enhance the capabilities of depth sensing systems. The embodiments of the invention provide a novel depth sensing system and method that uses a fixed coded aperture (mask) to significantly reduce the number of required sensors, compared to conventional systems. As an advantage, the mask can be cheaply printed on any transparent substrate.

[0019] The embodiments also provide a modeling and reconstruction framework based on model-based compressed sensing, which characterizes a large variety of depth sensing systems. Experimental data indicate that it is possible to reduce the number of sensors by more than 85% with a negligible effect on the sensing quality.

[0020] The main contributions of the invention include:

• • [0021] a real-time compressive hardware architecture for light-based time-of-flight (TOF) depth sensors;
  • [0022] a very general system model that characterizes a large number of depth sensing system architectures; and
  • [0023] a model-based compressive sensing (CS) reconstruction method that exploits a low variation (TV) of depth scenes to improve the reconstruction.

[0024] In contrast with the prior art, the embodiments use a small number of sensors and a fixed coded aperture. The aperture is a very inexpensive component compared to a spatial light modulator (SLM). Because the code does not change, one embodiment of the invention only requires a single pulse transmission per frame—that is a depth map in a sequence of depth maps—which may or may not carry a code. Thus, frame-rates equivalent to much more expensive single-pulse systems can be achieved at a significantly lower cost.

[0030] The embodiments of our invention provide a method and system for reconstructing a scene using optical sensing and a fixed coded aperature.

[0031] Optical Sensing

[0032] System Architecture

[0033] As shown in FIG. 1, a stationary light source 110, e.g., a laser, transmits an optical pulse for each frame to be acquired of a scene 100. A lens turns the pulse into a wide beam pulse 111. By using a wide beam pulse, the scene does not need to be scanned. The pulse can be, for example, Gaussian shaped 113 in time. Alternatively, a coded sequence of pulses 114 can be used.

[0034] A reflected pulse 112, from each scene point 102 that is a reflector, passes through a lens 120, which disperses the reflected pulse to a set of sensors 140 arranged in a plane with the light source. This way the total distance is simply 2d. The dispersed pulse passes through a fixed code aperture 130 in the form of a mask in an optical path between the light source and the scene. The mask can be cheaply printed on a transparent material, e.g., an acetate sheet. The mask spatially modulates the dispersed pulse.

[0035] In one embodiment, the arrangement of the sensors can be random. In another embodiment, the sensors can be arranged as a uniform array. One advantage of these embodiments, particularly when the cost of optical time-of-flight (TOF) sensors is significant, is that that the number of sensors required is reduced by about 85% when compared with conventional systems, while maining the sensing quality.

[0036] The output 301 of the set of sensors is processed by a sparse reconstruction method 300 in a processer 150 to obtain a reconstructed scene in a form of a depth map 101 using a scene model. The processor can also control 151 the timing of the pulses transmitted by the light source 110.

[0037] FIG. 2 shows an alternative structure of the system. In this case, the mask 130 is in an optical path between the light source and the scene 100. In this embodiment, the scene is scanned 115 by a nonstationary light source 110, and the set has only one sensor 140. As the light source scans the scene, the source transmits to the scene, through the coded aperture, one pulse at each position in the scanning path 115. Generally, the fixed coded aperture in an optical path between the light source and the set of sensors.

[0038] In addition, or alternatively, the sensors can be nonstationary while scanning the scene. An alternative embodiment could combine the two architectures and use both a scanning laser with or without a fixed coded aperture in front of it, combined with multiple sensors behind another fixed coded aperture on the return path of the reflection, with the same or different code.

[0039] Reconstruction Method

[0040] FIG. 3 shows the reconstruction method according to embodiments of the invention. For one embodiment (CoSaMP), the pseudo code for steps 1-10 are shown in FIG. 4. For an alternative embodiment (IHT), the pseudo code for steps 1-9 are shown in FIG. 5. The steps iterate until the estimated signal converges to the depth map 101. The variables used in the pseudo code are described herein.

[0041] Initially, the signal estimate is set 310 (step 1) to zero. The input data are in the form of a matrix of acquired signals by the set of sensors 140. Residuals that describe how the input data disagree with the signal estimate are determined 320 (steps 4-6) from signal the estimate and input data (signals).

[0042] Next, 330, the signal estimate is updated using the current signal estimate, the signal model, and the residuals. Depending on the reconstruction procedure, the model might or might not be used in this step.

[0043] Next, in step 340, the signal estimate is updated by enforcing the scene model 210. The process is repeated from 320 until convergence 350. In one embodiment, the scene model imposes a constraint that there is at most a single reflector along any direction.

[0044] Embodiments of the method rely on a K-term truncation step, which can take various forms in different embodiments. For example, the truncation step can be a simple K-sparse approximation, as known in the art, or it can be a selection of the largest in magnitude coefficient, related to a single reflector, along every spatial direction, as described in U.S. patent application Ser. No. 13/405,504. Additionally or alternatively a graph-cuts formulation can be used as described below.

[0045] Time-of-Flight

[0046] The reconstruction for the depth map is based on time-of-flight (TOF) principles because the sensing of reflected pulse is delayed from the time the pulse is transmitted due to the distance 2d. The delay corresponds to the distance to, or depth of each point 102 in the scene.

[0047] Conventional systems must carefully separate and sense the reflected pulses from each direction to determine the TOF and assign the delay to the correct place in the depth map 101.

[0048] In contrast, we use compressive sensing (CS) principles and intentionally mix the the reflected pulse in all directions to reduce the sampling burden. Therefore, we rely on computational methods to separate the mixture and assign the correct depth value to each point in the scene.

[0049] The fixed coded aperture can be implemented using a mask that blocks light according to, e.g., a binary pattern. In one embodiment the pattern is random. A coded aperture where the pattern has continuous values in a range [0,1] is also possible. The combination of the lens and the coded aperture mixes reflected pulses received from each direction and projects the mixed pulses to the set of sensors. Each sensor receives a different mixture, determined by the aperture code. Thus, the signal acquired by each sensor is a mixture of the reflected pulses, as received by the sensor. These mixtures are used to reconstruct the scene.

[0050] In the alternative embodiment, the single sensor receives a different mixtures as the laser with the coded aperture scans the scene. At every point of the scan path of the laser, the scene is illuminated differently, and therefore each reflection acquired by the sensor is a different mixture of the scene.

[0051] Signal Model

[0052] To model the system, we start with a single reflecting scene point at distance d from the sensing plane, and we assume there is no coded aperture. For convenience, we assume the light source 110 is coplanar with the sensors, although this assumption can be trivially relaxed. Thus, the distance the pulse travels from the light source to the reflecting scene point and back to the sensors is 2d, corresponding to a pulse delay, i.e., TOF, τ=2d/c, where c is the speed of light in the medium. Using this correspondence, time is equivalent to distance from the sensor plane, i.e., depth. Herein, we use distance and depth interchangeably.

[0053] We consider a 3-D scene comprising two spatial directions, transverse to the sensor plane, and one depth direction, which we refer to as time or delay, and is perpendicular to the sensor plane. Thus, the scene to be acquire can be represented by a function s_x,y,t, where s represents the reflectivity of the scene at point (x,y,t). If there is no reflection from a point, then the corresponding reflectivity is zero.

[0054] We assume Lambertian surfaces with no transparency, which implies that, for any (x,y) pair, there is only one depth t that has non-zero reflectivity. In other words, if there is a reflection from (x,y,t), then there is no reflector in-between that location and the sensing plane, i.e., in front of the reflector. Furthermore, the light does not reach any reflectors behind that point for the same (x,y) pair, and, therefore, their reflectivity is also zero. This constraint, in the context of coherent sensing systems, as a depth constraint on s. A valid scene signal s should satisfy this depth constraint.

[0055] Given a depth map d_x,y, representing the depth of the scene at coordinates (x,y), and a reflectivity (albedo) map α_x,y for the same coordinates, the scene is can be represented by

s_x,y,t=α_x,yδ_t-2ddx,y^/c,  (2)

where δ_t is the Dirac impulse function. For any scene satisfying the depth constraint, it is trivial to extract depth and albedo maps. In other words, a 3D scene corresponds to a 2D depth map and a 2D albedo map. Similarly, a 2D (one space and one depth dimension) scene corresponds to a 1D depth map and a 1D albedo map.

[0056] Acquisition Model

[0057] We use p_t to denote the transmitted pulse 111 or pulse sequence to be reflected by the scene 100. Assuming that a pinhole aperture is present, then the received signal for a location (x,y) is

{circumflex over (r)}_x,y,t=a_x,yp_t-2dx,y/c=p_t{circle around (*)}_tS_x,y,t)  (3)

where {circle around (*)}_t denotes a linear convolution along the time direction. The coded aperture 130, with spatial code m_x,y, introduces a mixing of the reflected pulses, which can be shown to be the spatial convolution of {circumflex over (r)}_x,y,t with the mask. Thus, the received light flux at the sensor plane is

{tilde over (r)}_x,y,t=m_x,y{circle around (*)}_x,y{circumflex over (r)}_x,y,t=m_x,y{circle around (*)}_x,yp_t{circle around (*)}_tS_x,y,t.   (4)

[0058] This signal is sensed by M sensors, indexed by m=1, . . . , M. Each sensor is positioned at location (x_m,y_m). Each sensor samples the light flux at the sensor plane and samples r_m,t={tilde over (r)}_xm,ym,t in time.

[0059] We discretize the scene to sε^Nx,Ny,Nt, where N_x, N_y is a number of spatial pixels to be acquired, or desired resolution of the system, and N_t is the number of time samples. The discretization is such that each reflector in the scene can be assumed to be approximately flat and parallel to the sensor plane over the area of a pixel, such that the depth constraint is preserved in the discrete representation. This discretization also implies a discretized 2D depth map dε^Nx,Ny, which also takes discrete values d_x,yε{1, . . . , N_t}.

[0060] Furthermore, the sampling time is at a rate higher than the Nyquist rate of the transmitted pulses. Similarly, we discretize the pulse, the coded aperture mask, and the received signals, such that the convolution with the pulse and the mask pattern can be expressed in discrete-time.

[0061] The linear combination of the discretized received signal rε^M,Nt′ at sensor m can be decomposed as a a sequence of simpler linear transformations

r=S(M(P(s)))=A(s)  (5)

where P, M and S denote, respectively, the linear transformations due to the pulse, the mask of the coded aperture, and the sampling of the optical field by the sensors. Their composition A is the forward linear transformation mapping of the scene s to the received signal r.

[0062] While the particular model described the embodiment in FIG. 1, i.e., when the laser is stationary, a very similar linear model can be easily derived that describes other embodiments, such as the embodiment in FIG. 2, in which the laser uses a coded aperture to scan the scene. The reconstruction method only requires that such a linear model exists, which can take a scene s and express the received signal r as a linear transformation, r=A(s), of that scene.

[0063] Model Implementation

[0064] An efficient computation using the linear model in equation (5) is paramount in modern iterative reconstruction methods. To this end, the operators P and M can be efficiently implemented in a discrete time using fast Fourier transform (FFT) based convolution methods.

[0065] Furthermore, S is trivial computationally because it simply selects the appropriate signals. The adjoint A*, necessary for most reconstruction methods, can be implemented by composing the adjoint of each operator in the reverse order, i.e., using A*(r)=P*(M*(S*(r))). This implementation can also exploit the separability of space and time operators in equation (5), and rearrange the linear combinations as S(M(P(s)))=P(S(M(s))) to simplify the complexity. Thus, the pulse convolution P is applied to an M×N_t-dimensional object, rather than a N_x×N_y×N_t dimensional object, while the complexity of S and M do not change.

[0066] An efficient implementation can use a depth map d and albedo a representation to significantly reduce memory requirements. The forward operator can be determined efficiently from such a representation, although the adjoint requires additional temporary memory to store the entire s.

[0067] Depth Reconstruction

[0068] Depth Scenes, Variation and Network Flows

[0069] To reconstruct a subsampled signal and reduce noise, CS exploits the structure of the acquired signal. In particular, depth maps tend to have a low variation (TV) norm. Signals with low TV norms are generally flat, with very few discontinuities and very few areas with small gradients.

[0070] The (l₁) TV norm of a discrete map dε^Nx×Ny is defined as

∥d∥_TV=∥∇_xd∥₁+∥∇_yd∥₁,  (6)

where ∇′_xd is a discrete gradient along direction x and ∥·∥₁ is an element-wise l₁ norm of a matrix, i.e.,

∥∇_xd∥₁=Σ_ny=1^NΣ_nxx=1^N-1|d_nx,ny−d_nx,ny|,  (7)

with the y direction similarly defined.

[0071] Given a scene satisfying the TV norm depth constraint, an additional constraint on the TV norm of the depth map is equivalent to a constraint on the support of the non-zeros in s. Specifically, two spatially adjacent non-zero coefficients of s, i.e., with respect to the n_x and n_y coordinates, should also have similar depths, i.e., n_t coordinates, except for very few discontinuities. However, this one-to-one correspondence requires that s satisfies the depth constraint. A general dense s does not. Thus, a projection to both the depth and the TV constraints is required for model-based CS methods, such as CoSaMP and IHT.

[0072] For a one-dimensional depth map problem, i.e., for a two-dimensional scene sε^Nx,Nt and dε^Nx, the projections described by Hegde et al., “Approximation-tolerant model-based compressive sensing,” Proc. ACM Symposium on Discrete Algorithms (SODA). SIAM, pp. 1544-1561, January 2014 provide a solution. There, 2D signals, represented as a matrix, are considered. Each row of the matrix has only S non-zeros, and the support of those non-zeros from row to row changes very little, according to a pre-determined constraint on the EMD between supports.

[0073] It can be shown that with S=1, the EMD constraint applied to s becomes a TV constraint on its support, i.e., the depth map. The projection onto the constraint can be determined using a sequence of dynamic programs solving a network flow problem. Unfortunately, that approach does not generalize to 2D depth maps, i.e., 3D scenes. For those, we use a different formulation.

[0074] One formulation that can be used for 3D scenes is described in U.S. patent application Ser. No. 13/405,504, in which the largest in magnitude coefficient along the depth direction is selected for each spatial position x, y. The K largest coefficients and their location are selected among those. One embodiment of the invention can use this formulation. However, this formulation does not enforce a low TV norm constraint in the 2D depth map (i.e., in the support of the non-zeros in the 3D scene).

[0075] A different formulation that can enforce this constraint is described next. This is a preferred embodiment of the invention as it enables further reduction of the number of sensors or the laser scanning time, compared to the embodiments using the other formulations.

[0076] 2D Variation and Graph Cuts

[0077] To generalize the projection to 3D objects, i.e., 2D depth maps, we use a graph cut formulation. For an undirected, weighted graph =(,), we consider the general graph cuts problem. That is, given a set of observations , the task is to assign each vertex ε a label Iε such that the joint labeling of all vertices l minimizes an energy function between labels and observations E(l,).

[0078] In our depth sensing method, we map each vertex to represent a spatial location =(n_x,n_y) of the scene, and each label to represents a discrete depth value l=d_nx,ny. Hence, the cardinality of sets and are N_x,N_y and N_t, respectively.

[0079] We also map the set of observations to the scene s. We express the energy function as a sum of unary and pairwise terms

$\begin{array}{cc}E\left(l,s\right)=-\underset{\underset{\mathrm{Unary}}{}}{\sum _{v\in }s_{v,l_v}^2}+\underset{\underset{\mathrm{Pairwise}}{}}{\sum _{u\in {}_v}^{v\in }\lambda l_v-l_u},=-\sum _{\left(n_x,n_y\right)\in }s_{n_x,n_y,d_{n_x,n_y}}^2+& \left(8\right)\\ \sum _{\left(n_x,n_y,\right)\in {}_{n_x,n_y}}^{\left(n_x,n_y\right)\in }\lambda d_{n_x,n_{y^{\square }}}-d_{n_{x^{\prime }},n_{y^{\prime }}},& \left(9\right)\end{array}$

where _nx,ny={(n_x+1,n_y),(n_x−1,n_y),(n_x,n_y+1),(n_x,n_y−1)} is the neighborhood of _nx,ny, i.e., contains all vertices that are directly adjacent to =(n_x,n_y) in the graph.

[0080] The unary term is a fidelity term, which uses the label, i.e., the depth value, to select the appropriate data point from the scene s and impose an l₂ data penalty. The pairwise term imposes a smoothness constraint between the label of v and the label of vertices in the neighborhood set . Thus, the pairwise term from equation (8) is the l₁ norm of the gradient of depth values, i.e., the TV norm of d. Analogous to Rudin-Osher-Fatemi total variation, the parameter λ represents a tradeoff between data fidelity and smoothness.

[0081] Solvers for minimizing equation (8) that are well known include alpha-expansion and alpha-beta swap, as well as Boolean approaches. In one embodiment, we use an alpha-expansion technique.

[0082] The truncation step (6, 8 in FIG. 4, and 8 in FIG. 5) of our model-based method incorporates a K-term truncation of s by first optimizing equation (8) to obtain a candidate depth map d, which corresponds to a candidate support set (n_x,n_y,d_nx,ny). From this candidate support set, the largest K components of s are selected to be used by an appropriate step of model-based CoSaMP or IHT.

[0083] In contrast, conventional truncation just selects the K largest components of the data s not constrained by the graph cuts solution. Compared to conventional sparsity, our model-based thresholding produces a scene s that satisfies the depth and TV constraints.

[0084] Scene Reconstruction

[0085] Given a projection, such as the one described above, it can be incorporated to a reconstruction algorithm to recover the actual depth map. Two alternative embodiments are demonstrated in FIGS. 4 and 5.

[0086] In FIG. 4, an iterative algorithm to perform reconstruction is described, based on model-based CoSaMP (compressive sampling matching pursuit). An iteration count l and a signal estimate {tilde over (s)}^l are maintained and updated throughout each iteration. Step 1 initializes them to 0. Step 3 increases the iteration count. Step 4 computes a residual u^l, which captures the disagreement of the signal estimate {tilde over (s)}^l-1 with the data r. A proxy g^l for the scene is computed in Step 5, which is also the gradient of a cost function at the current scene estimate {tilde over (s)}^l-1. Step 6 uses a model-based truncation function such as the one described above to determine the support of the proxy, truncated to K coefficients according to the model, denoted supp(trunc(g^l,K)). Note that alternative embodiments could use a different truncation length such as 2K or K/2. This support is united with supp({tilde over (s)}^l-1) the support of the current signal estimate to form a support candidate ^l. A new temporary signal estimate b^l is formed in Step 7 by inverting the system restricted over the support candidate. The temporary signal estimate is truncated in step 8 to update the signal estimate {tilde over (s)}^l. Steps 3-8 are iterated until the converged of the algorithm or for a finite number of iterations. Convergence can be measured, for example, by the magnitude of the residual computed in step 4, or by the amount of change in the signal estimate {tilde over (s)}^l from iteration to iteration. The final signal estimate after convergence, {tilde over (s)}^l, is returned at Step 10.

[0087] In FIG. 5, an alternative iterative algorithm to perform reconstruction is described, based on model-based IHT (iterative hard thresholding). An iteration count l and a signal estimate {tilde over (s)}^l are maintained and updated throughout each iteration. Step 1 initializes them to 0. Step 3 increases the iteration count. Step 4 computes a residual u^l, which captures the disagreement of the signal estimate {tilde over (s)}^l-1 with the data r. A proxy g^l for the scene is computed in Step 5, which is also the gradient of a cost function at the current scene estimate {tilde over (s)}^l-1. An alternative temporary scene estimate {circle around (s)}^l is formed in Step 6, by following the gradient with a step size τ. The temporary signal estimate is truncated in step 7 to update the signal estimate {tilde over (s)}^l. Steps 3-8 are iterated until the converged of the algorithm or for a finite number of iterations. Convergence can be measured, for example, by the magnitude of the residual computed in step 4, or by the amount of change in the signal estimate {tilde over (s)}^l from iteration to iteration. The final signal estimate after convergence, {tilde over (s)}^l, is returned at Step 9.

[0088] Although the invention has been described by way of examples of preferred embodiments, it is to be understood that various other adaptations and modifications may be made within the spirit and scope of the invention. Therefore, it is the object of the appended claims to cover all such variations and modifications as come within the true spirit and scope of the invention.