Improving the image reconstruction quality of compressed ultrafast photography via an augmented Lagrangian algorithm

Capturing ultrafast dynamic scenes at high imaging speed has been a long-term dream of scientists, as it can enable the discovery of new physical phenomena and the development of new optical imaging technologies. It has been shown that charge-coupled device (CCD) or complementary metal oxide semiconductor (CMOS)-based imaging techniques provide well-established tools to measure dynamic scenes. Using this imaging method, the dynamic scene at each moment is recorded in turn, and thus the acquisition rate is limited by the on-chip storage and electronic readout speed. So far, the maximal frame rate of a CCD or CMOS is about in the order of 10⁷ frames per second (fps) [1]. It is a great challenge to further increase the imaging speed under the current technical conditions. Recently, a compressed ultrafast photography (CUP) technique employed a new imaging strategy to break this technical limitation [2–6], which recorded the whole dynamic scene and then recovered it by a computational imaging method. Now, the imaging speed of CUP can be up to 10¹³ fps; it can increase by six orders compared to that of a CCD or CMOS. The number of reconstructed images depends strongly on the mathematical processing algorithms and the temporal sweeping velocity. So far, CUP has been successfully applied to measuring some typical transient optical events, such as laser pulse reflection and refraction [2], photon racing in two media [2], photonic Mach cones [3], and spatially modulated pulsed laser spots [4].

In CUP, the dynamical scene was encoded with random codes and decoded using a two-step iterative shrinkage/thresholding (TwIST) algorithm [7]. However, the selection of the penalty parameter in the TwIST algorithm will significantly affect the image reconstruction quality, which is unfavorable for determining the unknown dynamic scenes. To overcome the drawback of the TwIST algorithm, here we report an augmented Lagrangian (AL) algorithm to reconstruct the ultrafast dynamic scenes for CUP. To verify the advantages of the AL algorithm, we numerically simulate a flashing Shepp–Logan (S–L) phantom and a moving circular spot in the space, and experimentally measure a collimated femtosecond laser obliquely illuminating a stripe pattern and the temporal evolution of a spatially modulated picosecond laser spot. Both the simulation and experimental results show that the AL algorithm is almost unaffected by the selection of the penalty parameter, and can obtain higher image reconstruction quality compared to the TwIST algorithm.

In CUP, as shown in figure 1, the ultrafast dynamic scene I(x,y,t) is measured by a CCD after a spatial encoding operator, C, a temporal shearing operator, S, and a spatiotemporal integration operator, T. Thus, the measured two-dimensional (2D) image E(x,y) can be mathematically formulated as [2]

$$\begin{eqnarray}&&E\left(x,y\right)=TSCI\left(x,y,t\right).\end{eqnarray} \tag{ 1 }$$

For convenience, we define O = TSC, thus equation (1) can be further written as

$$\begin{eqnarray}&&E\left(x,y\right)=OI\left(x,y,t\right).\end{eqnarray} \tag{ 2 }$$

To reconstruct the original ultrafast dynamic scene I(x,y,t), it needs to inversely solve equation (2). A common method is to employ a convex total variation (TV) model [8–11], which can be expressed as

$$\begin{eqnarray}&&\left\{\begin{array}{l}\mathop{Min}\limits_{I(x,y,t)}{\rm{\Phi }}\left(I\left(x,y,t\right)\right)\\ s.t.E\left(x,y\right) \mbox{-} OI\left(x,y,t\right)=0\end{array}\right.,\end{eqnarray} \tag{ 3 }$$

where Φ(I(x, y, t)) is the form of TV for I(x, y , t). Usually, the two constrained problems in equation (3) can be optimized by a Lagrangian function that is defined as [12]

$$\begin{eqnarray}&&L\left(I\left(x,y,t\right),\lambda \right)={\rm{\Phi }}\left(I\left(x,y,t\right)\right) \mbox{-} \lambda \left(E\left(x,y\right) \mbox{-} OI\left(x,y,t\right)\right),\end{eqnarray} \tag{ 4 }$$

where λ is the Lagrange multiplier matrix. Here, E(x,y)-OI(x, y, t) is a matrix, but it will be arranged as a vector in operation, and the same operation is performed for λ. Thus, λ(E(x, y)-OI(x,y,t)) is a scalar. Only when λ approaches the unique right matrix due to the linear independence constraint qualification (LICQ), I(x, y, t) is the solution of equation (3) by minimizing L(I(x, y, t), λ) [12]. Once L(I(x, y, t),ζ) achieves the minimal value, the derivative of P(I(x,y,t),ζ) to I(x,y,t) should be zero, which can be written as

$$\begin{eqnarray}&&{\rm{\nabla }}L\left(I\left(x,y,t\right),\lambda \right)={\rm{\nabla }}{\rm{\Phi }}\left(I\left(x,y,t\right)\right)+\lambda O=0.\end{eqnarray} \tag{ 5 }$$

To deal with the big data problem in equation (4), here a computational method of quadratic penalty function is used and written as [12]

$$\begin{eqnarray}&&P\left(I\left(x,y,t\right),\zeta \right)={\rm{\Phi }}\left(I\left(x,y,t\right)\right)+\displaystyle \frac{\zeta }{2}{\parallel E\left(x,y\right) \mbox{-} OI\left(x,y,t\right)\parallel }_{2}^{2},\end{eqnarray} \tag{ 6 }$$

where ∣∣·∣∣ is a l₂ norm, which is utilized for vectors, and ζ is the penalty parameter with ζ > 0. In the same way, the derivative of P(I(x,y,t),ζ) to I(x,y,t) is also zero, and it can be formulated as

$$\begin{eqnarray}&&\,\begin{array}{r}\begin{array}{l}{\rm{\nabla }}P\left(I\left(x,y,t\right),\zeta \right)={\rm{\nabla }}{\rm{\Phi }}\left(I\left(x,y,t\right)\right)\\ \,-\zeta \left(E\left(x,y\right)-OI\left(x,y,t\right)\right)O=0.\end{array}\end{array}\,\end{eqnarray} \tag{ 7 }$$

In previous studies, the image reconstruction method for CUP was based on the quadratic penalty function, such as TwIST [2–4] and the fast iterative shrinkage/thresholding algorithm (FISTA) [5]. Comparing equations (7) and (5), the Lagrange multiplier λ can be deduced as

$$\begin{eqnarray}&&\lambda =-\zeta \left(E\left(x,y\right)-OI\left(x,y,t\right)\right).\end{eqnarray} \tag{ 8 }$$

By rearranging equation (8), we have

$$\begin{eqnarray}&&E\left(x,y\right)-OI\left(x,y,t\right)=-\displaystyle \frac{1}{\zeta }\lambda .\end{eqnarray} \tag{ 9 }$$

We can use the value of E(x, y)-OI(x, y, t) to illustrate the feasibility of the image reconstruction, and the smaller value corresponds to the greater feasibility. It is easy to see from equation (9) that the small value ζ will bring great infeasibility. Increasing ζ can improve the feasibility, but it will cause ill conditions for the quadratic penalty function [12]. Therefore, the selection of ζ will greatly affect the image reconstruction quality. The TwIST algorithm in previous studies was faced with the problem of image reconstruction instability [2–4]. To overcome this problem, here we introduce an AL function that is written as [10, 13–15]

$$\begin{eqnarray}&&\begin{array}{l}{L}_{a}\left(I\left(x,y,t\right),\gamma ,\zeta \right)={\rm{\Phi }}\left(\left(I\left(x,y,t\right)\right)-\gamma \left(E\left(x,y\right)\right.\right.\\ \,-\left.OI\left(x,y,t\right)\right)+\displaystyle \frac{\zeta }{2}{\parallel E\left(x,y\right)-OI\left(x,y,t\right)\parallel }_{2}^{2},\end{array}\end{eqnarray} \tag{ 10 }$$

where γ is a variable Lagrange multiplier matrix. Similarly, the derivative of L_a(I(x, y, t), γ, ζ) to I(x, y, t) will be zero when L_a(I(x, y, t), γ, ζ) approaches the minimal value, and it is given by

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\nabla }}{L}_{a}\left(I\left(x,y,t\right),\gamma ,\zeta \right)={\rm{\nabla }}{\rm{\Phi }}\left(I\left(x,y,t\right)\right)\\ \,+\left(\gamma -\zeta \left(E\left(x,y\right)-OI\left(x,y,t\right)\right)\right)O=0.\end{array}\end{eqnarray} \tag{ 11 }$$

By comparing equations (11) and (5), the Lagrange multiplier λ can be calculated as [12]

$$\begin{eqnarray}&&\lambda =\gamma -\zeta \left(E\left(x,y\right)-OI\left(x,y,t\right)\right).\end{eqnarray} \tag{ 12 }$$

Equation (12) can be further formulated as

$$\begin{eqnarray}&&E\left(x,y\right)-OI\left(x,y,t\right)=-\displaystyle \frac{1}{\zeta }(\lambda -\gamma ).\end{eqnarray} \tag{ 13 }$$

One can see from equation (13) that the value of E(x, y)-OI(x, y, t) depends on ζ and γ. By optimizing γ in each iteration, it is easy to minimize E(x, y)-OI(x, y, t). Therefore, the AL-based image reconstruction method can be independent of the selection of ζ. Moreover, simultaneously controlling the two parameters ζ and γ can make the term E(x,y)-OI(x,y,t) easier to converge, and accordingly improves the image reconstruction accuracy. Obviously, compared to the previous TwIST algorithm, the proposed AL algorithm is less dependent on the selection of ζ and can obtain higher image reconstruction quality.

Zoom In Zoom Out Reset image size

In the image reconstruction, we initially set ΔLa⁰ = 1, γ⁰ = 0, tolerance ρ (0 < ρ < 10⁻³), penalty parameter ζ (ζ > 0), and start the point I⁰. Here, we adopt the AL algorithm framework as follows.

While ${\rm{\Delta }}{{L}_{a}}^{i-1}$ > ρ do

Find the minimizer Ii of ${{L}_{a}}^{i}$ by minimizing equation (10) from Ii-1;

Compute ΔLai = $\left({\parallel {{L}_{a}}^{i}\parallel }_{2}-{\left|\left|{{L}_{a}}^{i-1}\right|\right|}_{2}\right)/{\parallel {{L}_{a}}^{i-1}\parallel }_{2};$

Update γi by γi = γi-1-ζ(E(x,y)-OI(x,y,t));

End do

In each iteration, we minimize L_a(I(x, y, t), γ, ζ) in equation (10) by utilizing a TVAL3 [13]. Moreover, it should be noted that γⁱ is a variable matrix and will be automatically updated.

We first numerically simulate a simple dynamic scene with a flashing 200-by-200 S–L phantom. In the dynamic scene, the intensity of the measured target varies, but the position remains unchanged. The ultrafast dynamic scene consists of six frames, and the odd frame is bright while the even one is dark. Here, we utilize a normalized correlation coefficient (CC) value to characterize the image reconstruction quality. In mathematics, the CC is used to describe the image similarity, which can be written as [16]

$$\begin{eqnarray}&&CC=\displaystyle \frac{\displaystyle {\sum }_{n=1}^{N}\left({a}_{n}-{\bar{a}}_{n}\right)\left({e}_{n}-{\bar{e}}_{n}\right)}{\sqrt{\displaystyle {\sum }_{n=1}^{N}{\left({a}_{n}-{\bar{a}}_{n}\right)}^{2}}\sqrt{\displaystyle {\sum }_{n=1}^{N}{\left({e}_{n}-{\bar{e}}_{n}\right)}^{2}}}\times 100 \% \end{eqnarray} \tag{ 14 }$$

where N is the number of total pixels in the image, and a_n and ${\unicode{x00101}}_{n}$ (or e_n and ${\unicode{x00113}}_{n}$) are a pixel value and mean of the reconstructed image (or base image), respectively. In general, the larger the CC value, the higher the image reconstruction quality. When the flashing S–L phantom is well reconstructed, the CC value for the bright frame is large while that for the dark frame is small. Figure 2(a) shows the base S–L phantom image for the bright (up panel) and dark (low panel) frames. Figures 2(f) and (g) present the CC values with the increase in the penalty parameter ζ by the AL (black line) and TwIST algorithms (red line) for the bright and dark frames, respectively. With the increase in ζ, the CC value for the bright (or dark) frame shows an increase (or decrease) followed by a decrease (or increase) process for both algorithms. However, the ζ value significantly affects the image reconstruction quality for the TwIST algorithm, but has little effect for the AL algorithm. Moreover, the maximal (or minimal) CC value for the bright (or dark) frame is 94.17% (or 11.67%) for the AL algorithm, while it is 94.08% (or 24.17%) for the TwIST algorithm. Obviously, compared to the TwIST algorithm, the AL algorithm can obtain higher image reconstruction quality. In addition, the CC value for the TwIST algorithm will experience a mutation phenomenon, but will not for the AL algorithm. For further comparison of the image reconstruction results by the two algorithms, we give the best and worst reconstructed images for the bright and dark frames by the AL and TwIST algorithms, as shown in figures 2(b)–(e). By comparing the two algorithms, the AL algorithm can obtain clearer details of the image information, especially in the region labelled with arrows. Furthermore, the dark frame also carries less information from the bright frame.

Zoom In Zoom Out Reset image size

To further illustrate that, compared to the TwIST algorithm, the AL algorithm has the greater advantage in the image reconstruction for CUP, we simulate another type of ultrafast dynamic scene: that is, a circular spot with the Gaussian intensity distribution moving from left to right. Here, the position of the measured target moves, but the intensity remains unchanged, which is the opposite of the dynamic scene above. In our simulation, the ultrafast dynamic scene contains 20 base images, as shown in figure 3(a). Figures 3(b) and (c) show the reconstructed images by the AL and TwIST algorithms, respectively. As shown in figure 3(b), these reconstructed circular spots by the AL algorithm keep the original shape well. However, those reconstructed by the TwIST algorithm have an obvious deformation, especially in the front and back of a few of images, as shown in figure 3(c). That is to say, the TwIST algorithm will affect the spatial structure of the reconstructed ultrafast dynamic scene, but the AL algorithm will not. Figures 3(d) and (e) present the tenth reconstructed images with ζ from 1 to 5 by the AL and TwIST algorithms, respectively. As expected, the reconstructed circular spots by the AL algorithm remain the same shape, which is independent of ζ, but those by the TwIST algorithm are deformed with the smaller value of ζ. Thus, if the ultrafast dynamic scene is unknown in advance, the reconstructed images cannot be correctly determined by the TwIST algorithm. It is obvious that the AL algorithm can solve the problem of image reconstruction instability in CUP.

Zoom In Zoom Out Reset image size

As shown above, we have proved from theoretical formulation and numerical simulations that the AL algorithm has more advantages in the image reconstruction for CUP than the TwIST algorithm. Next, we further compare the two image reconstruction methods in experiments. We measure such an ultrafast dynamic scene with a collimated femtosecond laser pulse obliquely illuminating a stripe pattern, and the experimental arrangement is shown in figure 4(a). A Ti:sapphire amplifier (Spectra-physics, Spitfire Ace-35F) is used to generate the femtosecond laser pulse with a pulse width of about 50 fs, a central wavelength of 800 nm, and repetition rate of 10 Hz. The output femtosecond laser pulse obliquely illuminates a transverse stripe pattern at an angle of ∼38° with respect to the surface normal. The wavefront movement on the stripe pattern is imaged via a camera lens and a 4f imaging system, and then a digital micromirror device (DMD) (Texas Instruments, DLP LightCrafter) encodes the image. Finally, the encoded dynamical scene is sent to a streak camera (Hamamatsu, C7700) for measurement by a beam splitter. In our experiment, the spatial resolution of the streak camera is 18 lp/mm, and the temporal resolution is 4 ps (i.e. the imaging speed of 2.5 × 10¹¹ fps). Figures 4(b) and (c) show the best reconstruction images by the AL and TwIST algorithms, respectively. As can be seen, both algorithms can recover the original dynamical scene of the wavefront movement, but the AL algorithm can better reflect the true situation. Additionally, compared to the TwIST algorithm, less background noise is involved for the AL algorithm.

Zoom In Zoom Out Reset image size

What is more, we also perform another experiment to further show the advantage of the AL-based image reconstruction method. Here, we measure the temporal evolution of a spatially modulated picosecond laser spot, and the experimental design is shown in figure 5(a). The output 50 fs laser pulse from the Ti:sapphire amplifier is broadened to 20 ps by a stretcher, and a thin wire is used to divide the laser spot into two components in space. The spatially modulated laser spot illuminates on a thin white paper, and a small fraction of photons can pass through the thin white paper. Thus, the temporal evolution of the laser spot can be measured by our CUP system. The best reconstruction images by the AL and TwIST algorithms are shown in figures 5(b) and (c), respectively. As expected, the AL algorithm can clearly distinguish the blocked part in the center of the laser spot, but the TwIST algorithm does not, especially around the time zero. Moreover, similar to the experimental result in figure 4, more background noise exists for the TwIST algorithm.

Zoom In Zoom Out Reset image size

In summary, we have developed an AL algorithm for CUP to reconstruct the ultrafast dynamic scenes. Our theoretical and experimental results showed that, compared to the TwIST algorithm, the AL algorithm is less dependent on the selection of the penalty parameter in the image reconstruction, and can obtain higher image reconstruction quality. The AL algorithm overcomes the limitation of image reconstruction instability demonstrated using the TwIST algorithm in previous studies, which is very meaningful for correctly identifying the unknown ultrafast dynamic scenes. Moreover, this study also presents a clear physical and mathematical insight for the advantages of the AL algorithm in the image reconstruction of CUP, which may further promote the future applications of the AL algorithm in the computational imaging area. In addition, considering the standard compressed sensing model in equation (3), some other algorithms can also be used for the image reconstruction of CUP, such as gradient projection for sparse reconstruction (GPSR) [17]. In future studies, an important aim for CUP is to optimize the algorithm to further improve the image reconstruction accuracy.

This work was partly supported by the National Natural Science Foundation of China (Nos. 91850202, 11774094, 11727810, 11804097 and 61720106009) and the Science and Technology Commission of Shanghai Municipality (Nos. 17ZR146900 and 16520721200).