## Progressive compressive imaging from Radon projections |

Optics Express, Vol. 20, Issue 4, pp. 4260-4271 (2012)

http://dx.doi.org/10.1364/OE.20.004260

Acrobat PDF (1060 KB)

### Abstract

In this work we propose a unique sampling scheme of Radon Projections and a non-linear reconstruction algorithm based on compressive sensing (CS) theory to implement a progressive compressive sampling imaging system. The progressive sampling scheme offers online control of the tradeoff between the compression and the quality of reconstruction. It avoids the need of a priori knowledge of the object sparsity that is usually required for CS design. In addition, the progressive data acquisition enables straightforward application of ordered-subsets algorithms which overcome computational constraints associated with the reconstruction of very large images. We present, to the best of our knowledge for the first time, a compressive imaging implementation of megapixel size images with a compression ratio of 20:1.

© 2012 OSA

## 1 Introduction

1. A. Stern and B. Javidi, “Random projections imaging with extended space-bandwidth product,” J. Disp. Technol. **3**(3), 315–320 (2007). [CrossRef]

7. A. Stern, O. Levi, and Y. Rivenson, “Optically compressed sensing by under sampling the polar Fourier plane,” J. Phys. Conf. Ser. **206**, 012019 (2010). [CrossRef]

1. A. Stern and B. Javidi, “Random projections imaging with extended space-bandwidth product,” J. Disp. Technol. **3**(3), 315–320 (2007). [CrossRef]

7. A. Stern, O. Levi, and Y. Rivenson, “Optically compressed sensing by under sampling the polar Fourier plane,” J. Phys. Conf. Ser. **206**, 012019 (2010). [CrossRef]

5. R. M. Willett, R. F. Marcia, and J. M. Nichols, “Compressed sensing for practical optical imaging systems: a tutorial,” Opt. Eng. **50**(7), 072601 (2011). [CrossRef]

6. A. Stern, “Compressed imaging system with linear sensors,” Opt. Lett. **32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

6. A. Stern, “Compressed imaging system with linear sensors,” Opt. Lett. **32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

8. E. J. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory **52**(2), 489–509 (2006). [CrossRef]

9. M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressed sensing MRI,” IEEE Signal Process. Mag. **25**(2), 72–82 (2008). [CrossRef]

6. A. Stern, “Compressed imaging system with linear sensors,” Opt. Lett. **32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

**32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

**32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

**32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

## 2 Progressive compressive sensing with fixed step angular sampling

### 2.1 Compressive sensing with optical Radon projections

**32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

*s*is the radial coordinate and

*g*(

*s*) is captured with a linear sensor which may be a one dimensional array of pixels. The linear sensor

**S**is aligned with the imaging axis

*y’*of the cylindrical lens

**L**and is located in the image plane of the system. In order to obtain the Radon projections at various angles, the lens

**L**and the sensor

**S**rotate concomitantly around the optical axis

*z*.

*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

**32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

7. A. Stern, O. Levi, and Y. Rivenson, “Optically compressed sensing by under sampling the polar Fourier plane,” J. Phys. Conf. Ser. **206**, 012019 (2010). [CrossRef]

5. R. M. Willett, R. F. Marcia, and J. M. Nichols, “Compressed sensing for practical optical imaging systems: a tutorial,” Opt. Eng. **50**(7), 072601 (2011). [CrossRef]

**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

^{12}) . 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

*Q*that will be sufficient for the reconstruction.

*q*

^{th}sample can be written as

**32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

*Q*is the number of projections. In fact only

*Q*, and since each object has different sparse representation, the number of projections needs to be adjusted by means of changing the angular step

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

*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

## 3 Reconstruction with ordered subsets

### 3.1 The Ordered Subsets reconstruction concept

13. H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imaging **13**(4), 601–609 (1994). [CrossRef] [PubMed]

*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

13. H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imaging **13**(4), 601–609 (1994). [CrossRef] [PubMed]

*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].

_{1}*M*being the subset order (the number of projections in a subset). The second subset is

_{t}*M*samples

_{t}*q*, acts on the sub-set

*G*. 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

_{q}*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

_{1}*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.

_{2}**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

## 4 Results

15. J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new twIst: two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Trans. Image Process. **16**(12), 2992–3004 (2007). [CrossRef] [PubMed]

### 4.1 Simulated experiment

*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:where

*R*is the reference image (Fig. 3(c)) and

*I*is the reconstructed image, both of size

### 4.2 Optical experiment

**32**(21), 3077–3079 (2007). [CrossRef] [PubMed]

**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.

*n*, so that the entire Radon transform will have approximately

*n*

^{2}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.

**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.

## 5 Conclusions

## Acknowledgment

## References and links

1. | A. Stern and B. Javidi, “Random projections imaging with extended space-bandwidth product,” J. Disp. Technol. |

2. | E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. |

3. | D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory |

4. | Y. Rivenson and A. Stern, “An efficient method for multi-dimensional compressive imaging,” Computational Optical Sensing and Imaging, COSI OSA Technical Digest (CD), paper CTuA4 (2009). |

5. | R. M. Willett, R. F. Marcia, and J. M. Nichols, “Compressed sensing for practical optical imaging systems: a tutorial,” Opt. Eng. |

6. | A. Stern, “Compressed imaging system with linear sensors,” Opt. Lett. |

7. | A. Stern, O. Levi, and Y. Rivenson, “Optically compressed sensing by under sampling the polar Fourier plane,” J. Phys. Conf. Ser. |

8. | E. J. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory |

9. | M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressed sensing MRI,” IEEE Signal Process. Mag. |

10. | H. Niederreiter, |

11. | M. Kleider, B. Rafaely, B. Weiss, and E. Bachmat, “Golden-Ratio sampling for scanning circular microphone arrays,” IEEE Trans. Audio, Speech, Lang. Process. |

12. | M. Livio, |

13. | H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imaging |

14. | H. Zaidi, |

15. | J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new twIst: two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Trans. Image Process. |

**OCIS Codes**

(110.0110) Imaging systems : Imaging systems

(110.1758) Imaging systems : Computational imaging

**ToC Category:**

Imaging Systems

**History**

Original Manuscript: December 20, 2011

Manuscript Accepted: January 27, 2012

Published: February 6, 2012

**Citation**

Sergei Evladov, Ofer Levi, and Adrian Stern, "Progressive compressive imaging from Radon projections," Opt. Express **20**, 4260-4271 (2012)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-4-4260

Sort: Year | Journal | Reset

### References

- A. Stern and B. Javidi, “Random projections imaging with extended space-bandwidth product,” J. Disp. Technol.3(3), 315–320 (2007). [CrossRef]
- E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag.25(2), 21–30 (2008). [CrossRef]
- D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory52(4), 1289–1306 (2006). [CrossRef]
- Y. Rivenson and A. Stern, “An efficient method for multi-dimensional compressive imaging,” Computational Optical Sensing and Imaging, COSI OSA Technical Digest (CD), paper CTuA4 (2009).
- R. M. Willett, R. F. Marcia, and J. M. Nichols, “Compressed sensing for practical optical imaging systems: a tutorial,” Opt. Eng.50(7), 072601 (2011). [CrossRef]
- A. Stern, “Compressed imaging system with linear sensors,” Opt. Lett.32(21), 3077–3079 (2007). [CrossRef] [PubMed]
- A. Stern, O. Levi, and Y. Rivenson, “Optically compressed sensing by under sampling the polar Fourier plane,” J. Phys. Conf. Ser.206, 012019 (2010). [CrossRef]
- E. J. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory52(2), 489–509 (2006). [CrossRef]
- M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressed sensing MRI,” IEEE Signal Process. Mag.25(2), 72–82 (2008). [CrossRef]
- H. Niederreiter, Uniform Distribution of Sequences (Dover Publications, 2006).
- M. Kleider, B. Rafaely, B. Weiss, and E. Bachmat, “Golden-Ratio sampling for scanning circular microphone arrays,” IEEE Trans. Audio, Speech, Lang. Process.18, 2091–2098 (2010).
- M. Livio, The Golden Ratio: The Story of Phi, the World's Most Astonishing Number (Broadway Books, 2003).
- H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Trans. Med. Imaging13(4), 601–609 (1994). [CrossRef] [PubMed]
- H. Zaidi, Quantitative Analysis in Nuclear Medicine Imaging (Springer, 2006).
- J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new twIst: two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Trans. Image Process.16(12), 2992–3004 (2007). [CrossRef] [PubMed]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.