Open Access
Add to BibSonomy
Abstract
Single-pixel imaging is a technique that can reconstruct an image of a scene by projecting a series of spatial patterns on an object and capturing the reflected light by a single photodetector. Since the introduction of the compressed sensing method, it has been possible to use random spatial patterns and reduce its number below the Nyquist-Shannon limit to form a good quality image but with lower spatial resolution. On the other hand, Hadamard pattern based methods can reconstruct large images by increasing the acquisition measurement time. Here, we propose an efficient strategy to order the Hadamard basis patterns from higher to lower relevance, and then to reconstruct an image at very low sampling rates of at least 8%. Our proposal is based on the construction of generalized basis vectors in two dimensions and then ordering in zigzag fashion. Simulation and experimental results show that the sampling rate, image quality and computational complexity of our method are competitive to the state of the art methods.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Imaging is commonly performed using 2D array of sensors such as digital cameras at the visible spectral region. In contrast, single-pixel imaging (SPI) uses only a single detector to retrieve an image of a scene [1,2]. The SPI technique can reconstruct an image by projecting spatially modulated patterns while measuring the correlated intensities using a photodetector. Among the various applications of SPI are the implementation of low-cost imaging systems at non-visible wavelengths, multispectral [3,4], hyperspectral imaging with applications to biomedical imaging [5], the reconstruction of three-dimensional video at a frame rate of 12 Hz [6], the development of secure imaging systems, and visual cryptography schemes for QR-codes [7,8].
SPI is related to the ghost imaging technique that reconstructs an image using the information of two correlated coherent beams that are sensed by a camera and a photodetector [9]. This technique has also been implemented using entangled photons where one of the photons interacts with the object and a photodetector, while the second one is sensed by a camera [10]. On the other hand, the fundamental theory of compressive sampling (CS) allows the development of single-pixel cameras that take advantage of the sparsity of the target and the interaction of random patterns with the scene [11,12]. CS belongs to the class of computational imaging methods where the image is reconstructed through an optimization process. In particular, in CS, an undetermined system of linear equations is solved by imposing the restriction of sparsity, that is, signals with relatively few terms. This sparsity/compression criterion allows an image to be reconstructed with a series of measurements below the Nyquist-Shannon sampling limit [13]. Recently, the sampling basis method proposed uses an orthogonal basis pattern to obtain the measurements, for example, an orthogonal basis formed with Hadamard matrices [14], Fourier [15,16] and wavelets [17]. In particular, Hadamard basis is suitable for experimental realization due to their binary nature of patterns that can be projected using DMD (Digital Micromirror Device) technology as spatial light modulator (SLM) [18].
In the Hadamard basis, a scan method is required to obtain all measurements described as an inner product between the patterns and the object. This last process takes a long time and, in principle, reconstructs an image perfectly. However, if there are time constraints, as in real-time applications, we must obtain the measurements by choosing the most relevant patterns. The first option to deal with this problem is to use the well-known orderings of the Hadamard matrix, such as the dyadic ordering, the Walsh-Paley Hadamard matrix, and the normalized Haar Transform [19]. On the other hand, there are several proposals that take advantage of the two-dimensional structure of Hadamard patterns to sort it in descending order of importance. In [20] proposed the Russian Doll (RD) method, which is based on ordering the basis patterns with respect to their significance when using natural images. The main idea of the RD method is to exploit the symmetry of the Hadamard matrix and assume that the fewer blocks a pattern contains, the more relevant it will be. A similar approach was proposed in the so-called Cake-Cutting (CC) [21], based on the assumption that the fewer connected regions (consisting of -1s or 1s) a basis pattern contains, the more significant measurement is obtained. In [22], the Origami Pattern Construction (OPC) method based on the CC method was proposed. In the OPC method, patterns are constructed by taking advantage of the symmetry and axial symmetry of the patterns. As in [23], other approaches propose to order the Hadamard basis patterns in ascending order with respect to their total variation and total wavelet transformed coefficients in ascending order to have the best performance. In [24], a general strategy was proposed to construct the Hadamard basis for SPI with multi-resolution imaging applications based on the Hadamard pipeline encoder (HPE). The HPE allows the construction of high-order basis patterns according to four basic rules and generates other orderings such as the RD. As a result, the PHE requires much less memory consumption for multi-resolution imaging applications and has better computational performance than traditional basis pattern construction.
The sorting methods aforementioned, such as RD, CC, and OPC, are based on determining the number of regions or blocks connected within a pattern. Counting the number of connected regions is a time-consuming process that increases with the dimension of the Hadamard basis. Therefore, it is important to have an efficient procedure for selecting the most significant patterns or, equivalently, ordering the Hadamard basis patterns from highest to lowest importance. Here, we propose an efficient method that orders the Hadamard basis patterns in computational complexity of $O(N^2)$ and a performance similar to the CC method considered as one of the best ordering method [21], where $N=2^n$ and $n$ is the dimension of the Hadamard matrix. Our proposal consists of two steps: (1) a two-dimensional representation of the Hadamard basis patterns of dimension $2^n$ is constructed using an initial ordering of the sub Hadamard matrices of dimension $2^{n/2}$, and (2) a zigzag traversal of the two-dimensional representation given in step (1) is made to select the patterns with the highest energy. In the first step, the initial ordering is obtained by sorting the basis patterns with respect to the sign change in ascending order. The second step is based on the fact that the energy distribution accumulates in the upper left corner of the two-dimensional representation of the basis patterns [25]. Here, instead of using the natural ordering of the Hadamard matrix when constructing its four-dimensional representation as in [25], a different ordering that improves the performance of the SPI method at very low sampling rates is used.
The organization of this article is as follows: in section 2, the theory and methods are given; in section 3, the simulated and experimental results are shown; and section 4 provides conclusions and discussion.
2. Theory and methods
SPI consists of projecting a series of patterns onto a scene and integrating the reflected light with a bucket detector i.e., using a photodetector. When a pattern is reflected by the SLM, then the photodetector integrates the interaction of the pattern with the scene to obtain a measurement or bucket detection. A measurement can be modeled as a linear process as follows:
where ${\mathbf f}\in \mathbb {R}^d$ is the desired image of the scene with $d$ the total number of pixels in the image, i.e. $d=\mbox {height}\times \mbox {width}$. The matrix $P=\left [ {\mathbf p}_1 \ \cdots \ {\mathbf p}_k \right ]^T\in \mathbb {R}^{k\times d}$ whose rows ${\mathbf p}_j\in \mathbb {R}^d$ are patterns for $j=1,\dots,k$, and ${\mathbf m}\in \mathbb {R}^k$ are the corresponding measurements. In the simplified model given in Eq. (1), the contribution of the dark current caused by the photodetector is omitted, which is justified when using the differential method [14,26]. The process of obtaining the ${\mathbf f}$ image given the set of $P$ patterns and ${\mathbf m}$ measurements is called restoration. There are several restoration algorithms for solving the given system in Eq. (1) that can be grouped into three categories [27]: (a) non-iterative methods based on solving the system in Eq. (1) by matrix inversion when possible, (b) linear iterative methods using the gradient descent optimization algorithm, the Poisson maximum likelihood method and alternating projection, and (c) nonlinear iterative methods using methods such as the TV regularization algorithm [28] and the NESTA algorithm [29].Figure 1 shows the experimental configuration for SPI that uses a DMD as SLM and a single photodetector. The DMD is used to filter the incident light on it and depending on the angle of each micromirror. Then, the backscatter light of the scene is measured by the photodetector. The DMD and photodetector are synchronized by the computer which also stores all the measured values. As can be seen, the DMD is limited to generating only binary patterns, so the ${\mathbf p}_j$ pattern series must be binary arrays. In [17], a splitting strategy was shown to project integer patterns ${\mathbf p}_j\in \mathbb {Z}^d$ using a DMD that is based on splitting a positive pattern into a negative and a positive components, then each component is decomposed as a weighted sum of binary patterns. In this research, we concentrate on binary patterns called Hadamard matrices and therefore do not required the splitting strategy, rather than the differential method being used.
Fig. 1. Schematic of the experimental configuration for SPI: the DMD is used to partially filter the incident light by individually moving each micromirror at an angle of $\pm 12^{\circ }$. The object (O) in the scene is illuminated by white LED diode and the detector (PD) measures the backscattered light collected by lens L1, L2, and M.
A square matrix $H_N$ of order $N$ with elements $-1$ and $+1$ is called a Hadamard matrix if the following statement is satisfied
where $I_N$ is the identity matrix of order $N$.It is known that if there exists a Hadamard matrix of order $N>2$, then $N \equiv 0 \ (\mbox {mod} \ 4)$. The latter statement is called the Hadamard conjecture which remains unproved. Hadamard matrices of order $N=2^n$ where $n$ is an integer are called Sylvester’s matrices. In this case, if $H$ is a Hadamard matrix of order $N$, then a matrix of order $2N$ is obtained as
For instance, the Hadamard matrices of order 2 and 4 are given by
In particular, the Hadamard transform has the property that it is a symmetric matrix, as well as its inverse. For any positive $n\geq 1$, each entry of $H_N=[\mbox {wal}_h(x,y)]$ with $x,y=0,\dots,N-1$, can be expressed as
with $x_i, y_i$ are the $i$-th coefficients in the binary expansion of $x,y$ using $n$ bits.The SPI based on Hadamard transform consists on constructing a Hadamard matrix whose $j$ rows are used as patterns, then the measurements are obtained as
where the restored image ${\mathbf f}$ has a dimension of a power of 2. For instance, if we want to restore an image of dimension $d=2^n\times 2^n=2^{2n}$ then the order of the Hadamard matrix is $N=2^{2n}$. The restored image after the projection of all basis patterns is theoretically obtained as where ${\mathbf m}$ is a vector with all $m_j$ measurements. To deal with negative $H_N$ values, the differential method is used in which two matrix patterns are defined as where $H_N^{+}$ and $H_N^{-}$ are binary patterns such that the matrix $H_N^{+} = (H_N+{\mathbf 1}_N)/2$ is zero in each negative entry of $H_N$, the matrix $H_N^{-}=({\mathbf 1}_N - H_N)/2$ is zero in each positive entry of $H_N$, and ${\mathbf 1}_N$ is a matrix of ones of order $N$. Therefore, the measurements ${\mathbf m}_{+}$ and ${\mathbf m}_{-}$ acquired using the patterns $H_N^{+}$ and $H_N^{-}$, respectively, are subtracted to obtain the measurement ${\mathbf m}={\mathbf m}_{+}-{\mathbf m}_{-}$.In practice, for a typical configuration like the one shown in Fig. 1, to obtain the ${\mathbf m}_{+}$ and ${\mathbf m}_{-}$ measurements, the DMD mirrors will adjust their position appropriately for each of the deployed Hadamard patterns $H_N^{+}$ and $H_N^{-}$, respectively. Once in position, the photodetector will record a voltage value corresponding to each pattern, which allows the measurement ${\mathbf m}$ to be obtained. This procedure is repeated for each Hadamard pattern $H_N$ that is used. Therefore, once all the measurements ${\mathbf m}$ obtained by the photodetector are available, they are used together with the $H_N$ patterns as inputs of a reconstruction algorithm to obtain the ${\mathbf f}$ image, as described in [30]. Section 3.3 shows the experimental results of our proposal.
2.1 Generalized basis vectors
Here we consider a generalization of basis vectors provided in [31]. Let $V_s^{(k)}$ be the set of one-dimensional basis vectors as follows:
If $|{\mathbf s}|=t$ is the length of the seed vector ${\mathbf s}$, then the set $V_s^{(k)}$ forms an orthogonal set of $2^k$ vectors of length $2^k t$. The construction of the set $V_s^{(k)}$ can be visualized as a binary tree of depth $k$. Figure 2 (top tree) shows an example of the set $V_s^{(k)}$ for $k=2$ using the seed vector ${\mathbf s}=[1]$. The nodes at level $i$ represent the vectors of $V_s^{(i)}$, and the leaves of the tree represent the four vectors of $V_s^{(2)}$. It can be seen that the vectors in $V_s^{(2)}$ coincide with the row vectors of the Hadamard matrix of size $4\times 4$. The branches of the tree are labelled with the sign +/- of the values $\alpha \in \{+1,-1\}$ used to create the vectors. The sequence $\alpha = \alpha _1,\dots,\alpha _k$ along the branches of the tree uniquely defines every vector in ${\mathbf v}\in V_s^{(k)}$, which will be called the $\alpha$-index of ${\mathbf v}$. It is said that two vectors ${\mathbf v}_i,{\mathbf v}_j\in V_s^{(k)}$ are $\alpha$-related if and only if the Hamming distance from their $\alpha$-index is one. In particular, if an ordered set of vectors ${\mathbf v}_1,\dots,{\mathbf v}_n\in V_s^{(k)}$ that are consecutively $\alpha$-related then it is called a gray code sequence (GCS).
Fig. 2. (top tree) one-dimensional orthogonal set $V_s^{(k)}$ for $k=2$ using the seed vector ${\mathbf s}=[1]$. (full picture) two-dimensional orthogonal set $V_s^{(k_1,k_2)}$ for $k_1,k_2=2$ using the seed vector ${\mathbf s}=[1]$. Below each 2D pattern is the corresponding $\alpha$-index which can be seen as a binary string $\rho$ of length $n/2$.
The two-dimensional version of the set $V_s^{(k)}$ can be defined as
2.2 Ordering of the Hadamard basis
In the SPI method based on the Hadamard transform, an image can be restored by projecting all basis vectors. However, if there are time constraints, such as in video restoration, then we must select the most significant basis vectors in order to achieve the restoration process. This is equivalent to ordering the rows of the Hadamard matrix in descending order of relevance and only using a percentage of the rows according to the stated sorting. To achieve this objective, it is necessary to define a criterion for assigning a numerical value of relevance to each row of the Hadamard matrix.
The proposed ordering method consists of two steps: (1) the 2D basis vectors are created using an initial ordering of the one-dimensional basis vectors used in the outer product. This is based on the fact that by changing the order of the basis vectors in the outer product a different indexing of the 2D patterns is obtained. (2) the 2D patterns are traversed in a zigzag fashion from the left upper pattern to the right bottom pattern. The zigzag traversal allows us to choose the patterns with highest energy or that contribute most to the image restoration [25]. In the following, we provide a detailed description of the proposed ordering method.
The 2D basis vectors defined in Eq. (11) give us an indexing of the Hadamard patterns or rows. For instance, for a Hadamard matrix $H_N$ with $N=2^n$ each row has length $N$ which can be seen as a pattern of size $2^{n/2}\times 2^{n/2}$ pixels. The matrix $H_N$ is equivalent to construct the basis set $V_{[1]}^{(k_1,k_2)}$ for $k_1,k_2=n/2$, where for any ${\mathbf v}\in V_{[1]}^{(k_1,k_2)}$ we have that ${\mathbf v}={\mathbf v}_1\times {\mathbf v}_2$ and ${\mathbf v}_1, {\mathbf v}_2\in V_{[1]}^{(n/2)}$. In other words, each 2D basis vector in the outer product of $V_{[1]}^{(n/2)}$ by itself, corresponds to a pattern of size $2^{n/2}\times 2^{n/2}$ pixels and also coincides with a row of $H_N$. Indexing 2D patterns in $V_{[1]}^{(n/2,n/2)}$ depends on the order of the basis vectors ${\mathbf v}\in V_{[1]}^{(n/2)}$ since ${\mathbf v}$ has associated a unique $\alpha$-index that can be seen as binary string $\rho$ of length $n/2$ where $\rho _i=0(1)$ if $\alpha _i=-1(+1)$ for $i=1,\dots,n/2$. Let $\sigma$ be a permutation of the set of integers $S=\{1,\dots,n/2\}$ and let ${\mathbf v}_{\sigma (i),\sigma (j)}$ be the indexing of the basis vectors in $V_{[1]}^{(n/2,n/2)}$ where ${\mathbf v}_{\sigma (i),\sigma (j)}={\mathbf v}_{\sigma (i)}\times {\mathbf v}_{\sigma (j)}$ with ${\mathbf v}_{\sigma (i)}, {\mathbf v}_{\sigma (j)}\in V_{[1]}^{(n/2)}$ and $i,j\in S$. Thus, the permutation $\sigma$ imposes an ordering of the basis vectors ${\mathbf v}_{\sigma (k)}$, and consequently of the 2D patterns in $V_{[1]}^{(n/2,n/2)}$, where the index $\sigma (k)$ corresponds to the decimal value of the corresponding binary string $\rho$ of ${\mathbf v}$. For example, Fig. 2 shows the $\alpha$-index of each pattern, and Fig. 3(a) shows the decimal values associated with each $\alpha$-index. Notice that the given indexing depends on the order of the one-dimensional basis vectors used to construct $V_{[1]}^{(n/2,n/2)}$.
Fig. 3. Example of an indexing and the zigzag traversal: (a) Indexing of the 2D patterns of the basis vectors in $V_{[1]}^{(n/2,n/2)}$ as constructed in Fig. 2. The arrows indicate the traversal order of the 2D patterns or zigzag scanning. (b) Sorting of patterns using the GCS+S method and the Cake-Cutting. The number above each pattern is its index and the number in parentheses is the number of blocks that the pattern contains.
In [25], it was shown that a zigzag scan of a four-dimensional rearrangement of the Hadamard basis patterns, allows to take advantage of the concentration of energy in the left upper corner that decreases vertical and horizontally. However, the authors only consider traversing the basis vectors in a fixed order using the natural ordering of the Hadamard matrix. Here, we found that changing the indexing of the basis vectors allows us to improve the sampling rate by using the most important vectors to restore the image. Figure 3(a) shows the zigzag scanning of the 2D basis vectors using a specific indexing, where the traversal starts at the left upper corner and finishes in the bottom right corner. We choose an indexing derived by permuting the one-dimensional vectors of $V_{[1]}^{(n/2)}$ in ascending order according to the sign change (or Sequency ordering) which is the best order when there is no knowledge about the target scene. The combination of the Sequency ordering in conjunction with the zigzag traversal will be called the GCS+S method, since the first one coincides with a GCS. Figure 3(b) shows the zigzag traversal with the proposed indexing, and a comparison with the CC method. In case of the GCS+S method, the patterns are listed in the order of traversal together with their corresponding index and the number of blocks it contains. In addition, it is shown the ordering obtained using the CC method of the Hadamard matrix of size $2^4\times 2^4$. Note that GCS+S and CC methods are similar with respect to the number of blocks; for instance, if we compare the two orderings, the number of blocks coincide in exactly 10 positions and differ in 6, and the difference in the number of blocks is at most 2.
In the next section, we show by numerical simulations that the GCS+S method has at least the same performance and better algorithm complexity than the CC method for image restoration.
3. Results and discussion
3.1 Performance measures and setup
This section shows the results of our proposed sorting method through numerical simulations and experimental validation. To measure performance, the root mean square error (RMSE) and the structural similarity index (SSIM) are used. The RMSE is calculated using the following equation:
where $I_0,I_1$ are 2D images of size $m\times n$ pixels. Generally, the RMSE measures the difference between a reference image $I_0$ and a target image $I_1$. In the context of SPI, the lower the RMSE value, the better the quality of the target or restored image.The SSIM($x$, $y$) is calculated as follows [32]:
The calculations were performed on a laptop with an Intel processor i7-8750H CPU at 2.20GHz, 6 cores and 16 GB of RAM. The programming language was Matlab R2017b under the Windows 10 operating system. The test images were selected from two datasets: the USC-SIPI image database [33] was used as proof of concept of our proposed method, for example, to show the results of image restoration, and the COCO database [34] consisting of 85783 images of more than 80 object categories of natural scenes, including 250000 people with key points. All images used were converted to gray scale when necessary and normalized to the range from 0 to 1, and re-scaled using bilinear interpolation. The source code of the implemented GCS+S algorithm can be downloaded at [https://bit.ly/35ekSPU]. Our contribution is expected to motivate further research and allow systematic comparisons with other proposals.
3.2 Numerical simulations
An important step in constructing the 2D representation of the Hadamard matrix as described in Section 2.2, is the selection of the permutation used for indexing the basis vectors. For example, to restore an image of size $64\times 64$ pixels, it is required to construct the set $V^{(6,6)}_{[1]}$ which is equivalent to the Hadamard matrix $H_{2^{12}}$. The number of elements ${\mathbf v}_{\sigma (i),\sigma (j)}\in V^{(6,6)}_{[1]}$ is $2^{12}$ and each ${\mathbf v}_{\sigma (i),\sigma (j)}$ is a 2D basis vector of size $64\times 64$ obtained as the outer product of two vectors as ${\mathbf v}_{\sigma (i),\sigma (j)}={\mathbf v}_{\sigma (i)}\times {\mathbf v}_{\sigma (j)}$ where ${\mathbf v}_{\sigma (i)},{\mathbf v}_{\sigma (j)}\in V^{(6)}_{[1]}$. The permutation $\sigma$ is a one-to-one mapping of the integers $\{1,\dots,2^6\}$ or equivalently of the index rows of the Hadamard matrix $H_{2^6}$. To justify our selection of the permutation $\sigma$ of the basis vectors in $V^{(6)}_{[1]}$, we compare several well known orderings of the Hadamard matrix $H_{2^6}$. After choosing a permutation, the ordering of the basis vectors in $V^{(6,6)}_{[1]}$ is obtained by the zigzag traversal.
Figure 4 shows the results of restoring images of size 64$\times$64 pixels using an image from the USC-SIPI database for various permutations. Figure 4 shows the RMSE of image restoration with respect to sampling rate from 10% to 50% in steps of 10%. Here, the sampling rate refers to the percentage of rows of the Hadamard matrix taken from the first row up to a given number of rows. Image restoration was obtained using Eq. (7) and the differential method given in Eq. (8), where the measurement values corresponding to the non-projected patterns are set to zero. The permutation of the well-known Hadamard matrices was used [19]: Cal-Sal, natural ordering, Walsh-Paley, the Sequency obtained by ordering in ascending order with respect to the sign change per row, and a random permutation. Although Fig. 4 shows a single performance example, in most cases using the USC-SIPI database with 39 images, the Sequency ordering obtains the best performance, followed by the Walsh-Paley ordering, while the other three matrices obtained the worst performance. As the figure shows, the image restoration results with sampling rates from 10% to 50% along with the computed RMSE and SSIM values with respect to the original image. As can be seen, the Sequency and Walsh-Paley orderings show better visual restoration results. It is important to remember that the different orderings or permutations are used to construct the 2D representation of the Hadamard matrix, in conjunction with the zigzag traversal. Figure 4 justifies why the Sequency ordering was chosen to construct the set $V^{(6,6)}_{[1]}$, which we call the GCS+S method. In simulation and experimental settings, a restored image is considered good if the RMSE is less than 0.1 and the SSIM value is greater than 0.6 [35]. It can be seen in Fig. 4, that the GCS+S method obtains a good restoration result from a sampling rate of 10%.
Fig. 4. Image restoration results using various sampling rates. Visual comparison of image restoration for permutation matrices with sampling rates from 10% to 50%, along with the RMSE and SSIM values with respect to the original image.
We made a comparison of our GCS+S sorting method against CC and the well-known Hadamard, Walsh-Paley and Sequency matrices. In the case of the Walsh-Paley and Sequency matrices, they were constructed using the usual definition as described in [19]. These last two matrices are considered as a performance reference and are different from the matrices built using the zigzag traversal. On the other hand, the CC method was chosen because is one the best algorithms to order the patterns or rows of the Hadamard matrix, and uses the least number of them to restore an image [21]. Figure 5 shows the results of the comparison using images in the COCO database and re-scaled to 64$\times$64 pixels. In Fig. 5(top) we show the restoration result for three images with sampling ratio from 10% to 40% in 10% steps. For each method, the restored images are shown along with their RMSE and SSIM values. It can be seen that the Walsh-Paley and Sequency orderings perform the worst, even at the 20% sample rate, the restored images have poor visual quality with SSIM values below 0.5. Similarly, the GCS+S method shows better similarity values that can be seen at low sampling rate. To show that the above results are not a particular case, Fig. 5(bottom) shows the average RMSE and SSIM values for a set of 150 images randomly taken from the COCO database. As expected, the Walsh-Paley and Sequency orderings report large error values and low similarity values. On the other hand, GCS+S has better mean RMSE values at low sampling rates of 5%, 10% and 20%. However, at the 5% sampling rate, the SSIM value of the GCS+S method does not reach the lower limit of 0.6, which is considered a good quality image; this limit is reached approximately at 8%. Furthermore, it can be observed in Fig. 5(bottom) that the variation of the RMSE and SSIM values decreases as the sampling rate increases.
Fig. 5. Comparison of GCS+S and CC methods, and the Walsh-Paley and Sequency orderings: (top) image restoration results for three images from sampling rates of 10% to 40% in 10% steps. For each restored image, the RMSE and SSIM values are shown. (bottom) Plots of mean RMSE and SSIM values with vertical error bars of equal length for 150 randomly chosen images from the COCO database. The horizontal dashed lines in the plots serve as reference values at 0.1 and 0.6 for the RMSE and SSIM curves, respectively.
An important aspect of SPI image restoration methods is that they must be tolerant to the noise sources at low sampling rates. The most common sources of noise in SPI are the dark current caused by the bucket detector and ambient light. In the former, it can be reduced significantly using the SPI differential method, that is, decomposing each Hadamard pattern into positive and negative input values as binary patterns [14]. In the second one, it is more difficult to model since the complexity of the scattering of the ambient light. Figure 6 shows the result of image restoration using in the COCO database and re-scaled to 64$\times$64 pixels where measurements were corrupted by white Gaussian noise at SNR=50dB. Here we consider a fixed amount of noise since the behavior of the sorting methods was similar. In Fig. 6(top) we show image restoration for three images using the GCS+S, CC, Walsh-Paley, and Sequency orderings. The Walsh-Paley and Sequency orderings have the worst results, and in the GCS+S and CC methods, the RMSE and SSIM values are acceptable and have good image restoration. However, in the case of GCS+S, it obtains better results at sampling rates of 10% and 20%. Similarly, Fig. 6(bottom) shows the mean values with error bars of length given as the standard deviation, for various sampling rates. It can be seen that the RMSE and SSIM values do not converge to zero and one, respectively, as in the case of ideal measurements, as in Fig. 5(bottom). However, the GCS+S and CC methods obtain acceptable restoration results at sampling rates above 8%. Note that the SSIM value has a greater variation than the RMSE value, which means that it is more difficult to obtain good image quality.
Fig. 6. Comparison of GCS+S and CC methods, and the Walsh-Paley and Sequency orderings under corrupted measurements by white Gaussian noise with SNR=50dB: (top) image restoration results for three images from 10% to 40% sampling rates in 10% steps. For each restored image, the RMSE and SSIM values are shown. (bottom) Plots of mean RMSE and SSIM values with vertical error bars of equal length for 150 randomly chosen images from the COCO database. The horizontal dashed lines in the plots serve as reference values at 0.1 and 0.6 for the RMSE and SSIM curves, respectively.
One of the important factors in SPI is the sparsity of the scene or image to be restored. In particular, the CS method can reduce the number of measurements by taking advantage of the sparsity of the image. In general, the sparsity decreases as the size of the image increases. To see the influence of the sparsity of the image to be restored, Fig. 7 shows the RMSE and SSIM values for images of size 128$\times$128 pixels using the COCO database. We use Total Variation minimization using the Augmented Lagrangian algorithm (TVL) to restore images from ideal and noisy measurements [28]. It is worth mentioning that the TVL algorithm solves the minimization problem with an objective function given by a TV regularization term with norm L1 subject to a system of linear equations.
Fig. 7. Comparison of the GCS+S, CC methods and the Walsh-Paley and Sequency orderings: (a)-(b) plots of mean RMSE and SSIM values with vertical error bars of equal length for 75 images of size 128$\times$128 pixels randomly chosen from the COCO database. (c)-(d) plots the mean RMSE and SSIM values as in (a)-(b) under corrupted measurements by white Gaussian noise with SNR=50dB. The horizontal dashed lines in the plots serve as reference values at 0.1 and 0.6 for the RMSE and SSIM curves, respectively.
Figures 7(a) and 7(b) show RMSE and SSIM values under ideal measurements using the CC, GCS+S, Walsh-Paley and Sequency orderings. Notice that Fig. 7(a) is similar to Fig. 5(bottom-left); nevertheless, in the former, the RMSE values do not converge to 0 at 100% sampling rate, and the standard deviation does not decrease as the sampling rate increases. On the other hand, the SSIM values as shown in Fig. 7(b) improve for all cases compared to Fig. 5(bottom-right). Similarly, the noisy measurements corrupted with white Gaussian noise with SNR=50dB, shown in Fig. 7(c), have a better convergence of the RMSE values compared to Fig. 6 (bottom-left). We also notice an improvement of the SSIM values shown in Fig. 7(d) compared to Fig. 6(bottom-right). It is important to remember that the TVL algorithm was initially introduced to image denoising problems and restoration [28,36]. In particular, the TVL algorithm takes advantage of the sparsity of the image, whose regularization term favors the reconstructed images to be sharper and preserve the edges or boundaries more accurately.
The simulation results in this section show that the GCS+S ordering method has a similar performance as the CC method when restoring images of size 64$\times$ 64 and 128$\times$128 pixels. The previous results were based on a statistical analysis of a set of randomly chosen images from the COCO database. Therefore, evaluation metrics, such as SSIM or RMSE values vary depending on the selected image.
3.3 Experimental results
This section shows results using experimental data of two objects of interest. The structured detection configuration is illustrated in Fig. 1. A white LED diode coupled to a multimode fiber was used as the illuminating source and a camera lens was used to image the object onto the plane of a DMD (VIALUX, Mod V – 9501, resolution 1920x1080 pixels). A collecting lens sends the light to a photodiode (Thorlabs, DET36A Si Biased Detector) placed at its focal distance. The signals were digitalized using National Instruments converter (NI USB–6003) connected to a computer. The generated patterns were loaded to the DMD and displayed at rate of 1K Hz. All experiments were performed under ambient illumination.
Figure 8 shows the restored 128x128 pixel resolution images for an USAF 1951 test chart pattern printed on acetate. As the figure shows, the CC and GCS+S methods obtain sharper images starting from a sampling rate of 10%. For the Sequency and Walsh-Paley Hadamard orderings, you can see blurred images with sampling rates of 10% and 20%. In particular, in the natural ordering appears the typical periodic silhouette ghosting due to the random structure of the matrix. The results in this experiment are consistent with the simulations in Section 3.2.
Fig. 8. Comparison results of image restoration using the USAF 1951 resolution test chart. The images were restored at a resolution of 128$\times$128 pixels using the TVL algorithm for sampling rates from 10% to 50%.
A second experiment is shown in Fig. 9 for INAOE’s institutional logo printed on acetate. Here, we observe similar results for the CC and GCS+S methods in which sharper images are obtained in both cases. Notice that for this specific object, the restored images with the GCS+S method are visually better, in terms of definition, than obtained by the CC method. On the other hand, the natural, Sequency, and Walsh-Paley Hadamard orderings behave similarly to those shown in Fig. 8. The experimental results in Fig. 8 and Fig. 9 demonstrate that the GCS+S method performs similarly to CC and, in specific cases, outperforms the latter.
Fig. 9. Comparison results of image restoration using INAOE’s institutional logo. The images were restored at a resolution of 128$\times$128 pixels using the TVL algorithm for sampling rates from 10% to 50%.
Table 1 shows a quantitative comparison using the RMSE and SSIM metrics for the USAF 1951 test chart pattern and the institutional logo. For each method, the RMSE and SSIM were calculated by taking as ground truth the restored image at a sampling rate of 100%. As this table shows, most of the RMSE values for the GCS+S and CC are less than the lower bound of 0.1, and the SSIM values are greater than the upper bound of 0.6. We observe an improvement of the SSIM values of the GCS+S method for the Logo compared to the CC method. The latter coincides with the experimental results shown in Fig. 9.
Table 1. Quantitative comparison for the USAF test chart pattern and INAOE’s institutional logo.
3.4 Algorithmic complexity of Hadamard sorting methods
The construction of a Hadamard matrix of order $2^n$ as indicated in Eq. (3) can be done in linear time, assuming there is enough RAM in our computer to store it. However, if you want a specific order of their rows, then the complexity of construction can become intractable. Here, we show that for an image of size $N\times N$, the GCS+S method can construct the Hadamard sorted matrix in a complexity of $O(N^2)$, in contrasts to the CC method whose complexity is of $O(N^4)$ where $N=2^n$. It is important to remember that the construction of the Hadamard matrix is done prior to the process of measuring and restoring the image, and it is stored in the main memory of the system. Algorithm 1 shows the steps of the GCS+S method that takes dimension $n$ as an input and returns a Hadamard matrix $H_{2^{2n}}$ used to restore images of size $2^n\times 2^n$ pixels. To know the complexity of Algorithm 1, we must determine which of its steps has the greatest complexity in terms of the main parameter $n$; this type of analysis corresponds to the worst case scenario, also known as the asymptotic big ‘$O$’ notation [37]. In step 1, the Hadamard basis $V_{[1]}^{(n)}$ can be constructed in $\log _2 N$ iterations by following Eqs. ;(9) and (10) since the number of vectors grows exponentially when concatenating lower-dimensional vectors. By definition, the set $V_{[1]}^{(n)}$ has $N$ vectors with $N$ entries each; therefore step 2 can be performed in $N+N\log N$ iterations which is the summation of the complexities of computing the sign change of each vector and ordering a vector of size $N$. In step 3, the construction of the two-dimensional basis set $V_{[1]}^{(n,n)}$ can be performed in $N^2$ iterations by following Eq. (11). On the other hand, the zigzag traversal can be performed in a maximum of $N^2$ iterations, so steps 4 and 5 have a complexity of $2N^2$. Consequently, the largest number of iterations of steps 1 to 5 is $N^2$ so the GCS+S algorithm has a complexity of $O(N^2)$.
Algorithm 1. The GCS+S algorithm
Algorithm 2 and Algorithm 3 show the steps of the CC method, which has the same input and output as Algorithm 1. Step 1 only requires a constant time to perform an assignment, and step 2 can be done in $\log _2(N)$ iterations by following Eq. (3). In steps 3 through 19, each row of the Hadamard matrix is reshaped into a pattern of size $2^n\times 2^n$ and the number of blocks it contains, either black or white, is counted. Here, a black (white) block is a region of adjacent cells in a pattern with values equal to -1 (1). The process of counting the number of blocks in a pattern is equivalent to a flood-fill algorithm [38] used to determine a bounded area connected to a starting pixel with a given attribute in a multi-dimensional array. Steps 5 through 17 describe a modified version of the flood-fill algorithm: it selects an initial pixel and inserts its adjacent pixels into a stack data structure, then iterates until the stack is not empty, repeating the same process with the stack pixels. Algorithm 3 updates the stack by inserting the adjacent pixels with the same attribute or color. It can be seen that the maximum number of iterations of the flood-fill algorithm is $N^2$, and the total number of iterations of all rows of the Hadamard matrix is $N^4$. Finally, in steps 20 and 21, the one-dimensional array $T$ containing the number of blocks per pattern is sorted, and a Hadamard matrix is constructed according to the new indexing of the array $T$, respectively. These last two steps can be performed at most in $N^2+N^2\log N^2$ iterations. Therefore, it can be seen that the modified flood-fill has the largest number of iterations, so the Cake-Cutting algorithm has a complexity of $O(N^4)$.
Algorithm 2. The Cake-Cutting algorithm
Algorithm 3. Update-adjoin-list
We have shown that the GCS+S method has better computational complexity than the CC method. Figure 10 shows a comparison of the experimental running time of Algorithm 1 and Algorithm 2. Figure 10 shows the runtime for image sizes of $2^n\times 2^n$ with $n=5,6,7$. We were unable to run any more simulations because our system does not have the memory capacity needed to store the Hadamard matrix for $n>7$. Note that for fixed $n$, the complexity of the GCS+S and CC methods is polynomial in $N$; but in general, the complexity is exponential in $n$. Therefore, a major problem in Hadamard-based SPI is the design of efficient sorting methods capable of optimally managing the main or temporary memory of the system. The HPE method mentioned in Section 1 has taken a step forward for memory management in constructing a Hadamard basis applied to multi-resolution imaging. Although bases built in [24] are used to reconstruct several images of different resolutions, the HPE method generates a basis like other methods such as RD, CC, and GCS+S. In addition, it is described how to use the HPE method to build the RD ordering. However, it is not possible to compare the execution time of the HPE method with the GCS+S since we do not have the source codes to carry out the simulations, and there is also no analysis of its complexity.
Fig. 10. Experimental running time of sorting methods GCS+S and CC shown in Algorithm 1 and Algorithm 2, respectively. The vertical axis gives the execution time on a logarithmic scale and the horizontal axis gives the input dimension $2^n$ for $n=5,6,7$.
4. Conclusions
We have proposed an efficient Hadamard basis sorting method for SPI based on a generalized two-dimensional orthogonal basis that can be traversed in a zigzag fashion to order the basis vectors from highest to lowest importance. In addition, the proposed GCS+S method has a better computational complexity of $O(N^2)$ than the well-known CC method, whose complexity is $O(N^4)$ where $N=2^n$. The latter is particularly important since the complexity of building and sorting the Hadamard basis grows exponentially with the dimension $n$. Therefore, future research will optimize memory consumption to perform image reconstruction using SPI with limited hardware resources as proposed in [24]. Furthermore, as explained in Section 2.2, we found that if the one-dimensional basis vectors $V_{[1]}^{(n/2)}$ form a GCS, together with a zigzag traversal of the derived two-dimensional basis vectors, then an efficient sorting method is obtained. In particular, the chosen GCS coincides with the traditional order of the Hadamard matrix based on the sign change or Sequency order. Therefore, we can pose the following research problem: find GCSs that improve the performance of the GCS+S method. This problem can be formulated by assigning a cost to each GCS and selecting the GCSs with a minimum cost according to a given optimality criterion. We note that this problem is difficult to solve since the number of GCSs grows exponentially with the dimension of the orthogonal basis [31].
Acknowledgments
L. López-García and A. García–Arellano express their gratitude to the National System of Research.
Disclosures
The authors declare no conflicts of interest.
Data Availability
The Data underlying the results presented in this paper are not publicly available at this time, but can be obtained from the authors at reasonable request.
References
1. G. M. Gibson, S. D. Johnson, and M. J. Padgett, “Single-pixel imaging 12 years on: a review,” Opt. Express 28(19), 28190–28208 (2020). [CrossRef]
2. T. Lu, Z. Qiu, Z. Zhang, and J. Zhong, “Comprehensive comparison of single-pixel imaging methods,” Opt. Lasers Eng. 134, 106301 (2020). [CrossRef]
3. Z. Li, J. Suo, X. Hu, C. Deng, J. Fan, and Q. Dai, “Efficient single-pixel multispectral imaging via non-mechanical spatio-spectral modulation,” Sci. Rep. 7(1), 41435 (2017). [CrossRef]
4. F. Rousset, N. Ducros, F. Peyrin, G. Valentini, C. D’Andrea, and A. Farina, “Time-resolved multispectral imaging based on an adaptive single-pixel camera,” Opt. Express 26(8), 10550–10558 (2018). [CrossRef]
5. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,” Proc. Natl. Acad. Sci. 109(26), E1679–E1687 (2012). [CrossRef]
6. M.-J. Sun, M. P. Edgar, G. M. Gibson, B. Sun, N. Radwell, R. Lamb, and M. J. Padgett, “Single-pixel three-dimensional imaging with time-based depth resolution,” Nat. Commun. 7(1), 12010–1723 (2016). [CrossRef]
7. Z. Zhang, S. Jiao, M. Yao, X. Li, and J. Zhong, “Secured single-pixel broadcast imaging,” Opt. Express 26(11), 14578–14591 (2018). [CrossRef]
8. S. Jiao, J. Feng, Y. Gao, T. Lei, and X. Yuan, “Visual cryptography in single-pixel imaging,” Opt. Express 28(5), 7301–7313 (2020). [CrossRef]
9. D. V. Strekalov, A. V. Sergienko, D. N. Klyshko, and Y. H. Shih, “Observation of two-photon “ghost” interference and diffraction,” Phys. Rev. Lett. 74(18), 3600–3603 (1995). [CrossRef]
10. S. Walborn, C. Monken, S. Pádua, and P. Souto Ribeiro, “Spatial correlations in parametric down-conversion,” Phys. Rep. 495(4-5), 87–139 (2010). [CrossRef]
11. D. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006).
12. E. J. Candes and T. Tao, “Near-optimal signal recovery from random projections: Universal encoding strategies?” IEEE Trans. Inf. Theory 52(12), 5406–5425 (2006). [CrossRef]
13. R. G. Baraniuk, “Compressive sensing [lecture notes],” IEEE Signal Process. Mag. 24(4), 118–121 (2007). [CrossRef]
14. Z. Zhang, X. Wang, G. Zheng, and J. Zhong, “Hadamard single-pixel imaging versus fourier single-pixel imaging,” Opt. Express 25(16), 19619–19639 (2017). [CrossRef]
15. Z. Zhang, X. Ma, and J. Zhong, “Single-pixel imaging by means of Fourier spectrum acquisition,” Nat. Commun. 6(1), 6225 (2015). [CrossRef]
16. Z. Zhang, X. Wang, and J. Zhong, “Fast Fourier single-pixel imaging using binary illumination,” Sci. Rep. 7(1), 12029 (2017). [CrossRef]
17. F. Rousset, N. Ducros, A. Farina, G. Valentini, C. D’Andrea, and F. Peyrin, “Adaptive basis scan by wavelet prediction for single-pixel imaging,” IEEE Trans. on Comput. Imaging 3(1), 36–46 (2017). [CrossRef]
18. M. Edgar, G. Gibson, and M. Padgett, “Principles and prospects for single-pixel imaging,” Nat. Photonics 13(1), 13–20 (2019). [CrossRef]
19. H. Sarukhanyan, K. Egiazarian, and J. Astola, Hadamard transforms (SPIE, 2011).
20. M.-J. Sun, M. Tong, M. Edgar, M. Padgett, and N. Radwell, “A Russian Dolls ordering of the Hadamard basis for compressive single-pixel imaging,” Sci. Rep. 7(1), 3464 (2017). [CrossRef]
21. W.-K. Yu, “Super sub-Nyquist single-pixel imaging by means of Cake-Cutting Hadamard basis sort,” Sensors 19(19), 4122 (2019). [CrossRef]
22. W.-K. Yu and Y.-M. Liu, “Single-pixel imaging with origami pattern construction,” Sensors 19(23), 5135 (2019). [CrossRef]
23. X. Yu, F. Yang, B. Gao, J. Ran, and X. Huang, “Deep compressive single pixel imaging by reordering Hadamard basis: A comparative study,” IEEE Access 8, 55773–55784 (2020). [CrossRef]
24. C. Zhou, X. Zhao, H. Huang, G. Wang, X. Wang, L. Song, and K. Xue, “Multi-resolution single-pixel imaging via hadamard pipeline coding,” Appl. Phys. B 126(10), 163 (2020). [CrossRef]
25. H. Ma, A. Sang, C. Zhou, X. An, and L. Song, “A zigzag scanning ordering of four-dimensional walsh basis for single-pixel imaging,” Opt. Commun. 443, 69–75 (2019). [CrossRef]
26. C. Watts, D. Shrekenhamer, J. Montoya, G. Lipworth, J. Hunt, T. Sleasman, D. Smith, and W. Padilla, “Terahertz compressive imaging with metamaterial spatial light modulators,” Nat. Photonics 8(8), 605–609 (2014). [CrossRef]
27. L. Bian, J. Suo, Q. Dai, and F. Chen, “Experimental comparison of single-pixel imaging algorithms,” J. Opt. Soc. Am. A 35(1), 78–87 (2018). [CrossRef]
28. W. Yin, H. Jiang, and Y. Zhang, “An efficient augmented lagrangian method with applications to total variation minimization,” Comput. Optim. Appl. 56(3), 507–530 (2013). [CrossRef]
29. S. Becker, J. Bobin, and E. Candés, “Nesta: A fast and accurate first-order method for sparse recovery,” SIAM J. Imaging Sci. 4(1), 1–39 (2011). [CrossRef]
30. T. Mizuno and T. Iwata, “Hadamard-transform fluorescence-lifetime imaging,” Opt. Express 24(8), 8202–8213 (2016). [CrossRef]
31. G. Ben-Artzi, H. Hel-Or, and Y. Hel-Or, “The gray-code filter kernels,” IEEE Trans. Pattern Anal. Mach. Intell. 29(3), 382–393 (2007). [CrossRef]
32. Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. on Image Process. 13(4), 600–612 (2004). [CrossRef]
33. “The USC-SIPI Image Database,” https://sipi.usc.edu/database/. Accessed: 2021-09-30.
34. T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick, “Microsoft coco: Common objects in context,” in Computer Vision – ECCV 2014, D. Fleet, T. Pajdla, B. Schiele, and T. Tuytelaars, eds. (Springer International Publishing, Cham, 2014), pp. 740–755.
35. M. Ma, Q. Sun, X. Gao, G. Wang, H. Deng, Y. Zhang, Q. Guan, and X. Zhong, “High-efficiency single-pixel imaging using discrete hartley transform,” AIP Adv. 11(7), 075211 (2021). [CrossRef]
36. M.-J. Sun, M. P. Edgar, D. B. Phillips, G. M. Gibson, and M. J. Padgett, “Improving the signal-to-noise ratio of single-pixel imaging using digital microscanning,” Opt. Express 24(10), 10476–10485 (2016). [CrossRef]
37. T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, Introduction to Algorithms, 3rd Edition (MIT, 2009).
38. J. D. Foley, A. van Dam, S. K. Feiner, and J. F. Hughes, Computer Graphics: Principles and Practice, 2nd Ed. (Addison-Wesley, 1990).
