Article Contents
Article Contents

# Edge-adaptive $\ell_2$ regularization image reconstruction from non-uniform Fourier data

Rick Archibald's work is sponsored by the Applied Mathematics Division of ASCR, DOE; in particular under the ACUMEN project (RA). Anne Gelb's work is supported in part by the grants NSF-DMS 1502640, NSF-DMS 1732434, AFOSR FA9550-18-1-0316 and AFOSR FA9550-15-1-0152

• Signals and images recovered from edge-sparsity based reconstruction methods may not truely be sparse in the edge domain, and often result in poor quality reconstruction. Iteratively reweighted methods provide some improvement in accuracy, but at the cost of extended runtime. This paper examines such methods when data are acquired as non-uniform Fourier samples, and then presents a new non-iterative weighted regularization method that first pre-processes the data to determine the precise locations of the non-zero values in the edge domain. Our new method is both accurate and efficient, and outperforms reweighted regularization methods in several numerical experiments.

Mathematics Subject Classification: Primary: 68U10, 65F22; Secondary: 42A10.

 Citation:

• Figure 1.  Non-uniform sampling ${ λ}_{\bf{k}}$ as in (5)

Figure 13.  $257\times257$ pixel Shepp-Logan phantom

Figure 2.  $f_1(x)$ (top) and $f_2(x)$ (bottom) reconstructed via Algorithm 1 using $m = 1$ from $257$ Fourier modes on $257$ grid points. The red reconstructions indicate recovery from noise-free data as in (1), while yellow reconstructions are recovered from Fourier coefficients with zero-mean complex Gaussian noise added as in (3). Here we use a signal to noise ratio (SNR) of $20$ dB. For $f_1(x)$, we used parameters $\rho = 1$, $\ell_{max} = 25$, and $\epsilon = 1.9$. For $f_2(x)$, we used parameters $\rho = 1$, $\ell_{max} = 25$, and $\epsilon = 2.9$

Figure 3.  (Left) $257\times257$ pixel function $f_3(x,y)$; (Middle) reconstruction via Algorithm 2 with $m = 2$ from $257\times257$ jittered Fourier coefficients on $257\times257$ grid points; (Right) pointwise error plotted on a logarithmic scale. The algorithm parameters used are $\rho = .01$, $\epsilon = .9$, and $\ell_{max} = 5$

Figure 4.  (Left) Ideal weighting matrix for the $y$-direction edges computed from (19) using the exact $257\times257$ two-dimensional function $f_3(x,y)$; (Right) final weighting matrix for the $y$-direction edges produced by Algorithm 2. The minimum weight of $0$ is indicated by black while maximum weight of $\frac{1}{\epsilon}\approx 1.11$ is indicated by white, while gray indicates a weight in between $0$ and $\frac1\epsilon$

Figure 5.  (Left) $f_1(x)$ and $f_2(x)$; (Middle) Jump function approximations using (37); (Right) Jump function approximations with additive noise starting from (3) with SNR $= 20$ dB. Here we use $257$ reconstruction points, $257$ jittered Fourier modes, and regularization parameter $\mu = 1$

Figure 6.  Edge detection (center and right) in the $x$ and $y$ directions using (38) for $f_3(x,y)$ (left). Here we use $257\times257$ jittered Fourier modes, $257\times257$ reconstruction points, and regularization parameter $\mu = 1$

Figure 17.  (Top) Reconstruction comparison of SAR vehicle data using (left) equation (44) and (right) Algorithm 6. (Bottom) A close up of the lower right tire

Figure 7.  (Left) Comparison of Algorithms 1 and 5 on $f_1(x)$ using PA order $m = 1$ given $257$ jittered Fourier samples reconstructed on $257$ grid points; (Right) corresponding pointwise errors. For parameters, we use $\rho = 1$, $\mu = 1$, $\lambda = 1$, $\epsilon = 1.9$, $\ell_{max} = 25$, and threshold $\tau = 1/257$

Figure 8.  Comparison of Algorithms 2 (left) and 6 (right) for reconstructing $f_3(x,y)$ from $257\times257$ noise-free Fourier modes on $257\times257$ grid points. For parameters, we use PA order $m = 2$ due to the piecewise quadratic nature of the function, $\rho = .01$, $\epsilon = .9$, and $\ell_{max} = 5$, $\mu = .1$, $\tau = .025$, and $\lambda = 1$

Figure 9.  Comparison of errors using Algorithms 2 (left) and 6 (right) for reconstructing $f_3(x,y)$, using the same parameters as Figure 8

Figure 10.  Pointwise error plot comparisons between Algorithm 1 and Algorithm 5 for $f_1(x)$ given $257$ jittered Fourier samples reconstructed on $257$ grid points. (Left) Parameters $m = 1$, $\ell_{max} = 25$, $\epsilon = 1.9$, $\tau = 1/257$, $\lambda = 1$, and vary $\rho = .01,.1,1,10,100$; (Right) Parameters $m = 1$, $\rho = 1$, $\ell_{max} = 25$, $\epsilon = 1.9$, $\mu = 1$, $\tau = 1/257$, and vary $\lambda = .01,.1,1,10,100$

Figure 11.  Comparison of Algorithms 1 and 5 on $f_2(x)$ using (top) $m = 1$ (middle) $m = 2$ and (bottom) $m = 3$ given $257$ jittered Fourier samples reconstructed on $257$ grid points. Left are the reconstructions, right are the corresponding pointwise errors. For parameters, we use $\rho = 1$, $\ell_{max} = 25$, $\epsilon = 2.9$, $\mu = 1$, $\tau = 1/257$, and $\lambda = 1$

Figure 12.  (Left) Comparisons of Algorithms 1 and 5 on $f_1(x)$ and $f_2(x)$ given $257$ noisy jittered Fourier samples. $SNR = 15$ for $f_1(x)$ and SNR $= 20$ for $f_2(x)$. In both cases we use PA order $m = 1$ and reconstruct on $257$ grid points; (Right) corresponding pointwise errors. For $f_1(x)$ we use parameters $\rho = 1$, $\epsilon = 1.9$, $\ell_{max} = 25$, $\mu = 1$, $\tau = 1/257$, and $\lambda = 1$. For $f_2$ we use parameters $\rho = 1$, $\epsilon = 2.9$, $\ell_{max} = 25$, $\mu = 1$, $\tau = 1/257$, and $\lambda = 1$

Figure 14.  (Top) Reconstruction on $257\times 257$ grid points of the Shepp-Logan phantom using Algorithms 2 (left) and 6 (right) with PA order $m = 1$ from $129^2$ Fourier coefficients randomly chosen from a grid of $257\times257$. (Bottom) respective pointwise errors. For parameters, we use $\rho = .01$, $\epsilon = .9$, $\ell_{max} = 5$, $\mu = .01$, $\tau = 0.1$, $\lambda = .1$. The relative error using Algorithm 2 was $RE = .7160$, while Algorithm 6 yielded $RE = .4873$

Figure 15.  (Top) Reconstruction on $257\times 257$ grid points of the Shepp-Logan phantom using Algorithms 2 (left) and 6 (right) with PA order $m = 1$ from $181^2$ Fourier coefficients randomly chosen from a grid of $257\times257$. (Bottom) respective pointwise errors. For parameters, we use $\rho = .01$, $\epsilon = .9$, $\ell_{max} = 5$, $\mu = .01$, $\tau = 0.1$, $\lambda = .1$. The relative error using Algorithm 2 was $RE = .5180$, while Algorithm 6 yielded $RE = 0.3259$

Figure 16.  (Top) Reconstruction on $257\times 257$ grid points of the Shepp-Logan phantom using Algorithms 2 (left) and 6 (right) with PA order $m = 1$ from $225^2$ Fourier coefficients randomly chosen from a grid of $257\times257$. (Bottom) respective pointwise errors. For parameters, we use $\rho = .01$, $\epsilon = .9$, $\ell_{max} = 5$, $\mu = .01$, $\tau = 0.1$, $\lambda = .1$. The relative error using Algorithm 2 was $RE = .3458$, while Algorithm 6 yielded $RE = .2930$

Figure 18.  Comparison of edge-adaptive $\ell_1$ and edge-adaptive $\ell_2$ methods on $f_2(x)$ using 257 Fourier coefficients. (left) error plot from reconstructions with noise-free data. (right) error plot from reconstruction with 20 dB zero-mean Gaussian noise added to the data

Table 1.  Run time comparison between Algorithms 2 and 6 for reconstructing $f_3(x,y)$. We used $\ell_{max} = 5$. The run time includes the time to perform Algorithm 4

 Data size Algorithm 2 Algorithm 6 $129\times129$ 4 mins 10 secs 5.6 secs $257\times257$ 13 mins 2 secs 22 secs $513\times513$ 49 mins 26 secs 1 min 33 secs $1025\times1025$ 3 hours 16 mins 6 mins 28 secs
•  [1] R. Archibald, A. Gelb and R. B. Platte, Image reconstruction from undersampled Fourier data using the polynomial annihilation transform, Journal of Scientific Computing, 67 (2016), 432-452.  doi: 10.1007/s10915-015-0088-2. [2] R. Archibald, A. Gelb and J. Yoon, Polynomial fitting for edge detection in irregularly sampled signals and images, SIAM Journal on Numerical Analysis, 43 (2005), 259-279.  doi: 10.1137/S0036142903435259. [3] M. S. Asif and J. Romberg, Fast and accurate algorithms for re-weighted $\ell_1$-norm minimization, IEEE Transactions on Signal Processing, 61 (2013), 5905-5916.  doi: 10.1109/TSP.2013.2279362. [4] E. J. Candès, J. Romberg and T. Tao, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Transactions on Information Theory, 52 (2006), 489-509.  doi: 10.1109/TIT.2005.862083. [5] E. J. Candes, J. K. Romberg and T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Communications on Pure and Applied Mathematics, 59 (2006), 1207-1223.  doi: 10.1002/cpa.20124. [6] E. J. Candes and T. Tao, Near-optimal signal recovery from random projections: Universal encoding strategies?, IEEE Transactions on Information Theory, 52 (2006), 5406-5425.  doi: 10.1109/TIT.2006.885507. [7] E. J. Candes, M. B. Wakin and S. P. Boyd, Enhancing sparsity by reweighted $\ell_1$ minimization, Journal of Fourier Analysis and Applications, 14 (2008), 877-905.  doi: 10.1007/s00041-008-9045-x. [8] T. Chan, A. Marquina and P. Mulet, High-order total variation-based image restoration, SIAM Journal on Scientific Computing, 22 (2000), 503-516.  doi: 10.1137/S1064827598344169. [9] R. Chartrand and W. Yin, Iteratively reweighted algorithms for compressive sensing, in Acoustics, Speech and Signal Processing, 2008. ICASSP 2008. IEEE International Conference on, IEEE, 2008, 3869-3872. [10] I. Daubechies, R. DeVore, M. Fornasier and C. S. Güntürk, Iteratively reweighted least squares minimization for sparse recovery, Communications on Pure and Applied Mathematics, 63 (2010), 1-38.  doi: 10.1002/cpa.20303. [11] D. L. Donoho, Compressed sensing, IEEE Transactions on Information Theory, 52 (2006), 1289-1306.  doi: 10.1109/TIT.2006.871582. [12] K. E. Dungan, C. Austin, J. Nehrbass and L. C. Potter, Civilian vehicle radar data domes, in Algorithms for Synthetic Aperture Radar Imagery XVII, International Society for Optics and Photonics, 7699 (2010), 76990P. doi: 10.1117/12.850151. [13] A. Dutt and V. Rokhlin, Fast Fourier transforms for nonequispaced data, SIAM Journal on Scientific Computing, 14 (1993), 1368-1393.  doi: 10.1137/0914081. [14] G. Franceschetti and  R. Lanari,  Synthetic Aperture Radar Processing, CRC press, 2018.  doi: 10.1201/9780203737484. [15] A. Gelb and D. Cates, Detection of edges in spectral data iii-refinement of the concentration method, Journal of Scientific Computing, 36 (2008), 1-43.  doi: 10.1007/s10915-007-9170-8. [16] A. Gelb and G. Song, A frame theoretic approach to the nonuniform fast Fourier transform, SIAM Journal on Numerical Analysis, 52 (2014), 1222-1242.  doi: 10.1137/13092160X. [17] A. Gelb and G. Song, Detecting edges from non-uniform Fourier data using Fourier frames, Journal of Scientific Computing, 71 (2017), 737-758.  doi: 10.1007/s10915-016-0320-8. [18] A. Gelb and E. Tadmor, Detection of edges in spectral data, Applied and Computational Harmonic Analysis, 7 (1999), 101-135.  doi: 10.1006/acha.1999.0262. [19] A. Gelb and E. Tadmor, Detection of edges in spectral data ii. nonlinear enhancement, SIAM Journal on Numerical Analysis, 38 (2000), 1389-1408.  doi: 10.1137/S0036142999359153. [20] A. Gelb and E. Tadmor, Adaptive edge detectors for piecewise smooth data based on the minmod limiter, Journal of Scientific Computing, 28 (2006), 279-306.  doi: 10.1007/s10915-006-9088-6. [21] T. Goldstein and S. Osher, The split bregman method for l1-regularized problems, SIAM Journal on Imaging Sciences, 2 (2009), 323-343.  doi: 10.1137/080725891. [22] G. H. Golub and C. F. Van Loan, Matrix Computations, Fourth edition. Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, 2013. [23] I. F. Gorodnitsky and B. D. Rao, Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm, IEEE Transactions on Signal Processing, 45 (1997), 600-616.  doi: 10.1109/78.558475. [24] L. Greengard and J.-Y. Lee, Accelerating the nonuniform fast Fourier transform, SIAM Review, 46 (2004), 443-454.  doi: 10.1137/S003614450343200X. [25] W. Guo and W. Yin, Edge guided reconstruction for compressive imaging, SIAM Journal on Imaging Sciences, 5 (2012), 809-834.  doi: 10.1137/110837309. [26] J.-Y. Lee and L. Greengard, The type 3 nonuniform fft and its applications, Journal of Computational Physics, 206 (2005), 1-5.  doi: 10.1016/j.jcp.2004.12.004. [27] H. Mansour and Ö. Yilmaz, Support driven reweighted? 1 minimization, in Acoustics, Speech and Signal Processing (ICASSP), 2012 IEEE International Conference on, IEEE, 2012, 3309-3312. [28] S. Osher, M. Burger, D. Goldfarb, J. Xu and W. Yin, An iterative regularization method for total variation-based image restoration, Multiscale Modeling & Simulation, 4 (2005), 460-489.  doi: 10.1137/040605412. [29] D. Rosenfeld, New approach to gridding using regularization and estimation theory, Magnetic Resonance in Medicine, 48 (2002), 193-202.  doi: 10.1002/mrm.10132. [30] L. I. Rudin, S. Osher and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena, 60 (1992), 259-268.  doi: 10.1016/0167-2789(92)90242-F. [31] T. Sanders, A. Gelb and R. B. Platte, Composite sar imaging using sequential joint sparsity, Journal of Computational Physics, 338 (2017), 357-370.  doi: 10.1016/j.jcp.2017.02.071. [32] T. Scarnati, Recent Techniques for Regularization in Partial Differential Equations and Imaging, PhD thesis, Arizona State University, 2018. [33] L. A. Shepp and B. F. Logan, The Fourier reconstruction of a head section, IEEE Transactions on Nuclear Science, 21 (1974), 21-43.  doi: 10.1109/TNS.1974.6499235. [34] W. Stefan, A. Viswanathan, A. Gelb and R. Renaut, Sparsity enforcing edge detection method for blurred and noisy Fourier data, Journal of Scientific Computing, 50 (2012), 536-556.  doi: 10.1007/s10915-011-9536-9. [35] Y. Wang and W. Yin, Sparse signal reconstruction via iterative support detection, SIAM Journal on Imaging Sciences, 3 (2010), 462-491.  doi: 10.1137/090772447. [36] D. Wipf and S. Nagarajan, Iterative reweighted $\ell_1$ and $\ell_2$ methods for finding sparse solutions, IEEE Journal of Selected Topics in Signal Processing, 4 (2010), 317-329.  doi: 10.1109/JSTSP.2010.2042413. [37] W. Yin, S. Osher, D. Goldfarb and J. Darbon, Bregman iterative algorithms for $\ell_1$-minimization with applications to compressed sensing, SIAM Journal on Imaging Sciences, 1 (2008), 143-168.  doi: 10.1137/070703983.

Figures(18)

Tables(1)