University of Southern California Dissertations and Theses

Multi-constrained inversion algorithms for microwave imaging

(USC Thesis Other)

Multi-constrained inversion algorithms for microwave imaging

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Multi-constrained Inversion Algorithms for Microwave Imaging by Pratik Shah A dissertation presented to the Faculty of the USC Graduate School In partial fulfillment of the requirements for the degree of Doctor of Philosophy (Ming Hsieh Department of Electrical Engineering) in the University of Southern California May 2018 Doctoral Committee: Professor Mahta Moghaddam, Chair Assistant Prof. Justin Haldar Research Assistant Prof. John Stang Professor Chunming Wang To my parents, Kosha ii A C K N O W L E D G M E N T S First of all, I would like to express my heartfelt gratitude to my faculty advisor and committee chair, Prof. Mahta Moghaddam for her years of support, encouragement, and motivation. I am thankful to her for providing an excellent environment for doing research. I would also like to thank my committee members, Prof. Justin Haldar, Prof. Antonio Ortega, Prof. John Stang, and Prof. Chunming Wang for their time, suggestions and numerous helpful comments. I learned and benefited a lot from discussions with many labmates of the Microwave and Systems, Sensors and Imaging Laboratory (MiXIL) and enjoyed every interaction with them. Special thanks to Uday Khankhoje, John Stang, Mark Haynes, Ruzbeh Akbar, Richard Chen, Agnelo Silva, and Guanbo Chen. I would also like to thank the entire MIXiL group for the fun, hours-long discussions, and friendship over the years. I am grateful to my close friends and family for encouragement. Many thanks to my friends at USC for making my life more balanced and enjoyable. I would also like to express my deep sense of gratitude to my parents. Their tremendous support, love, and encouragement had let me pursue what I am interested in. I am very fortunate to have my wife, Kosha along with me during my Ph.D. journey. I am thankful for her sacrifice, emotional support, and making my school years more pleasant with her sense of humor, laugh, and love. Last but not the least; I am very thankful to the ‘Brahman’ who blessed me with the zeal to learn and work hard. iii TABLE OF CONTENTS Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Microwave Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.1 Small scale examples . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.5 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.6 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2D Microwave Imaging 13 2 Joint Norm Constrained Imaging . . . . . . . . . . . . . . . . . . . . . . . 14 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Basic Theory and Inverse Problem Formulation . . . . . . . . . . . . 15 2.2.1 Inverse scattering problem . . . . . . . . . . . . . . . . . . . 15 2.2.2 Forward problem . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Applicability of Compressive Sensing . . . . . . . . . . . . . . . . . 16 2.4 Joint Norm Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4.1 Update of theL1 norm constraint . . . . . . . . . . . . . . 19 2.4.2 Minimize the residue using projected gradient method . . . . 21 2.4.3 Projection on a feasible set . . . . . . . . . . . . . . . . . . . 22 2.4.4 Update of the parameter 2 . . . . . . . . . . . . . . . . . . . 23 2.4.5 Undoing the effect of regularization . . . . . . . . . . . . . . 24 iv 2.5 Numerical Analysis and Results . . . . . . . . . . . . . . . . . . . . 24 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3 Shape Constrained Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3 Level Set Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3.1 Fast level set method . . . . . . . . . . . . . . . . . . . . . . 45 3.3.2 regularization . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4 Multi-material Formulation . . . . . . . . . . . . . . . . . . . . . . . 48 3.5 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5.1 Numerical approximation . . . . . . . . . . . . . . . . . . . 50 3.5.2 Selecting parameters . . . . . . . . . . . . . . . . . . . . . . 50 3.5.3 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5.4 Speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.6 Analysis and Results . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.6.1 Regularization effect . . . . . . . . . . . . . . . . . . . . . . 52 3.6.2 Single material . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.6.3 Robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.6.4 Comparative assessment . . . . . . . . . . . . . . . . . . . . 62 3.6.5 Multiple materials . . . . . . . . . . . . . . . . . . . . . . . 63 3.6.6 Breast phantoms . . . . . . . . . . . . . . . . . . . . . . . . 65 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4 Matrix Norm Constrained Imaging . . . . . . . . . . . . . . . . . . . . . . 68 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2 Inverse Problem Formulation . . . . . . . . . . . . . . . . . . . . . . 69 4.3 Cost Functions and Implementation . . . . . . . . . . . . . . . . . . 69 4.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.5 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . 72 3D Microwave Imaging 73 5 Combined Shape and Contrast Constrained Imaging . . . . . . . . . . . . 74 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . 76 5.3 Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.3.1 Shape estimation . . . . . . . . . . . . . . . . . . . . . . . . 78 5.3.2 Contrast estimation . . . . . . . . . . . . . . . . . . . . . . . 79 5.3.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.4 Implementation Details . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.4.1 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.4.2 Numerical approximation . . . . . . . . . . . . . . . . . . . 84 v 5.4.3 Selecting parameters . . . . . . . . . . . . . . . . . . . . . . 84 5.4.4 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.5 Analysis and Results . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.5.1 Comparative assessment . . . . . . . . . . . . . . . . . . . . 99 5.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Compressive Sensing For Microwave Imaging 102 6 Theoretical Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.2 RIPless Probabilistic Compressive sensing . . . . . . . . . . . . . . . 103 6.3 Generalized RIPless CS . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.3.1 Some useful derivations . . . . . . . . . . . . . . . . . . . . 107 6.3.2 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 6.5 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 7 Inverse Scattering Problems . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.2 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 7.3 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 7.4 Significance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 7.5 Application for Contrast Recovery . . . . . . . . . . . . . . . . . . . 127 7.6 3D Green’s Function . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 8 Concluding Remarks and Future Work . . . . . . . . . . . . . . . . . . . . 139 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 List of Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 vi LIST OF FIGURES 1.1 A graphical illustration of one of the equations for three unknowns points in 1D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 A graphical illustration of one of the equations for three unknowns points in 1D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1 Relative permittivity distribution for (a) actual object, (b) using just L2 norm with positivity constraint, (c) using joint norm with positivity con- straint, and (d) using joint norm with permittivity bound of [1, 3.5]. Imag- inary part of permittivity for (e) actual object, (f) using justL2 norm with negativity constraint, (g) using joint norm with negativity constraint, and (h) using joint norm being zero. Imaginary part is shown asIm. . . . . . 26 2.2 Horizontal cross section view from the center of the cylinder; the plot of real part of permittivity versus image points. Blue dotted line corresponds to actual profile, red line corresponds to joint norm and black line corre- sponds toL2 norm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3 Performance at various SNR values. . . . . . . . . . . . . . . . . . . . . . 29 2.4 Plot for performance of the reconstruction versus the ratio of the number of measurements to total number of unknowns. . . . . . . . . . . . . . . . . 30 2.5 Distribution of real part of permittivity for (a) original configuration, (b) justL2 norm recovery, (c) joint norm recovery. Distribution of the negative of imaginary part of permittivity for (d) original configuration, (e) justL2 norm recovery, (f) joint norm recovery. . . . . . . . . . . . . . . . . . . . 31 2.6 Original geometry for (a) rectangular cylinder, (b) triangular cylinder. Dis- tribution of relative permittivity recovered using the joint norm for (c) rect- angular shape object, (d) triangular shape object. The colormap shows the relative permittivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.7 Cross sectional view of relative permittivity retrieval for (a) rectangular cylinder and (b) triangular cylinder. . . . . . . . . . . . . . . . . . . . . . 34 2.8 Complex permittivity retrieval for rectangular cylinder. (a) real part of per- mittivity, (b) negative of imaginary part of permittivity. . . . . . . . . . . . 34 2.9 (a) Actual chili shape object. The reconstructed distribution of (b) real part of permittivity, (c) negative of imaginary part of permittivity. The colormap shows the relative permittivity. . . . . . . . . . . . . . . . . . . . . . . . . 35 vii 2.10 (a) original image for three objects configuration, (b) original image for two void holes in the box. Distribution of the recovered relative permittivity for (c) three objects and (d) void holes in a box. The colormap shows the relative permittivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.11 (a) large object, (b) permittivity 3, (c) permittivity 3.5, (d) permittivity 4, and (e) permittivity 6. The boundary of the actual object is overlaid (black line). The colormap shows the relative permittivity. . . . . . . . . . . . . . 37 2.12 (a) square dielectric object at center; the image reprinted from Fig. 6(d) of [1] (true values), (b) reconstruction results using TV; the image reprinted from Fig. 6(g) of [1], and (c) recovery using the joint norm. The colormap shows the permittivity contrast. . . . . . . . . . . . . . . . . . . . . . . . 39 3.1 Effect of regularization: (a) original geometry, recovered objects when (b) no regularization, (c) ==15. . . . . . . . . . . . . . . . . . . . . . . . 52 3.2 Single material configuration. Top row : actual geometry; center row : recovered geometry using fast level set method; bottom row : difference from the center row to the top row. The color bar indicates the contrast. . . . . . . . . . . . . . . . 53 3.3 For two ellipses case : (a) initial level set function, recovered objects after (b) BIM 1, (c) BIM 2, (d) BIM 3, (e) BIM 4, (f) BIM 5, and (g) BIM 6. The colormap shows the contrast. . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4 High contrast and fine resolution : (a) original configuration, (b) recovered, and (c) difference. The colormap shows the contrast. . . . . . . . . . . . . 57 3.5 Noise analysis: reconstructed images for SNR a) 25 dB, b) 20 dB, c) 15 dB, d) 10 dB, and e) 5 dB . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.6 Resolution analysis: two ellipses ( = 2). Reconstructed shape for the minimum distance (a) 1:1, (b) 0:9, (c) 0:5, (d) 0:3, (e) 0:1, and (f) 0:1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.7 Prior uncertainty: recovered objects when the prior contrast is (a) 20% lower than the actual contrast and (b) 20% higher than the actual contrast. 61 3.8 Three objects inL = 5 (a) reconstruction results using IMSMR-LS; the image reprinted from Fig. 5(a) of [2] , (b) reconstruction results using Bare-LS; the image reprinted from Fig. 5(c) of [2] (reprint permission to be acquired), and (c) recovery using our method. The colormap shows the contrast. The red contours show the boundary of the original configuration. 62 3.9 Three multiple materials cases - two circles, an object in a ring, and three objects. Each row corresponds to one of the cases. The left column shows the original configuration. The center column shows the recovered results using our method. The right column shows the difference from recovered to original. The colormap represents the contrast. . . . . . . . . . . . . . . 64 viii 3.10 Breast phantoms : Top row shows a sample from Class 1- Mostly Fatty. Bottom two rows show samples from Class 3- Heterogeneously Dense. Left column is the original relative permittivity distribution. Right column is the recovered using our method after 5 BIM iterations. The colormap shows the relative permittivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.1 Contrast source (Z) for clustered points: From left actual matrix represen- tation, mix norm recovery, and joint sparse and low rank recovery. Top row shows recovery for the real part and bottom row for imaginary part. . . . . 71 4.2 Contrast recovery. Geometry configuration - top row: clustered points, cen- ter row : distributed points, bottom row : large object. From left actual object, mix norm recovery, and joint recovery. . . . . . . . . . . . . . . . 72 5.1 Imaging chamber with 36 patch antennas. . . . . . . . . . . . . . . . . . . 84 5.2 Cube. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view at y = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. . . . . . . . . . . . 87 5.3 Vertical cross section view from the center of the domain; the plot of rela- tive permittivity versus points on z-axis. Black dotted line corresponds to actual profile and blue line corresponds to joint recovery. . . . . . . . . . . 88 5.4 Convergence plot. Relative error is plotted from the first iteration to the last iteration of the method in terms of the error at each level set iteration. . . . 89 5.5 Hollow sphere. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view aty = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. . . . . . . . . 91 5.6 Two spheres. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view aty = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. . . . . . . . . 92 5.7 Cross sectional view along z axis for the permittivity of (a) hollow sphere, (b) two spheres, (c) complex sphere, (e) two objects, and (f) three objects. (d) shows the cross section view along z axis for the conductivity of the complex sphere. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.8 Complex Sphere (real part). left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view aty = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. 95 ix 5.9 Complex Sphere (imaginary part). left column : actual geometry, right col- umn: joint contrast and shape recovery. The top row shows a 2D view at y = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates conduc- tivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.10 Two Objects. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view aty = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. . . . . . . . . 97 5.11 Three Objects. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view aty = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. . . . . . . . . 98 5.12 Comparative analysis. Top row : recovery by the proposed joint method, center row : recovery by the level set method with only bound constraint for the contrast estimation, bottom row : recovery by a traditional method. The left column shows the a 2D view aty = 0, the center column shows orthogonal slice view and the right column shows the max projection along the y-axis. The color bar indicates relative permittivity. . . . . . . . . . . . 100 6.1 the maximum off-diagonal element ofA A, as a function of number of measurements for fixed number of unknowns. . . . . . . . . . . . . . . . . 113 6.2 Average behavior of the maximum off-diagonal element of A A, as a function of number of measurements for fixed number of unknowns. . . . . 114 6.3 Comparative analysis. variation of inverse of sparsity to the number of un- knowns for three set of a number of measurements as a fraction of number of unknowns. The fractions are 0:2, 0:3, and 0:4. . . . . . . . . . . . . . . 115 6.4 Effect of random measurements on. . . . . . . . . . . . . . . . . . . . . 116 7.1 the maximum off-diagonal element of F , as a function of number of measurements for fixed number of unknowns. . . . . . . . . . . . . . . . . 118 7.2 One of the columns ofF , for different number of measurements. . . . . . 119 7.3 Bessel functions of different order versus the arguments . . . . . . . . . . 122 7.4 Bessel functions versus the order for different arguments . . . . . . . . . . 122 7.5 The effect on diagonal variations for different measurement radius,b for a given domain radius,a. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 7.6 The points that are failed to satisfy the condition onv T . . . . . . . . . . . . 125 7.7 error in the sparse recovery using CoSaMP. . . . . . . . . . . . . . . . . . 126 7.8 case 1: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 x 7.9 case 2: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.10 two points: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 7.11 small object: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 7.12 two high contrast: The top row shows the real part (in the left) and imagi- nary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configura- tion in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.13 high contrast: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the orig- inal contrast configuration in the left and recovered contrast configuration in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 7.14 zoom-in of 7.13: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the orig- inal contrast configuration in the left and recovered contrast configuration in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 7.15 square object: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the orig- inal contrast configuration in the left and recovered contrast configuration in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 7.16 off-center square object : The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 xi LIST OF TABLES 2.1 key symbols used in the article. . . . . . . . . . . . . . . . . . . . . . . . 21 2.2 quantitative measurement of the reconstructed images. . . . . . . . . . . . 32 3.1 quantitative measurement of the reconstructed images. . . . . . . . . . . . 55 3.2 Robustness Analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.3 Analysis of uncertainty in prior . . . . . . . . . . . . . . . . . . . . . . . 60 4.1 recovered the matrix rank and the mix-norm using two proposed methods . 71 5.1 Quantitative measurement of the reconstructed images. . . . . . . . . . . . 90 7.1 Quantitative measurement of the reconstructed images. . . . . . . . . . . . 138 xii ABSTRACT Imaging in the microwave regime of electromagnetic waves is a common tool for med- ical imaging, through-the-wall imaging, buried object detection, locating oil-gas in the Earth subsurface and many other industrial applications. The aim of the imaging is to estimate or map the electric properties such as permittivity and conductivity noninva- sively. Many developments have taken place in both hardware development and in the image reconstruction aspect of microwave imaging. However, still, the reconstructed images have poor resolution and often low estimation accuracy of the electrical proper- ties. The primary reason is the challenge of solving a highly non-linear inverse problem accurately. To address it, the first part of this thesis presents a detailed investigation on the effect of various constraints on the electrical properties of a 2D microwave imaging problem. In various simulation setups, different vector norms (` 1 and` 2 ), matrix norms (mix-norm and nuclear norm), bound constraints on permittivity and conductivity, as well as shape constraints have been applied separately or in some combination as con- straints. The ` 1 , ` 2 , and bound constraint together can produce reconstructed images with better estimation of the electrical properties than conventional constraints when the scatterers are not considered strong. Shape constraint via a level set function can have an almost perfect estimation of the electrical properties if the properties are known a priori but their location in the imaging domain is unknown. The matrix norm-based constraints can recover sparse strong scatters in a small imaging domain with a very low computational burden. xiii In the second part of the thesis, an experiment-ready 3D imaging setup is considered in simulation mode and a method is developed to solve two subproblems rather than solving in a single problem: shape estimation and the properties (contrast) estimation. The modular framework enables the freedom to choose any arbitrary constraints on the contrast in the cost function. Here, we demonstrate the applicability of a total-variation (TV) constraint for the contrast cost function. The reconstructed images show accurate estimates of object location, shape, and size, while recovering the contrast very well even with limited information. In addition to investigating the different constraints, the applicability of compressive sensing to inverse scattering problems from a theoretical point of view is also inves- tigated. The findings of the investigation are applied to derive new Lemmas, which broaden the applicability of compressive sensing theory from isotropic measurements to measurements that have a theoretical limit on maximally possible independent num- ber of measurements such as in inverse scattering problems. Then, the Lemmas are applied for a single source formulation in a 2D homogeneous domain, which show the possibility of reconstructing the image uniquely for the first-time for very sparse objects. xiv CHAPTER 1 Introduction 1.1 Background and Motivation The need of estimating material’s properties such as permittivity, conductivity, density, and compressibility arise in many sensing and imaging applications. If the estimation is carried out noninvasively, a large investigation area can be scanned quickly and can be monitored regularly. Non-accessible objects can be analyzed or visualized. The estimation is usually carried out by launching electromagnetic or acoustic waves into an investigation domain and then the resulting scattered field distribution is sampled with receiving antennas to estimate the properties by running a mathematical optimization or solving an inverse problem usually referred inverse scattering problems. One of the major applications of estimating permittivity or conductivity is in oil- gas industry, one of the crucial industrial sectors. The properties can be used to detect, locate, and quantify the amount of oil and gas in the earth subsurface because these properties of hydrocarbon, a principal ingredient of oil, are different than those of water. The measurements can be collected at the surface or through a single borehole or multi- borehole and then the inverse problem is solved. However, there are some challenges. The number of measurements is usually very limited. The inverse problem is non-linear and it is still an open problem to solve it correctly. We face similar challenges also when the density is estimated. The properties can also be used in estimation of flow parameters in the multiphase flows of oil, gas, and water. The estimation is very critical for maintaining the quality of products[3]. Even though a direct flow measurement is available, it is not reliable because it can estimate the flow for only single phase and separating the multiphase flow into three single phase flows may not be accurate. Also, sometimes the direct flow measurements are inaccessible. In these scenarios, indirect measurements through electromagnetic waves or acoustic waves are appropriate as the 1 properties of oil, gas, and water have different values. The another valuable application, where the properties can be used is to detect and localize the hidden objects, behind the wall or inside a room. These scenarios come up in military applications, security checking, law enforcement operations, or even in suspect-hostage situations. The problem formulation in these applications is termed as through-the-wall imaging. Initial approaches were focused on solving the approximated problem in a linear space so that images can be generated quickly and the problem can be solved in a real-time [4, 5, 6]. A recent review article, [7], summaries the different imaging schemes and highlights that the solving a full problem is needed to have better focused images. However, currently available methods for the full problem are not fast enough. For the near-surface geophysical exploration, archaeological, mine explorations, and environmental applications, objects are buried in the lossy earth at unknown depth. To cover larger depth, lower frequency waves are used but it imposes challenges on the in- version algorithm as it has to have an ability to produce super-resolution images beyond subwavelength so that objects can be located properly [8]. Multi-frequency approach can be used to increase the number of measurements to have better posed inverse prob- lem. [9] determined the location of buried objects by spectral analysis of surface waves. A grain storage monitoring system is also an evolving application where the per- mittivity estimation can be helpful because the moisture (water content) has different permittivity than that of grains. For wheat, an experimental validated system is already developed to test the moisture estimation performance on full-scale data [10]. The com- pute time and memory requirement to produce images and under estimation of permit- tivity are two major challenges that need to be addressed soon. In addition to above industrial applications, the electrical properties can also be used to perform a nondestructive evaluation in the nuclear industry [11] and inspect under- ground metal water pipelines. The noninvasive inspection of the underground pipes can help in identifying water leakage and hence reduce the cost of servicing underground water and sewage lines, which are rising rapidly in the UK [12]. Density and velocity of the acoustic waves (surface waves) have applications in seis- mology, especially in locating the source of seismic events and investigating the crust [13, 14]. The velocity estimation can also be useful in mitigating the hazards to build- ings by earthquakes [15, 16]. Full wave solutions of density and velocity are used in detecting shallow cavities [17]. However, the acoustic problems are large in terms of 2 wavelength, which required a longer time to solve the problem. Also, sometimes low density or velocity contrast makes harder to use acoustic waves for low-mid SNR sys- tems. Estimating the properties can also be useful in medical imaging. Tissue’s electric properties is determined by its molecular constituents, ion concentration and mobility, concentration of free water and bound water, and tissue temperature. Among them, water content is a major contributor, therefore, there is a natural contrast between water- rich tissues versus water deprive tissues such as bone and fat. This contrast can be used to identify the fractures in bones, soft tissue injuries, sports-related injuries in limbs, and abnormalities in liver. A malignant tissue can also have different properties because of higher cellular density and their hyper activity. Cancers in breast, brain, and colon have been investigating by many researchers. Among them, breast cancer has gained a lot of attractions. Breast cancer occurs in about 1 in 8 women (about 12%) in the United States [18] and is the second most common cause of death in US women [19]. The most significant action in fighting against breast cancer is to diagnose the shape and type of the tumor at early stages, before it alters into a malignant tumor. There is, therefore, a growing and urgent need for early stage detection and treatments of this cancer. Currently, it has been implemented by performing the routine check-up by mammography. A mammo- gram is the 2D projection of X-ray absorption of breast issue. Ultrasound and MRI are also used for the detection, whereas biopsy remains the gold standard. To estimate the effectiveness of different imaging modality in detecting the breast cancer a large scale analysis has been carried out in [20]. The mammography has 68% sensitivity (ability to detect malignant tumor) and 75% specificity (ability to avoid benign tumor to be de- tected as malignant). For the diagnosis by a ultrasound system alone, the sensitivity increases to 83% but very poor specificity (35%). MRI has even better sensitivity (94%) but the lowest specificity (26%) for an imaging system operating alone for the diagno- sis. A very high false positive rate or very low specificity requires unnecessary biopsies, which apart from being invasive can cause psychological problems. The overall bet- ter accuracy of 70% in the detection is achieved, when different modalities are used together. There are also health concerns related to ionizing radiation exposure with X- rays. These limitations and concerns have motivated research in alternative biomedical imaging techniques. 3 1.2 Microwave Imaging Frequencies in the microwave regime can be of interest in most of the above applications because microwaves can penetrate deep inside objects and can generate images with spatial resolutions that are suitable for many applications. For biomedical applications, microwaves are harmless, can penetrate inside a tissue, and have shown good sensitivity to benign and malignant lesions [21]. Based on their expected sensitivity, low cost, and non-ionizing nature, there has been substantial interest in developing microwave imaging systems for breast cancer detection. It should be noted that there are two main ways microwaves are used for biomedi- cal imaging: 1) radar based imaging [22], and 2) tomographic imaging. In radar based imaging, signal processing techniques are applied directly on the received signal simi- lar to conventional radars. The techniques indicate locations of strong scatterers, which can correspond to the contrast between healthy tissue and cancerous tissue. Sometimes, the radar based imaging is also classified as confocal imaging. Here, the image is re- constructed qualitatively by focusing the reflections [23] or by a space-time beamformer [24], which distinguish reflections from clutter and noise, taking into account frequency- dependent propagation effects. On the hand, tomography imaging is used to estimate the electrical properties at each pixel/voxel by a numerical optimization. We focus on the second type of microwave imaging. An extensive review on this topic [25] covers the description of different approaches as well as current fundamental challenges. The microwave imaging formulation can be derived from Maxwell’s equations [26] and can be represented either in differential equations or by an integral equation. We consider the microwave imaging, the scattering problem, is governed by the following electric field volume integral equation, expressed for a heterogeneous, isotropic, non- magnetic medium in a regionD, with measurement points on a surfaceS that encloses D, as E( r) = E in ( r) +k 2 b Z D G( r; r 0 )( r 0 ) E( r 0 )dv 0 ; (1.1) whereD is the imaging region, r 0 is the imaging domain position vector, and r is the measurement point position vector, E( r) is the total field, and E i ( r) is the incident field. The quantity G( r; r 0 ) is the Dyadic Green’s function, k b is the wavenumber of the background medium with lossless permittivity, and( r 0 ) is the dielectric contrast in 4 terms of the permittivity and conductivity contrast, defined as ( r 0 ) = [ r ( r 0 )= b 1]j[(( r 0 ) b )=( b 0 !)]; (1.2) where b is the lossless relative permittivity of a background medium, r denotes the relative permittivity, and b is the background conductivity. If we define the scattered field, E s ( r), as the difference between the total field, E( r), and the incident field, E i ( r), (1.1), the scattered field can be calculated outside the imaging domain as E s ( r) =k 2 b Z D G( r; r 0 )( r 0 ) E( r 0 )dv 0 : (1.3) The imaging problem deals with determining the dielectric contrast, ( r 0 ), of an unknown medium, given the measurement of scattered fields, E s ( r). Due to the pres- ence of the electric field within the integral, which itself is a function of the dielectric contrast, Eq. (1.3) is non-linear. We can consider( r 0 ) and E( r 0 ) together as a single unknown. This can reduce the problem into linear problem. However, the number of measurements increases an order lower than the number of unknowns. For example, the receiving antennas in 2D can lie in total ’length’ of 2r at a radiusr, while the number of unknowns within the measurement radius can lie in the ’length’ ofr 2 . For 3D, we get the unknowns proportional to r 3 and the number of measurements proportional to r 2 . This analysis is done with an assumption that the antenna can be as small as the pixel/voxel size, which may not possible. Therefore, for a relatively larger domain, r, the number of measurements will be too few. One of the ways, we can overcome the lower number of measurements is to use more than one source but thenE( r 0 ) changes with the source, which can increase the number of unknowns at the same rate. The inverse scattering problems exist for a long time and there are different ways, the non-linear is solved. We will discuss them in more details during the review of the previous work. 1.2.1 Small scale examples The equations are hyperbolic equations and to understand the behavior, we consider a very small scale experiment, where we can derive everything analytically. Let’s consider infinitely long objects in two directions of 3D and both are placed next to each other. This configuration reduces the problem in 1D and has only two unknowns. Let’s denote 5 their position by 1 and 2 and two receivers’ positions bya andb. The scattered field at receivera can be written in a scalar form from (1.3) as E s (a) =k 2 b 2 X j=1 G(a;j)(j)E(j): (1.4) Here, both (j) and E(j) are unknowns. First, we estimate E(j) as a function of using the forward solver equation, (1.1), analytically. For point 1, the forward equation can be written in a scalar form as E(1) =E in (1) +k 2 b 2 X j=1 G(1;j)(j)E(j): (1.5) For both points in a matrix form, " 1k 2 b G(1; 1)(1) k 2 b G(1; 2)(2) k 2 b G(2; 1)(1) 1k 2 b G(2; 2)(2) #" E(1) E(2) # = " E in (1) E in (2) # : (1.6) The matrix on the left side can be inverted analytically so we get closed form solution forE(j). Plugging that value in (1.4) gives " E s (a) E s (b) # = k 2 b " G(a; 1) G(a; 2) G(b; 1) G(b; 2) # " (1k 2 b G(2; 2)(2))E in (1)(1)k 2 b G(1; 2)(2)E in (2)(1) k 2 b G(2; 1)(1)E in (1)(2) + (1k 2 b G(1; 1)(1))E in (2)(2) # (1.7) where is the determinant of the matrix. if 6= 0, it has four terms-(1),(2), (1)(2), and constant. The matrix vector product in the right hand side can produce all terms but the constant. Rearranging terms, normalizing both equations separately so that the constant term = 1, and representing the factors by i , we get 1 (1) + 2 (1)(2) + 3 (2) + 4 = 0: (1.8) If we translate the origin such that 2 (1) (2) + 4 = 0, it is a standard equation of a rectangular hyperbola. The two rectangular hyperbolas can intersect and give two points: one in the first quadrant and the other is in third quadrant for the transformed coordinates. Non uniqueness! However, if we put put the physical constraint such as 6 permittivity has to be greater than 1, we can achieve a unique solution. For one of the equations, an equation is plotted for a set of values of(1) and(2) and is shown in Fig. 1.1. The hyperbolic behavior is clearly visible. Figure 1.1: A graphical illustration of one of the equations for three unknowns points in 1D. If we expand the domain to be of size 3 and perform the same analysis above (in- verting a 3 3 matrix), we get the equation of the following form: 1 (1) + 2 (2) + 3 (3) + 4 (1)(2) + 5 (2)(3) + 6 (3)(1) + 7 (1)(2)(3) + 8 = 0: (1.9) If we keep one of the variables constant, the form of the equation is again of a rectangular hyperbola. An exemplary shape is shown in Fig. 1.2. For a different measurement, there will also a rectangular hyperbola but with a separate set of parameters. For three mea- surements taken in three different directions, the intersection of three shapes is a point in each quadrant. This is an assumption at this point but can be rationalized by not- ing these are independent measurements and ‘hyperbolas’ are symmetric and behaving monotonically as we vary the variable. Similar to the two-points example, constraining to be physically valid gives us a unique solution. 7 Figure 1.2: A graphical illustration of one of the equations for three unknowns points in 1D. The method described above can not be applied to a larger set of unknowns because finding an inverse of a matrix analytically is not an efficient way and may not be practi- cally feasible. Forn unknowns, the inversion of the matrix by following the above steps can lead to estimaten(2 n 1) + 2 n2 (n 2 n) terms. This is a very large number. For example, n = 64, we get 9:88 10 21 terms. Also, the number of terms in the equa- tion grows exponentially as 2 n . Estimating the coefficient of each term will be a huge task and solving a equation in this form is very sensitive to a perturbation error, which can occur in noisy scenario. Therefore, to estimate the electric unknowns, the iterative approaches are more useful, which will be discussed in the next section. 1.3 Previous Work There has been substantial work done in developing algorithms to solve the nonlinear problems, in addition to the methods developed for an approximated nonlinear problem such Born or Rytov approximation [26]. We give you a short review here, where we cluster the methods based on what variables are being solved. The first cluster is based on Born series, where contrast and fields are explicit un- knowns. Usually, the first order approximation is used [27]. Here, the non-linear prob- lem is iteratively solved by estimating the total electric field and dielectric contrast in 8 an alternating manner during the inversion process. The inversion procedure is referred as Born Iterative Method. If the Green’s function is also updated, the inversion pro- cess is called Distorted Born Iterative Method (DBIM)[28]. Mathematically, DBIM is equivalent of Newton-Kantorovich method [29].The electric field estimation process is generally referred as a forward solver. Some standard forward solvers rely on one of the following methods- Finite-Difference Time-Domain Method, Finite Element Method, and Method of Moments. The inversion is usually carried out by the conjugate gradient method with Tikhonov regularization[30], or by Singular Value Decomposition [31]. In the second cluster, the contrast and the field are considered together as a sin- gle unknown referred as contrast-source. Hence, the unknowns are contrast-source and contrast. The problem becomes linear in contrast-source as unknown but this formula- tion has known to have a nontrivial null space. To overcome this problem, the authors in [32] presented a method, where minimum norm solution was found first and then contrast was estimated in the orthogonal complement basis of the solution using the forward equation. In [33], rather than finding the minimum norm solution, iteratively contrast and contrast-source were estimated. Recently, on a similar line, subspace-based methods were developed and extended to be applicable for moderately strong scatterer [34]. A different cost function is considered in Contrast source Inversion method (CSI) [35]. Here, the forward problem is considered as a constraint in the inversion cost func- tion. This single cost function is solved by conjugate gradient methods. Later, various variations of CSI are proposed to handle inhomogeneous background [36], differential formulation, and better contrast recovery[37]. CSI and DBIM have similar performance in contrast recovery [38]. In the third and the last cluster, contrast is only the explicit unknown. Global ap- proaches such as a genetic algorithm [39], a particle swarm algorithm [40], and an evo- lutionary algorithm[41] are applied for microwave imaging. When the new update on the contrast is generated by the global approach, a forward solver has to be run to esti- mate the scattered field. The forward solver has to be run for many iterations and hence the approach is very slow. A learning-based method is also proposed in [42] but it has a very restrictive use. Shape-based methods, either through explicit representation [43] or by an implicit method, [44, 45] have also been investigated. The methods can perform better but they are usually very slow. Developments also have been made in the understanding of the inverse scattering problems theoretically. The uniqueness of an inverse scattering solution is first proved 9 for 3D geometry. If all measurement directions, all aspect angles, and independent po- larizations are used, the refractive index is uniquely determined by the far-field pattern of the electric field represented in a 3D vectorial form [46]. For 2D geometry, it was eventually shown globally in [47]. As a next step, we ask how many samples of the scattered fields need to be collected? we generally assume that as more and more re- ceivers are used, the scattered fields can be represented more correctly. However, after the number of receivers reach to the Nyquist number, proportional to the domain ra- dius in an electrical dimension (/ (k b a),a radius of the domain) for 2D geometry, the difference in the scattered fields is indistinguishable[48]. This is sometimes referred as the degree of freedom of the scattered fields. In other words, we do not need infinite measurements to reconstruct the scattered field precisely. The analysis can be extended to multiple incident angles and 3D geometry [49]. If the plane waves are incident fields, the number of required sources is equal to the Nyquist number. 1.4 Objectives Microwave imaging targeted for medical applications are developed and eventually transformed to clinically viable devices [21, 50]. However, still the reconstructed im- ages have very poor resolution and lower estimation of contrast because there are still opened challenges in developing the hardware as well as unsolved mathematical prob- lems. These problems arise from the underlying highly non-linear inverse problem that we need to solve in the microwave-imaging problem to capture multiple scattering ef- fects. Also, there is a theoretically limit on the number of independent measurements that can be acquired. The number of independent measurements is far less than the number of unknowns present, which makes the system highly under-determined. The imaging system needs improvement in hardware design, efficient forward problem solver, and a better inversion algorithm. This thesis is focused on the last part. To be more specific, the objectives of this thesis are as follows. 1. Development of algorithms under multiple constraints, which can (a) improve contrast recovery. (b) Use prior information efficiently, if it is available. (c) solves robustly and in a reasonable amount of time. 10 2. Detailed investigation of an applicability of compressive sensing to inverse scat- tering problems and changes in theorems if required. Given the non-linearity of the problem, the work presented here has explored and ana- lyzed the efficiency of different constraints on solving the problem. 1.5 Thesis Overview This dissertation contains eight chapters, including this one. Chapters 2,3, and 4 are part of 2D microwave imaging. We use the Born Iterative Method based approach for chapters 2 and 3, and contrast-source as unknown for chapter 4. Each chapter has a different set of constraints. We do this to understand the effect of different constraints on inversion results. Chapter 2 has` 1 norm,` 2 norm, and box constraints, which help in producing better contrast recovery with very ’clean’ background. Chapter 3 has shape based constraints, which uses prior information to produce close to perfect shape recov- ery for a discrete set of contrasts. We use Born-type iterations to speed up the inversion process by reducing the number of a forward solver run. Chapter 4 has matrix norm based constraints, which does not need BIM type iterations or an inversion of a big ma- trix and can produce reasonably good results for a small size problem very fast. We compare the performance of joint sparsity constraints to nuclear norm plus sparsity con- straints. The second part of the thesis, Part II, is focused on combining some of the constraints derived in Part I and applying on an experimental setup in simulation mode for a 3D problem. We use S-parameter based measurements on a cubic imaging cavity rather than the scatted electric fields in free space. The Green function is numerically com- puted, which includes the effect of the cavity. Chapter 5 uses shape based, box, and Total Variational (` 1 in a gradient domain) constraints and recovers contrast and shape very well for a highly under-determined system with a limited number of measurements in a reasonably less time. The last part of the thesis investigates the applicability of compressive sensing to inverse scattering problems from a theoretical point of view. The findings of the investigation are applied to derive new Lemmas, which broaden the applicability of compressive sens- ing theory from isotropic measurements to measurements that have a theoretical limit on maximally possible independent number of measurements such as in inverse scattering problems. In particular, chapter 6 takes a different stand relying on matrix properties 11 rather than just on the number of measurements. We derive the conditions for unique recovery for` 1 minimization problem in terms of the magnitude of the element of the system matrix. Then, these conditions are validated on measurements taken in Fourier space. In chapter 7, we derive the uniqueness condition for a single source formulation in a homogeneous background in 2D and we show that it is possible to recover contrast uniquely for very sparse objects. Finally, Chapter 8 summarizes and concludes the thesis and gives several areas of future work. 1.6 Contributions Specific contributions of this dissertation include: 1. Chapter 2: Application of compressive sensing to nonlinear microwave imaging through Born Iterative Method (BIM) with joint norm constraints 2. Chapter 3: Reduction in the computational burden of a level-set method based shape recovery for microwave imaging by solving the problem in the BIM frame- work. Higher values of contrast are handled by mutli-level and multiple level-set functions. 3. Chapter 4: Direct contrast recovery under two different sets of matrix norm con- straints using contrast-source as unknown and without solving the forward prob- lem explicitly or as a constraint during the inversion. 4. Chapter 5: Development of an inversion method as a two-step recovery in an experiment-ready setup of 3D microwave imaging. Shape and contrast are esti- mated iteratively for improved object recovery performance. 5. Chapter 6: Derivation of two Lemmas for unique recovery in the compressive sensing framework so that the theory can be applicable to a limited number of measurements scenarios. 6. Chapter 7: Application of the derived Lemmas to inverse scattering problems. First-time demonstration of unique recovery for very sparse objects in the contrast- source formulation for a single source. 12 2D Microwave Imaging 13 CHAPTER 2 Joint Norm Constrained Imaging 2.1 Introduction A common practice to solve the non-linear imaging problem is to use local optimization strategies and linearize the problem iteratively, solving for the parameters and the field at each iteration [51] [28]. Global optimization methods have also been considered but are generally less popular due to large computational costs [41]. In local optimization methods, a linear approximation of the inverse problem can be described by an unconstrained optimization problem with the objective: min x kyAxk, where y denotes measurements of the scattered field, A is the system matrix (which is a function of the internal fields), and x denotes the unknown parameter of interest, such as the dielectric contrast. Generally, the system is under-determined as there are fewer measurements than the number of unknowns. In addition, the linear problem is ill-posed and ill-conditioned. To obtain a meaningful solution out of a multiplicity of solutions, various regularization schemes have used; generally inserted as a constraint or as a priori information. Regularization restricts the search space and thus reduces the degree of freedom in the solution. It also lessens the effect of noise, which is particularly important since the matrixA is usually ill-conditioned. L2 norm based regularizations are commonly used ([52] [53] [54] [55], among many); the Tikhonov method,L curve, and singular value decomposition (SVD) being examples of such methods [56] [57]. This type of regularization is known to produce an over-smooth solution [52], blur the edges/boundaries of the objects, overestimate the size of the object [58], and re-scale the solution [59]. Since the recent advent of compressive sensing [60],L1 norm based methods have gained popularity, showing promising results when the signal (x, in our case) is sparse in some basis, or is compressible [61] [62] and the system matrix obeys the restricted 14 isometric property (RIP). In such cases, the solution has sharp edges and preserves dis- continuities [58]. Work on compressive sensing in inverse scattering problems started with point-like targets [63, 64, 65, 66, 67, 68], and moved on to small-size objects of smaller size [69, 70, 71, 72]. In the latter, the inverse scattering problem was formu- lated as a Bayesian compressive sensing problem, and was shown to work without the need to satisfying RIP. However their approach did not consider complex-valued matri- ces, thus obstructing the use of the co-locations sparsity of real and imaginary values of the unknown vector. Other studies on compressive sensing in inverse scattering include total variation (TV) based approach and sparsifying the domain using a wavelet trans- form [73]. A recent article [1] based on this approach recovers the dielectric contrast within the Born approximation, although if the object has a smoothly varying profile of contrast, the assumption on sparsity fails. In this paper, an innovative joint norm based approach—based on the use of L1 andL2 norms as regularization—is proposed and is applied to a 2D microwave tomo- graphic imaging problem. We utilize the Born Iterative Method (BIM) to recover con- trasts higher than those allowed under the Born approximation. The proposed method estimates the regularization parameters and updates them as and when required. The approach exhibits the features of both norms–sharp at edges and smooth everywhere else. The main contributions of this work are (a) the derivation of a joint norm based inverse scattering strategy that, unlike similar state of art methods[70, 71, 72, 1] , effec- tively handles complex valued matrices and vectors for contrasts higher than the Born approximation and objects larger than point-like objects, (b) the introduction of anL2 norm regularizer to stabilize the electromagnetic imaging problems, when solve under the compressive sensing framework, and (c) an extension to the spectral gradient method by applying bound constraints and incorporatingL2 norm regularization. 2.2 Basic Theory and Inverse Problem Formulation 2.2.1 Inverse scattering problem The inverse scattering problem is governed by (1.3). In what follows, we assume az- independent two dimensional (2D) scatterer, and the electric field to be polarized along thez-direction, corresponding to a transverse magnetic (TM) polarization. 15 The inverse problem deals with determining the dielectric contrast, (r); r2 D, of an unknown medium, given some observations of the scattered field, E s (r); r 2 S. Due to the presence of the electric field under the integral sign, Eq. (5.1) is non- linear, which can be linearized by using an estimate of the electric field. This leads to a matrix equation of the form y = Ax that can be solved in an iterative manner; the Born iterative method (BIM) [74] is one such method. Here, y 2 C m , contains m measurements of the scattered field, x2 C n is the dielectric contrast for a domain discretized into n unknown values of contrast (in the pixel basis), and each element ofA2 C mn contains the product of the Green’s function of the background and the estimated field. The linearized inverse problem is still ill-posed and ill-conditioned, and requires regularization for successful solution. 2.2.2 Forward problem At every iteration of the BIM, we compute the electric field based on the current estimate of the dielectric contrast. For solving the forward electromagnetic scattering problem, we use a two-dimensional, vector element based finite element method (FEM) [75]. In this implementation of the FEM, we employ first-order Whitney edge elements as basis functions for the total electric field, and a first-order absorbing boundary condition to terminate the computational domain. This results in a sparse set of equations is solved very efficiently using a direct solver. 2.3 Applicability of Compressive Sensing In the compressive sensing framework, an accurate reconstruction of a sparse signalx can be obtained in the following reconstruction problem: minkxk 1 s:tkAxyk n ; (2.1) if the RIP holds, where n bounds the amount the noise in the data. The RIP is defined as [76]: given a matrixA and s 2 (0; 1), if the following relation holds (1 s )kx s k 2 2 kAx s k 2 2 (1 + s )kx s k 2 2 (2.2) 16 for alls-sparse vectorsx s ,A is said to obey the RIP to orders with an restricted isometry constant s . Even though the RIP is not a necessary condition, we are investigating the prop- erty here because we found in our numerical experiments that applying a method for L1 minimization directly in the BIM framework produces unstable and/or very sparse results. Hence, we are asking a question that can the RIP property be used to get some meaningful results? Re-arranging the terms in (2.2), we can express the RIP as: s max := max kAxk 2 2 kxk 2 2 (1 + s ) s min := min kAxk 2 2 kxk 2 2 (1 s ): (2.3) Given the vectorx s isssparse,s min is the square of the minimum singular value of allms submatrices [62] ofA. It is known that in inverse scattering problems, the sin- gular values ofA rapidly decay after a certain threshold, so the minimum singular value is close to zero [49]. As the number of non-zero entries inx increase (i.e.s increases) the minimum singular value decreases and for somes, the second condition in (2.3) is violated. We note in passing that in related literature on sparse targets (i.e. small values ofs), RIP may hold if thes th singular value (when the singular values are arranged in decreasing order) is in a range similar to the maximum singular value. 2.4 Joint Norm Formulation If we can modify the matrixA in such a manner that the minimum singular value of all submatrices can satisfy (2.3), the RIP can be satisfied and we can get anL1 minimized recovery of that modified system. One of the ways to achieve through a simple change is by using Tikhonov regularization. A Tikhonov regularized cost function is written as: J1(x) =kAxyk 2 2 + 2 kxk 2 2 ; J1(x) =kAxyk 2 2 ; (2.4) where we use the notation:A := h A p 2 I i T and y := h y 0 i T , and 2 is a regulariza- tion parameter. 17 Tikhonov regularization filters out smaller singular values, which can be shown through Singular Value Decomposition. The boost in minimum singular value is p 2 , so for an appropriate value of 2 , RIP can hold for the matrixA. We understand that satisfying the RIP forA may not seem logical but it can be argued that satisfying the RIP forA can be thought of as satisfying the RIP for an approximated ‘good ’ part of A. In the context of the ill-posed nature of the problem and imposing the L1 regular- ization, an additional regularization is needed to have bounded inverse operation, which the Tikhonov regularization provides and it also brings numerical stability during the inversion process. The updated conditions of (2.3) for the new matrix become: 2 +s max (1 + s ) 2 +s min (1 s ): (2.5) There is an upper bound [76] of p 2 1 on the value of s for the exact recovery of the reconstruction problem in (2.1). For any s < p 2 1, if 2 is selected within [1 s s min ; 1 + s s max ], (2.5) is satisfied and thus RIP holds. As an example, choose s =kxk 0 for the unknown vector, x, of the geometry in Example 1 of Section 2.5. The average values (over many sub-matrices ofA) ofs max ands min converge to 0:55 and 2:110 6 , respectively. Clearly, the RIP is not satisfied for thiss min . For the matrixA, if we choose s = p 2 1, the range for 2 is [0:59; 0:86]. Picking a value of 2 as 0:65, we get s A max = 1:20 and s A min = 0:65, and these values satisfy the conditions in (2.3). In this article, we justify the need of the joint norm and extend our previous work [77] with additional a priori information by means of placing box constraints on the unknown permittivity. For example, we assume that the relative permittivity of materials of interest must be greater than that of vacuum. Real part of the contrast is a relative measure with respect to the background and the background can have a value higher than that of the object. Thus, by considering two extreme values for background, namely lowest (vacuum) and highest ( b;max ), the range of permittivity contrasts for all materials we get is (1 + 1= b;max ; max 1]. Similarly the range for conductivity contrast is [ max = 0 !; 0]. Here, max and max represent the maximum relative permittivity and the maximum conductivity of the domain, respectively. Thus, our cost function becomes J(x) =kAxyk 2 2 + 1 kxk 1 ; x2X: (2.6) 18 where 1 is a regularization parameter of theL1 constraint, andX is the feasible set for the unknown vectorx. It should be noted that the 1 and 2 parameters are not related in any way, and are assumed to take any values inR + . Our solution strategy is essentially to estimate one parameter by minimizing the residue, while holding the other parameter constant. Our approach to solve the optimization problem is based on the ability to efficiently solve Lasso problems such as (2.6) (where a least squares solution is sought with a con- straint on theL1 norm of the solution) using the spectral projected gradient method [78] with range constraints and to update the 2 parameter if and when necessary. This prob- lem when solved as a Lasso problem is similar to [79] [80]. However, those approaches can be thought of as a special case of our generalized approach (no box constraints and 2 = 0). Let’s definekxk 1 =. The value of is related to 1 as follows. As the value of increases, the value of 1 should be decreased to keep the contribution of the L1 norm term to J(x) at a minimum. Thus, for appropriate choices of the parameters, solving the problem using or using 1 is equivalent. We consider the first approach (fixing kxk 1 =) to solve our problem. Hence, the cost function becomes J(x) =kAxyk 2 2 s.t. kxk 1 ; x2X: (2.7) The overall approach we use for solving this optimization problem is outlined in Algorithm 1 and elaborated further below. The key symbols are listed in Table 2.1. In the algorithm, , , , and r are error thresholds and l m is the maximum number of iterations. The algorithm has three key steps, which are iterated until a convergence criterion is satisfied; (i) updating theL1 norm constraint,, (ii) minimizing the residue using the projected gradient method by deploying the joint norm regularization, and (iii) estimating the parameter 2 . 2.4.1 Update of theL1 norm constraint Assume that at a given iteration we have an estimate of x for a given , but that the residue still needs reduction, i.e., the convergence criterion has not been satisfied. In such a case, has to be updated, generally to a larger value. We follow the method described in [79], where is expressed as a differentiable, convex, and strictly decreas- 19 Algorithm 1 Inverse scattering algorithm procedure JOINT NORM REGULARIZED ALGORITHM(A,x,y,, 2 , ,,, r ,l m ) set initialr 0 Axy set initialg 0 A H r 0 set initiall 0 setf l kr 0 k 2 2 define ^ f l = f=f l define ^ x l = x=x l X is the feasible set while (f l kl<l m ) do if ^ f l k ^ f l r kr l k 2 then % update + end if if ^ f l &k ^ x l k 1 then % update 2 estimate 2 updateA end if estimate % Barzilai-Borwein step length g l A H r l x l+1 X (x l g l ) % residue min. r l+1 Ax l+1 y setf l+1 kr l+1 k 2 2 l l + 1 end while end procedure 20 Table 2.1: key symbols used in the article. x vector of unknown dielectric contrasts X feasible set for the unknown dielectric contrast A new matrix generated by appendingL2 norm based constraint to the system matrixA y updated the measurement vector 1 regularization parameter forL1 norm based constraint 2 regularization parameter forL2 norm based constraint L1 norm ofx c newx in the gradient direction b upper or lower bound ofx, decided based on the sign ofc Lagrangian parameter for based constraint med normalized median value of the relative permittivity S the proportion of the number of measurements to total number of unknowns C o percentage of pixels correctly classified as background C i percentage of pixels correctly classified as the object ing function 1 (). The function represents an L1 constrained minimization problem (min x kAxyk 2 s.t. kxk 1 ) in a dual norm formulation [81]. Mathematically, the function is written for a variable,z2C n as follows; 1 () = inf x kAxyk 2 + sup z (x H zkzk 1 ) (2.8) For the optimum value of x := ~ x at a given , this function (called ~ 1 ()) becomes kA ~ xyk 2 . This monotonically decreasing function can be solved using Newton’s root finding method [82]. At each iteration, we update our previously estimated with an increment, = 1 () 0 1 () = ~ 1 () 2 kA H (A ~ xy)k 1 : 2.4.2 Minimize the residue using projected gradient method Given the parameters, and 2 , the given optimization problem becomes minimize x kAxyk 2 2 s.t. kxk 1 ; x2X: (2.9) 21 We rely on a projected subgradient method to get a feasible solution [78]. Projected subgradient methods minimize a function in the variablex subject to a constraint,x2 X, by generating a sequencex l via x l+1 = X (x l g l ) (2.10) where g l is the subgradient (gradient for non differentiable function (L1 projection) [83]) of the function, is the step length and X (x) is projection of x on X. Here, the projection also encompasses theL1 constraint. This minimization has two steps to perform: (i) minimize the residue using the gradient method, and (ii) project the estimate from the previous on the feasible set. For step (i), at every iteration, the functionalkAxyk 2 is minimized by moving in the opposite direction of the gradient and estimating the step length (computed using the method in [80]). For step (ii), the unconstrained update c l+1 :=x l g l has to be projected, which is discussed subsequently. 2.4.3 Projection on a feasible set The projection has to be performed in a way that the estimatedx l+1 satisfies the L1 norm constraint and is within the bounds. We want to estimatex l+1 from the current gradient estimatec l+1 by solving the following optimization problem. min x kcxk 2 2 s.t. kxk 1 ; x i;min x i x i;max ; i = 1;:::;n: (2.11) where the iteration number has been omitted for better readability, x i;min 2 C and x i;max 2 C are bound constraints on x i , which defines the feasible set X. Both are assigned based on a priori information. Different constraints for eachx i can incorporate more priors than constant bounds. The priors could include, apart from obvious range restriction onx i , positions and electric properties of the background medium. It can also include solutions of other methods. We solve this problem by formulating in the Lagrangian form and considering the 1 st order Karush?Kuhn?Tucker (KKT) conditions. Each element ofx, x i is estimated using the following relations, where b i is one of the bounds, chosen based on sign of c i , and is the Lagrangian parameter corresponding to the L1 constraint (refer to the 22 Appendix for details). x i = 8 > > > < > > > : 0 ifjc i j b i ifjc i j>b i + jc i j if<jc i jb i +: (2.12) 2.4.4 Update of the parameter 2 As mentioned earlier, the value of 2 has to be chosen in an interval as per Equation (2.5). However, this interval is not known a priori for all values of the sparsity parameter, s, and moreover, the value ofs itself is unknown. So rather than fixings or using the same interval for all cases, we obtain the interval from the current estimate of x, and then select one value from it. There are various methods in the literature which describe the estimation of 2 and also solve the Euclidean norm problem [52] [84] [85]. We use a method described in [86] because it can easily be altered to estimate 2 within a given interval. Mead [86] extended Rao’s fundamental theorem [87] to L2 regularized cost func- tions, stating that when model covariance matrix, C, is available, we may choose the data (x in our case) inverse covariance matrixC 1 x such that the cost function follows a 2 distribution with meanm , as closely as possible. If we assumeC 1 x = 2 I n then our cost function in (2.7) (without theL1 norm constraint) for the current estimate ofx can be equivalently written as 2 ( 2 ) = (Axy) H C 1 (Axy) + 2 kxk 2 2 . To satisfy Mead’s theorem, the cost function should stay within the following range: m p 2m z =2 < 2 ( 2 )<m + p 2m z =2 ; (2.13) wherez =2 is thez-value for a standard normal distribution having (1) confidence interval. If the tolerance related to confidence interval parameter is low, the above inequality approximately equals to the following equation after some algebra: 2 ( 2 ) =r H ((L 1 ) H U diag( 1 1 + 2 i 2 )U H L 1 )rm 0: (2.14) where L is Cholesky factorization matrix of C, U is the matrix corresponding to the SVD ofL 1 A asUV H , and i ; i = 1;:::;m are the singular values ofL 1 A. The equation is monotonically increasing as 2 varies from 0 to1. It is solved using 23 the Newton root finding algorithm, where we search for 2 within the current estimate of the interval. If the root exists and is positive, it is unique [86]. There are two other possibilities: (i) 2 (0) > 0, which effectively suggests that the equation does not need any regularization and we set 2 = 0. (ii) 2 (1) < 0, which indicates that no root exists. Here, we choose not to assign any value to 2 and use the previously estimated value. The advantage of this approach is that model parameters come directly from the ex- periment, thus avoiding any tuning parameters. It also considers statistical information available to define a physically meaningful regularization. 2.4.5 Undoing the effect of regularization AnL2 norm based solution scales the solution and anL1 norm based solution translates the solution [59]. Since the parameters are updated throughout the minimization, we consider the final values of the parameters as the scale/translation factors and use them to undo the effect. The scaling factor is not known a priori and is not easy to determine, so we apply the same factor which is known for the orthonormal case, i.e., (1 + 2 ). The translating factor is, which is the same as in (2.12). Generally, is higher than kxk 1 and therefore is zero. Therefore, we don’t apply any correcting factor for theL1 regularization. 2.5 Numerical Analysis and Results We performed several experiments to validate the proposed method. In the reported numerical results, scattered field data has been obtained by using the FEM. In order to avoid inverse crime, the mesh used for solving the forward problem has been chosen finer than that used for the inverse problem. The simulation configuration is as follows: The transmitting and receiving antennas are located at 20 apart from each other on the circumference of the circular domain of radius 1.5 (free space wavelength). The operating wavelength is 14:28cm and the domain is meshed at=100 to maintain the error in field calculation within bound. The object domain used for calculating scattered field is a circular domain of radius 1. The domain is meshed at =17 for the “measured” scattered field integral and at =50 for the forward solver. This gives an image size of 34 34 pixels inscribing the circular domain and thus the number of unknowns (only in the circular domain) for 24 the inverse problem has 2 912 = 1824 pixels. The total number of measurements is 18 18 = 324, making the inverse problem under-determined. When we estimate 2 , the model covariance matrixC is considered to be the identity matrix (because noise is zero). The bound forx is uniform for all pixels. Initially the value for is set to zero. The starting point for vectorx is the minimum energy solution (A H y) and 2 is estimated using this residue. Each image is formed with 7 iterations of BIM. An iteration consists of one run of the forward solver and up to 150 gradient iterations. The computation time to run a single forward solver is approximately 60 secs for code written in C++ and average time to run inverse problem is approximately 80 secs for code written in Matlab on a Linux desktop with a Xeon processor and 24GB RAM. Example 1 : circular cylinder : Here, the object under test is a homogeneous circular cylinder located at the center. The cylinder is 0:8 in diameter and has a relative permittivity of 2:5. The background is air. The reconstructed cylinders are shown in Fig. 2.1. The actual object and its real and imaginary permittivities are shown in Figs. 2.1 (a) and (e) respectively. 25 Figure 2.1: Relative permittivity distribution for (a) actual object, (b) using just L2 norm with positivity constraint, (c) using joint norm with positivity constraint, and (d) using joint norm with permittivity bound of [1, 3.5]. Imaginary part of permittivity for (e) actual object, (f) using justL2 norm with negativity constraint, (g) using joint norm with negativity constraint, and (h) using joint norm being zero. Imaginary part is shown asIm. 26 We compare the reconstruction performance by the joint norm to the performance by using only L2 norm. As discussed in Section 2.4.4, any method can be incorporated for 2 estimation. To make a fair comparison, we have used the same method as has been used in our joint norm, i.e. set 1 = 0, and apply positivity bound constraints for the real part ofx and negativity constraints for the imaginary part ofx. Figs. 2.1 (b) and (f) show the reconstruction of real and imaginary parts of the permittivity, respectively, us- ing justL2 norm. As expected and can be seen, the object in the images are smoothened and also the imaginary part has non zero values. Figs. 2.1 (c) and (g) show the recon- struction of the permittivity using joint norm with the same constraints as using justL2 norm. Figures clearly show the improvement in the sparsity of the background. The energy leakage in the imaginary part can be understood by the facts that y is not in the range ofA in the linearized BIM approach and the system is a highly under-determined system. If we have a priori information that the domain is lossless, bound on permittivity can be set to [1 +j0; 3:5 +j0]. The corresponding images are shown in Figs. 2.1 (d) and (h). The real part image has recovered higher contrast and image of the imaginary part remains as-is. In all of the cases, the object is detected at the correct location and its shape has been recovered with sharp boundaries. The peak relative permittivity at last iteration is approximately 2.72. Fig. 2.2 shows the cross section (horizontal) view of the actual and reconstructed permittivity using joint and justL2 norm. OnlyL2 norm based reconstruction underestimates the contrast and has non-zero values in background. The joint norm method estimates a reasonably good contrast and also produces a background that is very close to zero. 27 location 0 10 20 30 40 Dielectric Constant 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 actual l2 norm jointnorm Figure 2.2: Horizontal cross section view from the center of the cylinder; the plot of real part of permittivity versus image points. Blue dotted line corresponds to actual profile, red line corresponds to joint norm and black line corresponds toL2 norm. To quantify our observations, we used the following performance metric: C o (C i ) is the ratio of the number of pixels from outside (inside) the object classified as back- ground (object) to total number of pixels outside (inside) the object. We have used two thresholds to determine the class of each pixel: (i) lower threshold (1% in our case): all pixels below this level are classified as background, (ii) upper threshold (we took this to be 30% of the contrast): all pixels above this value are classified as object. The metric is applied to both permittivity images separately. Table 2.2 shows the performance values for all reconstructed objects. In the ideal case, C o andC i are 1. As can be seen from Table 2.2, the reconstruction of imaginary part of permittivity using joint norm is ideal. The contrast leakage into the background for L2 norm based approach is captured by a very small value ofC o . The improvement in contrast estimation using joint norm is evidenced by the high value ofC i . The effect of noise on the quality of reconstruction has also been analyzed. We added 28 white Gaussian noise with zero mean. Fig. 2.3 reports the errors (whole domain) for sev- eral values of SNR (defined with respect to the scattered field as in [1]). The quantity med represents the normalized median value of the relative permittivity. Mathemati- cally, med = median( l r ) o r , where l r is the estimated relative permittivity atl object points and o r is the average value of the original relative permittivity of the object. As can be seen, contrast recovery is almost constant after 10 dB SNR, whereas C i reaches a maximum value of 1 for the no-noise case. C o keeps its value around 90% for all SNR values. Overall, the contrast is very well recovered for all SNRs and the object is also localized correctly showing that noise at these levels does not have a major impact on the performance. SNR(dB) 0 10 20 30 40 50 Normalized value 0.6 0.7 0.8 0.9 1 C o C i @ med Figure 2.3: Performance at various SNR values. We have also analyzed the effect of the number of measurements on the reconstruc- tion quality. Fig. 2.4 shows the plot forC i , C o and med values for various values of S. The quantity S is the proportion of the number of measurements to total number of unknowns. As the number of measurements increase, the quality of reconstruction becomes better. However, we do not always have the luxury of having a large number of measurements in practice. Still, we still have reasonable accuracy forS = 8%. 29 S(%) 6 8 10 12 14 16 18 Normalized value 0.7 0.75 0.8 0.85 0.9 0.95 1 C o C i @ med Figure 2.4: Plot for performance of the reconstruction versus the ratio of the number of measurements to total number of unknowns. Next, using the same geometric parameters, we have formed an image of a cylin- der having the imaginary part of permittivity -1.5 as shown in Fig. 2.5(a) and (d). The complex permittivity recovery is shown in Fig. 2.5(c) and (f). The object is detected with good estimated conductivity contrast. The reconstruction using just L2 norm is also shown in Fig. 2.5(b) and (e). We have used the same constraints for both cases : positivity for real part and negativity for imaginary part. The conductivity contrast is underestimated for justL2 norm. The real part is zero everywhere for joint norm, which is significantly different than that forL2 norm. MinimumL1 norm based constraint has played a role to make it sparse. 30 y/ 6 1 0 -1 (a) 1 0 x/ 6 1 1.1 1.2 -1 0 r y/ 6 1 0 -1 (d) 1 0 x/ 6 -1 0.5 1 0 1.5 - 0 r " y/ 6 1 0 -1 (e) 1 0 x/ 6 1.5 -1 0 0.5 1 - 0 r " y/ 6 1 0 -1 (b) 1 0 x/ 6 1 1.1 1.2 -1 0 r y/ 6 1 0 -1 (f) 1 0 x/ 6 1.5 0 -1 1 0.5 - 0 r " y/ 6 1 0 -1 (c) 1 0 x/ 6 1 1.1 1.2 -1 0 r Figure 2.5: Distribution of real part of permittivity for (a) original configuration, (b) justL2 norm recovery, (c) joint norm recovery. Distribution of the negative of imaginary part of permittivity for (d) original configuration, (e) justL2 norm recovery, (f) joint norm recovery. 31 Table 2.2: quantitative measurement of the reconstructed images. Object L2 norm joint norm C o C i C o C i C o C i C o C i real real imag. imag. real real imag. imag. circular 2:5 0:27 0:93 0:44 0 0:87 1 1 1 circular j1:5 0:64 0:67 0:40 0:59 1 1 0:77 0:91 rectangular 2.5 0.07 0.62 0.07 0 0.90 0.97 1 1 rectangular 1:6j0:8 0.52 0.75 0.24 1 0.84 0.96 0.91 1 triangular 2.5 0.05 0.92 0.1 0 0.87 1 1 1 chilli 1:6 j0:7, 2:1j0:3 0.38 0.47 0.56 0.82 0.76 0.74 0.85 0.98 three objects 2.5, 2.1, 1.9 0.52 0.97 0.20 0 0.89 0.93 1 1 rings in a box 2.5 0.45 0.70 0.22 0.1 0.85 0.95 1 1 big triangular 2.5 0.01 0.5 0.32 0 0.57 0.91 1 1 high contrast 1 3 0.01 0.97 0.65 0 0.81 1 1 1 high contrast 2 3.5 0.01 0.37 0.24 0 0.71 0.50 1 1 high contrast 3 4 0.01 0 0.24 0 0.59 0.1 1 1 high contrast 4 6 0.01 0 0.24 0 1 0.03 1 1 Example 2 : other shape cylinders : The object is a homogeneous rectangular cylinder located at the center of the domain. The rectangle is 0:50:8 (shown in Fig. 2.6(a) ) and has a relative permittivity of 2.5. The reconstructed distribution is shown in Fig. 2.6(c). As can be seen, there is a good agreement with the actual configuration. Corresponding values forC o andC i are reported in Table 2.2. The bound is [1+j0; 3:5+ j0]. The cross sectional view of the recovery for joint norm and L2 norm is shown in Fig. 2.7(a). The edge of the rectangular cylinder is detected at the correct positions for joint norm, whereas L2 norm based approach fails to localize the object. Next, when the object is not loss-less and has complex permittivity (1:6j0:5), joint norm based approach is still able to recover both parts correctly. The reconstruction quality is shown in Fig. 2.8 and the quantitative values are reported in Table 2.2. Here, the bound is [1j1:5; 3:5 +j0]. 32 Figure 2.6: Original geometry for (a) rectangular cylinder, (b) triangular cylinder. Dis- tribution of relative permittivity recovered using the joint norm for (c) rectangular shape object, (d) triangular shape object. The colormap shows the relative permittivity. Next, we have considered a cylinder having the shape of a triangle, which is also located at the center (Fig. 2.6(b)). The equilateral object has a side-length of 0:866 and has a relative permittivity of 2.5. The reconstruction after 7 BIM iterations is shown in Fig. 2.6(d). The triangle’s edges are well localized. The recovered peak relative permittivity is 2.62. The bound is also the same as the bound used in recovering only relative permittivity for the circular cylinder. The localized quality shown in Fig. 2.7(b) is also supported by having high values for performance metric (Table 2.2). 33 Figure 2.7: Cross sectional view of relative permittivity retrieval for (a) rectangular cylinder and (b) triangular cylinder. Figure 2.8: Complex permittivity retrieval for rectangular cylinder. (a) real part of per- mittivity, (b) negative of imaginary part of permittivity. Example 3 : chili shape object : In this example, the object is composed of three el- lipses placed close to each other to generate an irregular ‘chilli’ shaped object as shown in Fig. 2.9(a). The complex permittivity is (1:6j0:7) at head and toe of the object and (2:1j0:3) at remaining places. Fig.2.9 shows reconstructed permittivity distribution. The peak values are 1.9 and -0.6 respectively for real and imaginary parts. The bound for permittivity is [1j1; 3:5 +j0]. 34 Figure 2.9: (a) Actual chili shape object. The reconstructed distribution of (b) real part of permittivity, (c) negative of imaginary part of permittivity. The colormap shows the relative permittivity. Example 4 : three objects : A more complicated object is considered here, consist- ing of three objects: circular, elliptical, and rectangular. Their sizes are 0:3 in diameter, 0:4 0:3 axial length, and 0:2 0:3 side length respectively. Their relative per- mittivities are 2.5, 1.9, and 2.1 respectively. The reconstruction distribution is shown in Fig. 2.10(c). As can be seen from Fig. 2.10(a), there is good agreement with the actual configuration. 35 Figure 2.10: (a) original image for three objects configuration, (b) original image for two void holes in the box. Distribution of the recovered relative permittivity for (c) three objects and (d) void holes in a box. The colormap shows the relative permittivity. Example 5 : two rings in the box : The box is characterized by relative permittivity = 2.5, length = 0:85 and width = 0:4. The circular void holes (background permittivity) of radius 0:15 and centers at (0;0:2) and (0; 0:2) respectively, are present inside the box cylinder. The reconstructed distribution is shown in Fig. 2.10(d). Comparing it with Fig. 2.10 (b), the proposed approach provides quite good results. In particular, the algorithm is able to detect the edges of box correctly and also the location and size of the circles. 36 Figure 2.11: (a) large object, (b) permittivity 3, (c) permittivity 3.5, (d) permittivity 4, and (e) permittivity 6. The boundary of the actual object is overlaid (black line). The colormap shows the relative permittivity. 37 Example 6 : limits : We use sparsity based constraints, so that if we do not apply any compressible transformations to convert non sparse objects into sparse objects, the recovery may not be accurate. The recovery of the non sparse object depends on the number of measurements, and the size of the object. The higher the number of measure- ments, the larger the size of the object that can be recovered. As an example we tried to recover the triangular cylinder of 1:4 side length (instead of 0:866) and kept the same number of the measurements (324) and unknowns (1824). The permittivity was set to 2.5. The recovered object is shown in Fig. 2.11(a). The algorithm was not able to recover the shape but contrast has recovered reasonably well. Since the object is big, the valid range for 2 in some iterations is not feasible so the recovery is as good as just usingL2 norm (see Table 2.2). There are two possible ways to address the problem; i) increase the number of measurements and reconstruction domain so that it becomes sim- ilar to previous examples, or ii) apply a transformation (wavelet, Fourier, total variation etc.) and solve the problem in that domain (in this case, the elements of the unknown vector,x, become the coefficients in the transformed domain). However, the bounds can not be known in the transformed domain and need further study. We attempted to recover a circular cylinder having a very high contrast (maximum permittivity = 6). The reconstruction results for permittivity values 3, 3.5, 4, and 6 are shown in Figs. 2.11 (b)-(e). It is a well known fact that the BIM based linearized ap- proach can not recover very high contrasts and multiple scattering effects can not be well incorporated in the first order approximation. Generally, the observed limit on ac- ceptable recovery of absolute value of permittivity is around 2.5. In our case also, the reconstruction algorithm fails to recover the shape and contrast of the object, when the contrast is higher in similar range. We can adopt DBIM or some other technique here but we are more interested to see if the sparsity constraint can recover higher contrast than other regularized methods. Example 7 : comparative assessment : To make a fair comparison with a state of art TV based compressive sensing method for inverse scattering problem, we selected one of their results as a reference. The method is valid only in the BA region [1] so we restricted our method to run only for one iteration. We keep the number of measure- ments, the geometry, and SNR same but the synthetic data is generated using FEM with the triangular mesh. Each side of the pixel in the reconstruction domain has the length same,=6, as used in the TV based method. The images of [1], reprinted in Fig. 2.12(a) and (b) are the original image and the image recovered using the TV based constraint. 38 The image reconstructed using our method is shown in Fig. 2.12(c). (a) (b) (c) x/ 6 -1 -0.5 0 0.5 1 y/ 6 -1 -0.5 0 0.5 1 0 0.02 0.04 0.06 0.08 0.1 Figure 2.12: (a) square dielectric object at center; the image reprinted from Fig. 6(d) of [1] (true values), (b) reconstruction results using TV; the image reprinted from Fig. 6(g) of [1], and (c) recovery using the joint norm. The colormap shows the permittivity contrast. As can be seen, the joint norm based method has more accurate and sparse recon- struction. Quantitatively, the error index tot , metric used in [1] to quantify the error, drops from 1:5% to 0:7%. The size of the object considered in this experiment is bigger than other experiments mentioned in the article. The reason for achieving better result for the large object is the use of larger number of measurements. We also tested the performance of the joint norm method for other examples in [1], but in the interest of 39 space, we are not putting all the graphs but just mentioning that the joint norm method performs better than the TV based method. It recovers the contrast much better, keeps the most of the background to zero and also has the better definition of edges. 2.6 Conclusion We presented a jointL1-L2 norm based regularization for inverse scattering. The pro- posed approach is based on sparsity and smoothness, which allows obtaining the so- lutions without over-smoothening at discontinuities. The regularizing parameters are estimated during the iterative process. The algorithm can also incorporate a priori infor- mation about the upper and lower bounds of the range of contrasts. The reconstructed images for various shapes showed very good shape recovery, localization, and contrast estimation up to a contrast of about 2.5. Some theoretical open issues still exist and need further consideration. In particular, bound constraints for large objects need to be inves- tigated, if the objects are represented in a sparsifying domain. Furthermore, recovery for higher contrasts must be analyzed. Future work includes evaluating the performance of the method on experimental data. 40 CHAPTER 3 Shape Constrained Imaging 3.1 Introduction The non-linear problem has been solved in several ways. Approaches such as the Born Iterative Method (BIM) and the Distorted Born Iterative Method (DBIM) solve the non- linear problem by applying linear approximations and updating them iteratively[26]. The linear approximations are to keep either the contrasts (complex permittivity of the object with reference to that of the background) or the electric fields constant, when we update the other. In practice, they perform well when the object is not too large compared to the wavelength or the contrast is not too high. The imaging problem is also ill-posed, which makes the choice of the regularization a salient factor [88]. A different formulation, Contrast Source Inversion (CSI), where the electric fields and the contrasts are formulated in a single cost function, has also similar performance [38]. Alternatively, the level set method has been applied for microwave imaging [45, 89, 90, 91, 2, 92, 93], which is a popular technique in image segmentation and image registration problems because of its implicit and simple representation of shapes [94, 95, 96]. For the microwave imaging problem, if accurate priors about dielectric properties are available, the problem can essentially be translated to detecting, locating objects, and estimating their shapes. In other words, it can be thought of as a contour identification or as a shape segmentation problem. It should be noted that even if the contrast values are known, conventional microwave imaging techniques such as BIM, DBIM, or CSI do not necessarily perform better. The published level set methods are computationally intensive. For example, for a single type of inhomogeneity (single material) the approach used in [90] takes around 4-7 hours on a single 64-bit AMD Opteron 246 processor running at 2 GHz to recover a 2D shape profile in the investigation domain of 10 cm radius in the frequency range 41 of 3-8 GHz. The main computation burden lies in the estimation of the electric field. This step is commonly known as a forward solver. The published work on level set methods for the imaging problem or for the inverse scattering problem are formulated using an adjoint field representation and calculation[97], where the forward solver is run twice at each iteration. Generally, the number of iterations is on the order of hundreds or sometimes thousands. Recently, a method for 3D level set has been proposed [98], where only a single run of the forward solver is needed at each iteration. It reduces the computational complexity by half; however, it still requires more than 50 iterations (50 runs of the forward solver) even for the setup where the initial guess is derived from the result of a DBIM method. In another approaches, the linear sampling method is integrated with the level set method to speed up the computation [99, 100]. We also proposed a method in [101], where we applied some approximations that allowed us to run the forward solver only a few times. However, it is defined only for one type of inhomogeneity. For more than one type of inhomogeneity (multi-material case), the level set methods are successfully applied in many areas [102, 103, 104, 105, 106]. For microwave imag- ing, the multi-material level set methods can be solved using adjoint-based approaches [107]. The level set functions are represented in either color, vector, or single format and within these different representations, color representation gives a better result in general[108]. However, all of them also suffer the same problem of being computa- tionally intensive. In addition, some of the approaches have shown better performance for realistic and complex problems, when the problem is restricted to a single type of geometry (breast phantom). However, here the adjoint-based method has to run until it converges, and this has to repeat for each level set function [109, 110]. In this article, we present a level set based method, which is not based on the adjoint formulation and requires very few forward solver iterations, resulting in much reduced computation time. We formulate the level set method as a gradient flow optimization, and borrow a concept from BIM to reduce the number of forward solver runs. We make the assumption that the contrasts are known with some uncertainty, although it is under- stood that the prior knowledge about contrasts in the domain may not be very accurate or easily available. We test this assumption later in this article. We then extend this formu- lation to a more generalized multi-material representation and show results on simulated data. The results also include numerical breast phantoms, where unlike previous work, all level set functions are updated together, so the required runs of the forward solver 42 are independent of the number of the level set functions. These developments and im- provements are important in solving a well known inverse scattering problem, which is known to have many minima, by leveraging the prior to reduce non-uniqueness and by formulating the level set method in an efficient way to speed up the computations. The novelty of the article lies in the development of an approach that requires very few forward solver runs, while producing results as good as or better than other available methods. Our approach can be thought of as the combination of BIM and level set method. This method can handle multiple materials and can generate accurate results with fewer forward solver runs and fewer algorithmic stages. The paper is organized as follows. In Section II the problem formulation is detailed. Section III describes the level set method for microwave imaging and explains key com- ponents of our approach. Section IV extends this to multiple type of inhomogeneities. Implementation details are described in Section V and evaluation of the algorithm on simulated data is performed in Section VI. Finally, conclusions and future directions are drawn in Section VII. 3.2 Theoretical Background The inverse scattering problem is governed by (1.3). In what follows, we assume az- independent two dimensional (2D) scatterer, and the electric field to be polarized along thez-direction, corresponding to a transverse magnetic (TM) polarization. The imag- ing problem deals with determining the dielectric contrast, (r 0 ); r 0 2 D, of an un- known medium, given some observations of the scattered field, E s (r); r2 S. Due to the presence of the electric field within the integral, which itself is a function of the dielectric contrast, Eq. (5.1) is non-linear. We use the linear operator Lf:g to repre- sent k 2 b R D G(r;r 0 ):dv 0 . Using this notation, the scattered field equation is written as E s =Lf;Eg. For a forward problem, we use the volume integral equation to solve the electric field,E(r 0 ), given the dielectric contrast. Mathematically, E(r) =E in (r) +k 2 b Z D G(r;r 0 )E(r 0 )(r 0 )dA 0 ; (3.1) wherer is defined everywhere. The integral equation is solved on a discrete grid by a standard method, e.g., methods in [111, 112, 113]. 43 3.3 Level Set Formulation The imaging problem is solved by minimizing the errorkE s Lf;Egk, while esti- mating the contrast and the field. Even though permittivity and conductivity are well documented for many materials [114], their effective use in the inversion process varies a lot. For example, they may be used to define an upper/lower bound of the search space or may be used as an average value of the material to guide the minimization process. In this article, we utilize this information as a contrast prior. We are aware that this as- sumption may not be very accurate, and later in the article, we examine the effect of an imprecise prior on object recovery. Given the contrast, the imaging problem essentially becomes estimating shapes of inhomogeneity inD. We consider a level set method for shape estimation. In a level set formulation, the contrast is parameterized using a higher dimensional continuous function, called a level set function . We first formulate for one type of inhomogeneity (single material) in the domainD. The contrast is represented by as (r 0 ) = i U((r 0 )l) + e (1U((r 0 )l)); (3.2) where i and e are prior information on the object and background contrast, respec- tively, andU(r 0 ) is the standard Heaviside function. The region occupied by the object is represented byl and the background is represented by<l, so =l traces the boundary of the object. For simplicity, we considerl = 0. Point should be noted that there is no constraint on connectivity of an object. The minimization of any cost function,J, is achieved in the form of the evolution of the level set function as @ @t + @J @ = 0;r = 0 onD; (3.3) where J in our case is J() =kE s Lf;Egk, called residual error, and D is the boundary of the domainD. The equation can be interpreted as follows. If the residual error is zero, the level set function remains constant during evolution. For a non-zero error, the level set function changes. The change can be set such that the error always de- creases. The change is normally set through a gradient flow. Formally differentiatingJ and the first order Taylor series approximation yields the Gateaux derivate of functional 44 J() as @J @ =Ref[R 0 ()] + R()g( i e )(): (3.4) Here,Re indicates real part of the argument,R is the residue (R =Lf;EgE s ) for the given estimation ofE using, ‘ + ’ indicates adjoint operation andR 0 denotes the derivative operation with respect to. Another way to interpret the evolution of the level set function is that the convergence or minimum residual error is achieved whenever @J @t decreases to zero forJ 0. Using the chain rule, we get @J @t = @J @ @ @ @ @t : (3.5) If @ @t is set to negative of the right side of (3.4), @J @t is always negative, i.e., J de- creases at every time step. Since J is greater than zero, eventually the residual error becomes zero. 3.3.1 Fast level set method Since the electric field is a function of object contrast, the above derivative is carried out by defining an adjoint problem. In the adjoint problem, the adjoint fields are defined and an equivalent forward problem is solved, which is computationally very intensive. This particular procedure needs two forward solver runs per iteration. We deviate from the standard procedure and solve the problem in the following way. Before we start, we make a note on the BIM method. BIM is a linearized iterative approach, where the fields and the contrasts are estimated alternatively in each iteration. When we update one, we keep the other constant. The contrast is generally updated using a gradient method (with some additional regularizations), which can be written as i+1 = i g l ; (3.6) whereg l is the gradient of the approximated linear system and is some constant de- noting step length along the gradient direction. Along the same lines, we formulate an iteration for the level set method as, i+1 = i tg; i+1 = threshold( i+1 ); (3.7) 45 where g is the Gateaux derivate given in (3.4). If we apply the BIM concept,R becomes a function of only in (3.4). In this case, the quantity in (3.4) becomes linear and can be solved easily and faster than the adjoint based problem. Let us use the notationA to representk 2 b R D G(r;r 0 )E(r 0 ):dv 0 , where the fields,E(r 0 ), are the current approximation of the electric fields. Using this concept and notation, (3.4) becomes @J @ =RefA + R()g( i e )(): (3.8) Even though approximatedg overlaps with the linear system’s gradientg l and differs by some weight factor, a major difference between these two is that the effect ofg on is not linear (‘threshold’ function). Apart from the non-linear effect, the prior in our method reduces chances of the solution getting trapped in local minima. The rationale behind this is as follows. The derivative of the residue,R 0 (), can be decomposed into two gradients as R 0 () =A + G @E @ ; (3.9) where is the Hadamard product and G is the Green’s function defined above. The derivative of the electric field with respect to the contrast, @E @ , can be derived from (3.1), which can be used to demonstrate that the small change in contrast at one point can affect the electric field in the entire domain (we refer this as a nonlocal effect.). This nonlocal effect is ignored in the BIM approximation, which is reasonable when the change in the total field is small. However, when the approximation does not hold, the incorrect gradient can lead the solution to local minima. Since in our approach l can take only a few distinct values, the effective gradient is different than that of the approximation and can compensate for the nonlocal effect eventually as the distinct values of l are correct. So our method can outperform a linearized approach. Summing up, our approach utilizes the best from both methods - speed and accuracy. The overall approach we use for solving the imaging problem is outlined as follows. The approach has two iterations. (1) outer iteration, which consists of a forward solver run and level set based optimization; we refer to this iteration as a BIM iteration; (2) level set iterations, which is the number of times the level set based optimization has run in a single BIM iteration. The BIM iterations are kept running until the method converges or the maximum number of BIM iterations is reached. 46 3.3.2 regularization During the level set evolution, irregularities develop for many reasons, which eventually destroy the stability of the level set evolution and create spurious objects. One of the popular remedies is reinitialization [96]. The problem with reinitialization is that it is not known when exactly the level set function has to be reinitialized. Too-frequent reinitialization makes the overall process very slow, and too infrequent initialization defeats its purpose. However, a recently developed method, called variational level set formulation [115] resolves this problem by having reinitialization as a regularization. Generally, reinitialization is carried out to make the level set function as the signed distance function (jrj = 1). The regularization term is defined such that it drives to keepjrj = 1. The definition, similar to [115] is R p () = Z p(jrj)dx; (3.10) where p is a potential energy, which is chosen such that the level set function stays as the signed distance function. The Gateaux derivative of the functionalR p () is @R p @ =r: p 0 (jrj) jrj r : (3.11) Whenjrj is close to zero, the derivative of R p () reaches1. To avoid this singularity behavior, the potential energy is defined similarly to [115] as p(s) = 8 < : 1 (2) 2 (1cos(2s)) ifs< 1 1 2 (s 1) 2 ifs 1: (3.12) We define our cost function as a weighted summation of these two terms, namely, residual energy and non uniform distance energy. Formally, it is written as E() =R p () +J(); (3.13) whereR p () is the regularization term,> 0 and are constant,J() is the difference between measured and estimated scattered fields. The value of helps in satisfying the necessary evolution conditions. The value of controls the effects of residual energy. The regularization is a function of only, and sometimes it is called the internal energy of the function, whereas the residual error is called the external energy of the function 47 . For the regularized cost function, equation (3.3) becomes @ @t + @E @ = 0. Substituting the derivatives, we get the level set evolution equation as @ @t =r: p 0 (jrj) jrj r RefA + R()g( i e )(): (3.14) 3.4 Multi-material Formulation We extend the inhomogeneity of the same type to N types (N-material). This can be achieved either by introducing additional level set functions or by more levels for each level set function. We use both approaches to generalize the formulation. The more general formulation is more useful when we are interested in solving a known geome- try, because certain features can be embedded more efficiently. For example in breast imaging, skin can be represented efficiently using an additional level of the level set function defined for fatty tissue, whereas fibroglandular tissue can be represented effi- ciently using an additional level set function. Mathematically, the formation is written as: @J @t = N 1 X j=0 @J @ @ @ j @ j @t ; (3.15) where the contrast is defined as (r 0 ) = N 2(1) N 2(2) :::N 2(N 1 ) X k=1 k N 1 Y j=1 U( j L kj & j <L kj+1 ); (3.16) where k are the contrast priors, N 1 is the number of level set functions, N 2(i) is the number of levels plus 1 for thei th level set function, andL is the permutation matrix of sizeN 2(1) N 2(2) :::N 2(N 1 ) N 1 . Each element,L kj , takes one of the values from all levels L a =fl 1 ;l 2 ;:::l N 2(j) 1 g. For a boundary element at the lower bound ofL a , the lower bound is1 and a boundary element at the upper bound ofL a , the upper bound is1. For example, 9 materials can be represented by 2 level set functions (N 1 = 2) and each has 2 levels (N 2(1) = N 2(2) = 3). Let’s choosel 1 = 0 andl 2 = l. Then, the contrast is 48 represented by (r 0 ) = 1 U( 1 >1 & 1 < 0)U( 2 >1 & 2 < 0) + 2 U( 1 >1 & 1 < 0)U( 2 0 & 2 <l) + 3 U( 1 >1 & 1 < 0)U( 2 l & 2 <1) + 4 U( 1 0 & 1 <l)U( 2 >1 & 2 < 0) + 5 U( 1 0 & 1 <l)U( 2 0 & 2 <l) + 6 U( 1 0 & 1 <l)U( 2 l & 2 <1) + 7 U( 1 l & 1 <1)U( 2 >1 & 2 < 0) + 8 U( 1 l & 1 <1)U( 2 0 & 2 <l) + 9 U( 1 l & 1 <1)U( 2 l & 2 <1): (3.17) In this article, we consider two special cases: 2material and 3material. For the 2material case, we consider two different settings N 1 = 1; N 2(1) = 3 and N 1 = 2; N 2(1) = N 2(2) = 2. If the number of materials is less than the number of priors required, multiple k can be of the same material, i.e., same value. For the 3material case, we considerN 1 = 2; N 2(1) =N 2(2) = 2. The contrast derivative for the 3material case can be written as, @ @ 1 =( 1 )f( 3 1 )U( 2 ) + ( 4 2 )U( 2 )g @ @ 2 =( 2 )f( 2 1 )U( 1 ) + ( 4 3 )U( 1 )g: (3.18) Applying the chain rule and the regularization, the following two level set functions are solved together: @ 1 @t =r: p 0 (jr 1 j) jr 1 j r 1 1 RefA + R() @ @ 1 g @ 2 @t =r: p 0 (jr 2 j) jr 2 j r 2 2 RefA + R() @ @ 2 g; (3.19) where 1 and 2 control the evolution of these two level set functions, respectively. The same concept can be applied to the remaining cases and in general to theNmaterial case. 49 3.5 Implementation Details 3.5.1 Numerical approximation For solving the evolution of the level set function, we use the central difference scheme to approximate all of the spatial partial derivatives, whereas we use the forward dif- ference scheme to approximate the temporal partial derivative. The approximation of (3.14) or (3.19) by the above difference scheme can be simply written as i+1 i t =L( i ); (3.20) where L( i ) is the approximation of the right-hand side of (3.14) or (3.19) and t is artificial time represents the evolution iterations. We choose and t such that the Courant-Friedrichs-Lewy(CFL) conditiont< 1=4 is satisfied [115]. 3.5.2 Selecting parameters We consider t = 1 and choose the remaining parameters accordingly. First, we de- cide the value of . Depending on the range, within which the gradient stays at the initial stage and the range of the level set function, the range of is decided and then by trial and error, the best value is considered. To avoid the longer run of the level set optimization, we restrict the level set function to lie within [2; 2]. This bound is not necessary but it helps in ‘generating’ and ‘clearing’ objects quickly. We select such thatt< 1=4 and<jj=2. The latter condition is to make sure that the regulariza- tion term does not dominate over the residual error during the minimization process. 3.5.3 Initialization In the fast level set method, the level set function has to be initialized at every BIM step. We can initialize the level set function once at the beginning of the process and for all successive BIM iterations, the initial level set function can continue from the level set function at the last step of the previous iteration. Since we have the regularization, which preserves the signed distance function, this option is viable. However, if the al- gorithm gets stuck in local minima, it may take a longer run to emerge so we choose to reinitialize the level set function at every iteration. We keep the same initialization for every BIM iteration. We tried four different initial configurations for the level set func- 50 tion; ‘concentric’- a map of the signed distance function as concentric circles, ‘zeros’- everywhere zero (all background), ‘binary’- some pixels as objects and the remaining as the background, and ‘many’ - many binary objects. Even though [115] suggests the use of binary as easy and reliable initialization, we found that it is not the case for this non-linear process. We consider ‘concentric’ and ‘zeros’ as our choices for the initial- ization. In our experience, both can produce the similar results with the proper choice of the parameters. We increase level set iterations by some number at every BIM iteration so that the level set method has enough iterations to converge. For multi-material, we can choose the same or two different initializations, if the second level set function is present. We initialize asA T y, whereA is created using the Born approximation. Then, we start the BIM iterations- the forward solver and the level set iterations. 3.5.4 Speed We may need more level set iterations than they are needed in other published level set methods. Does that compensate for a lower number of forward solver run? The answer is No for all practical problems. The reason is as follows. In our case, the level set optimization updates the geometry by matrix-vector product, whereas the forward solver generally needs a matrix inversion and/or very finer computational grid. Since the matrix-vector product has one order lower computational complexity than the matrix inversion, running longer run of the level set optimization does not adversely increase the computational complexity. 3.6 Analysis and Results We perform several experiments to validate the presented method. In the reported nu- merical results, the scattered field data are obtained by using (3.1). The simulation is performed using the following configuration unless otherwise specified. The T/R an- tennas are located on the circumference of the circular domain of radius 25 cm at 20 o apart from each other and the investigation domain is within a square of 36 cm length. The background medium is considered to be as air/vacuum. The wavelength in air is 15 cm and pixel size in both directions is 0:75 cm, which gives image size of 48 48 pixels. The current setup is clearly underdetermined (2304 unknowns against 324 mea- surements). The synthetic data is generated at 0:60 cm resolution to avoid the inverse 51 crime. 3.6.1 Regularization effect We illustrate the importance of the regularization for a case having noisy measurements (SNR= 20 dB). The objective is to recover the inhomogeneity shown in Fig. 3.1(a). We minimize (3.14) two times, where all parameters are kept the same for both times except. For the recovery shown in Fig. 3.1(b), it is set to 0, where for the recovery in Fig. 3.1(c), it is set to=15. The regularization prevents the appearance of the ghost objects/pixels and keeps the level set function as the signed distance function for better recovery. (a) -0.2 0 0.2 x(m) -0.2 0 0.2 y(m) 0 0.2 0.4 0.6 0.8 (b) -0.2 0 0.2 x(m) -0.2 0 0.2 y(m) 0 0.2 0.4 0.6 0.8 (c) -0.2 0 0.2 x(m) -0.2 0 0.2 y(m) 0 0.2 0.4 0.6 0.8 Figure 3.1: Effect of regularization: (a) original geometry, recovered objects when (b) no regularization, (c) ==15. 52 3.6.2 Single material The first experiment is on objects with the same type of inhomogeneity. The geometry of the problem and all four cases are shown in Fig. 3.2. In all simulations, the permittivity of the objects is 2. (a) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (e) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (i) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) -1 -0.5 0 0.5 1 (b) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (f) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (j) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) -1 -0.5 0 0.5 1 (c) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (g) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (k) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) -1 -0.5 0 0.5 1 (d) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (h) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (l) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) -1 -0.5 0 0.5 1 Figure 3.2: Single material configuration. Top row : actual geometry; center row : recovered geometry using fast level set method; bottom row : difference from the center row to the top row. The color bar indicates the contrast. In the first case, an equilateral triangle has side length of 21cm and is located at the center. In the second case, two same size ellipses of 18 cm major axis length and 7:5 cm minor axis length are located on x- axis at 6:75 cm either side of the origin. The third 53 case is a circle in a ring. The ring has 10:5 cm outer radius and 6:75cm inner radius. The circle is located at the center of the ring and has 2:25 cm radius. The entire geometry is shifted by 3:75 cm in both directions. In the fourth case, a star shape object is located at center has a side length of 21 cm. These configurations are shown in Fig. 3.2(a)- (d). The results shown here are generated by calling the forward solver only 15 times or until it converges, and 60 level set iterations for the first call of the forward solver. The iterations are increased by 10 for every successive forward solver run. The level set evolution equation’s parameters are set as = 0:1 and = =15. The respective reconstruction images are shown in Fig. 3.2(e)-(h), which show that the objects are located almost perfectly and the shapes are recovered very well. For the first three cases, the residual error is below the threshold just after six forward solver runs. The differences from the recovered images to the original images are shown in Fig. 3.2(i)- (l), which show the majority of the error at the boundary within a layer of single pixel size. The major reason for the error is the different grid size during synthesizing and reconstruction. For all cases, the level set function is initialized as shown in Fig. 3.3(a). The same initialization is used at every initial step of successive BIM iterations. The reconstructed results for the second case after every BIM iteration are shown in Fig. 3.3(b)-(g). As can be seen, from the third iteration onwards, the shape of the object is coming out for this case. To quantify our observations, we use the following performance metric: C o (C i ) is the ratio of the number of pixels from outside (inside) the object classified as background (object) to total number of pixels outside (inside) the object, where i is the index and its maximum value equals to the number of the materials present in the domain. The values inC i s are sorted according to the contrast from the lowest to the highest. If more than one object has the same contrast, we consider just one object (onlyC 1 ). In the ideal case,C o and allC i s are 1. As can be seen from Table 3.1, the reconstructions using the level set method are exceptionally well as the accuracies are higher than 90%. We also consider the percentage error, which is defined as, = 1 NM v u u t N X i=1 M X j=1 rec r (i;j) ori r (i;j) rec r (i;j) 2 100%; (3.21) where rec r and ori r are complex permittivity of a recovered object and an original object 54 Table 3.1: quantitative measurement of the reconstructed images. Object Accuracy Error (%) C o C 1 C 2 C 3 Triangle 2 0:98 0:99 - - 0:22 Two Ellipses 2 0:99 0:97 - - 0:34 Circle in Ring 2 0:91 0:98 - - 0:32 Star 2 0:98 0:99 - - 0:23 High Contrast 7 0:97 0:96 - - 1:56 Two Circles 1:5; 1:8; 2 0:96 0:83 0:88 0:96 0:31 Object in Ring 1:5; 2 0:95 1 0:76 - 0:27 Three Objects 1:5; 1:8; 2 0:96 0:97 0:97 0:61 0:34 Phantom 1 many (C) 0:91 0:96 1 - 16:70 Phantom 2 many (C) 0:95 0:76 0:92 - 16:39 Phantom 3 many (C) 0:97 0:70 0:95 - 19:43 respectively. Next, we consider the same shape as we had for one of the cases in the first exper- iment, but the permittivity of the object is increased to 7. Correspondingly, the pixel size is decreased to 0:5 cm and 0:3 cm for the inversion and synthesizing, respectively. The wavelength is increased to 30 cm and the number of measurements is increased to 900. All the parameters are kept the same for the level set function except. is set to 1:2. The recovered shape after the convergence (5 BIM iterations) is shown in Fig. 3.4(b). The recovery is very similar to the low contrast and coarse resolution image. The majority of the errors are at the boundaries and are within a layer of 1-2 pixels (see Fig. 3.4 (c)). The accuracy is still higher than 95%. As can be seen in Table 3.1, the error,, is higher than 1% because the contrast can take only one of the two values and if the value is incorrectly recovered, the percentage error increases with increase in the contrast. 55 (a) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) -0.5 0 0.5 1 1.5 2 (b) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 (c) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 (d) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 (e) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 (f) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 (g) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 Figure 3.3: For two ellipses case : (a) initial level set function, recovered objects after (b) BIM 1, (c) BIM 2, (d) BIM 3, (e) BIM 4, (f) BIM 5, and (g) BIM 6. The colormap shows the contrast. 56 (a) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 2 4 6 (b) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 1 2 3 4 5 6 (c) -0.1 0 0.1 x(m) -0.15 -0.1 -0.05 0 0.05 0.1 0.15 y(m) -6 -4 -2 0 2 4 6 Figure 3.4: High contrast and fine resolution : (a) original configuration, (b) recovered, and (c) difference. The colormap shows the contrast. 3.6.3 Robustness Next, we evaluate the performance of the method under noise in a sequence of exper- iments from SNR = 5 dB to SNR = 25 dB. All the parameters are kept the same for the level set function as they are for the first experiment. The reconstruction results are shown in Fig. 3.5. The method can recover shapes and size very well for SNR above 15 dB. Even for the low SNRs, the method detects the number of objects and positions correctly. The corresponding values of andCs are reported in Table 3.2. 57 Figure 3.5: Noise analysis: reconstructed images for SNR a) 25 dB, b) 20 dB, c) 15 dB, d) 10 dB, and e) 5 dB In order to evaluate the resolution efficiency of our method, the distance between the center of the two ellipses is varied in a sequence of experiments from 1:6 to 0:4, or in other words, the minimum distance between the two ellipses is varied from 1:1 to0:1. To avoid any error because of the geometry mismatched during synthesizing and reconstruction, we commit inverse crime here. The reconstruction results shown in Fig. 3.6 are using the parameters of the first experiment. Our approach can distinguish objects, which are 0:3 apart for 0:05 pixel size. Even for the distances smaller than this, the method can still recover objects very well as expected. A point should be noted that the scattered field is undersampled and all scattered field information [49] is not measured. The corresponding values of andCs are reported in Table 3.2. The better recovery is expressed correctly by higher values of performance metrics. 58 Table 3.2: Robustness Analysis. SNR Accuracy Error (%) Resolution Accuracy Error (%) C o C 1 (distance) C o C 1 25 dB 0:97 0:99 0:34 1:1 1:00 0:99 0:09 20 dB 0:97 0:99 0:34 0:9 1:00 0:99 0:08 15 dB 0:96 0:96 0:37 0:5 1:00 0:98 0:12 10 dB 0:95 0:93 0:42 0:3 0:98 0:94 0:29 5 dB 0:94 0:84 0:48 0:1 0:98 0:96 0:25 0:1 0:99 0:95 0:23 (a) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (b) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (c) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (d) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (e) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (f) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 Figure 3.6: Resolution analysis: two ellipses ( = 2). Reconstructed shape for the minimum distance (a) 1:1, (b) 0:9, (c) 0:5, (d) 0:3, (e) 0:1, and (f)0:1 . 59 Table 3.3: Analysis of uncertainty in prior object 20% lower 20% higher Accuracy Error (%) Accuracy Error (%) C o C 1 C o C 1 Triangle 0:96 0:84 0:30 0:99 0:94 0:21 Two Ellipses 0:93 0:93 0:41 0:98 0:90 0:33 Circle in Ring 0:94 0:96 0:38 0:98 0:78 0:37 Star 0:95 0:99 0:34 1:00 0:93 0:2 The next experiment aims to analyze the dependency of the recovery accuracy on the accuracy of the prior information for the cases shown in Fig. 3.2. We initialize the prior by introducing 20% error. More specifically, the permittivity of the objects is assumed to be 1:8 and 2:2, while the scattered fields are measured for permittivity 2. We apply our approach with no change in any parameters and the reconstructions are shown in Fig. 3.7. The results show the robustness against the incorrect prior. As can be seen, the quality of results turns out to be similar to that with the exact knowledge of the permittivity as in Fig. 3.2. The corresponding values of andCs are reported in Table 3.3. 60 (a) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 (b) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 (c) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 (d) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 (e) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 (f) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 (g) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 (h) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.5 1 Figure 3.7: Prior uncertainty: recovered objects when the prior contrast is (a) 20% lower than the actual contrast and (b) 20% higher than the actual contrast. 61 3.6.4 Comparative assessment To make a fair comparison with a state-of-the-art 2-D level set method for microwave imaging, a sample is selected as a reference. The geometry is matched to the setup in [2] and the same SNR is used. We increase the number of measurements to match with that in [2]. The method is multiregion with a variable scale. Therefore, to make a fair comparison, we choose an example where the method is compared to Bare-LS (single-region level set method). The images of [2], reprinted in Fig. 3.8(a) and (b) are the image recovered using IMSMR-LS and the image recovered using the Bare-LS. The image reconstructed using our method is shown in Fig. 3.8(c). (c) -2 -1 0 1 2 x/ 6 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 y/ 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Figure 3.8: Three objects inL = 5 (a) reconstruction results using IMSMR-LS; the image reprinted from Fig. 5(a) of [2] , (b) reconstruction results using Bare-LS; the image reprinted from Fig. 5(c) of [2] (reprint permission to be acquired), and (c) re- covery using our method. The colormap shows the contrast. The red contours show the boundary of the original configuration. 62 As can be seen, the fast level set method has better recovery than Bare-LS and sim- ilar to IMSMRA-LS. Quantitatively, the total number of iterations are 302 and 396 for Bare-LS and for IMSMRA-LS, respectively. For those methods, forward problem has to be solved two times for every iteration. To evaluate the computational complexity, we use the quantity defined in [2],f pos , which is the product of the cube of total number of unknowns, number of iterations, and number of forward solver runs per iteration.f pos is 9:44 10 12 and 4:63 10 11 , for Bare-LS and for IMSMRA-LS, respectively. f pos for the fast level set method is computed as 1:22 10 11 , which is lower by more than a fac- tor of 3 to IMSMRA-LS and 70 to Bare-LS. The reason is that even though we solve at full resolution, the number of forward solver runs is only 10. The error metrics reported in [2] are [3:4 10 2 ;; 3:2 10 2 ] and [2:8 10 3 ; 1:7 10 2 ; 4:5 10 2 ], respectively for the square object, the ring object, and the cross object. The same metric for our method is [6 10 3 ; 1:92 10 2 ; 7 10 3 ], which is similar in numbers to IMSMRA-LS. The values forC o ,C i , and are 0:74, 0:98, and 0:24% respectively. 3.6.5 Multiple materials The next experiment deals with the multiple materials scenario. We change the pixel size during inversion to 0:88 cm (=17), a non-integer multiplier of the pixel size during synthesis. We consider three sample cases. In the first case, two circles have the same radius of 6 cm and the center is shifted at [3 cm; 3 cm] and [3 cm; 3 cm] respectively. The second case is an object in a ring, which is located at the center and the ring has 3 cm width and has the outer radius of 10:5 cm. The object has 7:5 cm radius. The third case is three objects; triangle, circle, and square. The triangle is equilateral, located at [7:5 cm; 6 cm] and has 12 cm length. The circle has 5:25 cm radius and is located at [2:25 cm; 8:25 cm]. The square has 10:5 cm length and is located at [7:5 cm;4:5 cm]. The permittivities of the scatterers are 1:5; 1:8; and 2 as shown in Fig. 3.9 (a), (d), and (g). 63 (a) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (b) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (c) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) -0.5 0 0.5 (d) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (e) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (f) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) -0.5 0 0.5 (g) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (h) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) 0 0.2 0.4 0.6 0.8 1 (i) -0.1 0 0.1 x(m) -0.1 0 0.1 y(m) -1 -0.5 0 0.5 Figure 3.9: Three multiple materials cases - two circles, an object in a ring, and three objects. Each row corresponds to one of the cases. The left column shows the original configuration. The center column shows the recovered results using our method. The right column shows the difference from recovered to original. The colormap represents the contrast. The circle in a ring case is represented using only the single level set function with l 1 = 0; l 2 = 0:2 because of its geometry. To avoid any influence of the initialization, we choose the ‘zeros’ as initialization. For remaining cases, two level set functions are used. The initialization of these level set functions is ‘concentric’ and ‘zeros’, respectively. The maximum BIM iterations are 25 and for the first iteration, the number of level set iterations is 30 and is increased by 10 for successive iterations. 1 is 1:8; 0:5; and 1:5 for the Fig. 3.9 (a), (d), and (g) cases, respectively. = 0:01 and 2 = 0:5 1 for the multiple level set functions cases. The recovery results are shown in Fig. 3.9. The BIM converges at 15; 20, and 23 iterations, respectively. As can been seen, our method 64 can recover the shape and size very well. It is also able to detect the contrast at proper locations. The spatial distribution of the error has the same behavior as it has for single material cases (see Fig. 3.9 (c), (f), and (i)). The quantitative numbers are reported in Table 3.1, where values for are similar to those for the single material cases. FewC i s have relatively lower values due to the lower number of pixels for that material. 3.6.6 Breast phantoms Next, we consider a more complex geometry, where we used MRI-derived breast phan- toms from a public database [116]. The repository has phantoms of 4 ACR classes. These are used to generate maps of dielectric constant and conductivity, which are then used to synthetically generate the microwave measurements in a coronal plane of 0:5 cm thickness for each coronal plane, while neglecting inter-slice interaction. The distribu- tion of the permittivity for a few cases are shown in Fig. 3.10 (a), (c), and (e). For these cases, conductivity has similar distribution but in the interest of saving space, we are not including all of the graphs. When the scattered fields are generated, we keep tissue variations including transition tissue as they are and also consider the contrast as a complex quantity. We consider two level set functions with 2-material andl = 0. The two materials - fatty tissue and fibroglandular tissue are in a homogeneous background. Even though the transition tissue is present but since it is not a physical tissue, we have not considered it in our prior. The permittivity variations for fatty tissue and fibroglan- dular tissue are from 3j1:5 to 9j6 and from 36j5 to 65j19, respectively. The complex permittivity for tissues as priors are assumed to be 3j1 and 49j10. The background medium is lossless and has permittivity 10. Both level set functions are initialized as concentric circles centered at the origin. We have 324 measurements and the image size is 3726. The results are shown in Fig. 3.10 (b), (d), and (f) after 5 BIM iterations, which depict the correct identification of tissues. The method can recover size and shape of the tissues reasonably well. Here, we add the variation in each tissue, which is linearly proportional to the gradient at each point, at the end of every BIM iteration. The quantitative numbers are reported in Table 3.1, which show good values for performance metrics. 65 (a) -0.05 0 0.05 x(m) -0.05 0 0.05 y(m) 10 20 30 40 50 (b) -0.05 0 0.05 x(m) -0.05 0 0.05 y(m) 10 20 30 40 50 (c) -0.05 0 0.05 x(m) -0.05 0 0.05 y(m) 10 20 30 40 50 (d) -0.05 0 0.05 x(m) -0.05 0 0.05 y(m) 10 20 30 40 50 (e) -0.05 0 0.05 x(m) -0.05 0 0.05 y(m) 10 20 30 40 50 (f) -0.05 0 0.05 x(m) -0.05 0 0.05 y(m) 10 20 30 40 50 Figure 3.10: Breast phantoms : Top row shows a sample from Class 1- Mostly Fatty. Bottom two rows show samples from Class 3- Heterogeneously Dense. Left column is the original relative permittivity distribution. Right column is the recovered using our method after 5 BIM iterations. The colormap shows the relative permittivity. 66 3.7 Conclusion A fast level set based method for 2D microwave imaging has been presented. Our ap- proach is based on integrating concepts from BIM in the level set method, where BIM reduces the computation time and the latter provides good accuracy in recovery. In ad- dition, the use of the signed distance-based regularization helps avoid irregularities in the level set evolution. From a methodological point, the main contributions of the article are: 1) formula- tion of the non-adjoint based approach for a level set method; 2) introducing an itera- tive linearization approach inspired by BIM into the level set method; 3) integration of signed distance function into the cost function of microwave imaging; and 4) extension of our method to the multi-material case, where all level set functions are updated to- gether, so the required runs of the forward solver are independent of the number of the level set functions. These changes allow us to reduce the number of forward solver runs resulting in much reduced computation time. The following are inferred from the numerical analysis performed in Section VI. 1. The method can recover shapes, locations, sizes, and number of objects very ac- curately in a small number of forward solver runs even from a relatively small number of measurements. The majority of the errors in the recovery are at bound- aries. The reason may be that the grids during the synthesis and the inversion are different. 2. The method is robust against noise for picking up object’s locations and it can provide super resolution for identifying two closely placed objects. 3. In the comparative analysis, the reconstruction errors are lower than Bare-LS and as good as IMSMRA-LS with less computational complexity. 4. It performs well for high contrast and multiple materials. For multiple materials, the method can locate contrasts correctly. 5. The results using the multi-material approach on realistic phantoms are also rea- sonably good. Future work includes evaluating the performance of the method on experimental data and extending it to 3D. 67 CHAPTER 4 Matrix Norm Constrained Imaging 4.1 Introduction Inverse scattering problems arise in many practical applications and are formulated us- ing the volume integral equation. Recently, compressive sensing (CS) has been applied to address some of the challenges of inverse scattering problems and has shown im- proved performance [117] especially if the investigation domain is sparse or compress- ible in some transformed domain. In this framework, most of the state-of-the-art CS based inverse scattering solvers have been formulated using` 1 minimization of some quantity. The underlying formula- tion can be contrast source[66], Born approximation [1], or Born Iterative Method [88] but the objective is always the minimization of` 1 norm. However, given that the inverse scattering problem is nonlinear in nature, usually an underdetermined problem, and also suffers from having a limited number of independent measurements, there is still not a formulation, which can satisfy the CS requirements and hence solve the inverse scatter- ing problem uniquely. The CS formulation requires the system to be linear, unknowns to be compressible or sparse, and the system matrix to follow certain mathematical con- ditions. The CS theory is also extended to matrices, where instead of minimizing the` 1 norm, a matrix rank is minimized and has been found more useful in certain cases [118]. In this article, we formulate the cost function differently in a way that the unknowns are in a matrix form and we provide a logical explanation for minimizing the rank for the inverse scattering problem. We compare the formulation to other well-known matrix norm minimization, i.e. mix norm. A set of numerical results is reported to show their potential and limitations. 68 4.2 Inverse Problem Formulation We assume a z-independent, heterogeneous, isotropic, and non-magnetic two dimen- sional (2D) scatterer, immersed in a homogeneous background medium. The scattering problem for this geometry can be written as a linear mapping using volume integral equation with contrast source (product of the contrast and total field at a given point) as unknowns. ForM scattered field measurements in vectorb, andN unknown contrast source locations, the problem in a matrix form is written asb = Gz. If the contrast is zero, contrast source is also zero so the sparsity in the contrast domain is preserved. 4.3 Cost Functions and Implementation Since the contrast source is a function of incident field, forM s sources, the unknowns z are different for each M s . The unknowns can be concatenated in a long vector and the sparsity constraint can be put on the vector, but this formulation does not use the fact that eachM s vector is sparse. To use this fact, minimizing mix-norm (k:k 1;2 ) of the M s N matrix,Z, is a good strategy because it enforces` 1 norm constraint on columns and` 2 norm constraint on rows, thus enforcing sparsity in the solution ( column). We can formulate the cost function as, minkZk 1;2 s:tkGZBk 2;2 <= n ; (4.1) whereB := [b 1 ; b 2 ; :::; b Ms ], and n is the noise in the system. The cost function can be minimized similar to any` 1 minimization method. We follow the method in [79]. Sparsity in each M s vector can also be considered by localizing non-zero entries in the vector to only a few positions. Mathematically, this can be achieved by putting a joint constraint using matrix rank and element wise sparsity. The latter creates a sparse matrix, whereas the former controls the sparse matrix to be a low rank matrix, which im- poses localization. Rank of a matrix can be minimized using convex relaxation proposed in [119], which is to replace rank by the nuclear norm (sum of the singular values). So our second cost function is written as, min (kZk +kZk 1;1 )s:tkGZBk 2;2 <= n ; (4.2) where is a constant to balance the rank and the sparsity constraint. For the joint constraint, we adopt an accelerated proximal gradient based method described in [120]. 69 Here, essentially each solution is first projected on a low rank space using singular value decomposition and then sparsity is enforced by setting elements that are less than a certain value to zero. Once the contrast source is estimated, the permittivity and conductivity contrasts are estimated by first evaluating the volume integral equation to estimate the total field uniquely, and then estimating the contrast asx l = P Ms i=1 (E i;l Z i;l ) P Ms i=1 (E i;l E i;l ) ,l2 [1;N]. 4.4 Numerical Results The algorithm has been tested with 2-D synthetic data for high contrast and sparse ob- jects. The following geometry is considered: 36 receivers are placed uniformly at a measurement radius of 4, where is the free space wavelength. The incident waves are assumed to be plane waves at a wavelength of 15 cm. The investigation region is assumed to be within a square of length 2. The domain is discretized at 8 . Thus, the system has 288 measurements and since it is the contrast source formulation, the number of unknowns is 196M s . In the first set of experiments, 8 transmitters are considered. Two types of objects are considered- distributed point targets and clustered points at the center. The contrast source recovery for clustered points shown in Fig.1 demonstrates the low rank nature and sparsity of the unknown matrix. As can be seen, for both mix-norm recovery and joint sparse and low rank recovery, the matrix is very sparse but for the latter, the re- covery matches the truth very well. The corresponding contrast recovery is shown in Fig. 2. Because of the perfect recovery ofZ, for cluster points, the object is recovered almost perfectly even for the given contrast as high as 20. The same can not be said for distributed points. For mix-norm based recovery, the results are completely reversed. The behavior is quite surprising but as seen from Table 1, the joint recovery is still not converged for distributed points and mix-norm recovery has a minimum point different than the actual object but at a lower value. 70 Table 4.1: recovered the matrix rank and the mix-norm using two proposed methods parameters Matrix Rank Mix-norm Method Actual mix-norm joint Actual mix-norm joint Object recovery recovery recovery recovery clustered 4 5 4 31.5 25 31.4 distributed 4 4 6 40.3 40.1 619 large 18 18 24 358.5 99.9 1400 In the next set of experiments, the number of sources is increased to 24 and relatively larger object with contrast = 6 is considered. The recovery is shown in Fig. 2. Mix norm based recovery is able to identify the shape of the object a little but could not get the contrast correctly. Joint recovery performs reasonably well. -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -1 0 1 -1 0 1 -1 0 1 Figure 4.1: Contrast source (Z) for clustered points: From left actual matrix represen- tation, mix norm recovery, and joint sparse and low rank recovery. Top row shows recovery for the real part and bottom row for imaginary part. 71 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 2 4 6 0 2 4 6 0 2 4 6 Figure 4.2: Contrast recovery. Geometry configuration - top row: clustered points, center row : distributed points, bottom row : large object. From left actual object, mix norm recovery, and joint recovery. 4.5 Conclusion and Future Work We demonstrated the feasibility of a matrix norm minimization based formulation for inverse scattering. The results show that it has the potential to be applicable for high contrasts but further analysis is needed. Future work includes improving on the method for low-rank recovery and further incorporation of the matrix structure into the cost function. 72 3D Microwave Imaging 73 CHAPTER 5 Combined Shape and Contrast Constrained Imaging 5.1 Introduction Level set methods (LSMs) have long been popular in image segmentation and im- age registration problems because of their implicit and simple representation of shapes [94, 95, 96]. For the same reason, they have also been applied to broader imaging problems[121, 97, 122, 123], including microwave imaging. A more detailed review on level set methods for imaging is given in [45] and the challenges are described in [107]; however, a brief review of the major developments and challenges to the application of LSM towards microwave imaging problems is given below. Traditionally, the imaging problem has been formulated using an adjoint field rep- resentation and calculation[97], where the forward solver is run twice at each iteration. Typically, the number of iterations is on the order of hundreds or sometimes thousands. Therefore, the method is slow and was initially limited to 2D formulations. In [124] and [89], LSMs were used to recover PEC objects in 2D. Next, the method in [89] was validated experimentally [90], and later the same group also reported that TM waves demonstrate better convergence behavior than TE [92]. Many approaches have been suggested to address the large number of iterations required, including a parallel implementation on the Message Passing Interface (MPI) technique [91] and a multi-resolution or variable grid size approach [2]. The linear sampling method has been used as an initial guess to reduce number of iterations [100, 99]. We also recently proposed a fast level set method [125], where we formulated the problem as Born Iterative Method (BIM) type iterations to reduce the number of forward solver runs. 74 For many of the above articles, the constitutive parameters or contrasts are known a priori, but if these are also to be estimated, the same PDE formulation is used as in the case of shape estimation. Recently, [93] presented the modified level set method to estimate the contrast for 2D microwave imaging. The contrast and shape are itera- tively estimated by level set evolution equations with some regularizations. However, the method is still based on adjoint representation, which requires a large number of forward solver runs. Currently, there are only a handful of articles for LSMs applied to microwave imag- ing in 3D, likely due to the large number of forward solver runs required by the tradi- tional adjoint-based formulation. For 3D, a frequency-hopping based approach is more common given that the number of measurements is very small compared to the number of unknowns. The method was developed for PEC objects [126] and also for dielectric objects having different contrast values [127]. For dielectric objects, the contrasts are known a priori and a different level set function is used to represent each contrast value. To our knowledge, there are currently only two other studies of 3D LSMs for microwave imaging, both of which focused on breast imaging: [128] considered 3D synthetic ho- mogeneous breast phantom and recovered the shape of fibroglandular tissue immersed in fatty tissue, [98] considered heterogeneous breast phantoms and estimated both the shape and contrast in 80 iterations with a single run of the forward solver per iteration, but used the final outcome of the Distorted Born Iterative Method (DBIM) as an initial guess. In this article, we extend the fast level set method to 3D, and incorporate both shape and contrast estimation by formulating the imaging problem as two separate subprob- lems: shape estimation and contrast estimation. This decoupling allows more complex methods for the contrast estimation at a reduced computational cost, as the contrast sub- problem is solved only within the estimated object’s locations. The first subproblem, shape estimation, is carried out by the level set formulation. For the contrast estimation, a total variation (TV) regularized approach, similar to those used in compressed sensing, is considered as a cost function. By bridging different inversion regularizations together and integrating level set methods within a BIM inversion scheme, we are able to achieve equal or better recovery at greatly reduced computational complexity. We demonstrate the method here on synthetic data generated with an experiment- ready forward model. Compared with the simple dipole sources and scattered field used in [98], we use a conformal FDTD technique [129] to model a cubic imaging cavity 75 with 36 transmitting and 36 receiving patch antennas at a single frequency and link the object property to the S-parameter measurements from the vector network analyzer with a full numerical waveport vector Green’s function [130]. We achieve acceptable shape and contrast recovery for less than 10 runs of the forward solver, a factor of at least 10 reduction from published methods. The paper is organized as follows. In Section II, the problem formulation is de- tailed. Section III describes our inversion approach for microwave imaging and derives key components of our approach. Implementation details are described in Section IV , and evaluation of the algorithm on simulated data is performed in Section V . Finally, conclusions and future directions are discussed in Section VI. 5.2 Theoretical Background In microwave imaging, the inverse scattering problem in a heterogeneous, isotropic, nonmagnetic medium is governed by the following electric volume integral equation [131], E s ( r) =k 2 b Z D G( r; r 0 )( r 0 ) E( r 0 )dv 0 ; (5.1) whereD is the imaging region, r 0 is the imaging domain position vector, and r is the measurement point position vector. The quantity E s ( r) is the scattered electric field which is the difference between the total field, E( r), and the incident field, E i ( r). The quantity G( r; r 0 ) is the Dyadic Green’s function, k b is the wavenumber of the back- ground medium with lossless permittivity, and( r 0 ) is the dielectric contrast in terms of the permittivity and conductivity contrast, defined as ( r 0 ) = [ r ( r 0 )= b 1]j[(( r 0 ) b )=( b 0 !)]; (5.2) where b is the lossless relative permittivity of a background medium, r denotes the relative permittivity, and b is the background conductivity. In inverse scattering experiment, it is often difficult to measure the scattered elec- tric field E s , instead, we use the vector network analyzer to measure the scattered S- parameterS s , which is a scalar instead of a vector. In addition, with the presence of the imaging cavity, the measurement background is often inhomogeneous, which does not have a Green’s function with an analytical form. Therefore, we developed a full numer- ical waveport vector Green’s function (WVGF) which directly links the object contrast 76 and the vector total electric field to the scattered S-parameters measurements [132]. The WVGF takes the form, Rl G( r; r 0 ) = j E inc ( r 0 ) p Z n ! 0 RR An e t n ( r) h t n ( r)d A p Z m (5.3) wherem andn represent the receiving and the transmitting waveport, respectively. The quantity E inc is the incident field in the imaging domain produced by exciting the re- ceiving portm, which is computed with a 3D conformal finite difference time domain method (CFDTD) described in [129]. The quantity Z m and Z n are the impedance of the receiving and transmitting port, respectively. The electric and magnetic mode field template of transmitting portn is represented by e t n and h t n . which are obtained with a 2D conformal finite difference frequency domain (CFDFD) method described in [129]. With the full numerical waveport vector Green’s function G( r; r 0 ), the S-parameter vol- ume integral equation can be written as S s ( r) =k 2 b Z D G( r; r 0 )( r 0 ) E( r 0 )dv 0 : (5.4) The imaging problem deals with determining the dielectric contrast, ( r 0 ), of an unknown medium, given the measurement of scattered S-parameters, S s ( r). Due to the presence of the electric field within the integral, which itself is a function of the dielectric contrast, Eq. (5.4) is non-linear. We use the linear operator Lf:g to repre- sentk 2 b R D G( r; r 0 ):dv 0 . Using this notation, the scattered S-parameter volume integral equation is written asS s =Lf; Eg. To solve the nonlinear inverse scattering problem, we use a linearized approach such as BIM, which iteratively solves the total electric field and dielectric contrast in an al- ternating manner during the inversion process. To solve the total electric field, E( r 0 ) inside the imaging domain, for a given or reconstructed dielectric contrast, we use con- formal FDTD method [129] which numerically excites and computes the fields within the imaging domain. This step is commonly known as the forward solution. 5.3 Inversion The imaging problem is solved by minimizing an objective function,k S s Lf; Egk, by estimating the dielectric contrast and hence directly or indirectly also estimating 77 the total electric field E. For N s transmitters and N r receivers, let the scattered S- parameters be represented as a vector,y, of length (N s N r ) 1. We follow the fast level set method [125] for the overall minimizing strategy. In the fast level set method, the objective function is minimized in a Born-iterative style by a level set function. The objective function is approximated by a linearized cost function, defined asky Ak 2 , to recover the shape of the object, while assuming that the domain has known contrast values. Here,A is the approximatedL estimated by the calculated total electric fields. In this article, we do not assume the contrast values are known a priori. We propose the following approach to solve for both contrast and shape in the fast level set framework. We estimate the dielectric contrast at every point in the domain by recovering the shape of objects and then recovering the contrast for each shape. Even though the contrast can be recovered simultaneously with the shape, we separate them to increase computational efficiency and to accommodate the introduction of constraints in the cost function for contrast estimation. The estimation details are explained in the following two subsections. 5.3.1 Shape estimation In a level set formulation, the contrast is parameterized using a higher dimensional con- tinuous function, called a level set function. We formulate the level set function for background versus object. The object may have more than one value or type of contrast and may be disconnected. The contrast is represented by as (r 0 ) = i U((r 0 )l) + e (1U((r 0 )l)); (5.5) where i and e are estimated average contrast for the object and background, respec- tively, andU(r 0 ) is the standard Heaviside function. The region occupied by the object is represented by>l and the background is represented byl, so =l traces the boundary of the object. For simplicity, we considerl = 0. The minimization of any cost function,J, is achieved in the form of the evolution of the level set function as @ @t + @J @ = 0;r = 0 onD; (5.6) where D is the boundary of the domain D and J in our case is J() = ky 78 A()k 2 , called residual error. Regularization To diminish any irregularities during the level set function evolution and hence avoid the need to reinitialize the level set function, we add distance-based regularization[115]. In addition, we expect the objects will have smooth boundaries, so we add mean-curvature- based regularization. The mean curvature of is formally defined asr: r jrj . We keep the mean-curvature-regularization active only at the object boundaries, so minimizing the regularizing term decreases the curvature of the boundaries and hence promotes the edge with less local variations. Thus, the objective function of the shape estimation is to minimize the residual error under the constraints from two regularizations. These all can be written in the form of the level set evolution equation as Rl @ @t = 1 RefA + R()g( i e )() (5.7) +r: p 0 (jrj) jrj r + 2 ()r: r jrj ; (5.8) whereRe indicates real part of the argument,R is the residue (R =Ay), ‘ + ’ indi- cates adjoint operation,() is the Dirac delta function,p is a potential energy chosen such that the level set function stays as the signed distance function, and 1 ; 2 ; are the parameters. For further details, please refer to [125]. 5.3.2 Contrast estimation In the previous step, we obtain an estimate of the position and shape of the objects. We then use this information in the contrast estimation step. The multi-view scattering equation can be written as y =A =A b b +A o o ; (5.9) where the subscript b corresponds to the background and o corresponds to the object. The notationA i ; i2fb;og represents a submatrix, generated by selecting the columns corresponds to either the background (i = b) or the object (i = o). The notation for 79 represents selecting corresponding indexes. For a homogeneous background, the back- ground contrast is always zero (see (5.2)) so the cost function for contrast estimation becomes, minkyA o o k 2 +R c (z); (5.10) whereR c is regularization as the subproblem may be still underdetermined depending on the size of the objects and number of measurements and it may have many solutions. It is clear at this point thatR c can be any regularization terms and can have more than one term. To demonstrate the feasibility, we use total variation based regularization, commonly used in compressed sensing literature, with the box constraints. Equivalently, it can be written as minkyA o o k 2 +TV ( o ) s:t: l o u; (5.11) where l and u are lower bound and upper bound vectors respectively and TV is the anisotropic total variation. Writing explicitly, the cost function is Rl min o kyA o o k 2 + jr x o j +jr y o j +jr z o j (5.12) s:t: l o u: (5.13) This cost function is difficult to solve because the gradient operator is not invertible. We apply the Split-Bregman framework [133] to solve the above problem because the framework converges quickly and regularization parameter remains constant throughout iterations. In Split-Bregman iteration, the ` 1 and ` 2 portions of the cost function are decoupled. The box constraint terms, which are absent in [133], are here part of the` 2 portion. In particular, the TV norm constraint is separated by an introduction of auxiliary variables and splitting the problem into solving two subproblems- one for unknowns and the other for the auxiliary variables. The equivalent problem: Rl min o;dx;dy;dz kyA o o k 2 +jd x j +jd y j +jd z j (5.14) +kd x r x o k 2 +kd y r y o k 2 (5.15) +kd z r z o k 2 s:t: l o u; (5.16) whered x ; d y ; d z are auxiliary variables and is a regularization parameter. Applying 80 Bregman iteration, we get Rl min o;dx;dy;dz kyA o o k 2 +jd x j +jd y j +jd z j (5.17) +kd x r x o b k x k 2 +kd y r y o b k y k 2 (5.18) +kd z r z o b k z k 2 s:t: l o u; (5.19) whereb k p ; p2fx;y;zg is chosen by Bregman-iteration asb k+1 p =b k p + (r p o d k+1 p ) [133] andk is an iteration number. Eq. (5.19) can now be split into four subproblems - minimizing with respect to o ; d x ; d y ; andd z respectively. For the last three subproblems, the cost function is sim- ilar so we only need to derive the solution once. 5.3.2.1 o subproblem The cost function for o subproblem is written as min o kyA o o k 2 + X p2fx;y;zg kd p r p b k p k 2 s:t: l o u (5.20) Rearranging terms, the equivalent cost function can be written as, min o kyB o k 2 s:t: l o u; (5.21) where we use notation :B := h A o ; p r x ; p r y ; p r z i T andy := h y; p (d x b k x ); p (d y b k y ); p (d z b k z ) i T . The equivalent problem is a standard optimization problem. We consider the projected gradient method to estimate o at every step. The estimated contrast at stepn is written as o n+1 = X ( o n l g l ); (5.22) where is the projection on the box andg l is the gradient. We use the method described in [134] to estimate optimal l by Barzilai-Borwein step along with backtracking line- search. The projection step is carried out by deriving and solving Karush Kuhn Tucker (KKT) conditions similar to our previous work[88]. It can be summarized by 81 stating that the value of the cost function is optimal when any element of o outside the box is projected on the closest surface or edge of the box. For unknowns in the complex domain, applying the box constraints on real and imaginary parts separately give the optimal value of the cost function. 5.3.2.2 d p subproblem Eachd p ; p2fx;y;zg subproblem can be written as min dp jd p j +kd p r p o b k p k 2 : (5.23) This is a standard ` 1 minimization problem. Since the unknown variable d p is sepa- rate elementwise, each element can be obtained by the shrinkage operator or the soft thresholding operator: d k+1 p = d k p jd k p j max(r p o +b k p ; 0): (5.24) The first ratio in the right hand side captures the sign or phase of the unknowns and the max operation estimates the magnitude. 5.3.3 Algorithm The inversion has three main steps: field estimation, shape estimation, and contrast es- timation. The overall approach we use for solving this optimization problem is outlined in Algorithm 1. 5.4 Implementation Details 5.4.1 Problem setup In the numerical experiment, we designed a microwave imaging chamber to simulate the measurement of S-parameters around the object. The imaging chamber has a size of 25 15 25 cm 3 . As shown in Fig. 5.1, on the 4 vertical sides of the imaging chamber, there are 72 rectangular patch antennas with a dimension of 8:4 33:6 mm 2 working at 810 MHz. The imaging chamber is filled with a lossless coupling fluid having a 82 Algorithm 2 Microwave imaging algorithm procedure INVERSION ALGORITHM(JOINT SHAPE AND CON- TRAST RECOVERY) initialization:, 0 , i whilekyAk 2 > do estimate the fields by conformal FDTD for n =1 to N ork n n1 k>tol do n+1 = n + tL( n ) end for getA o and o form = 1 toM do fork = 1 toK do k o = argminkyB o k 2 ; s:t:l o u d k+1 p = d k p jd k p j max(r p k o +b k p ; 0) end for b m+1 p =b m p + (r p k o d k+1 p ) end for i = mean( k o ) end while end procedure 83 permittivity of 10. Of the 72 antennas, 36 are as transmitters and the remaining 36 are used as receivers. The imaging region,D, has the size of 19:2 16 19:2 cm 3 . The 3D CFDTD method [129] is utilized to model the imaging chamber and compute the electric and magnetic field within. Figure 5.1: Imaging chamber with 36 patch antennas. 5.4.2 Numerical approximation For solving the evolution of the level set function, we use the central difference scheme to approximate all of the spatial partial derivatives, whereas we use the forward differ- ence scheme to approximate the temporal partial derivative. The approximation of (5.8) by the above difference scheme can be simply written as i+1 i t =L( i ); (5.25) whereL( i ) is the approximation of the right-hand side of (5.8) and t is an artificial time that represents the level-set evolution iterations. We choose the parameters such that the Courant-Friedrichs-Lewy(CFL) condition is satisfied [115]. 5.4.3 Selecting parameters The parameters for shape estimation and contrast estimation are chosen independently. For the level set method, we consider t = 1 and choose the remaining parameters 84 accordingly. First, we select the value of 1 by initially estimating its range from the gradient of the residual error at the first BIM step and then a value is chosen by trial and error that produces a low residual error. To reduce number of the level-set iterations, we restrict the level set function to lie within [2; 2]. This bound is not necessary, but it helps in ‘generating’ and ‘clearing’ objects quickly. We select such thatt < 1=4 and <j 1 j=2. The latter condition is to make sure that the regularization term does not dominate over the residual error during the minimization process. The value for 2 is selected in a similar way. For contrast estimation, the value of is set such that the contribution of the associated term is less than 50% of the contribution from the residual error. The value is set such that the ratio, = is lower than the minimum spatial contrast difference we are interested in detecting so that the contrast difference is preserved in (5.24). 5.4.4 Initialization As pointed out in the Algorithm I, three quantities need initialization: , which is the unknown contrast of the domainD in a vectorized form; 0 , which is an initialization of the level set function; and i , which is the average contrast value of the object used during the shape estimation. For , we initialize as zeros, which is equivalent to starting the inversion process with the Born approximation. Since the level set function needs initialization at each BIM iteration, we initialize 0 as zeros everywhere (all background) in the first BIM run, and from each BIM iteration onward, we can either start from the converged level set function from the last step or we can reinitialize as zeros again in order to avoid getting stuck in local minima [125]. We choose the latter initialization. For the value of i , we pick few sample values uniformly within the given bound (l u), run the level set routine, and select the i which has produced the lowest residual error. This step is only for the initialization step. For all successive iterations, the mean value of the contrast estimated by the contrast step is used as i . Even though we do contrast sam- pling within the bound during shape estimation of the first BIM iteration, the additional computational burden is less than half of the computational burden of a single forward solver run. 85 5.5 Analysis and Results We perform several experiments to validate the presented method. In the reported nu- merical results, the scattered field data are obtained by using (5.4). The average voxel size in all directions is 4 mm, which gives an image size of 484048 voxels. The cur- rent setup is clearly underdetermined (92; 160 unknowns against 1296 measurements). Each image is formed with 6 iterations of BIM. An iteration consists of one run of the forward solver, up to 75 iterations for shape estimation, and M = 2 and K = 7 for contrast estimation (see Algorithm I). The computation time to run a single forward solver is approximately 30 mins on an Nvidia K20 GPU and the average times to run the shape estimation and the contrast estimation are 25 secs and 40 secs respectively for code written in Matlab on a Linux desktop with a Dual Xeon processor and 24GB RAM. The level set evolution equation’s parameters are set as t = 1, = 10 4 , 2 = 10 4 , and 1 = 50=kyk, wherekyk is on the order of ten thousands. The contrast estimation parameters are set as = 0:005, = 6 , andA andy are normalized bykyk. We set a uniform bound for object’s voxels,l = 0 andu = 5 for all the following cases unless otherwise specified. Example 1 : cube : Here, the object under test is a homogeneous cube located at the center. The cube is of 6:8 cm length and has relative permittivity of 15. The reconstructed cube is shown in the right column of Fig. 5.2. The max projection view 5.2(f) shows that the object’s shape and contrast are recovered accurately. 86 Figure 5.2: Cube. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view at y = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. 87 Fig. 5.3 shows the cross section (vertical) view at x = 0 plane (see Fig. 5.1 for axis reference) andy = 0 plane of the actual and reconstructed permittivity using joint shape and contrast recovery. The proposed method estimates the contrast very close to the actual contrast and localizes the object correctly. -10 -5 0 5 10 z (cm) 10 11 12 13 14 15 relative permittivity Figure 5.3: Vertical cross section view from the center of the domain; the plot of relative permittivity versus points on z-axis. Black dotted line corresponds to actual profile and blue line corresponds to joint recovery. The convergence of the algorithm is shown by plotting the residual error vs iteration number. The residual error is normalized bykyk to produce a relative error term. The relative error at every level set iteration is shown in Fig. 5.4. The larger jumps in the figure correspond to the update in the total field and henceA. The plot shows that error is decreasing and at the later iterations the error converges to the similar value for different As. 88 0 100 200 300 level set iteration # -60 -40 -20 0 relative error (dB) Figure 5.4: Convergence plot. Relative error is plotted from the first iteration to the last iteration of the method in terms of the error at each level set iteration. To quantify our observations, we use the Intersection over Union (IoU) performance metric for shape estimation and relative error metric to jointly characterize shape plus contrast estimation. IoU, which is defined as the ratio of the volume of overlap between recovery and original to the volume of the union of the two. In the perfect shape re- covery, IoU is 1. The second metric, relative error, which also considers the contrast, is defined as, = 1 N x N y N z v u u t Nx X i=1 Ny X j=1 Nz X k=1 rec r (i;j;k) ori r (i;j;k) rec r (i;j;k) 2 100%; (5.26) where rec r and ori r are complex permittivity of a recovered object and an original object respectively. The calculated values for all cases are reported in Table 5.1. As can be seen, the quantitative shape error for the cube case is 17% and the relative error is 0:01%. Example 2 : hollow sphere : The object consists of two concentric spheres located at the center of the domain. The outer sphere has a 4:8 cm radius and a relative per- mittivity of 15, whereas the inner sphere has a 2:4 cm radius and has no contrast with 89 Table 5.1: Quantitative measurement of the reconstructed images. Object r IoU Error () (%) Cube 15 0 0:83 0:012 Hollow Sphere 15 0 0:81 0:016 Two Spheres 15; 20 0 0:79 0:013 Complex Sphere 15 0:35; 0:5 0:82 0:013 Two Objects 20; 22 0 0:79 0:026 Three Objects 17; 17 20 0 0:72 0:02 respect to the background (Fig. 5.5). In terms of wavelength, the inside sphere has a radius approximately equal to a fifth of a wavelength in the background medium and the width of the outer sphere is on the same order. The final object estimation is shown in the right column of Fig. 5.5. As can be seen, there is a good agreement with the ac- tual configuration. Corresponding metrics are reported in Table 5.1. The cross-sectional view of the recovery along z axis is shown in Fig. 5.7(a). The surfaces of the spheres are detected at the correct positions with accurate estimation of the contrast. Example 3 : two spheres : Here, we consider the same shape as in the previous example except the object is not at the center and has different permittivity. The center of the object is at (1:6;0:8;0:8) cm and the relative permittivity is 15 for the outer sphere and 20 for the inner sphere. The original object and the recovered object are shown in Fig. 5.6. As can be seen, the method is able to detect the presence of the inner sphere. Corresponding metrics are reported in Table 5.1 with similar values as those for previous cases. The cross-sectional view of the recovery along z axis is shown in Fig. 5.7(b). 90 Figure 5.5: Hollow sphere. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view at y = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. 91 Figure 5.6: Two spheres. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view aty = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. 92 Example 4 : complex sphere: Next, we again consider two concentric spheres at an off-center position but their permittivity is a complex number. The center of the object is at (2:8;1:2;1:6) cm. The inner sphere has relative permittivity of 15 and conduc- tivity of 0:35S=m and the outer sphere has relative permittivity of 15 and conductivity of 0:5S=m. The bound is the same for real part of the contrast as for the previous cases. For the imaginary part, we set the maximum absolute value of the contrast to be 3, re- sulting in bound of [03i; 50i]. The recovered object after 6 BIM iterations is shown in Fig. 5.8 for the relative permittivity and in Fig. 5.9 for the conductivity. The cross- sectional views are shown in Fig. 5.7(c) and (d), respectively. From Table 5.1, the error, , and IoU have values similar to previous cases. Example 5 : two objects : In this example, there are two objects: a cube and a sphere located at (2:4; 1:6; 2:4) cm and (3:6;3:2;3:6) cm, respectively, as can be seen from the projection view in Fig. 5.10(e). The cube has the relative permittivity of 22 and has a length of 5:2 cm. The sphere has relative permittivity of 20 and a 3:2 cm radius. The right column of Fig.5.10 shows the reconstructed permittivity distribution. The method is able to recover the shape of the objects really well and the estimated contrasts are very close to the actual contrast (see Fig. 5.7(e)). Example 6 : three objects : A more complicated case is considered here, where there are three objects: a sphere, a cube, and a prism. Their sizes are 6:4 cm in diameter, 4:4 cm length, and 6 cm side length and 3:2 cm height, respectively. Their relative per- mittivity is 17, 17, and 20, respectively. The estimated permittivity distribution is shown in Fig. 5.11. As can be seen from Fig. 5.11 and Fig. 5.7(f), there is good agreement with the actual configuration. Even though the estimated contrast is lower than the ac- tual contrast, the method is able to retrieve their relative contrast values correctly. 93 -10 -5 0 5 10 z (cm) 10 12 14 16 relative permittivity (a) -10 -5 0 5 10 z (cm) 10 12 14 16 18 20 22 relative permittivity (b) -10 -5 0 5 10 z (cm) 0 0.1 0.2 0.3 0.4 0.5 relative permittivity (c) -10 -5 0 5 10 z (cm) 0 0.1 0.2 0.3 0.4 0.5 conductivity (d) -10 -5 0 5 10 z (cm) 10 12 14 16 18 20 22 relative permittivity (e) -10 -5 0 5 10 z (cm) 10 12 14 16 18 20 relative permittivity (f) Figure 5.7: Cross sectional view along z axis for the permittivity of (a) hollow sphere, (b) two spheres, (c) complex sphere, (e) two objects, and (f) three objects. (d) shows the cross section view along z axis for the conductivity of the complex sphere. 94 Figure 5.8: Complex Sphere (real part). left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view aty = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. 95 Figure 5.9: Complex Sphere (imaginary part). left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view aty = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates conductivity. 96 Figure 5.10: Two Objects. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view at y = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. 97 Figure 5.11: Three Objects. left column : actual geometry, right column: joint contrast and shape recovery. The top row shows a 2D view at y = 0, the center row shows orthogonal slice view and the bottom row shows the max projection along the y-axis. The color bar indicates relative permittivity. 98 5.5.1 Comparative assessment To do a comparative assessment, we consider two scenarios. (a) How much improve- ment can the joint shape and contrast estimation have over a traditional method? (b) how much improvement can the more complicated contrast estimation provide? For the first question, we consider a TV regularized method (` 1 minimization) with the box constraints in the BIM framework. To make a fair comparison, we pick the same split- Bregman method as a way to solve the TV regularized cost function. So effectively, the traditional method is run on the entire grid of unknowns rather than restricted grid of the object voxels. The parameters of the method are tuned in a similar way as we did for the contrast estimation. The method is tested on all examples. The TV method can recover simple objects very well but fails to recover more complicated objects. To demonstrate the performance degradation, the object used in example 3, two spheres, is considered here. The recovery using this method is shown in the last row of Fig. 5.12. The method is able to localize the object correctly but not the shape and contrast. The background is not perfectly zero as in the case of the level set based method but we consider the contrast value lower than 10% of the lowest object contrast as the background. The re- covery by the proposed joint shape and contrast recovery is also shown in the first row of Fig. 5.12 for ease of the comparison. The joint method has clearly better contrast and shape recovery. It localizes the object very well. For the second question, we consider is zero and solve the relevant problem in the same manner. This can be considered as running the standard level set method with the explicit bound constraint for the contrast estimation. However, in our case, the param- eter used in the cost function of the contrast estimation is adaptive. The recovery using only bound constraint on the same example, two spheres, is shown in the center row of the Fig. 5.12, which clearly depicts the importance of the additional constraint. For a homogeneous contrast inside an object, this method does not preserve the homogeneity of the contrast. Quantitatively, the error, , is 0:013%, 0:016%, 0:017% for the proposed method, the solely TV-based method, and the level set method with only bound constraint for con- trast estimation, respectively. The IoU values are 0:79, 0:75, and 0:81, respectively. To show the improvement more clearly in raw numbers, we also compute the distance from the original object in an ` 2 sense, where the euclidean distance between the relative permittivity of the recovered object and the relative permittivity of the actual object is computed. The distance for each method is 16, 22:1, and 23:9, respectively. 99 Figure 5.12: Comparative analysis. Top row : recovery by the proposed joint method, center row : recovery by the level set method with only bound constraint for the contrast estimation, bottom row : recovery by a traditional method. The left column shows the a 2D view aty = 0, the center column shows orthogonal slice view and the right column shows the max projection along the y-axis. The color bar indicates relative permittivity. 5.6 Conclusion A level set function based contrast and shape estimation method for 3D microwave imaging has been presented. Our approach is based on estimating the shape by the fast level set method and the contrast by the total variation regularized cost function. The 100 integrated approach reduces the computation time and provides good accuracy in recov- ery. From a methodological perspective, the main contributions of the article are: 1) formulation of an end to end method for contrast and shape estimation; 2) introducing a modular framework, where based on problem requirements, the cost function can be de- fined for better performance; 3) demonstration of the use of the TV based regularization in the level set framework; 4) modification of the Split-Bregman method to include box constraints. These changes allow us to recover contrast in the fast level set framework without compromising the recovery or increasing the number of forward solver runs; and 5), implementing the inversion algorithm with an experiment-ready forward model that enables the direct usage of scattered S-parameter data from a VNA-based inverse scattering measurement system, such as our experimental system in [135]. The following are inferred from the numerical analysis performed in Section V . 1. The method can recover shapes, locations, sizes, and a number of objects very accurately in a small number of forward solver runs even from a relatively small number of measurements. 2. The method estimates the contrast very well and have relative error less than 0:03%. 3. In the comparative analysis, the reconstruction errors are lower than both- a tradi- tional method and a level set method. 4. It can handle different contrast very well. Future work includes extending the method to multi-phase level set method and evaluating the performance on numerical biological phantoms. 101 Compressive Sensing For Microwave Imaging 102 CHAPTER 6 Theoretical Derivation 6.1 Introduction Inverse scattering problems arise in many practical applications whereby the distribu- tion of key physical parameters are estimated based on the measured samples of the scattered field. Medical imaging, diffraction tomography, buried object detection, non- destructive testing, and industrial imaging are few examples where an inverse problem needs to be solved. The electromagnetic inverse problem is usually formulated using the volume integral equation. The equation is non-linear in general, since the electric field is a function of medium’s properties and the internal fields of the object appear in the integral. One of the approaches is to combine contrast and field in a single variable known as a contrast source. The formulation based on this notation is referred as source type integral equation (STIE). 6.2 RIPless Probabilistic Compressive sensing Compressive sensing theory provides a theoretical base to show that in certain linear systems, randomly selectedm measurements can recovers-sparsen dimensional signal by optimizing the cost function and selecting the solution, which has the minimuml1 norm. The required number of measurements for a given linear system can be estimated by a property called the restricted isometry property (RIP) [76]. However, computation- ally the RIP can be extremely difficult to prove and for real-world signals, which are either corrupt by noise or not exactly sparse but compressible, the RIP can not be apply directly. The later problem is addressed in [136] and proves that nearly sparse signal can be recovered faithfully fromm (wherem is smaller thann) random measurements even 103 in noisy setting. We will elaborate the faithful recovery in the following paragraphs. Letx be as near sparse vector fixed but unknown, and we have the values ofm mea- surements in a vectory, which is related tox by a linear operatorA. We are interested inx, which satisfy the following equation: minkxk 1 s:tAx =y: (6.1) The proof of the faithful recovery uses the concept of a dual certificate and proves in a probabilistic manner. In the dual certificate framework, one formulates the dual prob- lem of the given primal problem, (6.1), and find a dual feasible point (dual variable), which certifies the optimality of the primal variable. It is shown that if one can exhibit a dual certificate, the recovery ofx is exact with high probability [137]. A relax version of this first introduced in [138] for matrices and later applied to vectors in [136]. This approach does not establish universal recovery results but has been proved for the fixed but otherwise arbitrary sparsitys andA obeying certain properties, namely isotropy and incoherence. We can say a measurement population/distribution (for example Fourier) obeys the isotropy property, if (E (aa) = I), where a is the column of A and is in- dependently sampled for each column from the distribution. The incoherence property requires thatA can not be sparse. The rows ofA have to be global. In calculation, we take the coherence parameter F as maxha; i 2 , where i is a domain, in whichx is sparse. The more incoherent the measurements are i.e, smaller F , the fewer samples are needed for accurate recovery. In [136], the isometry property is relaxed and is needed only locally - only for the columns of A, where x is nonzero. Let T of s cardinality be the set of locations of nonzero withinx. We need the following property with high probability: kA T A T Ik 2 ; (6.2) whereI is the identity matrix,A T denotes is a submatrix ofA havings columns and is constant and should be less than 1=2 for noiseless measurements and less than 1=4 for noisy measurements. Similarly, incoherence is defined for off-support columns as none of the columns ofA outside the support ofx should be well approximated by any vector sharing the support ofx. Mathematically, max i2T ckA T A fig k 2 1: (6.3) 104 If the given matrix, A, follows the above mentioned two properties and there exists v2R N in the row space ofA following the properties, kv T sgn(x T )k 2 1 4 and kv T ck 1 1 4 ; (6.4) then using vector algebra, one can show thatx is the unique minimizer to (6.1). Next, it is shown that such a dual variable, v, can exist with higher probability. If M satisfies the constraint: MC s F s log(N); (6.5) then with probability at least 1e w , for any scalarw 1, the constraints onA hold and a vector v with the properties required exists. Hence, x is the unique minimizer. The proof is using the Matrix Bernstein and Vector Bernstein inequalities in the golfing scheme. The proof can be found in [136]. The isotropy property has been relaxed in [139] to anisotropic measurement dis- tribution. The distribution has to be complete in that sense thatE (aa) is invertible. Under this condition and following the above strategy, the authors in [139] are able to show the exact recovery for anisotropic measurements. This formulation generalize the theoretical applicability of sparse recovery to any invertible linear system. There are also other anisotropic compressive sensing articles but they are based on RIP based for- mulation, which provides a higher sampling rate (more measurements). For this RIPless anisotropic formulation, the condition on the number of measurements for successful re- covery isMC s s F s log(N), where s is the condition number ofE (a i a i ); i2T . This formulation has a simple condition onM, similar to isotropic case. The anisotropic effect has been captured by a single quantity, condition number, which provides logical extension to isotropic case, however they are still some open questions. 1. The constant in the condition equation is very high, greater than 18000, which makes its applicability only to a larger system. As the authors mentioned in [139] that the bounds are not optimal, is it possible to tighten up the constant?. 2. There is an underlying assumption that as the number of measurements increases, the information content in the measurements increases. For anisotropic cases, this is not always true. For example, in inverse scattering experiments, the scattered 105 fields are sampled at receiver positions. We expect that we need many more re- ceivers to represent the scattered field correctly as the fields can vary a lot locally but it has been shown that there is an upper bound on the number of required receivers [140]. Beyond that number, the measurements do not add any addition information. In this scenario, it is not the best to relate the number of measure- ments to the sparsity. Is it possible to formulate the condition differently? In this article, we will focus on answering the second question. The condition, we get may not have sharp bound as with most of the compressive sensing literature, but it can help in the development of the generalized theory for compressive sensing. 6.3 Generalized RIPless CS We start by understanding the isotropy condition for anisotropic measurements. Let U V 0 be Singular Value Decomposition of the matrix A, then A A can be decom- posed asA A = V U 0 U V 0 = V 2 V 0 . Putting the decomposition in the isotropy condition, we getkIA Ak 2 =kV V 0 V 2 V 0 k 2 =kV (I 2 )V 0 k 2 =k(I 2 )k 2 . This metric, the largest deviation of a singular value from 1, shows that the isotropy con- dition estimates a deviation of the given matrix from an orthogonal matrix. The same conclusion can be made for the local isometry condition. However, since the local isom- etry condition is defined for a submatrix ofA, the chances of its being satisfied are high even for ill-posed problems. The division by p m on both sides of (6.1) are implied in the presented notation. Even then, the singular values of A T may not be close to 1 for an ill-posed problem. They can be brought closer by different ways because (6.1) is invariant under scaling. To rescale, one of the ways is to multiplyA T by a factor such that all its singular values are closer to 1. Or equivalently, we can introduce a constantc 1 in the local isometry condition such that singular values are closer toc. We assume that (6.1) is normalized such that singular values for anyA T is greater than 1. For the setup described above, we will havekcIA T A T k condition for local isometry. is considered to be 0:5 in noiseless recovery in [136] but here we keep it less than 1. In singular values, the condition is written askcIA T A T k =k(cI 2 T )k. It is clear that for a perfectly orthogonal system, settingc = 1 can give the same condition as it is in [136]. Now, we have to prove that solution of the optimization problem (6.1) even for the 106 bounded number of measurements equals the unknown vector x. We use dual certificate approach as used before to first show the exact recovery and then show the existence of dual certificate. 6.3.1 Some useful derivations Since the additional measurement can not bring new information, we can not useM to relate it with the sparsity and the number of unknowns. We have to use available infor- mation in the matrixA. The matrix norms,kAk inftt andkAk 1 operate at the individual element level rather than using all elements ofA as it is the case in the standardkAk 2 . We will explore these norms in the derivation. For the reference, the norm definitions are as follows. For a matrixA2C MN , kAk 1 = max 1pM N X q=1 jA pq j kAk 1 = max 1qN M X p=1 jA pq j: (6.6) We will also need bounds fork(A T A T ) 1 k 2 , which can be derived as follows: It is assumed that (A T A T ) 1 exists. We know thatkak 2 kbk 2 kabk 2 . So given k(cIA T A T )k 2 , we get kcIk 2 kA T A T k 2 : ckA T A T k 2 : (6.7) Using the norm relation between a matrix and its inverse, we getk(A T A T ) 1 k 2 1 c , provided thatc>. 6.3.2 Main results Now, we will prove that` 1 minimizer of (6.1) is the unique solution ofAx =y. Lemma 1(inexact duality) Letx2C N be assparse vector, letT =supp(x). Assume thatk(A T A T ) 1 k 2 C A and max i2T ckA T A fig k ` 2 C At . Suppose there exists v2 C n in the row space of A obeyingkv T c sgn(x T )k ` 2 1=4 andkv T ck `1 C v such thatc>C v + C A C At 4 , then 107 x is the unique` 1 minimizer to (6.1). The proof is highly borrowed from [136]. The proof is by contradiction. Let’s assume that ^ x = x +h be the solution to (6.1), which states thatk^ xk 1 kxk 1 and Ah = 0 since bothx and ^ x are feasible. We will show that the null spaceh is trivial. k^ xk 1 can be expanded in the setsT andT c as, ck^ xk 1 =ckx T +h T k 1 +ckh T ck 1 ckx T k 1 +hcsgn(x T );h T i +ckh T ck 1 : (6.8) Letv =A w be a dual vector, we have hcsgn(x T );h T i =hcsgn(x T )v T ;h T i +hv T ;h T i =hcsgn(x T )v T ;h T ihv T c;h T ci; (6.9) where we used thathv;hi = hw;Ahi = 0 andhv T ;h T i = hv;hihv T c;h T ci = hv T c;h T ci. Applying the Cauchy-Schwarz inequality and then using the properties ofv, we get jhcsgn(x T );h T ijjhcsgn(x T )v T ;h T ij +jhv T c;h T cij kcsgn(x T )v T k 2 kh T k 2 +kv T ck 1 kh T ck 1 1 4 kh T k 2 +C v kh T ck 1 : (6.10) Plugging the lower bound of the above relation in (6.8), we get ck^ xk 1 ckxk 1 1 4 kh T k 2 + (cC v )kh T ck 1 : (6.11) Next, we will find the bound forkh T k with respect tokh T ck 1 . Since (A T A T ) has full rank,h T = (A T A T ) 1 A T A T h T =(A T A T ) 1 A T A T ch T c, where we use the fact that Ah = 0 for the last step. Using the bound forA T , we getkh T k 2 C A kA T A T ch T ck 2 . Next, kA T A T ch T ck 2 X i2T c kA T A fig k 2 jh i j max i2T ckA T A fig k 2 kh T ck 1 C At kh T ck 1 (6.12) 108 In conclusion,kh T k 2 C A C At kh T ck 1 , so from (6.11), k^ xk 1 kxk 1 + (1 C v c C A C At 4c )kh T ck 1 : (6.13) Sincec > 1 andc > C v + C A C At 4 , this impliesh T c = 0. Also Since (A T A T ) has full rank,A T A T h T = 0 only ifh T = 0. Soh has to be zero, which proves the lemma. In [136],C A = 2,C At = 1, andC v = 1 4 . Since, the isotropy condition is not valid in general, we need little flexibility in setting up the dual variable. The proof require an indirect relation betweenc and. During the derivation, it is required that 1 c C A soc has to greater than + 1 C A . In lemma 1, we assumed certain properties of a dual certificate,v. In the next lemma we will prove that a dual certificatev as described in Lemma 1 can indeed be constructed. Lemma 2 (existence of a dual certificate) Letx2C N be assparse vector. Let 0 < < 0:5, d > 0, andC At >= 0:75. Ifs 1 + ( d )= and can also satisfy other two constraints, then the constraints onA in Lemma 1 hold and a vectorv with the properties required for Lemma 1 exists. Proof LetB := cIA A. We use matrix elements based conditions to hold the constraints onA. SinceA is normalized, we can get 0 < < 0:5 such thatjB pq j < ; 8p6= q. The condition states that each off-diagonal is bounded above by . For the diagonal elements, the bound is required only for certainB pp ; p2T . If there exist 0< d such thatjB pp j< d ; 8p2T , we getkB T k 1 < d + (s 1). If the sparsity inx is restricted such that s 1 + ( d )=; (6.14) thenkcIA T A T k 2 =kB T k 2 p kB T k 1 kB T k 1 d + (s 1). For the second condition onA, max p2T ckA T A fpg k 2 = max p2T c s X l2T (A l A p ) 2 = max p2T c s X l2T B 2 lp p sC At : (6.15) 109 The last inequality can be derived as follows. It says thats C At 2 , so we need to show that the maximum allowable sparsity in (6.14) can satisfy this constraint. Starting with the constraint equation and plugging the value ofs, we get 1 + ( d ) = ( + d ) C At 2 ; ( + d ) C At ; C At + d : (6.16) since is less than 0:5 andC At > 0:75, C At 1, which proves the inequality. To get the dual certificate,v, we use golfing scheme introduced in [138]. The basic idea is that each iteration brings us closer to the target (v T closer toc times sign ofx) just as in golf, where each shot would bring the ball closer to the hole. Let’s partition the matrixA intol row blocks (not necessary sequentially) and represent thei th partition by A i having sizem i so thatA = [A 1 ;A 2 ;:::A l ] T and P l i=1 m i =m.A i;T is the restriction ofA i to the columns inT . In golfing scheme,v is defined inductively. We chose to define it as v i = m m i A i A i;T (sgn(x T ) 1 c v i1;T ) +v i1 ; i = 1;::;l; (6.17) where m m i is renormalization factor,v 0 = 0 andv is set tov l . Letq i = sgn(x T ) 1 c v i;T , we get the following relations after some algebra. q i = (I m cm i A i;T A i;T )q i1 = i Y k=1 (I m cm i A k;T A k;T )sgn(x T ) (6.18) v :=v l = l X i=1 m m i A i A i;T q i1 : (6.19) SincekcIA T A T k 2 , we claim thatkcI m m i A i;T A i;T k 2 is smaller than 1. We use i to represent the upper bound,kcI m m i A i;T A i;T k 2 i . The claim implies that the norm ofq i decreases asi increases andv T is getting closer toc sgnx T . 110 Now, we start proving the constraints onv. For the first constraint, kc sgn(x T )v T k 2 =ckq l k 2 cksgn(x T )k 2 l Y k=1 k(I m cm i A k;T A k;T )k 2 ksgn(x T )k 2 l Y k=1 k(cI m m i A k;T A k;T )k 2 p s l Y k=1 k 1 4 : (6.20) Where the last inequality comes by setting l =blog i ( 1 4 p s) c and i is defined as l i 1 2 ::: k :: l . Since each k is less than 1, given the sufficient number of measurements, we can have the required number of blocks,l, and the inequality satisfies. However, we may not have those many measurements. At max,l can be equal toM. If the inequality can not be satisfied for the givenm, we can reduce the sparsitys such that the inequality holds. The sparsitys can be bounded by s 1 16 2m i : (6.21) v T c can be written as P l i=1 m m i A i;T cA i;T q i1 . Therefore, kv T ck 1 l X i=1 k m m i A i;T cA i;T q i1 k 1 l X i=1 k m m i A i;T cA i;T k 1 kq i1 k 1 l X i=1 s i kq i1 k 1 s l X i=1 i 1 c ksgn(x T )k 1 l Y k=1 k(cI m m i A k;T A k;T )k 1 s l X i=1 i 1 c i1 Y k=1 k ; (6.22) where we used submultiplicative property of induced norm in the second line; defin- tion ofq from (6.18), Holder’s inequality, and symmetric matrix(kAk 1 =kAk 1 ) in the fourth line; and i is the upper bound for off-diagonal elements for thei th block. 111 The rightmost term can be made lower thanC v by choosing the proper value ofs as it is shown below. s c ( 1 + 2 1 + 3 1 2 +::: + l 1 2 ::: l1 ) s c m (1 + m + 2 m +::: + l1 m ) = s c m 1 l m 1 m ; (6.23) where m is the maximum of i ;8i and m is the maximum of k ;8k. To make kv T ck 1 C v , we need s< C v c m 1 m 1 l m : (6.24) Is sparsitys can satisfy (6.14), (6.21), and (6.24), then the constraints onA andv hold. It should be clear that the constraints appearing here are highly likely optimal. 6.4 Discussion The conditions here do not give universal results in the sense that for a given set of measurements, the exact recovery is applicable up to certain non-zeros in x. Beyond this sparsity level, there is a chance of getting a nontrivial null space. If the bound on the conditions are tight, there would be a nontrivial null space for the next nonzero right after the number of non-zeros inx is equal to the estimateds. The constraint on sparsity is essentially converting a rank-deficient system to a ’full-rank’ (for the submatrix) system and if we choose an optimization approach to solve (6.1) such that it finds at maxs sparse solution rather than inverting the full matrix, the recovery would be exact. The derivation in this article does not rely on the idea of randomness, which is the central idea in the compressive sensing. The reason is that the analysis starts from the given matrix. It is highly likely that random measurements can produce lower and d but the condition is not necessary. We will later in this article, test the effect of random measurements. For a sparse recovery, the usual practice is to relate the number of measurements to the number of unknowns and usually it is on the order of s log n, which is considered to be sub-optimal if not optimal. However, in the derivation above, we don’t have n in any of the constraints directly. Does it mean the irrelevance of n to sparse recovery? No, actually the effect of n is encoded in the values of and d . For example, if the 112 measurements are sampled from the DFT matrix, the encoding can be clear from the element of the matrixe i2kn i =n . The relation with an example will be demonstrated in the next section. For random sensing, the scaling of the matrix does the trick. It changes and s . 6.5 Analysis In this section, we test the variation in as we vary the number of measurements and the number of unknowns. We consider one the mostly used sensing matrix in compressed sensing, Discrete Fourier Transform (DFT). The measurements are randomly selected from the DFT matrix. In the first test, we keep n constant(n = 512) and vary the number of measurements from 0:033n ton. The columns of the matrix,A is normalized by 1 p m . The maximum off-diagonal element ofA A, is plotted against m n in Fig. 6.1. 0 0.2 0.4 0.6 0.8 1 m/n 0 0.1 0.2 0.3 0.4 0.5 0.6 zeta Figure 6.1: the maximum off-diagonal element ofA A, as a function of number of measurements for fixed number of unknowns. It shows the monotonically decreasing behavior of. Since is inversely proposal 113 to the sparsity,s, sparsity increases asm increases. The plot here is one of the instances of the randomly selected measurements. If we run 50 random sets of measurements, the curve becomes more smooth as shown in Fig. 6.2. 0 0.2 0.4 0.6 0.8 1 m/n 0 0.1 0.2 0.3 0.4 0.5 0.6 zeta Figure 6.2: Average behavior of the maximum off-diagonal element of A A, as a function of number of measurements for fixed number of unknowns. Next, we consider the relation between and number of unknowns. Since number of measurements can change the value of zeta and hence, can potentially change the relation of andn, we choose the number of measurements as a fraction ofn and take three different fraction-0:2n; 0:3n; 0:4n. The changes in as function ofm andn is shown in Fig. 6.3. To make a comparison, we usem/slog(n) for RIPless compressive sensing. We normalized every graph by its first element. As can been seen, both have similar behavior. Asn increases, number of measurement also increases in this setup so the sparsity also increases as expected. Even though, there is no explicit log factor in our formulation, it seems that may have embedded it with a different scaling factor. 114 0 500 1000 1500 2000 2500 # of unknowns 0 0.2 0.4 0.6 0.8 1 1/s fixed measurement to unknown ratio element CS (0.2) element CS (0.3) element CS(0.4) RIPless CS(0.2) RIPless CS (0.3) RIPless CS (0.4) Figure 6.3: Comparative analysis. variation of inverse of sparsity to the number of un- knowns for three set of a number of measurements as a fraction of number of unknowns. The fractions are 0:2, 0:3, and 0:4. Next, we evaluate the effect of random measurements on. The experiment setup is the same as we have in Fig. 6.1 but here we also consider a different set of measurements generated by uniformly sampled the rows of the DFT matrix. The values for different ms are plotted in Fig. 6.4. As can be seen, the is consistently higher for uniform samples, which is in harmony with the concept of compressive sensing. 115 0 0.2 0.4 0.6 0.8 1 m/n 0 0.2 0.4 0.6 0.8 1 1.2 zeta randomly sampled uniform sampled Figure 6.4: Effect of random measurements on. 6.6 Conclusion A theory of compressive sensing for a unique recovery is investigated. The theory is applicable only when the measurements are acquired in isotropic setting. To extend the applicability, a dual certificate based lemmas has been proposed. The proof for the Lemmas has been derived putting sparsity constraints through the values of matrix elements rather than on a number of measurements. The Lemmas are tested for a sensing matrix as the DFT matrix. The analysis shows that these new Lemmas have similar dependency on the number of unknowns and on the number of measurements as they have in the published literature of compressive sensing. Thus, the matrix element based conditions have shown theoretical consistency with the compressive sensing philosophy for measurement in Fourier domain, which is isotropic in nature. In the next Chapter, we consider anisotropic measurement setup, which has known to have an upper bound on the number of measurements. The setup we are considering is usually called as inverse scattering problems. 116 CHAPTER 7 Inverse Scattering Problems 7.1 Background We assume a z-independent, heterogeneous, isotropic, and non-magnetic two dimen- sional (2D) scatterer, immersed in a homogeneous background medium with the lossless wavenumberk b . The scattered field for this geometry is written as E s (r) =k 2 b Z D G(r;r 0 )J(r 0 )dv 0 ; r2S; r 0 2D; (7.1) where E s (r) is the scattered electric field, J(r 0 ) is the current source, and G(r;r 0 ) is the Green’s function. In the inverse problem, dielectric constant and conductivity are included inJ(r 0 ) in the above equation. If the contrast is zero,J(r 0 ) is also zero so the sparsity in the contrast domain is preserved by default. The above formulation is linear with respect to the current source and can be written for M scattered field measure- ments in vectorb, andN unknown current source locations in a matrix form such that b = CGz, whereC = j 2 k b aJ 1 (k b a),G = H (2) 0 (k b jrr 0 j),J is the Bessel function of 1 st order,H (2) 0 is the 0 th order Hankel function of the second kind, anda is the equiv- alent radius of a pixel. An inverse scattering problem can be written in the compressive sensing framework in noiseless condition as, minkzk 1 s:tGz =b: (7.2) The system matrix G certainly does not follow isotropy condition and has a theoret- ical limit on how long the measurement vector should be to fully represent scattered field[140] or in other words, the information in the scattered fields can be represented by a fix length vector. Any additional measurement can be estimated accurately from 117 the given set of the measurements. Hence, the rows ofG are linearly dependent beyond certain number of the measurements. Let’s consider a typical setup and evaluate the value of in the same manner as we did for the measurement in Fourier domain. We use the notationF :=G H G. The simulation is performed using the following configuration. The wavelength in air is 15 cm. The background medium is considered to be as air/vacuum. The T/R antennas are located on the circumference of the circular domain of radius 5 uniformly apart from each other and the investigation domain is within a square of 6 length. The pixel size in both directions is 10 , which gives image size of 60 60 pixels. In the first test, we vary the number of measurements from 5 to 50. The columns of the matrix,G is normalized by 1 p m after initially rescalingG by a constant factor for all different measurements. The maximum off-diagonal element ofF , is plotted against m in Fig. 7.1.It shows that decreases up to a certain point and then it remains the con- stant, which suggests that the measurement matrix becomes the same and the additional measurements do not add any information. 0 10 20 30 40 50 # of measurements M 0.55 0.6 0.65 0.7 0.75 zeta Figure 7.1: the maximum off-diagonal element of F , as a function of number of measurements for fixed number of unknowns. The similar effect can be seen if we visualize one of the rows of F . For better 118 visualization, the investigation domain is reduced to 3 length, which gives total 900 pixels. For one randomly selected column,F is plotted against the indexes for different measurements in Fig. 7.2. AS can been seen, for total measurements 30 and 60, there is no difference for all elements. The off-diagonal oscillation diminishes asM increases. Figure 7.2: One of the columns ofF , for different number of measurements. 119 The above plots show that the matrix element based Compressive sensing can cap- ture the band limited nature of the measurements. In the next section, the same result on F is derived analytically. 7.2 Derivation Let’s takeM measurements in the homogeneous background. For the contrast source formulation as described above, the sensing matrixG in 2D can be represented by Han- kel function. The element ofF ,F pq is expressed as F pq =G H T G T (p;q) = 1 M M X l=1 H (2) 0 (k b jr l r 0 p j) H H (2) 0 (k b jr l r 0 q j): (7.3) Let’s modify this by applying some approximations and theorems so we can analyze. Assuming the fields are measured in far field region, Hankel function is approximated byH (2) 0 (x) = q 2 x e j(x=4) . Hence, (7.3) becomes F pq = 2 k b M M X l=1 1 q jr l r 0 q jjr l r 0 p j e jk b (jr l r 0 p jjr l r 0 q j) : (7.4) For the far field, the distance can be approximated for magnitude as follows. jrr 0 j = p r 2 +r 02 2rr 0 cosr q jr l r 0 q jjr l r 0 p jr l (7.5) For the phase term, the approximation is jrr 0 j = p r 2 +r 02 2rr 0 cosrr 0 cos jr l r 0 p jjr l r 0 q jr 0 q cos( lq )r 0 p cos( lp ): (7.6) 120 Hence, (7.4) becomes F pq 2 k b M M X l=1 1 r l e jk b (r 0 q cos( lq )r 0 p cos( lp )) = 2 k b M M X l=1 1 r l e jk b (r 0 q cos( lq )) e jk b (r 0 p cos( lp )) : (7.7) Using the Bessel function property,e jzcos = P 1 n=1 j n J n (z)e jn , we get, F pq 2 k b M M X l=1 1 r l 1 X s=1 j s J s (k b r 0 q )e js lq 1 X t=1 j t J t (k b r 0 p )e jt( lp ) : (7.8) Changing the order of summation ande jt = (1) t , we get F pq 2 k b M 1 X s=1 j s J s (k b r 0 q ) 1 X t=1 (1) t j t J t (k b r 0 p ) M X l=1 1 r l e j(s lq t lp ) : (7.9) Scattered fields can be optimally sampled using uniform sampling scheme [48], which can be used to decompose the angle as lp = 2m l M p ; m l 2f0; 1;::;Mg. If the samples are on a constant radius,r l =r =b, we get F pq 2 k b Mb 1 X s=1 j s J s (k b r 0 q ) 1 X t=1 (1) t j t J t (k b r 0 p )e j(tpsq) M X l=1 e j 2m l M (st) : (7.10) It can be further simplified as F pq 2 k b b 1 X s=1 j s J s (k b r 0 q ) 1 X v;t=1 (1) t j t J t (k b r 0 p )e j(tpsq) [st +vM]: 2 k b b 1 X v;t=1 j tvM J tvM (k b r 0 q )(1) t j t J t (k b r 0 p )e j(tp(tvM)q) 2 k b b 1 X v;t=1 j vM J tvM (k b r 0 q )J t (k b r 0 p )e j(tp(tvM)q) 2 k b b 1 X v;t=1 J tvM (k b r 0 q )J t (k b r 0 p )e j = (t p (tvM) q vM=2): (7.11) 121 where we use the delta function identity in the first equation andj 2t (1) t = 1. Some of the Bessel functions with different orders are shown as a function of argument(k b r) in Fig. 7.3 and as a function of the order for different arguments is shown in Fig. 7.4. Figure 7.3: Bessel functions of different order versus the arguments Figure 7.4: Bessel functions versus the order for different arguments The figures show the oscillating and damping behavior for the argument and de- creasing to zero for the order. In other words, Bessel functions are band-limited. The Bessel functionJ n (x) 0 ifjnj > 2djxje [141]. Applying this condition on the last equation, we get 122 k b r 0 m <t<dk b r 0 m e tdk b r 0 m e M <v< t +dk b r 0 m e M ; (7.12) wherer 0 m is the far away object’s point from the center. AsM increases, t+dk b r 0 m e M be- comes smaller and afterM >t +dk b r 0 m e, the ratio becomes less than 1 andv takes only 0 value. Plugging this value ofv in the last equation of (7.11) gives F pq 2 k b b dk b r 0 m e X t=dk b r 0 m e J t (k b r 0 q )J t (k b r 0 p )e jt(pq) ; (7.13) which is M independent. In other words, remains the same for all M > 2dk b r 0 m e. The optimal value ofM obtained here matches with the previously derived values [49] [141]. Moreover, here we are able to link it to the element of the matrix F , which is used in estimating the sparsity in a signal that can be recovered exactly. In addition, this can also be helpful in gradient based methods as in those methods, the update on the unknowns are driven by the product of a matrix and its transpose. For the diagonal element, F pp becomes constant ( 2 k b b ). The constant value of the diagonal elements can help in choosing the valuec easily. The approximation error we have in estimation of theF pp is not major. The typical variation in the diagonal elements is shown in Fig. 7.5. a o - a i 1 1.5 2 2.5 3 diagonal variations 0 0.05 0.1 0.15 0.2 0.25 Figure 7.5: The effect on diagonal variations for different measurement radius,b for a given domain radius,a. 123 Here, the approximation error is defined as the maximum absolute difference be- tween the diagonal elements,b is the measurement radius anda is a given domain radius. The radii are defined in the wavelength scale. 7.3 Analysis In this section, we will show that the conditions in Lemmas hold for the inverse scatter- ing problem. Since the matrix,G does not have a property that can be used to validate the Lemma 2 universally, we consider the following setup to validate the conditions for a given setup. The wavelength in air is 15 m. The background medium is considered to be as air/vacuum. The 100 receiving antennas are randomly located on the circumference of a circular do- main of radius 3 and the investigation domain is within a square of 3 length. The pixel size in both directions is 4 , which gives image size of 12 12 pixels. Even though we may not need many receivers in practice as the current source has to be physically realistic, here we are considering more number of measurements because we can test only for any random current source and we use the golfing scheme in the derivation of v. The scheme, which divides measurements into small blocks, may need more mea- surements as we have not chosen carefully neither the measurements for each sub block nor the number of blocks for the optimality of the number of measurements. Since the conditions ons are not highly likely optimal, we estimated the expected sparse signal that can be recovered exactly usings < 1 + d , which turns out to be roughly 1:7. So we have considereds = 2 and tested all sparsity conditions on all of the pos- sible combinations. The possible combinations for G H T G T is 144 2 . We scaled the G matrix properly and setc = 1:5 and estimatedc A , c A t, and. The maximum value is 2:24; 0:85; 1:05, respectively. Using these values, we can see that the conditions onG hold and the relation betweenc and holds. Next, we have to validate the conditions onv. The last equation in 6.20 and 6.22 are used to validate the conditions. We initially started with a lower number of measurements but the conditions were not satisfied. The points, where the conditions failed for one of the lower measurements setup is shown in Fig. 7.6. Each axis represents a linearized index (point). As can be seen, the failure points are very close to the diagonal of the plot and it seems that they have a repeating pattern. When a different set of measurements but the same number is used, some of the points got disappeared. So it is highly likely that for some configuration ofm,m i , and 124 l, it is possible to satisfy the conditions. There are many permutations that we can try but rather than checking for different permutations, we chose to go higher in number of measurements. For the setup specified above, all points are lower thanC v forv T c and lower than l i forv T . Their maximum value is 1:55 and 0:065, respectively. 0 20 40 60 80 100 first point index 0 20 40 60 80 100 second point index Figure 7.6: The points that are failed to satisfy the condition onv T . As pointed out earlier, we would like to reiterate that this numerical exercise only validate the exact recovery for this particular configuration for the sparsitys = 2. We also tested many other setups, where we vary the measurement radius, pixel spacing, and domain length. For each of the setup, all the conditions did hold for sparsitys = 2. We also perform some numerical analysis, where the contrast source is recovered for different sparse signals. The simulation setup is the same as in the above analysis. For the unknowns, 50 samples are created as follows. The sample is configured by the value of the contrast source, the type of the contrast source- same vs different, and their 125 location. Five different contrast source values are picked arbitrarily; the different value is created by lowering the contrast source by a random factor; and five set of locations are considered- all centered, all in top left side of the domain, all in bottom right side of the domain and the other two random. Compressive Sampling Matched Pursuit (CoSaMP) algorithm [142] is used to recover the unknown contrast source. During each run, we pass the true sparsity value as a parameter because the performance of CoSaMP is sensitive to the value of sparsity. The performance is evaluated by the ` 2 error defined between original contrast source and estimated contrast source. The error fors = 1 tos = 6 is shown in Fig. 7.7. 1 2 3 4 5 6 sparsity(s) 10 20 30 40 50 sample # 0 1 2 3 4 5 Figure 7.7: error in the sparse recovery using CoSaMP. The samples are first ordered by location, then by value, and then by type. As can be seen, there is no error for any samples fors< 2. Fors = 3, the error comes for the points, which are randomly located in the domain but it is for few samples. Fors = 4; 5, 126 the error rate for randomly located points increases. It is interesting to see that there is no error fors< 5 when the points are clustered anywhere in the domain. Fors> 6, we see error even for the clustered points. 7.4 Significance For the contrast source as unknowns such as in the equation (7.1), it is well accepted that the non-trivial solution exist. In addition, the minimum` 2 norm solution produced, for example, by the conjugate gradient method is not the appropriate physical solution. So to overcome this problem, more than one incidence angle is considered and then the equation has been used in two different flavor. In the first approach, the minimum norm solution is obtained and then the solution is projected on a subspace that is defined by the V olume Integral Equation. A recent article on this approach is [34]. In the other approach, the V olume Integral Equation is used as a constraint (sometimes referred as the Helmholtz constraint) to the equation (7.1). Contrast Source Inversion[35] is a standard example. In both of the above approaches, the inversion process becomes computationally more intense as a large size square matrix has to be inverted iteratively. However, here, using the lemmas above, it can be shown that it is possible to recover the very sparse objects exactly. This approach does not need an inversion of a large square matrix and can be solved in a short time, which is desirable in scenarios such as through-the-wall imaging, medical therapy/monitoring or aerial surveillance. In these scenarios, we expect that the investigation domain can be sparse either in the scattering formulation or in a differential scattering formulation. The number of sources is limited and fast detection or identification is one of the requirements. 7.5 Application for Contrast Recovery In the scenarios described above, we can estimate the domain electric properties in a very simple two steps process. We assume the simulation setup is noiseless. Once we retrieve the contrast source, we can estimate the electric fieldE(r) within the domain using the following equation: E(r) =E in (r) + n X r 0 =1 C F G F (r;r 0 )z(r 0 ); r2D; (7.14) 127 where G F is the same Green’s function and C F is a constant. The superscript in G F is used to denote that this function takes both inputs from the domain. Point should be noted thatG F needs only one time computation and there is no matrix inversion. The estimatedE(r) is then used to calculate the contrast,x, as follows. x i = z i E H i kE i k 2 ; 8i2f1;::;ng: (7.15) We use the following configuration for the experiments. Plane waves are considered as the source with 0 o as an incidence angle. The wavelength in air is 15 m. The back- ground medium is considered to be as air/vacuum. The receiving antennas are located at equidistant on the circumference of a circular domain of radius 4 and the investigation domain is within a square of 3 length. The number of receivers is 30. For the first set of experiments, the pixel size in both directions is kept at 10 , which gives image size of 30 30 pixels. There are two nonzero randomly located pixels and both have the same value. Two cases are considered: contrast is 2 and 100, respectively. The contrast source recovery and also the estimated contrast are shown in Fig. 7.8 and Fig. 7.9. As can be seen, the recovery is very accurate. 128 0 200 400 600 800 1000 linear index -3 -2 -1 0 1 2 3 Re(z) true estimated 0 200 400 600 800 1000 linear index -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 -Im(z) true estimated -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 0.5 1 1.5 2 -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 0.5 1 1.5 2 Figure 7.8: case 1: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. 129 0 500 1000 linear index 0 0.5 1 1.5 2 2.5 3 Re(z) true estimated 0 500 1000 linear index -6 -5 -4 -3 -2 -1 0 -Im(z) true estimated -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 20 40 60 80 100 -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 20 40 60 80 100 Figure 7.9: case 2: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. To quantify our observations, we use the relative maximum error and relative eu- clidean error to characterize contrast estimation. The relative maximum error, n is de- fined as the ratio of the maximum absolute difference between the true contrast source and estimated contrast source to the maximum absolute value of the true contrast source. The second metric, relative error is defined as, = 1 N x N y v u u t Nx X i=1 Ny X j=1 rec r (i;j) ori r (i;j) rec r (i;j) 2 100%; (7.16) where rec r and ori r are complex permittivity of a recovered object and an original object respectively. The metric numbers are reported in Table 7.1, which shows that errors are 130 very low and on the order of 10 11 . For the case having contrast of 100, the grid size was not sufficiently small to dis- cretize the integration accurately. However, we used it to check if there is effects of the contrast value on the contrast recovery for a fixed size problem. We tested on many other cases and the results are very accurate. In the following set of experiments, we will make sure that the grid size is sufficiently small. First case has two random nonzero points in the domain with the contrast of 5. Second case has four nonzero clustered points at the center of the domain and has uni- form contrast of 3. The contrast estimation in 7.10 and 7.11, respectively shows that the object’s position and contrast are recovered accurately. 0 200 400 600 800 1000 linear index -5 0 5 10 Re(z) true estimated 0 200 400 600 800 1000 linear index -6 -4 -2 0 2 4 -Im(z) true estimated -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 1 2 3 4 5 -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 1 2 3 4 5 Figure 7.10: two points: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. 131 0 200 400 600 800 1000 linear index 0 0.5 1 1.5 2 2.5 3 Re(z) true estimated 0 200 400 600 800 1000 linear index -0.5 0 0.5 1 1.5 2 2.5 -Im(z) true estimated -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 0.5 1 1.5 2 2.5 3 -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 0.5 1 1.5 2 2.5 3 Figure 7.11: small object: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. Next, we consider the contrast of 20 for two non zeros points located randomly and the contrast of 25 for two clustered points at the center. The pixel size in both directions is kept at 25 , which gives image size of 75 75 pixels. There are more than 5000 unknowns and only 30 measurements are taken but since we know that the domain has at most two nonzeros,the contrast source recovery and contrast recovery are accurate as shown in Fig. 7.12 and Fig. 7.13. A more detailed or zoomed-in image for Fig. 7.13 is shown in Fig. 7.14, which shows the accurate contrast recovery even for two points having different contrast, which are next to each other. The quantitative numbers for two high contrast case are similar to other cases but for high contrast, the numbers are higher but still well less than noticeable change in contrast recovery. 132 0 2000 4000 6000 linear index -40 -20 0 20 40 Re(z) true estimated 0 2000 4000 6000 linear index -30 -20 -10 0 10 -Im(z) true estimated 0 2 10 20 2 y/ 0 x/ 0 -2 -2 0 2 10 20 2 y/ 0 x/ 0 -2 -2 Figure 7.12: two high contrast: The top row shows the real part (in the left) and imag- inary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. 133 0 2000 4000 6000 linear index 0 5 10 15 20 25 Re(z) true estimated 0 2000 4000 6000 linear index 0 10 20 30 40 -Im(z) true estimated -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 5 10 15 20 25 -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 5 10 15 20 25 Figure 7.13: high contrast: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. 134 2725 2730 2735 2740 2745 linear index 0 5 10 15 20 25 Re(z) true estimated 2725 2730 2735 2740 2745 linear index 0 10 20 30 40 -Im(z) true estimated -0.15 -0.1 -0.05 0 0.05 x/ -0.15 -0.1 -0.05 0 0.05 y/ 0 5 10 15 20 25 -0.15 -0.1 -0.05 0 0.05 x/ -0.15 -0.1 -0.05 0 0.05 y/ 0 5 10 15 20 25 Figure 7.14: zoom-in of 7.13: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. Since the theoretical analysis is carried out in noiseless scenarios, we committed in- verse crime in all above cases to validate the analysis by simulations. What will happen if we avoid the inverse crime? Before we go into inversion, we tested many cases, where we interpolate both the contrast and electric field and compute the estimated scattered field. This field is then compared to the simulated scattered field through the least square error. Since the field inside the object can vary even though the contrast is constant, the error can vary a lot. It depends on various factors such as a grid-size difference between imaging and inversion grid, interpolation technique for the field and the contrast, and if both grids are aligned or not. We consider only the cases, where the scattered field error is less than 5% of the synthesized scattered field. We consider a square object having contrast of 4 located at the center in one setup and in the top-right side in the other setup. 135 For the first case, the pixel size during imaging in both direction is 36 and for the second case, it is 30 . The inversion pixel is 12 for the first case and 15 for the second. The size of the square is fixed such that the object in the inversion has 4 nonzero elements. The recovery results are shown in Fig. 7.15 and 7.16. As can be seen, the object location, size, and contrast are recovered accurately with the low quantitative errors 7.1. 0 500 1000 1500 linear index 0 0.5 1 1.5 2 Re(z) true estimated 0 500 1000 1500 linear index -1 0 1 2 3 4 5 -Im(z) true estimated -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 1 2 3 4 -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 1 2 3 4 Figure 7.15: square object: The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. 136 0 500 1000 1500 2000 2500 linear index -3 -2.5 -2 -1.5 -1 -0.5 0 Re(z) true estimated 0 500 1000 1500 2000 2500 linear index -5 -4 -3 -2 -1 0 1 -Im(z) true estimated -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 1 2 3 4 -1 0 1 x/ -1.5 -1 -0.5 0 0.5 1 1.5 y/ 0 1 2 3 4 Figure 7.16: off-center square object : The top row shows the real part (in the left) and imaginary part (in the right) of the contrast source and the bottom row shows the original contrast configuration in the left and recovered contrast configuration in the right. 7.6 3D Green’s Function To derive the conditions for 3D, we start the analysis from a homogeneous background medium. The Green’s function can be written asG(r;r 0 )/ e j(jrr 0 j) jrr 0 j . The function is very similar to (7.4). In fact, the difference up to a constant factor is that for 3D, the equation has 1 r variation rather than 1 p r . We can follow the same derivation and can reach to the same conclusion. 137 Table 7.1: Quantitative measurement of the reconstructed images. Object r m Error () (%) low linear 2 7:68 10 8 9:62 10 9 high linear 100 2:2 10 8 3:98 10 8 two points 5 2:55 10 2 1:19 10 2 small object 3 1:54 10 2 4:57 10 3 two high contrast 20 7:65 10 2 1:84 10 2 high contrast 25; 12:5 0:29 7:45 10 3 square object 4 0:14 4:56 10 2 off-center square object 4 0.12 2:14 10 2 7.7 Conclusion An application of compressive sensing theory to inverse scattering problems has been discussed. We have considered a typical 2D imaging setup, where the background is homogeneous. The contrast-source formulation for a single source has used as an inver- sion process and the conditions for the unique recovery of the inversion has been tested. The following are inferred. 1. It is possible to have a unique recovery for the very sparse objects for the imaging setup considered here. The finding deviates from the conventional understanding of the existence of a non-unique solution in the single-source setup. 2. The finding were numerically validated, which shows that the contrast can be recovered very well and have relative error less than 0:03%. 3. The recovery process is very simple, and does not need any forward problem’s matrix inversion. Hence, it has a very low computational complexity. 4. If modeling error is within the bound for the cases where ’inverse crime’ is avoided, the performance for contrast recovery matches to the performance of the cases, where ’inverse crime’ is committed. Future work includes testing the method for diverse system matrices and applying it to noisy scenarios. 138 CHAPTER 8 Concluding Remarks and Future Work Here we summarize the thesis and note the key points of the thesis. For more detailed conclusions, please refer the respective chapter. In the thesis, various constraints have been applied especially for the 2D microwave imaging problem. We tested them on different objects and achieved a better recon- struction of the images. After seeing the results, it can be said that each constraint can improve a different aspect of the reconstructed image. Depending on requirements of an application of interest, any one of the methods can be useful. However, to solve the non-linear problem in general, we need to combine different constraints in a meaningful way to solve that multi-constrained cost function. We demonstrated the applicability and a solution approach for one of the multi-constrained cost function in 3D microwave imaging problem and were able to achieve improved reconstructed images. We generalized the applicability of the compressive sensing theory to more ‘general’ system matrices. Using this development, we showed that a unique recovery is possible for a single source 2D imaging problem formulated in the contrast-source formulation. For the future work, it will be interesting to test the theorems in various imaging setups (hence, different system matrices) and to be able to tighten the conditions for the unique recovery. As an extension of the work presented in this thesis on the constraints is to apply the 3D microwave imaging formulation to publicly available breast phantom database. Even for the linear case, the performance can be improved, if multi-frequency data is available and an additional constraint, such as wavelet-based constraint, can be incorporated in the contrast estimation cost function. This may be an involved problem since wavelets are not localized as total variation, therefore applying the wavelet transform on some part of the 3D domain will be an interesting problem to solve. Another direction that can be pursued is to use information available in the ‘data’ (such 139 as breast phantom database). For a linear problem such as computed tomography or MRI, deep learning based imaging methods are already published in the last year or two. One of the examples is learning artifacts from the available data. Once the arti- facts are learned using a deep learning model, then this learning is used to remove the artifacts generated during a new imaging scan. For most of the published articles, the common theme is to get a low resolution image or a well defined initial guess using a conventional method and then apply a learning based constraint to remove the arti- facts or reduce the noise. However, for microwave imaging, getting the correct initial guess will be a challenge because of the nonlinearity of the problem. It will be of great interest for the microwave imaging community if the learning based constraint can be applied similarly to the other imaging problems, which then can be used to improve the resolution, reduce computational time, or solve a strong scatterer scenarios. 140 Appendix Projection on a` 1 ball and a box Before starting to solve (2.11), let’s consider the optimization problem without box constraints and also make all vectors real. minimize x kcxk 2 2 s.t. kxk 1 : (1) We note that ifkck 1 then, the optimization problem is trivial and the solution is x =c. Therefore, from now on we assume thatkck 1 , which means that the optimal solution is on the boundary of the L1 ball, i.e.,kxk 1 = . Lemma 3 from Duchi et al. [143] states that each element of the optimal solution ~ x shares the same sign as each element of the given vectorc. In other words, ~ x andc are in the same orthant (quadrant for two-dimensional ). Therefore, solving the absolute value problem is equivalent to solving a sign value problem. Eq. (1) can be expressed as minimize x kjcjxk 2 2 s.t. kxk 1 ; x i 0 (2) and solution of (1) is ~ x =sgn(c) T x. Our optimization problem is an extension of the above and it has been shown that Lemma 3 of [143] can also be applicable for a real vector’s projection on a constrained box [144], provided necessary geometric transformation has been applied. If both bounds are in the same orthant, we can apply the geometric transformation such that the bound close to zero is shifted to the origin. We also update by jx min j 1 if bounds are positive or byjx max j 1 , if they are negative. When both bounds are in the different orthant, we select the upper bound based on the sign ofc i so that ~ x i is close toc i . The upper bound isx i;min ifc i < 0 or it isx i;max ifc i > 0. For all cases the 141 lower bound is always zero. Let’s call the upper bound as b i ; the equivalent problem then becomes minimize x kjcjxk 2 2 s.t. kxk 1 ; 0x i b i i = 1;:::;n: (3) The vector x is estimated by expressing the above optimization problem in La- grangian form L(x;;;) = 1 2 kjcjxk 2 2 +( n X i=1 x i ) T x + T (xb); (4) where2 R + and;2 R n + are Lagrangian parameters. Differentiating with respect tox i and comparing it with zero gives optimality condition, @L @x i =jc i j +x i + i + i = 0 (5) The complementary slackness conditions for parameters [81] are i x i = 0; i (x i b i ) = 0 (6) To satisfy all of the conditions,we get following relations; • When 0 < x i < b i , i = 0 and i = 0. This givesjc i j +x i + = 0. ) x i = jc i j. • Whenx i = 0, i = 0. This givesjc i j + i = 0. Since i 0 we getjc i j • Whenx i = b i , i = 0. This givesjc i j +b i + + i = 0. Since i 0 we get jc i jb i + • Ifjc i j>b i and i > 0, the distance betweenjc i j andx i is minimum whenx i =b i . This is only true if i > 0 because in generaljc i j>b i can be projected inside box also. Summarizing the relations in known variables gives the unknownx i as in Equation (2.12). The only unknown is estimated similar to an unbounded problem and is detailed in [143]. The above solution for box constraints is derived assuming thatkck 1 . Whenkck 1 , we still have to ensure thatx is within bounds, unlike the unbounded 142 case. The solution, whenkck 1 and box constraints are active, is exactly the same as (2.12) except that now = 0. For inverse scattering, the system of equations are complex valued. The concept used for the real vector, which stated that the distance between the original and pro- jected vectors is the lowest, if they are in the same orthant, can extend to the complex domain. Complex numbers have a phase instead of a simple sign. We define sign func- tion for a complex numberz assgn(z) = z kzk ;8z2 C (we consider positive sign for z = 0). When we apply the sign function to a complex number, the function projects the vector on the unit circle of a complex domain (phase of a vector). Therefore, we can use the same projection method as we have applied to the real vector to a complex vector. Though the method can work very well for this polar representation of complex-valued vectors, we can alternatively use the Cartesian (real and imaginary) representation of the same vectors. We get additional sparsity if any of the two (real or imaginary) com- ponents is zero. It should be noted that the entire problem can be solved in the complex domain, while still performing projections using real and imaginary expansion of the complex vector (c2C n ). 143 BIBLIOGRAPHY [1] G. Oliveri, N. Anselmi, and A. Massa, “Compressive sensing imaging of non- sparse 2d scatterers by a total-variation approach within the born approximation,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 10, pp. 5157–5170, 2014. [2] M. Benedetti, D. Lesselier, M. Lambert, and A. Massa, “Multiple-shape recon- struction by means of multiregion level sets,” IEEE Transactions on Geoscience and Remote Sensing, vol. 48, no. 5, pp. 2330–2342, 2010. [3] Z. Wu, H. McCann, L. Davis, J. Hu, A. Fontes, and C. Xie, “Microwave- tomographic system for oil-and gas-multiphase-flow imaging,” Measurement Sci- ence and Technology, vol. 20, no. 10, p. 104026, 2009. [4] L.-P. Song, C. Yu, and Q. H. Liu, “Through-wall imaging (twi) by radar: 2-d tomographic results and analyses,” IEEE Transactions on Geoscience and Remote Sensing, vol. 43, no. 12, pp. 2793–2798, 2005. [5] R. Persico, R. Bernini, and F. Soldovieri, “The role of the measurement config- uration in inverse scattering from buried objects under the born approximation,” IEEE Transactions on Antennas and Propagation, vol. 53, no. 6, pp. 1875–1887, 2005. [6] F. Soldovieri and R. Solimene, “Through-wall imaging via a linear inverse scat- tering algorithm,” IEEE Geoscience and Remote Sensing Letters, vol. 4, no. 4, pp. 513–517, 2007. [7] P. Nkwari, S. Sinha, and H. C. Ferreira, “Through-the-wall radar imaging: a review,” IETE Technical Review, pp. 1–9, 2017. [8] T. J. Cui, A. A. Aydiner, W. C. Chew, D. L. Wright, and D. V . Smith, “Three- dimensional imaging of buried objects in very lossy earth by inversion of vetem data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 41, no. 10, pp. 2197–2210, 2003. [9] V . Ganji, N. Gucunski, and A. Maher, “Detection of underground obstacles by sasw methodnumerical aspects,” Journal of geotechnical and geoenvironmental engineering, vol. 123, no. 3, pp. 212–219, 1997. 144 [10] C. Gilmore, M. Asefi, J. Paliwal, and J. LoVetri, “Industrial scale electromagnetic grain bin monitoring,” Computers and Electronics in Agriculture, vol. 136, pp. 210–220, 2017. [11] S. Ganapathy, W. Wu, and B. Schmult, “Analysis and design considerations for a real-time system for non-destructive evaluation in the nuclear industry,” Ultra- sonics, vol. 20, no. 6, pp. 249–256, 1982. [12] L. Crocco, G. Prisco, F. Soldovieri, and N. Cassidy, “Advanced forward model- ing and tomographic inversion for leaking water pipes monitoring,” in Advanced Ground Penetrating Radar, 2007 4th International Workshop on. IEEE, 2007, pp. 127–131. [13] H. K. Thio, X. Song, C. K. Saikia, D. V . Helmberger, and B. B. Woods, “Seis- mic source and structure estimation in the western mediterranean using a sparse broadband network,” Journal of Geophysical Research: Solid Earth, vol. 104, no. B1, pp. 845–861, 1999. [14] G. Ekstr¨ om, “Global detection and location of seismic sources by using surface waves,” Bulletin of the Seismological Society of America, vol. 96, no. 4A, pp. 1201–1212, 2006. [15] S. Nazarian, I. Stokoe, H. Kenneth, and W. Hudson, Use of spectral analysis of surface waves method for determination of moduli and thicknesses of pavement systems, 1983, no. 930. [16] C. B. Park, R. D. Miller, and J. Xia, “Multichannel analysis of surface waves,” Geophysics, vol. 64, no. 3, pp. 800–808, 1999. [17] C. G´ elis, D. Leparoux, J. Virieux, A. Bitri, S. Operto, and G. Grandjean, “Numer- ical modeling of surface waves over shallow cavities,” Journal of Environmental & Engineering Geophysics, vol. 10, no. 2, pp. 111–121, 2005. [18] C. DeSantis, J. Ma, L. Bryan, and A. Jemal, “Breast cancer statistics, 2013,” CA: a cancer journal for clinicians, vol. 64, no. 1, pp. 52–62, 2014. [19] R. A. Smith, D. Manassaram-Baptiste, D. Brooks, M. Doroshenk, S. Fedewa, D. Saslow, O. W. Brawley, and R. Wender, “Cancer screening in the united states, 2015: a review of current american cancer society guidelines and current issues in cancer screening,” CA: a cancer journal for clinicians, vol. 65, no. 1, pp. 30–54, 2015. [20] W. A. Berg, L. Gutierrez, M. S. NessAiver, W. B. Carter, M. Bhargavan, R. S. Lewis, and O. B. Ioffe, “Diagnostic accuracy of mammography, clinical examina- tion, us, and mr imaging in preoperative assessment of breast cancer,” Radiology, vol. 233, no. 3, pp. 830–849, 2004. 145 [21] P. M. Meaney, M. W. Fanning, T. Raynolds, C. J. Fox, Q. Fang, C. A. Kogel, S. P. Poplack, and K. D. Paulsen, “Initial clinical experience with microwave breast imaging in women with normal mammography,” Academic radiology, vol. 14, no. 2, pp. 207–218, 2007. [22] E. C. Fear, X. Li, S. C. Hagness, and M. A. Stuchly, “Confocal microwave imag- ing for breast cancer detection: Localization of tumors in three dimensions,” IEEE Transactions on Biomedical Engineering, vol. 49, no. 8, pp. 812–822, 2002. [23] E. Fear and M. Stuchly, “Microwave detection of breast cancer,” IEEE Trans- actions on Microwave Theory and Techniques, vol. 48, no. 11, pp. 1854–1863, 2000. [24] E. J. Bond, X. Li, S. C. Hagness, and B. D. Van Veen, “Microwave imaging via space-time beamforming for early detection of breast cancer,” IEEE Transactions on Antennas and Propagation, vol. 51, no. 8, pp. 1690–1705, 2003. [25] N. K. Nikolova, “Microwave biomedical imaging,” Wiley Encyclopedia of Elec- trical and Electronics Engineering, 2014. [26] W. C. Chew, Waves and fields in inhomogeneous media. IEEE press New York, 1995, vol. 522. [27] M. Moghaddam and W. C. Chew, “Study of some practical issues in inversion with the born iterative method using time-domain data,” IEEE Transactions on Antennas and Propagation, vol. 41, no. 2, pp. 177–184, 1993. [28] W. Chew and Y . Wang, “Reconstruction of two-dimensional permittivity distri- bution using the distorted born iterative method,” Medical Imaging, IEEE Trans- actions on, vol. 9, no. 2, pp. 218–225, 1990. [29] R. F. Remis and P. Van den Berg, “On the equivalence of the newton-kantorovich and distorted born methods,” Inverse Problems, vol. 16, no. 1, p. L1, 2000. [30] M. Haynes, S. Clarkson, and M. Moghaddam, “Electromagnetic inverse scat- tering algorithm and experiment using absolute source characterization,” IEEE Transactions on Antennas and Propagation, vol. 60, no. 4, pp. 1854–1867, 2012. [31] J. D. Shea, P. Kosmas, S. C. Hagness, and B. D. Van Veen, “Three-dimensional microwave imaging of realistic numerical breast phantoms via a multiple- frequency inverse scattering technique,” Medical physics, vol. 37, no. 8, pp. 4210–4226, 2010. [32] T. M. Habashy, M. L. Oristaglio, and T. Adrianus, “Simultaneous nonlinear re- construction of two-dimensional permittivity and conductivity,” Radio Science, vol. 29, no. 04, pp. 1101–1118, 1994. 146 [33] P. M. Van den Berg and K. F. Haak, “Profile inversion by error reduction in the source type integral equations,” Wavefields and Reciprocity, pp. 87–98, 1996. [34] Y . Zhong, M. Lambert, D. Lesselier, and X. Chen, “A new integral equation method to solve highly nonlinear inverse scattering problems,” IEEE Transac- tions on Antennas and Propagation, vol. 64, no. 5, pp. 1788–1799, 2016. [35] P. M. Van Den Berg and R. E. Kleinman, “A contrast source inversion method,” Inverse problems, vol. 13, no. 6, p. 1607, 1997. [36] A. Zakaria, C. Gilmore, and J. LoVetri, “Finite-element contrast source inversion method for microwave imaging,” Inverse Problems, vol. 26, no. 11, p. 115010, 2010. [37] A. Abubakar, P. M. Van den Berg, and J. J. Mallorqui, “Imaging of biomedical data using a multiplicative regularized contrast source inversion method,” IEEE Transactions on Microwave Theory and Techniques, vol. 50, no. 7, pp. 1761– 1771, 2002. [38] C. Gilmore, P. Mojabi, and J. LoVetri, “Comparison of an enhanced distorted born iterative method and the multiplicative-regularized contrast source inversion method,” Antennas and Propagation, IEEE Transactions on, vol. 57, no. 8, pp. 2341–2351, 2009. [39] M. Pastorino, A. Massa, and S. Caorsi, “A microwave inverse scattering technique for image reconstruction based on a genetic algorithm,” IEEE transactions on Instrumentation and Measurement, vol. 49, no. 3, pp. 573–578, 2000. [40] M. Donelli, G. Franceschini, A. Martini, and A. Massa, “An integrated multiscal- ing strategy based on a particle swarm algorithm for inverse scattering problems,” IEEE Transactions on Geoscience and Remote Sensing, vol. 44, no. 2, pp. 298– 312, 2006. [41] P. Rocca, M. Benedetti, M. Donelli, D. Franceschini, and A. Massa, “Evolution- ary optimization as applied to inverse scattering problems,” Inverse Problems, vol. 25, no. 12, p. 123003, 2009. [42] I. T. Rekanos, “Neural-network-based inverse-scattering technique for online mi- crowave medical imaging,” IEEE transactions on magnetics, vol. 38, no. 2, pp. 1061–1064, 2002. [43] M. El-Shenawee and E. L. Miller, “Spherical harmonics microwave algorithm for shape and location reconstruction of breast cancer tumor,” IEEE transactions on medical imaging, vol. 25, no. 10, pp. 1258–1271, 2006. 147 [44] O. Dorn, E. L. Miller, and C. M. Rappaport, “A shape reconstruction method for electromagnetic tomography using adjoint fields and level sets,” Inverse prob- lems, vol. 16, no. 5, p. 1119, 2000. [45] O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” Inverse Problems, vol. 22, no. 4, p. R67, 2006. [46] D. Colton and L. P¨ aiv¨ arinta, “The uniqueness of a solution to an inverse scat- tering problem for electromagnetic waves,” Archive for rational mechanics and analysis, vol. 119, no. 1, pp. 59–70, 1992. [47] A. I. Nachman, “Global uniqueness for a two-dimensional inverse boundary value problem,” Annals of Mathematics, pp. 71–96, 1996. [48] O. M. Bucci and G. Franceschetti, “On the degrees of freedom of scattered fields,” IEEE transactions on Antennas and Propagation, vol. 37, no. 7, pp. 918–926, 1989. [49] O. Bucci and T. Isernia, “Electromagnetic inverse scattering: Retrievable infor- mation and measurement strategies,” Radio Science, vol. 32, no. 6, pp. 2123– 2137, 1997. [50] M. Klemm, I. J. Craddock, J. A. Leendertz, A. Preece, and R. Benjamin, “Radar- based breast cancer detection using a hemispherical antenna arrayexperimental results,” IEEE transactions on antennas and propagation, vol. 57, no. 6, pp. 1692–1704, 2009. [51] Y . Wang and W. Chew, “An iterative solution of the two-dimensional electromag- netic inverse scattering problem,” International Journal of Imaging Systems and Technology, vol. 1, no. 1, pp. 100–108, 1989. [52] A. Tarantola, Inverse problem theory and methods for model parameter estima- tion. siam, 2005. [53] S. Sibisi, “Regularization and inverse problems,” in Maximum entropy and Bayesian methods. Springer, 1989, pp. 389–396. [54] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of inverse problems. Springer, 1996, vol. 375. [55] L. Tenorio, “Statistical regularization of inverse problems,” SIAM review, vol. 43, no. 2, pp. 347–366, 2001. [56] R. A. Renaut, I. Hnˇ etynkov´ a, and J. Mead, “Regularization parameter estimation for large-scale tikhonov regularization using a priori information,” Computational Statistics & Data Analysis, vol. 54, no. 12, pp. 3430–3445, 2010. 148 [57] L. Reichel and H. Sadok, “A new l-curve for ill-posed problems,” Journal of Computational and Applied Mathematics, vol. 219, no. 2, pp. 493–508, 2008. [58] C. Estatico, M. Pastorino, and A. Randazzo, “A novel microwave imaging ap- proach based on regularization in banach spaces,” Antennas and Propagation, IEEE Transactions on, vol. 60, no. 7, pp. 3373–3381, 2012. [59] R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society. Series B (Methodological), pp. 267–288, 1996. [60] E. Candes, R. Baraniuk, J. Romberg et al., Compressive Sensing Resources, http://www.dsp.rice.edu/cs. [Online]. Available: http://www.dsp. rice.edu/cs [61] R. G. Baraniuk, “Compressive sensing,” IEEE Signal Process. Mag., vol. 24, no. 4, 2007. [62] D. L. Donoho, “Compressed sensing,” Information Theory, IEEE Transactions on, vol. 52, no. 4, pp. 1289–1306, 2006. [63] L. S´ egui, “Sur un probl` eme inverse en diffraction d’ondes-identification des per- mittivit´ es complexes d’un mat´ eriau ` a partir de donn´ ees en champ proche.” 2000. [64] E. A. Marengo, “Compressive sensing and signal subspace methods for inverse scattering including multiple scattering,” in Proc. IEEE Int. Workshop Comput. Adv. Multi-Sensor Adapt. Process, 2008, pp. 7–11. [65] E. Marengo, R. Hernandez, Y . Citron, F. Gruber, M. Zambrano, and H. Lev- Ari, “Compressive sensing for inverse scattering,” XXIX URSI General Assembly, Chicago, Illinois, 2008. [66] G. Oliveri, P. Rocca, and A. Massa, “A bayesian-compressive-sampling-based in- version for imaging sparse scatterers,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 10, pp. 3993–4006, 2011. [67] L. Pan, X. Chen, and S. P. Yeo, “A compressive-sensing-based phaseless imaging method for point-like dielectric objects,” IEEE Trans. Antennas Propag, vol. 60, no. 11, p. 5472, 2012. [68] L. Poli, G. Oliveri, and A. Massa, “Microwave imaging within the first-order born approximation by means of the contrast-field bayesian compressive sens- ing,” IEEE Trans. Antennas Propag, vol. 60, no. 6, pp. 2865–2879, 2012. [69] L. Poli, G. Oliveri, P. Rocca, and A. Massa, “Bayesian compressive sensing ap- proaches for the reconstruction of two-dimensional sparse scatterers under te il- luminations,” Geoscience and Remote Sensing, IEEE Transactions on, vol. 51, no. 5, pp. 2920–2936, 2013. 149 [70] L. Poli, G. Oliveri, F. Viani, and A. Massa, “Mt–bcs-based microwave imaging approach through minimum-norm current expansion,” Antennas and Propaga- tion, IEEE Transactions on, vol. 61, no. 9, pp. 4722–4732, 2013. [71] F. Viani, L. Poli, G. Oliveri, F. Robol, and A. Massa, “Sparse scatterers imaging through approximated multitask compressive sensing strategies,” Microwave and Optical Technology Letters, vol. 55, no. 7, pp. 1553–1558, 2013. [72] L. Poli, G. Oliveri, and A. Massa, “Imaging sparse metallic cylinders through a local shape function bayesian compressive sensing approach,” JOSA A, vol. 30, no. 6, pp. 1261–1272, 2013. [73] L. Guo and A. M. Abbosh, “Microwave imaging of nonsparse domains using born iterative method with wavelet transform and block sparse bayesian learning,” IEEE Trans. Antennas Propag, vol. 63, no. 11, pp. 4877–4888, 2015. [74] Y . Wang and W. Chew, “An iterative solution of the two-dimensional electromag- netic inverse scattering problem,” International Journal of Imaging Systems and Technology, vol. 1, no. 1, pp. 100–108, 1989. [75] U. Khankhoje, J. van Zyl, and T. Cwik, “Computation of radar scattering from heterogeneous rough soil using the finite-element method,” IEEE Transactions on Geoscience and Remote Sensing, vol. 51, no. 6, pp. 3461–3469, Jun. 2013. [76] E. J. Cand` es and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2, pp. 21–30, 2008. [77] P. Shah, U. K. Khankhoje, and M. Moghaddam, “Joint l1-l2 regularization for in- verse scattering,” in Antennas and Propagation Society International Symposium (APSURSI), 2014 IEEE. IEEE, 2014, pp. 868–869. [78] D. P. Bertsekas, Nonlinear programming. Athena Scientific, 1999. [79] E. Van Den Berg and M. P. Friedlander, “Probing the pareto frontier for basis pursuit solutions,” SIAM Journal on Scientific Computing, vol. 31, no. 2, pp. 890–912, 2008. [80] E. G. Birgin, J. M. Mart´ ınez, and M. Raydan, “Nonmonotone spectral projected gradient methods on convex sets,” SIAM Journal on Optimization, vol. 10, no. 4, pp. 1196–1211, 2000. [81] S. Boyd and L. Vandenberghe, Convex optimization. Cambridge university press, 2009. [82] W. H. Press, Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007. 150 [83] A. P. Ruszczy´ nski, Nonlinear optimization. Princeton university press, 2006, vol. 13. [84] P. C. Hansen, “Regularization, gsvd and truncatedgsvd,” BIT Numerical Mathe- matics, vol. 29, no. 3, pp. 491–504, 1989. [85] J. Chung, J. G. Nagy, and D. P. OLeary, “A weighted gcv method for lanczos hybrid regularization,” Electronic Transactions on Numerical Analysis, vol. 28, pp. 149–167, 2008. [86] J. L. Mead and R. A. Renaut, “A newton root-finding algorithm for estimating the regularization parameter for solving ill-conditioned least squares problems,” Inverse Problems, vol. 25, no. 2, p. 025002, 2009. [87] C. R. Rao, Linear statistical inference and its applications. John Wiley & Sons, 2009, vol. 22. [88] P. Shah, U. K. Khankhoje, and M. Moghaddam, “Inverse scattering using a joint norm-based regularization,” IEEE Transactions on Antennas and Propagation, vol. 64, no. 4, pp. 1373–1384, 2016. [89] M. R. Hajihashemi and M. El-Shenawee, “Shape reconstruction using the level set method for microwave applications,” IEEE Antennas wireless propag. lett., vol. 7, no. 4, pp. 92–96, 2008. [90] D. A. Woten, M. R. Hajihashemi, A. M. Hassan, and M. El-Shenawee, “Experi- mental microwave validation of level set reconstruction algorithm,” IEEE Trans. Antennas Propag.,, vol. 58, no. 1, pp. 230–233, 2010. [91] M. R. Hajihashemi and M. El-Shenawee, “High performance computing for the level-set reconstruction algorithm,” Journal of Parallel and Distributed Comput- ing, vol. 70, no. 6, pp. 671–679, 2010. [92] ——, “TE versus TM for the shape reconstruction of 2-d pec targets using the level-set algorithm,” IEEE Transactions on Geoscience and Remote Sensing, vol. 48, no. 3, pp. 1159–1168, 2010. [93] M. R. Eskandari, M. Dehmollaian, and R. Safian, “Simultaneous microwave imaging and parameter estimation using modified level-set method,” IEEE Trans- actions on Antennas and Propagation, vol. 64, no. 8, pp. 3554–3564, 2016. [94] J. A. Sethian, “A fast marching level set method for monotonically advancing fronts,” Proceedings of the National Academy of Sciences, vol. 93, no. 4, pp. 1591–1595, 1996. 151 [95] F. Santosa, “A level-set approach for inverse problems involving obstacles,” ESAIM: Control, Optimisation and Calculus of Variations, vol. 1, pp. 17–33, 1996. [96] S. O. R. Fedkiw, Level set methods and dynamic implicit surfaces. Springer, 2003. [97] M. El-Shenawee, O. Dorn, and M. Moscoso, “An adjoint-field technique for shape reconstruction of 3-d penetrable object immersed in lossy medium,” IEEE Transactions on Antennas and Propagation, vol. 57, no. 2, pp. 520–534, 2009. [98] T. J. Colgan, S. C. Hagness, and B. D. Van Veen, “A 3-d level set method for mi- crowave breast imaging,” IEEE Transactions on Biomedical Engineering, vol. 62, no. 10, pp. 2526–2534, 2015. [99] M. R. Eskandari, R. Safian, and M. Dehmollaian, “Three-dimensional near-field microwave imaging using hybrid linear sampling and level set methods in a medium with compact support,” IEEE Transactions on Antennas and Propaga- tion, vol. 62, no. 10, pp. 5117–5125, 2014. [100] A. M. Hassan, T. C. Bowman, and M. El-Shenawee, “Efficient microwave imag- ing algorithm based on hybridization of the linear sampling and level set meth- ods,” IEEE Trans. Antennas Propag.,, vol. 61, no. 7, pp. 3765–3773, 2013. [101] P. Shah and M. Moghaddam, “Fast level set based method for high contrast mi- crowave imaging,” in PIERS Proceedings. PIERS, 2015, pp. 253–257. [102] L. A. Vese and T. F. Chan, “A multiphase level set framework for image segmen- tation using the mumford and shah model,” International journal of computer vision, vol. 50, no. 3, pp. 271–293, 2002. [103] Y . Wang, Z. Luo, Z. Kang, and N. Zhang, “A multi-material level set-based topol- ogy and shape optimization method,” Computer Methods in Applied Mechanics and Engineering, vol. 283, pp. 1570–1586, 2015. [104] M. Cui, H. Chen, and J. Zhou, “A level-set based multi-material topology opti- mization method using a reaction diffusion equation,” Computer-Aided Design, vol. 73, pp. 41–52, 2016. [105] S. Zhou and M. Y . Wang, “Multimaterial structural topology optimization with a generalized cahn–hilliard model of multiphase transition,” Structural and Multi- disciplinary Optimization, vol. 33, no. 2, p. 89, 2007. [106] T. Gao and W. Zhang, “A mass constraint formulation for structural topology op- timization with multiphase materials,” International Journal for Numerical Meth- ods in Engineering, vol. 88, no. 8, pp. 774–796, 2011. 152 [107] O. Dorn and D. Lesselier, “Level set methods for inverse scatteringsome recent developments,” Inverse Problems, vol. 25, no. 12, p. 125001, 2009. [108] A. Litman, “Reconstruction by level sets of n-ary scattering obstacles,” Inverse Problems, vol. 21, no. 6, p. S131, 2005. [109] N. Irishina, M. Moscoso, and O. Dorn, “Microwave imaging for early breast cancer detection using a shape-based strategy,” IEEE Transactions on Biomedical Engineering, vol. 56, no. 4, pp. 1143–1153, 2009. [110] N. Irishina, D. ´ Alvarez, O. Dorn, and M. Moscoso, “Structural level set inversion for microwave breast screening,” Inverse Problems, vol. 26, no. 3, p. 035015, 2010. [111] J. Richmond, “Scattering by a dielectric cylinder of arbitrary cross section shape,” IEEE Transactions on Antennas and Propagation, vol. 13, no. 3, pp. 334–341, 1965. [112] L. Sun and W. Chew, “A novel formulation of the volume integral equation for electromagnetic scattering,” Waves in Random and Complex Media, vol. 19, no. 1, pp. 162–180, 2009. [113] J.-M. Jin, Theory and computation of electromagnetic fields. John Wiley & Sons, 2011. [114] F. A. Duck, Physical properties of tissues: a comprehensive reference book. Academic Press, 1990. [115] C. Li, C. Xu, C. Gui, and M. D. Fox, “Distance regularized level set evolution and its application to image segmentation,” IEEE Trans. Imag. Proc.,, vol. 19, no. 12, pp. 3243–3254, 2010. [116] M. J. Burfeindt, T. J. Colgan, R. O. Mays, J. D. Shea, N. Behdad, B. D. Van Veen, and S. C. Hagness, “Mri-derived 3-d-printed breast phantom for microwave breast imaging validation,” IEEE antennas and wireless propagation letters, vol. 11, pp. 1610–1613, 2012. [117] A. Massa, P. Rocca, and G. Oliveri, “Compressive sensing in electromagnetics-a review,” IEEE Antennas and Propagation Magazine, vol. 57, no. 1, pp. 224–238, 2015. [118] E. J. Cand` es and B. Recht, “Exact matrix completion via convex optimization,” Foundations of Computational Mathematics, vol. 9, no. 6, pp. 717–772, 2009. [119] Y . Peng, A. Ganesh, J. Wright, W. Xu, and Y . Ma, “Rasl: Robust alignment by sparse and low-rank decomposition for linearly correlated images,” IEEE Trans- actions on Pattern Analysis and Machine Intelligence, vol. 34, no. 11, pp. 2233– 2246, 2012. 153 [120] P. Shah and M. D. Gupta, “Simultaneous registration and segmentation by l1 min- imization,” in International Workshop on Machine Learning in Medical Imaging. Springer, 2012, pp. 128–135. [121] M. Schweiger, O. Dorn, and S. R. Arridge, “3-d shape and contrast reconstruction in optical tomography with level sets,” in Journal of Physics: Conference Series, vol. 124, no. 1. IOP Publishing, 2008, p. 012043. [122] M. Schweiger, O. Dorn, A. Zacharopoulos, I. Nissil¨ a, and S. Arridge, “3d level set reconstruction of model and experimental data in diffuse optical tomography,” Optics express, vol. 18, no. 1, pp. 150–164, 2010. [123] K. Prieto and O. Dorn, “Sparsity and level set regularization for diffuse optical tomography using a transport model in 2d,” Inverse Problems, vol. 33, no. 1, p. 014001, 2016. [124] R. Ferray´ e, J.-Y . Dauvignac, and C. Pichot, “An inverse scattering method based on contour deformations by means of a level set method using frequency hopping technique,” IEEE Transactions on Antennas and Propagation, vol. 51, no. 5, pp. 1100–1113, 2003. [125] P. Shah and M. Moghaddam, “Multi-material recovery using fast level set based formulation for microwave imaging (under-review),” IEEE Transactions on An- tennas and Propagation, vol. 64, no. 4, pp. 1373–1384, 2017. [126] M. R. Hajihashemi and M. El-Shenawee, “Inverse scattering of three-dimensional pec objects using the level-set method,” Progress In Electromagnetics Research, vol. 116, pp. 23–47, 2011. [127] ——, “Level set algorithm for shape reconstruction of non-overlapping three- dimensional penetrable targets,” IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 1, pp. 75–86, 2012. [128] H. N. Patel and D. K. Ghodgaonkar, “3d level set based optimization of inverse scattering problem for microwave breast imaging,” in Microwave Bio Conference (IMBIOC), 2017 First IEEE MTT-S International. IEEE, 2017, pp. 1–4. [129] G. Chen, J. Stang, and M. Moghaddam, “A conformal fdtd method with accurate waveport excitation and s-parameter extraction,” IEEE Transactions on Antennas and Propagation, vol. 64, no. 10, pp. 4504–4509, 2016. [130] M. Haynes and M. Moghaddam, “Vector green’s function fors-parameter mea- surements of the electromagnetic volume integral equation,” IEEE Transactions on Antennas and Propagation, vol. 60, no. 3, pp. 1400–1413, 2012. [131] C. A. Balanis, Advanced engineering electromagnetics. Wiley New York, 1989, vol. 20. 154 [132] G. Chen, J. Stang, and M. Moghaddam, “Numerical vector green’s function for s-parameter measurement with waveport excitation,” IEEE Transactions on An- tennas and Propagation, 2017. [133] T. Goldstein and S. Osher, “The split bregman method for l1-regularized prob- lems,” SIAM journal on imaging sciences, vol. 2, no. 2, pp. 323–343, 2009. [134] E. van den Berg and M. P. Friedlander, “Probing the pareto frontier for basis pursuit solutions,” SIAM Journal on Scientific Computing, vol. 31, no. 2, pp. 890–912, 2008. [Online]. Available: http://link.aip.org/link/?SCE/31/890 [135] M. Haynes, J. Stang, and M. Moghaddam, “Real-time microwave imaging of differential temperature for thermal therapy monitoring,” IEEE Transactions on Biomedical Engineering, vol. 61, no. 6, pp. 1787–1797, 2014. [136] E. J. Candes and Y . Plan, “A probabilistic and ripless theory of compressed sens- ing,” IEEE Trans. Inf. Theory, vol. 57, no. 11, pp. 7235–7254, 2011. [137] E. J. Candes and T. Tao, “Near-optimal signal recovery from random projections: Universal encoding strategies?” IEEE Trans. Inf. Theory, vol. 52, no. 12, pp. 5406–5425, 2006. [138] D. Gross, “Recovering low-rank matrices from few coefficients in any basis,” Information Theory, IEEE Transactions on, vol. 57, no. 3, pp. 1548–1566, 2011. [139] R. Kueng and D. Gross, “Ripless compressed sensing from anisotropic measure- ments,” Linear Algebra and its Applications, vol. 441, pp. 110–123, 2014. [140] O. Bucci and T. Isernia, “Electromagnetic inverse scattering: Retrievable infor- mation and measurement strategies,” Radio Science, vol. 32, no. 6, pp. 2123– 2137, 1997. [141] U. K. Khankhoje and K. Shah, “Spatial bandlimitedness of scattered electromag- netic fields,” arXiv preprint arXiv:1505.00886, 2015. [142] D. Needell and J. A. Tropp, “Cosamp: Iterative signal recovery from incomplete and inaccurate samples,” Applied and Computational Harmonic Analysis, vol. 26, no. 3, pp. 301–321, 2009. [143] J. Duchi, S. Shalev-Shwartz, Y . Singer, and T. Chandra, “Efficient projections onto the l 1-ball for learning in high dimensions,” in Proceedings of the 25th international conference on Machine learning. ACM, 2008, pp. 272–279. [144] M. D. Gupta, S. Kumar, and J. Xiao, “L1 projections with box constraints,” arXiv preprint arXiv:1010.0141, 2010. 155 List of Publications 1. P. Shah and M. Moghaddam; ”RIPless Compressive Sensing: A Different View for Ill-posed Problems” (in preparation) 2. P. Shah, G. Chen, J. Stang, and M. Moghaddam; ”A Total Variation Regularized 3D Level Set Method for Joint Contrast and Shape Recovery in Microwave Imag- ing,” (to be submitted to Inverse Problems) 3. P. Shah and M. Moghaddam; ”Multi-Material Recovery Using Fast Level Set Based Formulation for Microwave Imaging,” (minor revision, IEEE Transactions on Antennas and Propagation) 4. G. Chen, P. Shah, J. Stang, S. Cho, M. Moghaddam, E. Leuthardt; ”A Microwave Imaging System and Real-Time Image Reconstruction Algorithm for Intraoper- ative 3D Monitoring of Thermal Ablation Therapies,” Radiological Society of North America (2017) 5. P. Shah and M. Moghaddam; ”Matrix Norm Based Method for Recovery of High Contrast and Sparse Objects in Microwave Imaging,” IEEE International Sym- posium on Antennas and Propagation and USNC-URSI Radio Science Meeting (2017) 6. P. Shah and M. Moghaddam; ”Super Resolution for Microwave Imaging: A Deep Learning Approach,” IEEE International Symposium on Antennas and Propaga- tion and USNC-URSI Radio Science Meeting (2017) 7. A. Batista, P. Shah, G. Chen, and J Stang; ”Three Dimensional Level Set Method For Microwave Imaging,” USNC-URSI National Radio Science Meeting (2016) 8. P. Shah and M. Moghaddam; ”Theoretical Derivation of RIP-less Compressive Sensing for Inverse Scattering,” IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting (2016) 9. P. Shah, U. Khankhoje, and M. Moghaddam; ”Inverse Scattering using a Joint L1-L2 based Regularization,” IEEE Transactions on Antennas and Propagation (2016) 156 10. P. Shah and M. Moghaddam; ”Fast Level Set Based Method for High Contrast Microwave Imaging,” Progress in Electromagnetics Research Symposium (2015) 11. R. Akbar, P. Shah, and M. Moghaddam; ”Joint Radar-Radiometer Multi-Parameter Estimation of Surface Soil Moisture and Roughness,” IEEE International Geo- science and Remote Sensing Symposium (2015) 12. M. Haynes, J. Stang, P. Shah, and M. Moghaddam; ”Real-Time 3D Microwave Thermal Imaging of a Microwave Ablation Probe,” IEEE International Sympo- sium on Antennas and Propagation and USNC-URSI Radio Science Meeting (2015) 13. P. Shah, U. Khankhoje, and M. Moghaddam; ”Joint L1-L2 Regularization For Inverse Scattering,” IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting (2014) 157

Abstract (if available)

Abstract Imaging in the microwave regime of electromagnetic waves is a common tool for medical imaging, through-the-wall imaging, buried object detection, locating oil-gas in the Earth subsurface and many other industrial applications. The aim of the imaging is to estimate or map the electric properties such as permittivity and conductivity noninvasively. Many developments have taken place in both hardware development and in the image reconstruction aspect of microwave imaging. However, still, the reconstructed images have poor resolution and often low estimation accuracy of the electrical properties. The primary reason is the challenge of solving a highly non-linear inverse problem accurately. To address it, the first part of this thesis presents a detailed investigation on the effect of various constraints on the electrical properties of a 2D microwave imaging problem. In various simulation setups, different vector norms ($\ell_1$ and $\ell_2$), matrix norms (mix-norm and nuclear norm), bound constraints on permittivity and conductivity, as well as shape constraints have been applied separately or in some combination as constraints. The $\ell_1$, $\ell_2$, and bound constraint together can produce reconstructed images with better estimation of the electrical properties than conventional constraints when the scatterers are not considered strong. Shape constraint via a level set function can have an almost perfect estimation of the electrical properties if the properties are known a priori but their location in the imaging domain is unknown. The matrix norm-based constraints can recover sparse strong scatters in a small imaging domain with a very low computational burden. ❧ In the second part of the thesis, an experiment-ready 3D imaging setup is considered in simulation mode and a method is developed to solve two subproblems rather than solving in a single problem: shape estimation and the properties (contrast) estimation. The modular framework enables the freedom to choose any arbitrary constraints on the contrast in the cost function. Here, we demonstrate the applicability of a total-variation (TV) constraint for the contrast cost function. The reconstructed images show accurate estimates of object location, shape, and size while recovering the contrast very well even with limited information. ❧ In addition to investigating the different constraints, the applicability of compressive sensing to inverse scattering problems from a theoretical point of view is also investigated. The findings of the investigation are applied to derive new Lemmas, which broaden the applicability of compressive sensing theory from isotropic measurements to measurements that have a theoretical limit on the maximally possible independent number of measurements such as in inverse scattering problems. Then, the Lemmas are applied for a single source formulation in a 2D homogeneous domain, which show the possibility of reconstructing the image uniquely for the first-time for very sparse objects.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Solution of inverse scattering problems via hybrid global and local optimization

PDF

Active and passive microwave remote sensing synergy for geophysical parameter estimation

PDF

Optimization methods and algorithms for constrained magnetic resonance imaging

PDF

On the theory and applications of structured bilinear inverse problems to sparse blind deconvolution, active target localization, and delay-Doppler estimation

PDF

Sparse feature learning for dynamic subsurface imaging

PDF

Improved brain dynamic contrast enhanced MRI using model-based reconstruction

PDF

Electromagnetic scattering models for satellite remote sensing of soil moisture using reflectometry from microwave signals of opportunity

PDF

Shift-invariant autoregressive reconstruction for MRI

PDF

Novel optimization tools for structured signals recovery: channels estimation and compressible signal recovery

PDF

Physics-based bistatic radar scattering model for vegetated terrains in support of soil moisture retrieval from signals of opportunity

PDF

Correction, coregistration and connectivity analysis of multi-contrast brain MRI

PDF

Signal processing for channel sounding: parameter estimation and calibration

PDF

Measuring functional connectivity of the brain

PDF

Interaction and topology in distributed multi-agent coordination

PDF

Machine learning methods for 2D/3D shape retrieval and classification

PDF

The role of rigid foundation assumption in two-dimensional soil-structure interaction (SSI)

PDF

Topics in algorithms for new classes of non-cooperative games

PDF

Grid-based Vlasov method for kinetic plasma simulations

PDF

Brain connectivity in epilepsy

PDF

Magnetic induction-based wireless body area network and its application toward human motion tracking

Asset Metadata

Creator Shah, Pratik (author)

Core Title Multi-constrained inversion algorithms for microwave imaging

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 01/31/2018

Defense Date 12/08/2017

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag compressive sensing,electromagnetics,inverse problems,inverse scattering,level set methods,microwave imaging,OAI-PMH Harvest

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Moghaddam, Mahta (committee chair), Haldar, Justin (committee member), Stang, John (committee member), Wang, Chunming (committee member)

Creator Email pratik0828@gmail.com,pratiks@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-467275

Unique identifier UC11268216

Identifier etd-ShahPratik-5985.pdf (filename),usctheses-c40-467275 (legacy record id)

Legacy Identifier etd-ShahPratik-5985.pdf

Dmrecord 467275

Document Type Dissertation

Rights Shah, Pratik

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA