Single-pixel imaging with origami pattern construction
Wen-Kai Yu Yi-Ming Liu






[1]Center for Quantum Technology Research, School of Physics, Beijing Institute of Technology, Beijing 100081, China. Correspondence and requests for materials should be addressed to W.-K.Y. (email: yuwenkai@bit.edu.cn)
[2]Department of Physics, Beijing Normal University, Beijing 100875, China




    Wen-Kai Yu^1,*, Yi-Ming Liu^2
    
=================================






Single-pixel compressive imaging can recover images from a small amount of measurements, offering many benefits especially for the scenes where the array detection is unavailable. However, the widely used random patterns fail to explore internal relations between the patterns and the image reconstruction. Here we propose a single-pixel imaging method based on origami pattern construction with a better imaging quality, but with less uncertainty of the pattern sequence. It can decrease the sampling ratio even to 0.5%, really realizing super sub-Nyquist sampling. The experimental realization of this approach is a big step forward toward the real-time compressive video applications.


Imaging is extremely important for acquiring the light field information of the target <cit.>. Single-pixel imaging <cit.>, as one of popular imaging technologies, can obtain images by using only a point detector without spatial resolution, thus it has extraordinary application prospects, especially at non-visible wavelengths. Except for point-by-point scanning, single-pixel imaging has also developed some imaging techniques based on Fourier <cit.>or Hadamard transforms <cit.> (i.e., complete orthogonal basis modulation, full sampling) and correlated/ghost imaging (GI) <cit.> (i.e., random modulation, oversampling). In recent ten years, with the rapid development of information theory, single-pixel imaging <cit.> based on compressed sensing <cit.> has come into being. This approach obtains images of sparse/compressive targets at very low sampling rates, breaking the Nyquist-Shannon sampling limit and bringing dawn of subsampling. However, compressive imaging sacrifices the measurement time and the computational complexity for acquiring spatial resolution, and generally uses random modulation. As a consequence, its sampling rate is still higher than 30% by a large probability, which has become the bottleneck of compressive imaging for realizing real-time practical applications, especially in large pixel-scale imaging. Recent studies have demonstrated that not all measurements contribute to the reconstructions of computational imaging, resulting in a technology named correspondence imaging (CI) <cit.>, where only the measurements well above or below the mean of detected values are used for reconstruction. Now that there exists such a phenomenon, we are thinking, why not just construct a universal deterministic (instead of random) measurement matrix where the modulated patterns with the most contributions to the reconstruction are preferentially displayed. In particular, an earlier research work <cit.> reported a ghost imaging method based on “Russian Dolls" (RD) ordering of the Hadamard basis, which provided us with inspiration.

In this Letter, we propose an origami model to construct deterministic patterns, decreasing the sampling ratio to 0.5%, which is ∼2 orders of magnitude lower than that of traditional compressive imaging. The row vectors stretched from each constructed pattern are orthogonal to each other, suitable for large-pixel-scale reconstruction. Therefore, it opens a door to the super subsampling applications for compressive video <cit.>.

Now let us begin with the definition of differential ghost imaging (DGI) <cit.>, one method that can greatly improve the quality of conventional GI. The bucket/single-pixel signal can be defined as S_B=∬ I_B(x_o,y_o)T(x_o,y_o)dx_ody_o, where I_B(x_o,y_o) and T(x_o,y_o) represents the intensity and transmission function at the spatial position (x_o,y_o) of the object x of p× q=N pixels, respectively. Similarly, the reference bucket signal can also be defined as S_R=∬ I_R(x,y)dxdy, where I_R(x,y) stands for the values at the spatial position (x,y) of the reference speckle patterns. By using a differential bucket signal S_Δ=S_B-⟨ S_B⟩/⟨ S_R⟩S_R instead of S_B in second-order correlation ⟨ S_B I_R(x,y)⟩, where ⟨⋯⟩ signifies the ensemble average, we can compute a differential ghost image U_DGI from M≫ N measurements via U_DGI=⟨ S_ΔI_R(x,y)⟩=⟨ S_B I_R(x,y)⟩-⟨ S_B⟩/⟨ S_R⟩⟨ S_R I_R(x,y)⟩. A high image quality can be guaranteed by DGI, but at the cost of oversampling.

Later, another method called CI <cit.> was reported, in which a positive (or negative) ghost image can be retrieved by only averaging some small fractions of the reference subset I_R^+ (or I_R^-), i.e., ⟨I_R_+⟩ (or ⟨I_R_-⟩), corresponding to {S_B^+|S_B≫⟨ S_B⟩} (or {S_B^-|S_B≪⟨ S_B⟩}). It can be seen that CI just need to calculate addition instead of second-order correlation, so CI takes less computation time than DGI. Although only a small number of patterns are involved in the calculation, the number of measurements is not decreased, and the conditional averaging of the reference patterns is made after oversampling. But it still can be concluded that not all the measured bucket/single-pixel intensities are necessary for reconstruction. Since single-pixel compressive imaging shares the same mathematical model with ghost imaging, this conclusion can be also applicable to compressive imaging. Obviously, the random patterns are not the best choices to minimize the number of measurements, because such a pattern has a high probability of transmitting or reflecting light at half of all pixels. It is really difficult for us to determine which pattern has a higher reconstruction contribution or produces a higher bucket intensity. However, we find that, if we encode on the spatial light modulator (SLM) with a pattern consisting of all-one, which can be treated as a connected domain, the detected intensity value will be the highest. If the pattern is divided equally into two connected domains with the white-black pixels 50%-50%, then the detected value will be the second highest with a high probability. On this basis, the problem can be partially deduced to the contributions of the piece number of connected domains (CDs) of patterns to the image quality. Here we propose a novel origami pattern construction method, which makes full use of symmetric reverse folding (reverse the values on the corresponding pixels at the symmetric positions) and the axial symmetry of the rescaled pattern. The steps of our origami pattern construction are as follows.

Step 1: Assume that the pixel-sizes of patterns are all p× q=n square matrices, with the same size as the object, where p, q (p=q) are all even numbers. Suppose there are n patterns in total, we divide them into n/4 groups, each of which with four patterns. The group number is denoted by i=1,2,3,⋯,n/4. Note, as an initial setting, the first pattern P_1,1 in the first group is a matrix with all ones. Then the second and the third patterns in this group P_1,2, P_1,3 are obtained by inversely folding the first one P_1,1 in this group up and down, left and right, respectively. That is, remain the upper (or the left) half unchanged while the lower (or right) half is axisymmetric with the upper (or left) half but with its values in the lower (or right) half taking the opposite. By this means, we realize black (-1) and white (+1) inverting. The fourth pattern P_1,4 in this group is formed by performing both up-down and left-right reverse folding operations (the sequence of operations is interchangeable) on the first pattern P_1,1. Similarly, for all the subsequent groups, the rest three patterns in the ith group, P_i,2, P_i,3 and P_i,4 will be generated from P_i,1 following the same processes.

Step 2: The first pattern P_i,1 in the ith group (also the 4(i-1)+1th pattern in the complete sequence) is built on the basis of the ith pattern in the current complete sequence. Take P_2,1 in the second group for example, we firstly scale both horizontal and vertical pixel dimensions of the second pattern in the complete sequence (i.e., P_1,2) to their 1/2 size, and place the scaled one on the upper left 1/4 part of a p× q square matrix with all zeros. Secondly, we perform the up and down, left and right axial symmetry about the midline lines on the vertical and horizontal axes of the matrix. After that, the pattern P_2,1 (also the 5th pattern in complete sequence) is generated. Then, repeat Step 1 to acquire the second to fourth patterns in the second group, P_2,2, P_2,3 and P_2,4 (also the 6th to 8th patterns in complete sequence). By analogy, we can create all n/4 groups. Fig. <ref>(a) gives a good illustration of Steps 1 and 2, and Fig. <ref>(b) shows the result obtained after Step 2.


Step 3: Adjust the pattern sequence to ensure that the number of CDs for each pattern in the ith group is in an incremental order. Here a CD is defined as an up-down-left-right connected area consisting of equal values (see the upper right corner of Fig. <ref>). The neighborhoods on both sides of the axis may cause the partition of CDs, which is worth of investigating. Let's first look at the pattern P_i,1 whose pixel values on both sides of the symmetry axis are the same, its contribution to the recovery is the largest in the ith group. Since the pattern P_i,4 needs double symmetric reverse folding, its number of CDs is largest (with lowest contribution) in the ith group. The patterns P_i,2 and P_i,3 are produced by only once symmetric reverse folding, whether to first perform up-down or left-right operation is approximately equivalent, thus they are particularly worth considering. As shown in Fig. <ref>, the patterns P_i,2 and P_i,3 with a group indication (ID) (note: not the absolute subscript) i=3, 9∼12, 15 need to be order reversed for the case of n=64, while those with a group ID i=3, 9∼12, 15, 33∼48, 51, 57∼60, 63 require adjustments for the case of n=256. This rule also exists in higher-order cases. For this, we set fourfold as a level. From top to bottom, the total n/4 groups (as a whole) can be catalogued into four parts: only the second part stays still, the groups in the third part all need to be adjusted, while the groups in the first and fourth parts also need to be adjusted but with the same group IDs, depending on the recursive layer. Then repeat the operation for the first part (as a new whole) until the number of groups in each new quarter part equals to 1. In the last layer, only the third part needs to be adjusted. Then, switch the order of the patterns P_i,2 and P_i,3 in these found groups.


If we perform left-right and up-down reverse folding to generate P_i,2 and P_i,3 in Step 1, respectively, then the ID positions should be adjusted accordingly in Step 3. After the above three steps, we will get the final pattern sequence. Fig. <ref> gives an example of the final pattern sequence for n=64. In order to deeply analyze the advantages of our method, we compare it with the RD method. It is interesting to find that these two approaches have similar numbers of CDs in low-order part (especially for the first 16 patterns, they are the same). Since their subsequent results are both derived from these first 16 patterns, it is advisable to set each 16 patterns as a comparison unit. Fig. <ref> shows the differences in their high-order parts (the last 64 patterns in 256-order sequence which also comprises the largest division of the RD method). As shown in Fig. <ref>(a), the RD ordering becomes very rough in this part, and has a lot of pattern pairs with the same number of CDs in each comparison unit. Because the maximum partition length of RD is n/4 while its minimum segmentation length is 4, which will inevitably incur too much uncertainty. Although the internal order exchange in these pattern pairs does not affect the RD sorting too much, it does bring hundreds of millions possibilities of finding the optimal pattern sequence to get the best imaging quality. While in our scheme, each four patterns are treated as a group, i.e., the partition length is always 4. As shown in Fig. <ref>(b), there are only four red marks in our method for the high-order part, with only 2^4 uncertainty. It makes the reconstructed images more accurate. Thus our method can be regarded as a better strategy.




Now, each pattern I_R can be sequentially reshaped into a row vector of 1× n, and then n such row vectors constitute a full-rank square matrix, which happens to be an orthogonal matrix. It is quite suitable for forming the measurement matrix A∈R^m× n of CS by just selecting the first m row vectors. Here we use a TVAL3 solver <cit.> to recover the image. For a fair comparison, we combine the RD ordering of the Hadamard basis with CS, instead of second-order correlation used in the original version <cit.>. Here we introduce an unitless performance quantitative measure, the root mean square error (RMSE), which is defined as RMSE=√(1/pq∑_i,j=1^p,q[Ũ(i,j)-U_o(i,j)]^2). It describes the square root of differences between the recovered image Ũ and the original image U_o for all pixels. Generally, the smaller the RMSE value, the better the quality of the image recovered. Since it is a pixel-by-pixel calculation, it cannot provide a very accurate characterization of performance. Therefore, in addition to the RMSE curves of CS, Russian dolls compressed sensing (RDCS) and origami pattern compressed sensing (OPCS), we also present their recovered results under different sampling ratios from 0.5%–8.5% with a 2% stepping increase, as shown in Fig. <ref>. It is obvious that the performance of our method outperforms those of CS and RDCS. Only 2.5% measurements are enough for OPCS to generate a nice reconstruction, and its limit sampling rate can reach 0.5%.

Our experimental realization is based on a conventional single-pixel camera, as shown in Fig. <ref>. The thermal light from a stabilized halogen tungsten lamp with a wavelength of 360∼2600 nm is collimated and attenuated via a beam expander and a neutral density filter (NDF), respectively. Here the NDF is used to attenuate the light to the ultra-weak light level. After that, the light passes through a negative 1951 USAF resolution test chart, and then it is imaged onto a digital micromirror device (DMD) via an imaging lens and a mirror. The applied DMD consists of 1024×768 micromirrors, each of which can be switched between two directions of ±12^∘, corresponding to 1 and 0. Since the values of our patterns are either 1 or -1, we need to divide each generated pattern I_R into a complementary matrix pair <cit.> Î_R=(I_R+1)/2 and Ǐ_R=1_N× N-Î_R, where I_R=Î_R-Ǐ_R, and 1_N× N stands for an all 1 matrix. Then, these complementary matrices are sequentially encoded on the DMD. A counter-type Hamamatsu H10682-210 PMT is placed on one of the reflection orientations to accordingly make a differential measurement of two adjacent total photon counts. By using OPCS algorithm, we can acquire a very satisfactory image quality, as shown in the illustration of Fig. <ref>.


In conclusion, a single-pixel compressive imaging method with origami pattern construction is proposed and demonstrated experimentally with differential modulation. In our scheme, following the fluctuation characteristic of measurement values, the patterns with the most contributions to the recovery should be preferentially displayed. By symmetric reverse folding, axial symmetry and partial pattern order adjustment, the generated deterministic patterns ensure a better image quality but with less uncertainty of the pattern sequence, compared with traditional CS and RDCS. The sampling ratio can be reduced to 0.5%, even about 2 orders of magnitude lower than that of traditional CS, actually realizing super sub-Nyquist sampling. Therefore, our scheme can be used to greatly improve real-time performance of single-pixel compressive video applications, especially where array detectors are unavailable.



§ 

This work is supported by the National Natural Science Foundation of China (61801022), the Beijing Natural Science Foundation (4184098), the National Key Research and Development Program of China (2016YFE0131500), the International Science and Technology Cooperation Special Project of Beijing Institute of Technology (GZ2018185101), and the Beijing Excellent Talents Cultivation Project - Youth Backbone Individual Project.



§ 

The authors warmly appreciate Shuo-Fei Wang and Ya-Xin Li for preparing Fig. <ref>.














 
  
  



99
YuOC2017sW.-K. Yu, A.-D. Xiong, X.-R. Yao, G.-J. Zhai, and Q. Zhao, Opt. Commun. 402, 413 (2017).

Edgar2019sM. P. Edgar, G. M. Gibson, and M. J. Padgett, Nat. Photonics 13, 13 (2019).
BSun2013sB. Sun, M. P. Edgar, R. Bowman, L. E. Vittert, S. Welsh, A. Bowman, and M. J. Padgett, Science 340, 844 (2013).
Radwell2014sN. Radwell, K. J. Mitchell, G. M. Gibson, M. P. Edgar, R. Bowman, and M. J. Padgett, Optica 1, 285 (2014).
Bian2016sL. Bian, J. Suo, G. Situ, Z. Li, J. Fan, F. Chen, and Q. Dai, Sci. Rep. 6, 24752 (2016).

Zhang2015sZ. Zhang, X. Ma, and J. Zhong, Nat. Commun. 6, 6225 (2015).

Huynh2016sN. Huynh, E. Zhang, M. Betcke, S. Arridge, P. Beard, and B. Cox, Optica 3, 26 (2016).
MJSunNC2016sM.-J. Sun, M. P. Edgar, G. M. Gibson, B. Sun, N. Radwell, R. Lamb, and M. J. Padgett, Nat. Commun. 7, 12010 (2016).
Phillips2017sD. B. Phillips, M.-J. Sun, J. M. Taylor, M. P. Edgar, S. M. Barnett, G. M. Gibson, and M. J. Padgett, Sci. Adv. 3, e1601782 (2017).

Shih1995sD. V. Strekalov, A. V. Sergienko, D. N. Klyshko, and Y. H. Shih, Phys. Rev. Lett. 74, 3600 (1995).
Shapiro2008sJ. H. Shapiro, Phys. Rev. A 78, 061802 (2008).
YuOE2014sW.-K. Yu, M.-F. Li, X.-R. Yao, X.-F. Liu, L.-A. Wu, and G.-J. Zhai, Opt. Express 22, 7133 (2014).
Gong2016sW. Gong, C. Zhao, H. Yu, M. Chen, W. Xu, and S. Han, Sci. Rep. 6, 26133 (2016).
Liu2017sH.-C. Liu and S. Zhang, Appl. Phys. Lett. 111, 031110 (2017).

Baraniuk2008sM. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F. Kelly, and R. G. Baraniuk, IEEE Signal Process. Mag. 25, 83 (2008).

Donoho2006sD. Donoho, IEEE Trans. Inform. Theory 52, 1289 (2006).
Candes2006sE. J. Candès, in Proc. Int. Cong. Math., European Mathematical Society, Madrid, Spain, 3, 1433 (2006).

Luo2011sL.-A. Wu and K.-H. Luo, in AIP Conf. Proc. 1384, 223 (2011).
Luo2012sK.-H. Luo, B. Huang, W.-M. Zheng, and L.-A. Wu, Chin. Phys. Lett. 29, 74216 (2012).
Shih2011sR. E. Meyers, K. S. Deacon, and Y. Shih, Appl. Phys. Lett. 100, 131114 (2012).
YuCPB2015sW.-K. Yu, X.-R. Yao, X.-F. Liu, L.-Z. Li, and G.-J. Zhai, Chin. Phys. B 24, 054203 (2015).
MJSunAO2015sM.-J. Sun, M.-F. Li, and L.-A. Wu, Appl. Opt. 54, 7494 (2015).

MJSunSR2017sM.-J. Sun, L.-T. Meng, M. P. Edgar, M. J. Padgett, and N. Radwell, Sci. Rep. 7, 3464 (2017).

Zheng2009sJ. Zheng and E. L. Jacobs, Opt. Eng. 48, 087006 (2009).
Yu2018sP. Paces, W.-K. Yu, and J. Klesa, in IEEE/AIAA 37th Digital Avionics Systems Conference, London, United Kingdom, 23 (2018).

Gatti2010sF. Ferri, D. Magatti, L. A. Lugiato, and A. Gatti, Phys. Rev. Lett. 104, 253603 (2010).

CBLi2010sC. B. Li, Masters of Science thesis, Rice University, (2010).

YuSR2014sW.-K. Yu, X.-F. Liu, X.-R. Yao, C. Wang, Y. Zhai, and G.-J. Zhai, Sci. Rep. 4, 5834 (2014).
YuAO2015sW.-K. Yu, X.-R. Yao, X.-F. Liu, L.-Z. Li, and G.-J. Zhai, Appl. Opt. 54, 363 (2015).
YuOC2016sW.-K. Yu, X.-R. Yao, X.-F. Liu, R.-M. Lan, L.-A. Wu, and G.-J. Zhai, Opt. Commun. 371, 105 (2016).

