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

Hosseinzadegan, Samar; Fhager, Andreas; Persson, Mikael; Meaney, Paul

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

International Journal of Antennas and Propagation

On this page

Abstract Introduction Conclusion Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Special Issue

Bioengineering Applications of Electromagnetic Wave Propagation

View this Special Issue

Research Article | Open Access

Volume 2019 | Article ID 9014969 | 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

Samar Hosseinzadegan,¹Andreas Fhager,¹Mikael Persson,¹and Paul Meaney^1,2

Guest Editor: Sandra Costanzo

Received26 Feb 2019

Revised06 May 2019

Accepted20 May 2019

Published17 Jul 2019

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.

(2) Incident Electric Field () for the Microwave Breast Imaging. The microwave imaging system at Chalmers University of Technology has a circular array of monopole antennas acting as both transmitters and receivers (Figure 1). The transmitting antenna is modeled via the electric field distribution caused by an electrical line source (ELS). In this formulation, the 2D ELS is written in the form ofwhere , , and are the current amplitude, operating frequency, and free-space permeability, respectively. The term is the distance of a dipole located at position from an antenna position .

2.2. Two-Dimensional Discrete Dipole Approximation as a System of Equations

The governing equation for the two-dimensional discrete dipole approximation is [36, 39]:Inserting (5) into (12) yieldsFor the 2D-DDA, is the scalar Green’s function for the 2D Helmholtz equation and describes the interaction between two dipoles located at and . The relationship is given bywhere is the zero-order Hankel function of the second kind. Discretizing the forward model zone into subdomains, (Figure 2), transforms Equation (13) into its matrix format ofOne problematic aspect of this formulation is that some terms on the diagonal can approach infinity when any one of the ’s goes to zero. This happens when the permittivity at a dipole is exactly the same as that of the background (see Equation (7)). One way to eliminate this problem is to multiply both sides of the equation by . However, this has the unintended consequence of placing the quantities in the off-diagonal terms of the matrix. As will be shown later, we exploit the fact that the off-diagonal terms are only functions of distance and the background permittivity which allows these portions of the matrix to become block Toeplitz and subsequently allows for important optimizations. In fact, for our imaging configuration, we regularly encounter instances where approaches zero. Because of this conflict, we assume a uniform background medium outside of the dipoles where the properties are set to . Simultaneously, for our actual measurement system, we utilize a uniform coupling bath everywhere outside the imaging domain, . Since our forward model zone is a uniform grid, the properties at all dipoles outside of the imaging domain (a circle in this case) are set to . Ideally would be set to to eliminate all reflections at the grid boundary. However, allowing the ’s to go to zero would make the diagonal terms go to infinity. To overcome this, we artificially set and to be slightly different. This essentially sets many diagonal terms to a large number. In addition, it also implies that there will be nonsignificant signal reflections at the outer boundary of the dipole domain. For the latter challenge, we are fortunate in that our imaging system only uses a highly lossy coupling bath ensuring that waves reflecting back into the DDA and imaging domains will be sufficiently attenuated such that they effectively have no impact on the forward solution. The primary trade-off in this situation is between the level of reflections that can be tolerated and whether the accuracy of the forward solution is adversely affected because the condition number of the matrix in (15) becomes too great. An analysis of this issue is presented in Section 3.1.

2.3. Computational Efficient Implementation

The task of computing the electric field distribution utilizing the DDA involves two primary steps: (a) solving the dipole moments (P) and (b) multiplying P by the matrix in (15) with the diagonal, , terms removed to compute the fields. Because the matrix in (15) is full, the computational order can be as high as when using the standard Gaussian elimination technique [47]. As gets large, this becomes prohibitive and alternatives such as iterative solvers become attractive. In these cases, because the time limiting step involves repetitive multiplication of a length vector by an matrix, the computational cost is , where is the number of iterations required for the solution to converge. For this implementation, we exploit the special nature of the matrix derived for the DDA which ultimately allows us to employ the FFT as part of the matrix-vector multiplication to dramatically reduce the computation time within each iteration.

Different iterative solvers have been suggested with most derived from the Krylov subspace methods [39, 48, 49]. Among these methods the conjugate gradient method (CG) is the most popular. In the derivation of the CG algorithm, symmetry and positive definiteness (SPD) are assumed; however, if one of these conditions is not satisfied for the system of equations, the algorithm does not converge to a solution. In our case, the matrix in (15) is symmetric since the off-diagonal terms are only a function of the background medium and the distances between dipoles and . Unfortunately, the matrix does not meet the requirement of positive definiteness. Alternatively, there are a class of Krylov subspace methods that have been suggested for these conditions [50–57]. For our problem, we have chosen the conjugate orthogonal conjugate gradient method (COCG); it requires one time matrix-vector multiplication per iteration, but similar alternatives such as the conjugate orthogonal conjugate residual (COCR) can also be used with equivalent computational costs [57, 58].

Examination of (15) provides important insights into how best to exploit the COCG method. By breaking the matrix in Equation (15) into its diagonal and off-diagonal components, the left-hand side of (15) can be rewritten asThe first important observation is that the matrix-vector multiplication in the first term can easily be reduced to an vector-vector multiplication. For the right-hand matrix, multiplication times the vector P is normally an operation. However, in this situation, the right-hand side matrix, , is block Toeplitz. That is, each of the (where ) portions of the matrix are themselves Toeplitz matrices (Figure 3). To complete the full matrix-vector multiplication of , it can be performed by multiplying the smaller Toeplitz matrices times the appropriate portion of and summing the results together afterwards.

Toeplitz matrices are unique in that they are not necessarily symmetric, but they have regular repetition of the individual coefficients, for example:For our case, the Toeplitz matrices are symmetric and are closely related to circulant matrices for which there are highly optimized means for performing matrix-vector multiplications (i.e., the convolution theorem). In this case, the symmetric Toeplitz matrix can be padded both in columns and rows (Figure 4) to produce a circulant matrix:where the upper left-hand side, matrix was the original block-Toeplitz matrix. By transforming the Toeplitz matrix to this form, several key observations can be made. First, since each row is identical to the next except for a modulo shift, only one row of the matrix needs to be computed and stored. This is a substantial savings in terms of memory requirements. Secondly, because the matrix is now circulant, the convolution theorem can be used to multiply it by the associated portion of the vector (in this case the vector also needs to be padded similarly to the matrix). Taking the inverse FFT of the product of the FFT’s of the first row of the matrix times the associated portion of the vector is equivalent to multiplying the whole matrix times the vector.

3. Computational Results

In this section, we discuss two different approaches to optimize the forward solver of the reconstruction algorithms. Additionally, the effects of background medium dielectric properties are investigated and its impact on the electric field distribution is presented.

3.1. Study of the Field Distributions as a Function of Permittivity Difference

As mentioned in Section 2.2, one challenge in using Equation (15) is the treatment of the diagonal of the system of equations coefficient matrix. As the terms tend to zero, implying that the electrical properties of background liquid are the same as those of the dipoles, the corresponding diagonal terms go to infinity. For medical microwave imaging systems, a coupling liquid reducing the contrast between skin tissues and its surrounding is used. Utilizing a coupling liquid that has identical dielectric properties to that of the outside results in a matrix with infinity values on its diagonal. To accommodate this challenge, we set the properties of the dipoles within the grid to something slightly different than that for the background. This problem motivates us to study the range of acceptable permittivity values for its surrounding. For this analysis, the permittivities of the background and coupling medium are denoted as and , respectively.

In this situation, we set the permittivity values for the grid region directly surrounding the imaging domain to and that for the surrounding background to a value slightly lower than that. We compute the electric field distributions for a single, large inclusion ( = 40, and = 1.0 S/m) and the bath properties of = 22 and = 1.0 S/m. The imaging domain consists of 6561 dipoles. The antenna array is located on a circle of diameter 15.2 cm and the target is positioned at the origin. In this situation, the background permittivity is varied from 21.9999 down to 21 where the difference between and is increased in increments of multiples 10 (i.e., = 21.9999, 21.999, 21.99, 21.9, and 21.0) and the conductivity is kept constant at = 1.0 S/m. Figure 5 shows plots of (a) the condition number of the matrix using the formulation in (15), (b) the condition number of the matrix using the formulation where (15) has been multiplied by , and (c) the relative error for the first formulation, all as a function of the difference in permittivity between the bath and background (log scale). The condition number is a useful metric because it provides a good measure of the digital accuracy of the inversion process utilizing a particular matrix [47]. Note that the condition number of the case where the matrix is multiplied by does not change appreciably since the values on the diagonal remain within a relatively tight bound. The field distributions (both amplitude and phase) for the cases of background permittivities of 21, 21.9, 21.99, 21.999, and 21.999 are shown in Figures 6(a)–6(e). For the most part, these distributions are quite similar. Except for , the amplitude plots all show a primary lobe towards 12:00 (clock-face orientation) with slight nulls to either side of it. The case does exhibit a slight bulge towards 12:00, but the nulls have been smoothed over. The phase distributions show roughly circular patterns with the sharp changes from red to blue designating phase wrapping as the distributions jump from to +. These distributions are roughly circular, with the cases for exhibiting a slight flattening in the 12:00 direction. The distribution is more circular for the entire region. From the visual examination in Figures 6(a)–6(e) and the error and condition number plots in Figure 5, it would appear that background permittivity values of or greater would be suitable. However, our preference is to keep the condition number to a modest level, so we choose a value of to restrain the condition number while simultaneously keeping the error low.

(a) 21.0

(b) 21.9

(c) 21.99

(d) 21.999

(e) 21.9999

3.2. Computational Efficiency with the FFT

We have implemented the 2D-DDA using both MATLAB and C++. The iterative and direct solutions for the given system of equations have also been calculated. Since the main concern on the choice of these algorithms stems from optimizing computation time to obtain vector , we have calculated and compared the computation times for each implementation and method with different numbers of dipoles. The computation times are divided into three categories based on the complexities , , and for the direct solvers, the COCG iterative algorithm and the FFT-COCG algorithm, respectively. The computation time for the COCG-FFT implemented in C++ is processed in two separate software packages, Armadillo and FFTW. Armadillo is an open source C++ library for linear algebra and scientific computing which provides a high level syntax and is balanced between speed and user friendliness [59, 60]. FFTW is a C subroutine library for computing the discrete Fourier transform (DFT) in one or more dimensions, of arbitrary input size, and of both real and complex data [37]. These two packages have been used in our implementation in this study. Table 1 shows the computation times for different solvers in MATLAB and C++. For these simulations, the forward model zone contains number of dipoles, where , and 6561. The computation times for these cases are given in Table 1. Table 1 shows that the MATLAB direct solver is more efficient compared to that of the C++ implementation using the Armadillo package. Generally, for the two smaller matrices, MATLAB is more efficient than C++. However, when using larger matrices with sufficient numbers of dipoles in the imaging domain () to ensure accuracy, the computation times for the direct and iterative solvers in MATLAB show that the time savings are on the order of a factor of two. Using the FFT-COCG in MATLAB, the time decreases by a factor of three. Our observations show that there is a significant difference in computation time between the Armadillo and FFTW packages using the built-in FFT algorithm.

3.3. two orders of magnitude faster than a conFurther Optimizations for Our Specific Microwave Imaging Problem

In this section, the optimal number of dipoles utilized in the imaging domain is investigated. In the previous section, it was convenient to sequentially double the number of dipoles to more easily identify the optimal number of dipoles. However, there may be more convenient intermediate values which have improved efficiency, especially with respect to the FFT.

For the block-Toeplitz matrices in the case with 6561 dipoles, the square Toeplitz matrices are 81 × 81. To convert these matrices to circulant ones, the size of the matrix needs to be nearly doubled in width and height along with duplication of coefficients to that in (18). In this case, the conversion means adding 79 columns with appropriate coefficients, before performing the actual FFT. In this form, the number of columns and the length of each column increases to 160.

The fast Fourier transformations of the first row of the matrix and the column vector also require zero-padding procedures to the next power of 2 to allow use of the FFT. In this way, a vector with the size of 160 would be converted to the one of size 256. One way to optimize these cases is to restrict our number of dipoles in a way that does not require significant zero-padding. One example of exploiting this observation is to reduce the overall grid size to dipoles. Without reducing the physical size of the grid, this would imply a node-to-node spacing of 6.3 mm instead of the previous 5.0 mm spacing, an increase of 25%. However, by decreasing the grid size by only (corresponding to spacings of 5.5mm and an overall grid size of 35.2cm x 35.2cm), the overall change to the accuracy is only minor, and we accomplish the goal of improved efficiency. Figure 7 shows the magnitude and phase distributions for the same case in Figure 6(c), the background permittivity of , except that the overall grid size is 35.2 (cm) × 35.2 (cm). While the physical boundary is closer to the antennas, the artifacts from the boundary still have minimal impact on the fields inside the array of antennas. In terms of efficiency, Table 2 shows the performance times for this grid size. For all solution techniques, this size has substantial time improvement over the previous 6561 case, especially when the FFT is exploited, by as mush as a factor of depending on method. This is a substantial improvement.

Moreover, for the iterative algorithms such as COCG, the number of iterations is not independent of prescribed accuracy tolerance levels. To reduce the computation time, we studied the effects of the error tolerance level on the computation time as well as the overall accuracy. Table 3 shows the corresponding computation times for tolerances of 1e-5 and 1e-3 for the 6561 dipole grid. For all implementations, there are significant time improvements, with the C++ times utilizing the FFT algorithm improving on the order of 35%.

4. Conclusion

We have implemented a version of the DDA for computing the forward solutions for a configuration used in our current 2D tomographic imaging system. 2D is intriguing in that it has already been successfully implemented for phantom and animal experiments along with considerable clinical studies. Even without the speed enhancements from the techniques presented here, the 2D approach is already considerably faster than any existing 3D inverse technique which positions it well with regard to being the foundation for a low cost and portable system that would be suitable for under-resourced settings.

We have previously demonstrated that the DDA is accurate for forward solution computation for our imaging algorithm. This study builds on that experience and formulates the problem in ways that facilitate dramatic speed optimization. Chief among these enhancements is the notion of utilizing conjugate gradient based iterative solvers in conjunction with breaking the core matrix into a simple diagonal matrix and a separate one that can be further broken into block-Toeplitz matrices. Through standard techniques, these smaller matrices can be easily expanded into circulant matrices for which the FFT can be employed to speed up the matrix/vector multiplications utilizing the convolution theorem.

Further enhancements including judicious selection of the grid size and analysis of appropriate error tolerance levels for the conjugate gradient type iterative process allow the forward solution time to decrease to the order of 0.1 second which is almost two orders of magnitude faster than a conventional COMSOL approach that is considered to be efficient. These computation time improvements along with significant reductions in memory usage because of the nature of circulant matrices make this an attractive approach for simplifying and greatly speeding up the image reconstruction process.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The article has been included in a thesis by the first author Samar Hosseinzadegan for the Licentiate Degree of Engineering at the Chalmers University of Technology, Göteborg, Sweden [61]. This work was sponsored by an NIH/NCI grant # R01-CA191227 and a Chalmers Foundation Excellence grant.

References

D. J. Brenner and E. J. Hall, “Computed tomography—an increasing source of radiation exposure,” The New England Journal of Medicine, vol. 357, no. 22, pp. 2277–2284, 2007.
View at: Publisher Site | Google Scholar
T. Jansson, J. E. Westlin, H. Ahlström, A. Lilja, B. Långström, and J. Bergh, “Positron emission tomography studies in patients with locally advanced and/or metastatic breast cancer: a method for early therapy evaluation?” Journal of Clinical Oncology, vol. 13, no. 6, pp. 1470–1477, 1995.
View at: Publisher Site | Google Scholar
D. M. Ikeda, D. R. Baker, and B. L. Daniel, “Magnetic resonance imaging of breast cancer: Clinical indications and breast MRI reporting system,” Journal of Magnetic Resonance Imaging, vol. 12, no. 6, pp. 975–983, 2000.
View at: Publisher Site | Google Scholar
S. G. Orel and M. D. Schnall, “MR imaging of the breast for the detection, diagnosis, and staging of breast cancer,” Radiology, vol. 220, no. 1, pp. 13–30, 2001.
View at: Publisher Site | Google Scholar
E. Porter, M. Coates, and M. Popović, “An early clinical study of time-domain microwave radar for breast health monitoring,” IEEE Transactions on Biomedical Engineering, vol. 63, no. 3, pp. 530–539, 2016.
View at: Publisher Site | Google Scholar
E. C. Fear, J. Bourqui, C. Curtis, D. Mew, B. Docktor, and C. Romano, “Microwave breast imaging with a monostatic radar-based system: a study of application to patients,” IEEE Transactions on Microwave Theory and Techniques, vol. 61, no. 5, pp. 2119–2128, 2013.
View at: Publisher Site | Google Scholar
A. W. Preece, I. Craddock, M. Shere, L. Jones, and H. L. Winton, “MARIA M4: Clinical evaluation of a prototype ultrawideband radar scanner for breast cancer detection,” Journal of Medical Imaging, vol. 3, no. 3, Article ID 033502, 2016.
View at: Google Scholar
H. Song, S. Sasada, T. Kadoya et al., “Detectability of breast tumor by a hand-held impulse-radar detector: performance evaluation and pilot clinical study,” Scientific Reports, vol. 7, Article ID 16353, 2017.
View at: Google Scholar
R. K. Amineh, M. Ravan, A. Khalatpour, and N. K. Nikolova, “Three-dimensional near-field microwave holography using reflected and transmitted signals,” IEEE Transactions on Antennas and Propagation, vol. 59, no. 12, pp. 4777–4789, 2011.
View at: Publisher Site | Google Scholar
D. Tajik, F. Foroutan, D. S. Shumakov, A. D. Pitcher, and N. K. Nikolova, “Real-time microwave imaging of a compressed breast phantom with planar scanning,” IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology, vol. 2, pp. 154–162, 2018.
View at: Google Scholar
R. A. Kruger, D. R. Reinecke, and G. A. Kruger, “Thermoacoustic computed tomography–technical considerations,” Medical Physics, vol. 26, no. 9, pp. 1832–1837, 1999.
View at: Publisher Site | Google Scholar
G. Ku and L. V. Wang, “Scanning microwave-induced thermoacoustic tomography: signal, resolution, and contrast,” Medical physics, vol. 28, pp. 4–10, 2001.
View at: Google Scholar
S. P. Poplack, T. D. Tosteson, W. A. Wells et al., “Electromagnetic breast imaging: results of a pilot study in women with abnormal mammograms,” Radiology, vol. 243, no. 2, pp. 350–359, 2007.
View at: Publisher Site | Google Scholar
P. M. Meaney, P. A. Kaufman, L. S. Muffly et al., “Microwave imaging for neoadjuvant chemotherapy monitoring: initial clinical experience,” Breast Cancer Research, vol. 15, article R35, 2013.
View at: Google Scholar
S. Y. Semenov, A. E. Bulyshev, A. Abubakar et al., “Microwave-tomographic imaging of the high dielectric-contrast objects using different image-reconstruction approaches,” IEEE Transactions on Microwave Theory and Techniques, vol. 53, no. 7, pp. 2284–2294, 2005.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar
P. Tournier, M. Bonazzoli, V. Dolean et al., “Numerical modeling and high-speed parallel computing: new perspectives on tomographic microwave imaging for brain stroke detection and monitoring,” IEEE Antennas and Propagation Magazine, vol. 59, no. 5, pp. 98–110, 2017.
View at: Publisher Site | Google Scholar
O. M. Bucci, L. Crocco, and R. Scapaticci, “On the optimal measurement configuration for magnetic nanoparticles-enhanced breast cancer microwave imaging,” IEEE Transactions on Biomedical Engineering, vol. 62, no. 2, pp. 407–414, 2015.
View at: Publisher Site | Google Scholar
Z. Miao and P. Kosmas, “Multiple-frequency DBIM-TwIST algorithm for microwave breast imaging,” IEEE Transactions on Antennas and Propagation, vol. 65, no. 5, pp. 2507–2516, 2017.
View at: Publisher Site | Google Scholar
I. Catapano, L. Crocco, M. D' Urso, and T. Isernia, “3D microwave imaging via preliminary support reconstruction: testing on the Fresnel 2008 database,” Inverse Problems, vol. 25, no. 2, Article ID 024002, 23 pages, 2009.
View at: Publisher Site | Google Scholar | MathSciNet
A. Fasoula, L. Duchesne, J. Gil Cano, P. Lawrence, G. Robin, and J. Bernard, “On-site validation of a microwave breast imaging system, before first patient study,” Diagnostics, vol. 8, no. 3, article 53, 2018.
View at: Publisher Site | Google Scholar
R. Scapaticci, P. Kosmas, and L. Crocco, “Wavelet-based regularization for robust microwave imaging in medical applications,” IEEE Transactions on Biomedical Engineering, vol. 62, no. 4, pp. 1195–1202, 2015.
View at: Publisher Site | Google Scholar
T. Rydholm, A. Fhager, M. Persson, and P. M. Meaney, “A first evaluation of the realistic supelec-breast phantom,” IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology, vol. 1, pp. 59–65, 2017.
View at: Google Scholar
T. Rydholm, A. Fhager, M. Persson, S. Geimer, and P. Meaney, “Effects of the plastic of the realistic GeePS-L2S-breast phantom,” Diagnostics, vol. 8, article 61, 2018.
View at: Publisher Site | Google Scholar
P. M. Meaney, M. W. Fanning, T. Raynolds et al., “Initial clinical experience with microwave breast imaging in women with normal mammography,” Academic Radiology, vol. 14, no. 2, pp. 207–218, 2007.
View at: Publisher Site | Google Scholar
P. M. Meaney, M. W. Fanning, D. Li, S. P. Poplack, and K. D. Paulsen, “A clinical prototype for active microwave imaging of the breast,” IEEE Transactions on Microwave Theory and Techniques, vol. 48, no. 1, pp. 1841–1853, 2000.
View at: Publisher Site | Google Scholar
S. Semenov, R. Svenson, A. Boulyshev et al., “Microwave tomography: two-dimensional system for biological imaging,” IEEE Transactions on Biomedical Engineering, vol. 43, no. 9, pp. 869–877, 1996.
View at: Publisher Site | Google Scholar
P. M. Meaney, F. Shubitidze, M. W. Fanning et al., “Surface wave multipath signals in near-field microwave imaging,” Journal of Biomedical Imaging, vol. 2012, Article ID 697253, 8 pages, 2012.
View at: Google Scholar
Q. Fang, P. M. Meaney, and K. D. Paulsen, “Singular value analysis of the Jacobian matrix in microwave image reconstruction,” IEEE Transactions on Antennas and Propagation, vol. 54, no. 8, pp. 2371–2380, 2006.
View at: Publisher Site | Google Scholar
P. M. Meaney, Q. Fang, T. Rubaek, E. Demidenko, and K. D. Paulsen, “Log transformation benefits parameter estimation in microwave tomographic imaging,” Medical Physics, vol. 34, no. 6, pp. 2014–2023, 2007.
View at: Publisher Site | Google Scholar
P. M. Meaney, S. D. Geimer, and K. D. Paulsen, “Two-step inversion in microwave imaging with a logarithmic transformation,” Medical Physics, vol. 44, no. 8, pp. 4239–4251, 2017.
View at: Publisher Site | Google Scholar
A. Fhager, M. Gustafsson, and S. Nordebo, “Image reconstruction in microwave tomography using a dielectric debye model,” IEEE Transactions on Biomedical Engineering, vol. 59, no. 1, pp. 156–166, 2012.
View at: Publisher Site | Google Scholar
Q. Fang, P. M. Meaney, and K. D. Paulsen, “Viable three-dimensional medical microwave tomography: theory and numerical experiments,” IEEE Transactions on Antennas and Propagation, vol. 58, no. 2, pp. 449–458, 2010.
View at: Publisher Site | Google Scholar | MathSciNet
P. Rocca, M. Benedetti, M. Donelli, D. Franceschini, and A. Massa, “Evolutionary optimization as applied to inverse scattering problems,” Inverse Problems, vol. 25, no. 12, Article ID 123003, 2009.
View at: Google Scholar
T. M. Grzegorczyk, P. M. Meaney, P. A. Kaufman, R. M. Diflorio-Alexander, and K. D. Paulsen, “Fast 3-D tomographic microwave imaging for breast cancer detection,” IEEE Transactions on Medical Imaging, vol. 31, no. 8, pp. 1584–1592, 2012.
View at: Publisher Site | Google Scholar
S. Hosseinzadegan, A. Fhager, M. Persson, and P. Meaney, “Application of two-dimensional discrete dipole approximation in simulating electric field of a microwave breast imaging system,” IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology, 2018.
View at: Google Scholar
M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” in Proceedings of the IEEE Special Issue on Program Generation, Optimization, and Platform Adaptation, vol. 93, pp. 216–231, 2005.
View at: Google Scholar
E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” The Astrophysical Journal, vol. 186, pp. 705–714, 1973.
View at: Publisher Site | Google Scholar
B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” Journal of the Optical Society of America A, vol. 11, pp. 1491–1499, 1994.
View at: Google Scholar
K. Skorupski, “Using the DDA (discrete dipole approximation) method in determining the extinction cross section of black carbon,” Metrology and Measurement Systems, vol. 22, no. 1, pp. 153–164, 2015.
View at: Publisher Site | Google Scholar
A. Lakhtakia, “Strong and weak forms of the method of moments and the coupled dipole method for scattering of time-harmonic electromagnetic fields,” International Journal of Modern Physics C, vol. 3, pp. 583–603, 1992.
View at: Publisher Site | Google Scholar
W. C. Chew, Waves and Fields in Inhomogeneous Media, IEEE press, 1995.
B. T. Draine and J. Goodman, “Beyond clausius-mossotti-wave propagation on a polarizable point lattice and the discrete dipole approximation,” The Astrophysical Journal, vol. 405, no. 2, pp. 685–697, 1993.
View at: Publisher Site | Google Scholar
F. Kahnert, “Numerical methods in electromagnetic scattering theory,” Journal of Quantitative Spectroscopy & Radiative Transfer, vol. 79-80, pp. 775–824, 2003.
View at: Publisher Site | Google Scholar
M. A. Yurkin, V. P. Maltsev, and A. G. Hoekstra, “Convergence of the discrete dipole approximation. I. theoretical analysis,” The Journal of the Optical Society of America A (JOSA A), vol. 23, pp. 2578–2591, 2006.
View at: Publisher Site | Google Scholar
A. Sihvola, “Peculiarities in the dielectric response of negative-permittivity scatterers,” Progress in Electromagnetics Research, vol. 66, pp. 191–198, 2006.
View at: Publisher Site | Google Scholar
G. H. Golub and C. F. Van Loan, Matrix Computations, vol. 3, JHU Press, 2012.
P. J. Flatau, “Improvements in the discrete-dipole approximation method of computing scattering and absorption,” Optics Expresss, vol. 22, no. 16, p. 1205, 1997.
View at: Publisher Site | Google Scholar
J. J. Goodman, P. J. Flatau, and B. T. Draine, “Application of fast-Fourier-transform techniques to the discrete-dipole approximation,” Optics Expresss, vol. 16, no. 15, pp. 1198–1200, 1991.
View at: Publisher Site | Google Scholar
P. Sonneveld, “CGS, a fast lanczos-type solver for nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 10, no. 1, pp. 36–52, 1989.
View at: Publisher Site | Google Scholar
H. A. van der Vorst, “BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 13, no. 2, pp. 631–644, 1992.
View at: Publisher Site | Google Scholar | MathSciNet
M. H. Gutknecht, “Variants of BICGSTAB for matrices with complex spectrum,” SIAM Journal on Scientific Computing, vol. 14, no. 5, pp. 1020–1033, 1993.
View at: Publisher Site | Google Scholar
Y. Saad and M. Schultz, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 7, no. 3, pp. 856–869, 1986.
View at: Publisher Site | Google Scholar
R. W. Freund, “A transpose-free quasi-minimal residual algorithm for non-hermitian linear systems,” SIAM Journal on Scientific Computing, vol. 14, no. 2, pp. 470–482, 1993.
View at: Publisher Site | Google Scholar
R. W. Freund, “Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices,” SIAM Journal on Scientific and Statistical Computing, vol. 13, no. 1, pp. 425–448, 1992.
View at: Publisher Site | Google Scholar
M. Clemens and T. Weiland, “Iterative methods for the solution of very large complex symmetric linear systems of equations in electrodynamics,” Front Range Scientific Computations, Inc, Lakewood, NJ, USA, 1996.
View at: Google Scholar
H. A. van der Vorst and J. B. M. Melissen, “Petrov-Galerkin type method for solving Ax = b, where A is symmetric complex,” IEEE Transactions on Magnetics, vol. 26, no. 2, pp. 706–708, 1990.
View at: Publisher Site | Google Scholar
T. Sogabe and S.-L. Zhang, “A COCR method for solving complex symmetric linear systems,” Journal of Computational and Applied Mathematics, vol. 199, no. 2, pp. 297–303, 2007.
View at: Publisher Site | Google Scholar
C. Sanderson and R. Curtin, “Armadillo: a template-based C++ library for linear algebra,” Journal of Open Source Software, 2016.
View at: Google Scholar
C. Sanderson and R Curtin, “A user-friendly hybrid sparse matrix class in C++,” in Proceedings of the International Congress on Mathematical Software, pp. 422–430, 2018.
View at: Google Scholar
S. Hosseinzadegan, A discrete dipole approximation forward solver for microwave breast imaging [Dissertation, thesis], Chalmers University of Technology, 2019.

Copyright

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.

PDF Download Citation

Download other formats

Order printed copies

Views

1363

Downloads

1402

Citations