A.

K-sparse and compressible signals

Let’s consider a real-valued, finite-length, one dimensional, discrete-time signal x, which can be viewed as a N×1 column vector in ℜ^N with elements x[n], with n = 1,2,…,N. Any signal in ℜ^N can be represented in terms of a basis of N×1 vectors For simplicity, let’s assume that the basis is orthonormal. Using the N×N basis matrix Ψ = {ψ₁,ψ₂,…,ψ_N} with the vectors {ψ_i} as columns, a signal x can be expressed as:

where s is the N×1 column vector of weighting coefficients s and x are equivalent representations of the signal with x in time or space domain and s in Ψ domain.

The signal x is K-sparse if it is a linear combination of only K basis vectors, which means that only K of the s_i coefficients in (1) are nonzero, while the remaining (N–K) coefficients are zero. In addition, the signal x is compressible if the representation in (1) has just a few large coefficients and many small coefficients, setting the basis of transform coding. Therefore, we can say that a signal x is sparse in the ψ domain if the coefficient sequence is supported on a small set, and compressible if the sequence is concentrated near a small set. In face of the typical data acquisition paradigm, huge amounts of data are collected only to be in large part discarded at the compression stage to facilitate storage and transmission. Imagine, for example, a digital camera that captures images with millions of sensors (pixels) but eventually encodes the image in just a few hundred kilobytes. Clearly, this is a tremendously wasteful process and suffers from three principal drawbacks. First, the initial number of samples N may be large, even if the desired K is small. Second, the set of all N transform coefficients {s_i} must be computed even though all but K of them will be discarded. Third, there is an overhead that is introduced by the encoding of the large coefficients locations.

B.

Recovering K-sparse signals

Following the work presented in [4], Candès and Tao developed a refined version of the Uniform Uncertainty Principle (UUP) [5], which has proved to be essential to the study of the general robustness of CS. This key notion was then named Restricted Isometry Property (RIP) and can be defined as follows:

For each integer K = 1,2,…, define the isometry constant δ_K of a measurement matrix A as the smallest number such that

holds for all K-sparse vectors x. Therefore, we can say that a matrix A obeys the RIP of order K if δ_K differs enough from one. When this condition is verified, A approximately preserves the Euclidean length of K-sparse signals, which in turn implies that K-sparse vectors cannot be in the null space of A. An alternative description of this property is to say that all subsets of K columns taken from A are in fact nearly orthogonal (they cannot be exactly orthogonal since we have more columns than rows).

Let’s imagine we want to acquire K-sparse signals making use of matrix A. Suppose that δ_2K is not close to one. This indicates that all pair-wise distances between K-sparse signals must be well preserved in the measurement space, which means that

is true for all K-sparse vectors x₁, x₂ [1, 6].

C.

Incoherence

Let’s now consider M < N linear measurements of x and a collection of test functions such that y[m] = ⟨x,φ_m⟩. By stacking the measurements y[m] into the M×1 vector y and the test functions as rows into an M × N sensing matrix Φ we can write

A condition related with RIP is incoherence, which requires that the rows of Φ (the measurement or sensing matrix) cannot represent the columns of Ψ in a sparse way (and vice-versa).

Incoherence extends the duality between time and frequency and expresses the idea that an object having a sparse representation in Ψ must be spread out in the domain in which it was acquired. This incoherence property is verified for many pairs of bases, including, for instance, delta spikes and sine waves of the Fourier basis, or the Fourier basis and wavelets. The coherence between the sensing basis Φ and the representation basis Ψ can be given by the following equation:

which, in simple words, is measuring the largest correlation between any two elements of Φ and Ψ. CS is essentially interested in low coherence pairs. For instance, for the previously referred delta spikes and sine waves (time-frequency) pair, μ(Φ,Ψ)=1, therefore, indicating maximal incoherence [1, 7].

A particular aspect of interest is that random matrices are largely incoherent with any fixed basis Ψ. This empowers the use of known fast transforms such as a Walsh, Hadamard, or Noiselet transform [8].

Furthermore, what is most remarkable about this concept is that it allows capturing information contained in a sparse signal in a very efficient way without trying to previously understand that signal.

D.

How compressive sensing works

Compressive sensing addresses the inefficiencies presented by the sample-then-compress framework by directly acquiring a compressed signal representation, avoiding the intermediate stage of acquiring N samples [5]. CS bypasses the sampling process and directly acquires a condensed representation y consisting of M linear measurements. Furthermore, The measurement process is non-adaptive in that Φ does not depend in any way on the signal x. The transformation from x to y is a dimensionality reduction and so loses information in general. In particular, since M < N, for a given y, there is an infinite number of x’ such that Φx’= y. The overwhelming capacity of CS is that Φ can be designed such that sparse/compressible x can be recovered exactly/approximately from measurements of y. To recover the signal x from the random measurements y, the traditional favorite method of least squares has been shown to fail with high probability. Instead, it has been demonstrated that using the l₁ optimization [4]

it is possible to exactly reconstruct K-sparse vectors and closely approximate compressible vectors stably with high probability using just M ≥O(K·log(N/K)) random measurements [2, 3].

A.

Principle of operation

In this work, we chose to implement an active illumination imaging LIDAR system, which compressed the spatial information through the incorporation of the random measurement codes to the light used to illuminate the scene, as described in the setup depicted in Fig. 2.

Observing the diagram presented in Fig. 2, the following principle of operation can be enunciated. The pulsed laser source (Alphalas Picosecond Pulse Diode Laser @ 670 nm) emits a beam, with the desired repetition rate, that will be spatially expanded by a beam expander (Edmund Optics Beam Expander 10x), in order to illuminate the DMD array in which the compressive codes are applied. Then, through a set of projection optics, the images representing the compressive codes applied in the DMD (LogicPD LightCommander) are projected into the field of view being acquired.

Fig. 2.

Diagram with the setup of the compressive-sensing based Imaging LIDAR System.

The light reflected from the scene being illuminated is then collected by a finder scope (Vixen 7x50 Finder Scope with a field-of-view of 7.2°) and a fiber collimator (Thorlabs Fiber Collimator) ultimately launches the light into an optical fiber that is connected to the single photon counter (PicoQuant tau-SPAD-50). The single photon counter is then connected to a photon counter board (Brilliant Instruments BI220) that is connected to the computer that is controlling the entire system. The software running on the control computer then gathers the timing histograms that will be used to calculate distances and reconstruct the images.

Fig. 3 illustrates the final assembly of the designed system in which it is observable the finder scope assembled on an L-shaped bracket with two rotational stages, so that its pitch and yaw can be adjusted. On the smaller end of the finder scope it is represented the collimator that is connected to the single photon counting module through an optical fiber. On a parallel direction to that of the finder scope, one can see the optical module comprising the projection optics and the DMD array, with the beam expander and the laser support using mechanical parts that have been designed and manufactured.

It should be said that the projection optics and the DMD array have been mounted based on a DMD projection kit produced by LogicPD, called LightCommander, in order to allow one to evaluate the feasibility of the proposed approach at this initial stage.

Fig. 3.

(left) View in lateral perspective of a 3D drawing of the prototype of the compressive sensing based Imaging LIDAR System. (right) Photo of the implemented compressive sensing based Imaging LIDAR System.

B.

Experimental procedure and results

For the acquisition of experimental data and testing of the system, our efforts relied on the development of the software so that the Imaging LIDAR System could illuminate the scene with the compressive random codes and acquire a timing histogram for each of the codes. These compressive random codes were based on either simple or permuted Hadamard codes. After the acquisition phase, we were then in conditions of reconstructing the ranged images using the timing histograms.

For the generation of the timing histograms the BI 220 photon counter board has been used. The BI220 board outputs the information necessary to produce a timing histogram every time that a predefined number of photons is counted. Therefore, in order to evidence the effect of using the different compressive random codes for each of the histograms a normalization process had to be implemented. For that, the number of counts of each timing histogram has been divided by the total acquisition time of that particular histogram. If this normalization had not been implemented, the effect of the spatial modulation of the light (either imposed by the compressive codes, or by the scene itself) would not be evident in the number of counts of the timing histograms, as the histograms would be considered complete only when the predefined number of photons was counted, and the time it took to complete each histogram would not impact the computation. The BI220 board did not have the possibility to perform measurements during a predefined amount of time, so an alternative strategy had to be implemented. For further clarification, imagine that the single-pixel camera was using a photodiode in forward biased mode instead of an avalanche photodiode operating as a single-photon detector and it was measuring the current passing at the photodiode with a fixed sampling interval. In that case, we could state that the current would vary with different compressive codes as the amount of light hitting the photodiode area would also vary with the compressive codes. This variation imposed by the compressive codes would then allow us to reconstruct the images being acquired [9].

In order to explain how the images are reconstructed based on the histograms, consider the example shown in Table 1.

Table 1.

Example to explain the reconstruction of ranged images from the timing histograms.

In Table 1 for each compressive code the corresponding histogram is registered. Each histogram contains in each bin a normalized value representing the ratio between the photons counted for that bin and the total acquisition time for that histogram. Then, to reconstruct an image representing the scene at a particular distance, the software used the compressive measurements associated with that particular distance as the input to the reconstruction algorithm. In the example provided, if one wanted to reconstruct an image of a scene at 5.0 m distance from the imaging LIDAR system, one would have to use the measurements of the column highlighted in red.

Some tests to characterize the timing information when the scene was being illuminated and the conditions of acquisition were kept unchanged were also performed. For that, we acquired 100 timing histograms (with 1000 bins) for a predefined amount of detected photons when the scene (white wall) was being totally illuminated (all-white code) and when the scene was being illuminated with a compressive random code. The situations in which the scene was illuminated using an all-white code or image, or when it was being illuminated with a compressive random code can be interpreted as a range finder mode operation or as a compressive imaging mode of operation, respectively. For these experiments, for each of the 100 timing histograms the number of counts of the histogram peak, its occurrence in time and the time required to perform all the measurements have been registered. Table 2 contains the results of these experiments and presents the average of the number of photon counts for each peak normalized by the time spent to acquire all the measurements for the corresponding histogram, the average distance and the average time consumed to acquire a histogram, plus the respective standard deviations (std).

From the analysis of Table 2, one can conclude that the operation mode (range finder [all-white code/image] or compressive imager [compressive random code]) did not affect significantly the distance values obtained from the timing information of the peak in each histogram.

From the information provided in Table 2 it is also obvious that the rate of detected photons along time was higher in the case in which the scene was being illuminated with an all-white image, as the average normalized photon counts have been always higher comparatively to the cases in which the scene was being illuminated with a compressive random code. This makes perfect sense since the amount of light hitting the scene when it is being illuminated with a compressive random code is lower than in the case when an all-white image is used for illumination.

In terms of measurement time, it is clear that it increases with the increase of the number of measurements associated with each histogram. This information is particularly important once it allows one to evaluate the amount of time required to acquire the desired number of histograms that will be used to reconstruct the images. In order to provide an estimate of the compromise that needs to be made, the ratio between the standard deviation and the corresponding average normalized photon counts was calculated. This allowed us to have an insight of the level of noise present in the corresponding measurements. Based on previous work [9], we can state it is possible to reconstruct images of similar quality in a compressive imaging manner when the amplitude of the noise is below 20% of the maximum amplitude of the signal without noise (SNR = 14.54 dB). Therefore, it is possible to infer from Table 2 that a histogram built using a number of measurements above 5000 shall address the stated desired conditions.

Table 2.

Results for the characterization of the measurements under different illumination conditions.

During operation, in order to ensure ideal conditions, one shall guarantee that the variation imposed by a compressive random code to the normalized photon count number is higher than the standard deviation associated with the current parameterization of the system, specifically the number of photon detections required for each histogram. If this is guaranteed we are in perfect conditions to reconstruct images from the resulting histograms as the effect of noise is strongly reduced.

Next, some results obtained with the system will be presented. To obtain these results, two modes of operation have been defined. One in which the system was operating as a range finder and another in which the system was operating as a 3D imager. During the tests performed the laser pulse repetition rate was set to 10 MHz.

For the case of range finding, the system was completely illuminating a flat target placed at a particular distance. Then, 400 blocks with a predefined amount of measurements were acquired and for each block a histogram of distance values was built (based on the time of flight of the photons that have been detected). Afterwards, to determine the distance to the flat target it was assumed that it corresponded to the distance value associated with the center of the Gaussian trace that best fit each histogram. In the end, 400 distance values have been collected. During these experiments two block sizes have been used – i.e., 5000 and 25000 measurements per block – and the number of bins used to build the histograms has also been varied – i.e., 1000 and 5000 bins per histogram –.

In order to evaluate how spread the 400 distance values were, they were distributed in histograms and were fit using a normal or Gaussian distribution. This time only 20 bins were used to distribute the distance values. The center of the fitting trace then corresponded to the average of the most frequent distance value.

The experiments previously described have been repeated with the target placed at three different distances (differing 1 cm from the immediately previous position) in order to assess the accuracy of the system. Table 3 presents the results obtained when the different parameterizations were used. From the values presented in Table 3, one can state that the system was able to detect the 1 cm displacement of the target. It can also be said that the resolution of the system is much better than 1 cm as the three histograms (in all the cases) are far from overlapping, i.e., their standard deviation is relatively reduced.

Table 3.

Average distance values obtained from the fittings for the different parameterizations.

To acquire 3D LIDAR images, the system illuminated the scene with the random binary compressive codes and for each code a distance histogram has been recorded. Each of these histograms had 1000 bins and was computed using 50000 measurements (detected photons). In the end, the images corresponding to the desired range of distance values were reconstructed and a datacube was created. For the imaging mode of operation, a scene consisting of a white flat wall and two boxes has been composed (see Fig. 4-left). The box that was lying closer to the wall had a thickness of 10 cm, while the other box had a thickness of 5 cm.

Fig. 4-center and Fig.4-right display a 27 x 29 pixels image and a 3D mesh, respectively, that has been reconstructed and that represent the distance information of the composed scene with color coding. Due to some reconstruction errors near the border of the image, caused by the non-uniformity of the illumination and by the acceptance cone of the collection optics, the 32 x 32 pixels image has been cropped to provide a better visual representation of its content. The lack of uniformity of the illumination derives from the projection optics and from the Gaussian nature of the laser beam profile.

Fig. 4.

(left) Photo of the scene, used to test the single-pixel Imaging LIDAR System, consisting of a flat wall and two boxes. The dashed red square has been added to indicate the location of the volume being imaged. Image reconstructed with 27 x 29 pixels (center) and respective 3D mesh (right) representing in color the distances (in meters) of the discrete targets that compose the scene displayed in the photo on the left.

Another test scene that was used for the 3D imaging mode of operation consisted of a futsal ball with 20 cm diameter (see Fig. 5). The results obtained for the test scene of Fig. 5 are displayed in Fig. 6 and were reconstructed using the measurements contained in 40 distance bins of 1024 histograms (with 1000 bins) that were acquired.

Fig. 5.

Photo of the scene, used to test the single-pixel Imaging LIDAR System, consisting of a futsal ball.

Fig. 6.

32 x 32 pixels image (right) and respective 3D mesh (left) representing in color the distances (in meters) of the ball surface in the scene displayed in Fig. 7.

Table 4 quantifies the main operational parameters of the implemented prototype. Considering the finder scope angular aperture of 3.6°, at 10 m the system’s field-of-view is defined by a circle with a diameter of 1.25 m.

Table 4.

Main operational parameters of the single-pixel imaging LIDAR system.