International Journal of Antennas and Propagation

Volume 2019, Article ID 9014969, 12 pages

https://doi.org/10.1155/2019/9014969

## A Discrete Dipole Approximation Solver Based on the COCG-FFT Algorithm and Its Application to Microwave Breast Imaging

^{1}Electrical Engineering Department, Chalmers University of Technology, 41296 Gothenburg, Sweden^{2}The Thayer School of Engineering, Dartmouth College, Hanover, NH 03755, USA

Correspondence should be addressed to Samar Hosseinzadegan; es.sremlahc@hramas

Received 26 February 2019; Revised 6 May 2019; Accepted 20 May 2019; Published 17 July 2019

Guest Editor: Sandra Costanzo

Copyright © 2019 Samar Hosseinzadegan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We introduce the discrete dipole approximation (DDA) for efficiently calculating the two-dimensional electric field distribution for our microwave tomographic breast imaging system. For iterative inverse problems such as microwave tomography, the forward field computation is the time limiting step. In this paper, the two-dimensional algorithm is derived and formulated such that the iterative conjugate orthogonal conjugate gradient (COCG) method can be used for efficiently solving the forward problem. We have also optimized the matrix-vector multiplication step by formulating the problem such that the nondiagonal portion of the matrix used to compute the dipole moments is block-Toeplitz. The computation costs for multiplying the block matrices times a vector can be dramatically accelerated by expanding each Toeplitz matrix to a circulant matrix for which the convolution theorem is applied for fast computation utilizing the fast Fourier transform (FFT). The results demonstrate that this formulation is accurate and efficient. In this work, the computation times for the direct solvers, the iterative solver (COCG), and the iterative solver using the fast Fourier transform (COCG-FFT) are compared with the best performance achieved using the iterative solver (COCG-FFT) in C++. Utilizing this formulation provides a computationally efficient building block for developing a low cost and fast breast imaging system to serve under-resourced populations.

#### 1. Introduction

The mortality rate due to the breast cancer in women worldwide has led numerous research groups to investigate early diagnosis programs. Along with other imaging methods such as X-ray computed tomography [1], positron emission tomography [2], and magnetic resonance imaging [3, 4], microwave imaging has been tested in multiple settings. In this context, microwave imaging is performed in four primary forms: radar, holography, thermoacoustic imaging, and tomography. The radar approaches have been studied in an array of simulation experiments and have advanced to several clinical tests [5–8]. Holography approaches have been primarily tested in simulation and phantom experiments [9, 10]. Thermoacoustic imaging work has been evaluated in several phantom experiments along with early clinical studies [11, 12]. In this report, we focus on microwave tomography (MWT) which has been tested in several breast imaging clinical trials and has provided relevant diagnostic information regarding diagnosis of cancer and monitoring of tumor progression during neoadjuvant chemotherapy [13, 14].

In spite of the demand and interest for microwave tomography, the computational costs of various algorithms have remained a primary obstacle in translation to real applications [15, 16]. The choice of three or two-dimensional imaging algorithms has considerable impact on the computation time and necessary hardware resources [17, 18]. Investigations into 3D imaging are considerably more common than that for 2D because of the superior measurement model match despite the associated costs [16, 19]. For many groups investigating microwave imaging, the amount of measurement data required is often a significant barrier to real 3D implementations [20, 21]. For most numerical techniques, solving the 3D imaging problem can be computationally expensive and requires use of multiprocessor computers working over many hours to even days to generate single images [16, 22]. While these 3D efforts are useful and necessary to advance the science of microwave imaging, the practical barriers to implementation, including measurement data costs and computation time, have greatly hindered its translation into the clinic and limited them primarily to simulation studies. Alternatively, 2D approaches have proved viable and have been demonstrated in numerous phantom and clinical studies [23–26]. In fact, largely due to continuing computer efficiency advances, 2D techniques are poised to be viable alternatives for conventional modalities in underresourced settings where cost and portability are significant concerns. In this context, reducing memory requirements and computation time for the 2D algorithm is an important concern.

2D imaging algorithms and system implementations are not new to the microwave imaging community. One notable example is the Semenov group when they were associated with the Carolinas Medical Center in Charlotte, North Carolina. In a relatively early work, they concluded that the reason their imaging technique did not work as well as desired was because of the mismatch between their 2D algorithm and the inherent 3D nature of the actual propagating fields [27]. From that experience, they surmised that they would achieve improved results by transitioning to 3D imaging to avoid the measurement/model mismatch. However, no proof of this being a general conclusion was given. Notwithstanding, significant technological advances have made 2D imaging feasible and attractive. These include: (1) the use of the lossy coupling medium to suppress unwanted multipath signal corruption [28], (2) the use of monopole antennas which can be positioned very close to the target and which we have demonstrated improving the overall images [29], and (3) the use of the log transform which improves the algorithm convergence behavior, eliminates the conversion to local minima, and makes a priori information unnecessary [30, 31]. As such, 2D imaging is poised for an expanded clinical role.

Microwave tomographic imaging algorithms require solving two problems—forward and inverse problems. The forward solver is by far the more computationally expensive part of an iterative image reconstruction algorithm and requires substantial attention to reduce its impact. Additionally, the need to perform many iterations to recover accurate images further motivates the requirement for improving computation time. The most common numerical methods used to solve the forward problem are the finite-difference time-domain (FDTD) [32], finite element method (FEM) [33], and volume integral equations such as the method of moments (MoM) [34]. Recently, the 3D discrete dipole approximation was used as a forward solver for 3D imaging and was found to be highly efficient regarding computational cost and accuracy [35]. We have also introduced the 2D DDA as a forward solver for 2D imaging [36]. As part of that paper, comparisons were performed between the 2D DDA results and those modeled using COMSOL Multiphysics along with comparisons with actual measurements. The agreement was quite good for both. Each numerical method provides important advantages that can be exploited depending on the circumstances. For instance, the finite element method is particularly well suited for situations where nonuniform shaped scatterers are involved. Likewise, the uniform grid formulation of the FDTD problem facilitates fast computation times. While the DDA can be utilized on both uniform and nonuniform grid configurations, as will be demonstrated in this paper, the computational advantages are most significant for the uniform setting. In addition, the associated uniform setting is optimal when the objects in the domain are mostly dielectric in nature. While it would be possible to incorporate metallic scatterers in the domain, for an accurate representation, extra dipoles would need to be deployed on the boundaries of the objects. This would immediately preclude the uniform grid representation. Conveniently, the imaging system developed at Dartmouth [26] essentially provides this feature. The array of monopole antennas is naturally made of metal; however, their radar cross section is sufficiently small that they only slightly perturb the field distribution when another antenna is radiating. In addition, the use of a very lossy coupling medium (glycerin and water mixtures) further dampens any perturbations from scattering off the antennas.

In this work, the two-dimensional discrete dipole approximation (2D-DDA) for calculating the electric field distribution for our microwave imaging system is proposed. The iterative solver for 2D-DDA has the potential to significantly improve the computational speed. A conjugate gradient based method, i.e., the conjugate orthogonal conjugate gradient method (COCG), is used for which the computational cost of the COCG is remarkably reduced when incorporating the fast Fourier transform (FFT). This is made possible because the coefficient matrix for the 2D-DDA is complex, symmetric, and block-Toeplitz after removal of the main diagonal and enables the possibility for employing the FFT after expansion of the block matrices to circulant form. The computation times for the direct and iterative solvers are calculated and have been investigated in this comparison using both MATLAB and C++ implementations. It is useful to compare performance both with an interpretive language such as MATLAB and a classic compiler-based code. While the interpretive code struggles computation time-wise with constructs such as loops, it contains highly optimized matrix operations which can often overcome such disadvantages. These examinations show that the computation time for the 2D-DDA is significantly decreased in the COCG-FFT approach and that the best performance is achieved in C++ using an open source C++ package, FFTW, for fast Fourier transform calculations [37].

#### 2. Derivations

In this section, we formulate the 2D-DDA and discuss possible computational efficiency for it as a forward solver of the reconstruction algorithms.

##### 2.1. The 2D-DDA for the Forward Problem

The three-dimensional discrete dipole approximation (3D-DDA) has been widely used for calculation of scattering and absorption properties caused by an external electromagnetic field [38–40]. In the volume integral equations for techniques such as the discrete dipole approximation, an arbitrary geometry is assumed to be a union of small volumes such that [41]. For the purpose of microwave tomography, the electromagnetic field distribution is expressed in the imaging domain with the forward solver, the discrete dipole approximation. The total electric field at a point in the imaging domain iswhere and are the incident and scattered electric fields, respectively. The term represents the electric field propagation due to a waveguide or an antenna for a homogeneous domain. The scattered electric field is the electric field caused by scatterers in the domain. The Helmholtz equation solution in form of (1) is usually written as [42]where is the dyadic Green’s function, is wavenumber at , and is the background wavenumber. The wavenumbers are expressed in terms of the material property and frequency as . For the DDA, the idea is to approximate the total electric field, on the right-hand side of (1) such that we assume that the forward model zone consists of multiple number of dipoles. Each dipole represents the macroscopic field at the position of the dipole.

On the macroscopic level, the total electric field is proportional to the polarization viawhere is the permittivity of a vacuum and is the electric susceptibility. The term in Equation (2) is related to such thatOn the microscopic level, the polarization field is related to dipole moments of individual molecules. Each individual molecule is affected by a local electric field and gets correspondingly polarized to exhibit dipole moment . The total contribution from all molecules in a unit volume is defined as the polarization P: Here is number of molecules per unit volume. The polarizability, , expresses the relationship between the local and macroscopic field. In multiple applications, this interaction has been modeled by the Clausius-Mossotti relation. In the following section, we discuss the polarizability for microwave breast imaging systems.

*(1) Constant Polarizability (**) for the 2D DDA*. The most common and well-known molecular polarizability, , for the 3D-DDA is the Clausius-Mossotti relationship [43] and is defined for a sphere aswhere and are the relative permittivity and volume of the sphere. The constant term has different formulations depending on the volume and concentration of the medium [41, 44, 45]. Previously, Grzegorczyk et al. [35] used the following 3D Clausius-Mossotti relation for the microwave imaging system at Dartmouth:where and are the complex permittivities of the background and inclusion, respectively.

Since our imaging system shown in Figure 1 consists of a tank filled with a liquid mixture of glycerin and water with different polarization behavior on the molecular level, we select a more general form of the Clausius-Mossotti relationship, the Maxwell-Garnett formula. The main advantage of the Maxwell-Garnett model is that it is also valid for composite media which is the normal situation for many applications including our microwave imaging system. The Maxwell-Garnett formula can be expressed asIn (8), denotes the composite medium consisting of media with permittivities, . The coefficients, , relate the volume of the inclusion to its concentration which is defined as where , , and are volume, mass, and the molecular weight of the inclusion, respectively. The general formula for the dimensionless in m-dimensions is [46]Based on (8) and (9), the two-dimensional in the form of the Maxwell-Garnett formula is expressed asThe concentration coefficients are modified based on the area of the inclusion for the 2D case.