1 Introduction

Compressive Sensing (CS) is a new and emerging field in the signal processing world; its main idea is direct acquisition of a signal in a compressed form [1–7]. Compressive Imaging (CI) [1, 4–7] is an optical implementation of CS where new optical designs replace conventional imaging systems with the aim of saving on sensor cost, time of acquisition or both.

Conventional digital compressing techniques use the redundancy feature of natural images. The compressed image is represented by small set of coefficients, compared to the number of pixels used for the acquisition. The theory of compressed sensing states that an image can be reconstructed with an overwhelming probability from fewer measurements than the conventional number of pixels provided it is compressible in some known basis (e.g. Fourier or wavelet) and an appropriate sensing mechanism is applied.

Two properties, which lie at the basis of the theory of CS, are the sparsity of the signal and the incoherence of the sensing mechanism. The sparsity property refers to the idea that the signal can be represented by a small set of significant coefficients. It implies that such a signal is compressible in the basis in which it is sparse. The sparsifying basis is denoted by$\Psi$. The incoherence property of the sensing part claims that in contrary to the sparse signal in a given basis, the measurements must be highly dense in that basis. Mathematical definition of the incoherence and its bound can be found in [2, 3]. In theory the best universal sensing operator$\Phi$ is one that performs uniformly random projections of the data [2, 3]. Full reconstruction is possible owing to optimization methods based on $l_1$norm minimization. The theoretical number of measurements is found to be directly related to both of these properties [2].

Optical implementation of random $\Phi$ involves several severe challenges [4, 5]. In [6] several compressive imaging (CI) architectures were proposed that compromise randomness for the benefit of optical implementation feasibility, yet exhibiting good performance. Those architectures are based on direct or indirect angular sub-sampling of the Fourier plane. One of the implementations proposed in [6] uses optical Radon projections of the object plane onto a line array sensor (vector sensor). Radon projections were used to demonstrate the CS concept in the seminal work of Candes et al. [8], while heuristics of their incoherence can be found in [9]. The principle of CI based on Radon projections [6] is briefly summarized in subsection 2.1. In this work we present an extension of the concept introduced in [6] by considering angular sampling that permits progressive acquisition of the image information. Whereas in [6] a uniform angular sampling is proposed, here we propose an alternative sampling scheme in which it is possible to progressively add information. This scheme is based on the irrational golden angle sampling step [10, 11].

Progressive sampling is beneficial when the sparsity of the signal is unknown in advance and therefore the number of samples required to be captured is not known prior to capturing. If the number of samples required is not known a-priori one would adopt a scheme that permits successive image improvement with respect to the numbers of samples. The sampling process can be stopped by the user once desired image quality has been achieved. Another benefit of a progressive sampling scheme is immunity to sudden stops during the sampling process; in such a case the partial samples suffice to provide a visual informative image. This point is demonstrated in sub-section 2.2. Another important advantage of the progressive sampling scheme proposed here is that it fits well with ordered-subset reconstruction algorithms, which as we show in Sec. 3 makes CI and reconstruction of large images feasible.

The structure of the paper is as follows. In Sec. 2 we review briefly the compressive sensing technique based on optical Radon projections presented in [6], and then explain the new angular sampling scheme for progressive compressive imaging. Section 3 presents the proposed ordered-sets (OS) estimation algorithm, while Sec. 4 shows the results of the computer simulations and real experiments done with the progressive sampling technique. Section 5 concludes and summarizes this work.

2 Progressive compressive sensing with fixed step angular sampling

2.1 Compressive sensing with optical Radon projections

The CI technique presented in [6] relies on capturing Radon projections of the image plane. Unlike Radon transform commonly used in medical tomography, the Radon transform used here is in the object plane, that is, perpendicular to the optical axis. If $f(x,y)$is the object plane, the Radon transform is given by

 

Fig. 1 The principle of optical Radon Projections. The detector S collects the ray sum of the object by means of the cylindrical lens.

Download Full Size | PDF

The point-spread function of the cylindrical lens is a line perpendicular to x’ axis, as demonstrated in Fig. 1 for points A, B and C. Therefore, all points of the image which lie on an axis perpendicular to y’ (e.g., A and C in Fig. 1) are summed on the sensor pixel according to the line integral formula (1). The height of the perpendicular line, y’, is the distance, s, in the line integral in Eq. (1). Figure 1 shows that the red point, C, and the green point, A, which lie on the same line are spread and overlapped in the image plane. This way, the central pixel of sensor S will capture the sum of all points, that lie on the x’ axis. The blue point, B, lies some distance away from the red and green points and will therefore contribute to the pixel proportional to the y’ coordinate of the point. The next projection is obtained by rotating the (x’-y’) plane by angle$\theta$ .

The original image can be reconstructed from Radon projections which have been captured with the angular sampling technique described above. By employing nonlinear sparsity-promoting algorithms, signals sampled with a compression ratio of approximately 1:10 were reconstructed in [6, 7].

In general, the reconstruction process involves minimization of the cost function $J(a)$in the form of [5]:

Technically, the forward operator $A$may be implemented explicitly by an appropriate transition matrix or implicitly, for instance as a function handle using Matlab’s Radon transform. Using A in the form of a matrix has the advantage of accessibility to the precise adjoint A^T_, which is necessary for iterative reconstruction algorithms. However, its disadvantage is its huge size when large images are considered. For example, if the matrix operator is used with megapixel size images, the number of coefficients in a typical A matrix is of the order O(10¹²) . Such matrices require memory of the order of terabytes and even in compact format available due the matrix sparsity, it requires gigabytes of memory. Yet, it is clear that in order to handle such a matrix large computational power is required. Fortunately, in Sec. 3 we show how this problem of dimensionality can be remedied significantly owing to an OS reconstruction algorithm available with the progressive sampling scheme proposed here.

2.2 Angular sampling for progressive compressive imaging

The natural and most straightforward angular sampling scheme is one that uses constant angular steps and that samples the object uniformly. Such a scheme requires a-priori knowledge of the number of projections Q that will be sufficient for the reconstruction.

If one samples the object uniformly, with fixed angular step, then these steps can be written in radians as$d\theta =2\pi \alpha$, with$\alpha$being a normalized step. Since the circle is being sampled repeatedly, starting at the origin$\theta =0$, the position of the q^th sample can be written as

In the uniform sampling with fixed step scenario described in [6], the angular step is set at$\alpha =1/Q$, where Q is the number of projections. In fact only $\pi$radians need to be scanned because of the symmetry of the Radon transform. The sparsity of the object dictates the number of projections Q, and since each object has different sparse representation, the number of projections needs to be adjusted by means of changing the angular step$\alpha$. easing the reconstruction effort. The reconstruction algorithm running time has been found to be approximately linear with the number of projections in advance. Therefore, if the angular sampling step a is not chosen sufficiently small to provide satisfactory reconstruction for a particular object, one may want to continue adding projections after the $\pi$ radians that have been scanned. However, if the angular steps are given by (3) with$\alpha =1/Q$, or more generally with $\alpha$ being any rational number, the new projections overlap the previous ones and therefore do not add any new structural information. Consequently the angular step needs to be an irrational number multiple of 2π, so that no overlapping between projections will occur. In addition, in order to fulfill low coherence of the samples, as required by CS theory, the angular samples should not be concentrated in a specific angular range. This implies that not only the projections must be all different, but also the angular separation should be as uniform as possible as shown in Fig. 2 .

 

Fig. 2 (a) 20 angular samples over 2π taken at uniform (+) and golden angle ($\odot$), (b) a scenario in which the sampling process is unintentionally stopped after 8 samples.

Download Full Size | PDF

Figure 2(a) (color online) demonstrates the two sampling schemes. The red cross marks relate to the uniform and fixed angular step with$\theta =2\pi /20$ radians; i.e., according to (3) with$\alpha =1/20$. The circle is being sampled uniformly by definition. The blue circles relate to the second sampling scheme, where the angular step is an irrational number known also as the “golden angle”${\theta }_{ga}\approx 137.45\text{[deg]}$, an angle subtended by the small arc (b) in Fig. 3 that obeys

 

Fig. 3 Arcs, a, and, b, are in a golden ratio relation. Golden angle ${\theta }_{ga}$is subtended by arc b.

Download Full Size | PDF

The ratio  between the two arcs is the golden ratio [12] which can be found as the solution of the quadratic equation:${\phi }^2-\phi -1=0$.

The golden ratio is well known to humanity since ancient times and was thought to be the divine sign, since it appears in nature in snail’s shells, the leaves of the sunflower and more [12]. For this reason, ancient artists, architects and philosophers used it extensively in their works. Here we choose to use the sampling step to be the golden angle ${\theta }_{ga}$because as shown in [11] it has the smallest discrepancy value. The discrepancy criterion [10] $D\left(\Theta /2\pi \right)$is used to measure how close the sequence $\Theta =\left[{\theta }_0,{\theta }_1,\mathrm{...},{\theta }_{Q-1}\right]$is to a uniform sampling and is found to be the closest for the golden angle sampling.

Figure 2 demonstrates that the discrepancy between the two sampling schemes is negligible. It can be seen that sampling with golden angle steps does not yield local aggregation of points and in fact it is almost as uniform as with angular steps of $\theta =2\pi /20$ .

Figure 2(b) shows angular samples obtained in a scenario in which the scanning process is suddenly stopped in the middle, so only 8 angular samples are captured. It can be seen that despite the small number of samples the angular sampling in steps of ${\theta }_{ga}$ spans the circle fairly uniformly, whereas sampling with a step of $\left(Q+\Delta Q\right)$$\alpha =1/20\text{[deg]}$naturally fails to cover the circle.

2.3 Advantages of angular sampling with golden angle step for progressive compressive imaging

Performing angular sampling with golden angle steps makes it possible to add information progressively. Since the samples do not overlap, the information addition is very simple compared to the regular uniform fixed step sampling. Progressive information addition may be useful when gradual quality improvement is acceptable. This may save time of the overall process, since the detail of interest may be identified at an early stage of acquisition, and the process may then be stopped completely or continued until the result is satisfactory. Consider for example the case in which, because of wrong assumption about object’s sparsity, Q angular samples are taken, leading to an unsatisfactory reconstructed image. In order to refine the image quality, with the sampling scheme using angular sampling steps of ${\theta }_{ga}$one just needs to continue sampling till the reconstruction is adequate. In contrast, with the uniform sampling scheme one needs to guess the number of additional samples required,$\Delta Q$ and perform a new scan taking additional new $\left(Q+\Delta Q\right)$samples with an angular step of$\theta \text{'}=2\pi /\left(Q+\Delta Q\right)$.

Progressive sampling is also useful also for transmission purposes; the receiver does not have to wait until all the samples are received and only then reconstruct the image, but can generate gradually improving images from partial data arriving continuously.

The second advantage of the progressive sampling scheme is robustness to sudden system scanning failure as demonstrated in Fig. 2(b). The case considered here is a sudden sampling process cessation, when only 40% of the required projections have been taken. It can be seen that while the red cross signs, denoting sampling step of $\theta =2\pi /20$ have completed sampling only half of the circle, the golden angle steps denoted by the blue signs have covered the whole circle. The meaning of this is that in spite of the malfunction of the system, the user still gets a reasonable representation, albeit of a degraded quality. We note that for illustrative purposes, in this example we did not account for the symmetry of the Radon transform and we demonstrate the point on a full circle. Obviously, the uniformity of golden samples holds also on part of the circle. Another important advantage of the non-overlapping property of the samples is the applicability of OS reconstruction algorithms, as explained in the following section.

3 Reconstruction with ordered subsets

3.1 The Ordered Subsets reconstruction concept

Reconstruction algorithms working on ordered subsets are probably best known in the realm of medical tomography (e.g., [13]). Ordered subsets or block-iterative methods have been applied together with a number of iterative algorithms, with Ordered Subsets Estimation Maximization (OS-EM) being the most widely used. The essential difference between the OS-EM algorithm and conventional algorithms (e.g., Maximum Likelihood Estimation Maximization, ML-EM) is the use of only a subset of projections for updating the estimate rather than using the entire set of measured projections. For OS-EM, the use of only part of the data during the updating process is termed a sub-iteration and the term single iteration usually refers to the use of all the data at once. With OS-EM the reconstruction proceeds by utilizing subsets of the projections, chosen in a specific order that attempts to maximize the new information being added in sub-iterations. The iteration process proceeds by using different projections in each subsequent subset until all projections are used. OS-EM was found to be much faster than the whole set of ML-EM with similar reconstruction quality results. It is important to note that the subset partition should be such that subsets are balanced, meaning that projections are as uniform distributed as possible within the subset [14].

3.2 Ordered sets obtained by golden angle sampling

Here we follow the main idea used within the statistical algorithm OS-EM. However, instead of performing Estimation Maximization (EM) on the ordered sets [13] we do an l₁ type minimization (2), as prescribed by the CS theory. The captured projections are divided into subsets and iteration proceeds by using different projections in each subsequent subset until all projections are used [14].

The uniformity property of samples taken with the golden angle step scheme becomes useful in the proposed OS algorithm. This property also holds on subsets of sequential samples. When the whole measurement set, acquired with the golden angle step scheme is divided into temporally captured subsets, the angular samples are inherently uniformly distributed in each subset, so the necessary condition of subset balance is automatically fulfilled. To put this more formally, let us denote by $G_1={\left\{g\left(s,m{\theta }_i\right)\right\}}_{{\theta }_i\in {\Theta }_1}$the first projection subset taken at angles ${\Theta }_1={\left\{m{\theta }_{ga}\right\}}_{m=1}^{M_t}$with ${\theta }_{ga}$being the golden angle and M_t being the subset order (the number of projections in a subset). The second subset is $G_2={\left\{g\left(s,m{\theta }_i\right)\right\}}_{{\theta }_i\in {\Theta }_2}$taken at the next M_t samples${\Theta }_1={\left\{m{\theta }_{ga}\right\}}_{m.=M_t+1}^{2M_t}$. Since the angular step ${\theta }_{ga}$is fixed for all the subsets, the second subset may be represented as a version of the first subset shifted by $M{\theta }_{ga}$phase or${\left\{g\left(s,(m{\theta }_{ga}+M_t{\theta }_{ga})\right)\right\}}_{m=1}^{M_t}={\left\{g\left(s,m{\theta }_{ga}\right)\right\}}_{m=M_t+1}^{M_{t+T}}$, thus the uniformity of the subsets is exactly the same. In other words,$D\left({\Theta }_1\right)=D\left({\Theta }_2\right)$, where the discrepancy function$D\left(\Theta /2\pi \right)$ represents the degree of the uniformity of the sampling in the subset relative to the uniform sampling [10].

The proposed OS algorithm works as follows. The algorithm stages are divided into inner sub-iterations and outer iterations. Each sub-iteration, q, acts on the sub-set G_q. One full iteration is defined when all the sub-sets have been optimized once. The first outer iteration sets the initial guess as zero and passes it to the first inner iteration, which feeds the set of measurements G₁ into the solver, which performs a CS optimization. The output of the first inner iteration is passed to the second inner iteration as an updated guess, while the set of the measurements is changed to G₂. This process continues until all the sub-sets have been optimized. After all the sets have passed the inner iteration, the minimization result is passed back to the outer iteration and is set as an initial guess for the following outer iteration. We found empirically that three to four outer iterations are sufficient, and there is no additional improvement with a larger number of outer iterations.

Besides the advantage of using OS in reducing the processing time, there is a significant advantage in terms of memory requirements. Recall that, as explained in section 2.1, the dimensions of the forward matrix A representing the set of Radon projections may be extremely large. However, when the measurements are divided into sets according to the golden angle sampling step, like in OS, the problem of$A\text{'s}$ dimensionality mentioned in section 2.1 can be avoided, and the computations may be performed on a regular PC. Recall that the golden angle sampled subsets are built in a progressive way, while containing partial information about the whole image and not a fraction of it. Thus a much smaller $A(G)$which can be created and stored separately from the other parts of is $A$facilitated.

4 Results

In order to demonstrate the progressive compressive imaging technique we performed simulations and optical experiments, where in both cases the image was reconstructed from a different number of projections acquired with the proposed sampling scheme. The reconstruction process is done using the TwIST solver [15] on the simulated or optically acquired measurements.

4.1 Simulated experiment

The simulation relates to the acquisition process. The purpose of the simulation is to verify the progressive improvement of the reconstruction with the increasing number of samples when a golden angle step, ${\theta }_{ga}$, is used. For comparison, we also simulated in parallel an experiment in which Radon projections were taken with a uniform sampling scheme; i.e., with steps of${\theta }_u=2\pi /Q$ radians. This latter auxiliary experiment refers to the case when the number of the samples that should be taken for good reconstruction quality is known in advance.

For the numerical test we used the synthetic image in inset Fig. 4(c) which is highly sparse in the “Haar” wavelet basis with sparsity K = 1.4% non-zero elements. In the simulation of the acquisition part, the Radon projections of the synthetic image are generated at an increasing number of angles by means of the forward projection matrix. The reconstruction part computes the approximations from the acquired projections, which are then compared to the original image in terms of PSNR:

 

Fig. 4 Reconstruction quality as function of number of Radon projections. (a) Reconstruction in the wavelet domain. (b) Reconstruction in the space domain, (c) Synthetic image used in simulation (1280 x 1280 pixels)

Download Full Size | PDF

Figure 4 shows the reconstruction PSNR, measured in [dB], versus the number of projections for the golden angle and uniform sampling scheme.

Figure 4(a) shows the results for signal reconstruction simulation in the space domain, where the vector $a$in Eq. (2) is the signal itself and$\Psi =I$. The gradual improvement in the reconstruction quality with increasing number of samples is clearly seen in Fig. 4(b). The golden angle scheme is best suited for this scenario since its advantage is in reducing acquisition time by adding new samples, in contrast to the uniform angle scheme, where re-sampling is required.

Figure 4(b) shows the reconstruction performed in the wavelet domain, i.e., with Ψ in (2) being Haar wavelet transform. It can be seen that the improvement of the reconstruction quality with increasing number of samples is faster in Fig. 4(b) than in Fig. 4(a) because the image is much sparser in the wavelet domain. The improvement rate is similar for both sampling schemes (with${\theta }_u\text{and}{\theta }_{ga}$), which implies that a similar number of samples is needed in order to achieve a given reconstruction quality.

4.2 Optical experiment

We built an optical setup following the incoherent CI approach outlined in [6] with the ability to automatically capture projections and with controllable and precise angular sampling. The scheme of the system is shown in Fig. 5 . Unlike the concept presented in Fig. 1, where both the lens, L, and the sensor, S, are rotated synchronously about the axis, the only rotating element in the system of Fig. 5 is a prism, which rotates an aerial image at the input of the cylindrical lens. The schematic drawing in Fig. 5 shows an object, which sends a bundle of rays into the beam splitter. Half of the rays are directed into the rotating prism, which reflects them back through the beam splitter and through the cylindrical lens onto the linear detector located in the image plane of the lens. The rotating right angle prism is mounted on the computer-controlled rotating stage with z being the axis of rotation. Each position of the prism enables a specific angle of the Radon projection.

 

Fig. 5 Schematic diagram of the automatic CI optical system. Radon projections are recorded with the linear detector aligned with the Y axis and located in the image plane of the system.

Download Full Size | PDF

As in the simulation part, we took several sets of projections from 1.5% to 5.5% of the nominal reconstruction number of samples. For an image of size $n\times n$we defined the nominal number of projections to be n, so that the entire Radon transform will have approximately n² point samples. For the quantitative evaluation of the reconstruction we compare the reconstructed images to the virtually perfect reference image (“golden standard”) taken to be a reconstruction made from 20% of the nominal number of the projections. This choice is reasonable since according to our experiments and simulations good reconstruction quality has been achieved with 20% of the projections. For instance, in Fig. 5 it can be seen that with 20% of the nominal projections (256 projections) the asymptotically best image is reached.

Figure 6 shows the results of reconstructions made from 20, 40, 50 and 70 projections, which are 1.56% (a), 3.13% (b), 4.3% (c) and 5.5% (d) of the nominal number of samples. It can be clearly seen that the visual quality improves successively with the increasing number of projections. It should be emphasized that the reconstruction was performed using the progressive scheme; for instance, image (b) is reconstructed from the same 20 projections as image (a) with addition of the next 20 new projections. Without prior knowledge of the required number of samples for good reconstruction, one would choose to acquire 256 projections (20%) in order to be on the safe side and would get the image (e). However, in the progressive sampling case one could stop sampling already after 70 projections (5.5%) to get image (d), which is virtually of the same quality. This way one would acquire a good version of the image nearly 3.5 times faster, save storage space and reduce the reconstruction effort. The reconstruction algorithm running time is approximately linear with the number of projections, which explains why in this example it is approximately 3.5 faster. This result confirms the simulation outcome and shows that the golden angle sampling scheme is, indeed, a feasible solution when no prior information about the object’s complexity is known and results in time saving, when compared to the conventional progressive sampling where re-sampling is necessary.

 

Fig. 6 (a-e) Progressive reconstructions from. 1.56%, 3.13%, 4.3%, 5.5% and 20% of the nominal samples needed for perfect reconstruction.

Download Full Size | PDF

As explained in section 2.3, the golden angle sampling scheme permits OS reconstructions with the explicit Radon forward operator even for large images. The advantage of using the precise forward matrix, A, instead of an implicit forward operator is demonstrated in Fig. 7 . Figure 7(a) shows the real life object of size 1280 x 1280 pixels. Figures 7(b)–7(d) show the reconstructions obtained from 70 projections, yielding a compression ratio of approximately 20:1. Figures 7(b) and 7(d) show the reconstructions obtained by applying the TwIST algorithm to the entire set of projections. The reconstruction is carried out with an implicit forward Radon operator implemented as a function handle referring to Matlab’s Radon transform function. The image in Fig. 7(b) is obtained from uniform samples and Fig. 7(c) from projections taken with golden angle sampling scheme. Figure 7(d) shows the reconstruction obtained from the OS reconstruction applied on the same number of projections captured following the golden angle sampling scheme (Sec. 2.2), where the number of sub-sets is 14 and each sub-set contains 5 projections. Here, the explicit Radon transition matrix, A, is appropriate for the subset. It can be seen that the visual quality of OS reconstruction in Fig. 7(d) is higher compared to the similar CS reconstruction without OS shown in Fig. 7(b) and Fig. 7(c). The OS result has less radial artifacts which are present in non-OS reconstructions. The difference between the reconstructions is attributed to the difference between the projecting operators. The transition matrix yields less numerical errors and permits a precise adjoint operator (simply the adjoint of the matrix A), and thus yields a stable and accurate reconstruction.

 

Fig. 7 (a) is the reference image, (b) is the reconstruction in the golden sampling scheme, image (c) is the reconstruction in the uniform sampling scheme, image (d) is OS reconstruction using the explicit Radon transition matrix. All images are cropped to 1024 x 1024 pixels.

Download Full Size | PDF

5 Conclusions

We have presented the progressive compressive imaging system, implementing the golden ratio based fixed sampling step. Progressive sampling makes it possible to increase the number of measurements, by adding new projections to the existing ones without re-sampling, while each measurement addition increases the quality of the previous reconstruction. This approach is useful when no prior knowledge about the required number of samples for good reconstruction exists. Another advantage of the proposed sampling scheme is its robustness to sudden scanning stop. Large images can be reconstructed thanks to the fact that non overlapping subsets of measurements are created from successive projections and are used in a step-wise way, decreasing the dimension of the projection matrix.

The simulation and laboratory results show that in spite of the pseudo-randomness of the sampling pattern, the reconstruction results display insignificant differences in quality compared to the standard uniform sampling scheme. We have demonstrated optical compression ratios of 20:1 with satisfying reconstruction quality. The progressive sampling scheme, upon which we rely, enables excellent tradeoff between the compression and the quality of reconstruction.

It is also shown that the angular sampling used for progressive CI fits well an ordered subset (OS) reconstruction approach. The OS algorithm accelerates the reconstruction and facilitates the processing of large images.