Distorted born iterative method
Proposal of measurement configuration
Consider the setting of transmitters and receivers in Fig. 1.
As shown in Fig. 1, the measurement configuration is arranged circularly around the target. The pressure signal from a transmitter is propagated, scattered and measured by receivers. In principle, transmitters and receivers can be distributed uniformly, randomly, or deterministically based on the measurement configuration. Besides, the number of transmitters and receivers depends on different scenarios, in line with practical requirements. It should be noticed that the number of transmitters and receivers should not be so large because the measurement system will be complex. In addition, a large amount of computation and memory requirement are needed to store and process the data. It should be also noticed that the locations of transmitters and receivers may be either different or identical. However, with current ultrasound transducer technology, one transducer can both transmit and receive ultrasound signals. In other words, the location of transmitter and receiver can be identical.
For the proposed measurement configuration, we assume that there are M transmitters and N receivers. M transmitters are arranged around the target to obtain complete information about the target at different angles. Transmitter-receiver procedure of ultrasound wave is performed as follows. Initially, only the first transmitter transmits ultrasound signal and all receivers (N) will simultaneously collect the scattered ultrasonic signal. As a result, we get the 1st set of measurements (i.e. N measurements) at the first position of the transmitter. The second transmitter then begins operating and all receivers simultaneously obtain scattering signal at the second position of the transmitter. At the end of this procedure, we obtain the 2nd set of measurements (i.e. 2 × N measurements). The process continues until the last transmitter (i.e. Mth) where we will obtain M sets of measurements (i.e. M × N measurements). With these measured values, full information about the target at different angles around the target are received.
The measured data will then be brought to DBIM to estimate the sound contrast. The change of the sound speed would be utilized to detect any tissue if exists. In this paper, the transmitters and receivers are working under the condition of homogeneous background medium, i.e. water, where there is a target (strange tumor). This target has sound contrast different from the background (denoted as ∆c) which is embedded in this medium. Because DBIM is used to detect strange tumors in early stage in this paper, sound constrast compared to the background medium is usually very small. Therefore, in our scenarios, which will be discussed in section Results below, value of ∆c will be selected to a few percents.
DBIM method
Because of the well-known properties, Bessel function [11] is usually used in the numerical simulation as a transmitted signal called Incident Wave whose frequency is f. Wavelength λ of this wave is calculated by
where c
0 is the sound speed in background medium. For homogeneous medium, the received signal at the receiver is the Incident Wave itself. In presence of tumors, the medium becomes inhomogeneous; when the incident wave hits the target, the following two situations may occur: if the target size is much larger than the wavelength of the incident wave, it is reflected; if the target size is less than or equal to the wavelength of the incident wave, it is scattered in all directions around the target. Notice that ultrasound frequency for clinical diagnosis is in the range from 20 KHz to 12 MHz. As a result, the wavelength ranges from 6.2 μm to 74.2 mm if sound speed in the background is 1484 m/s. Born iterative method is used to find out the linear relationship of the scattered pressure difference and the sound contrast difference. The key of this method is that the scattering signal is considered to be very small in comparison with the incident signal. This is entirely consistent with the real requirements which need to detect strange tumors in early stage. Therefore, we consider to reconstruct targets with very small sound contrast, i.e. scattering signal is very small.
The wave equation of the scheme discussed above can be expressed by
$$ \mathrm{p}\left(\overrightarrow{r}\right)={\mathrm{p}}^{\mathrm{inc}}\left(\overrightarrow{r}\right)+{\mathrm{p}}^{\mathrm{sc}}\left(\overrightarrow{r}\right), $$
(1)
where \( {\mathrm{p}}^{\mathrm{sc}}\left(\overrightarrow{r}\right) \), \( {\mathrm{p}}^{\mathrm{inc}}\left(\overrightarrow{r}\right) \), and \( \mathrm{p}\left(\overrightarrow{r}\right) \) are the scattered, incident, and total signals respectively. It can be seen that the known data arethe total signal and the incident signal. However, what we concern about is the reconstruction of the unknown target T (r) from the obtained data. This is the inverse problem.
Consider the wave numbers of the background and target mediums to be k
0 and k (r) respectively. According to [12], inhomogeneous differential equation has the form
$$ \left({\mathit{\nabla}}^2+{k}_0^2(r)\right) p(r)=- O(r) p(r) $$
Green function is an effective method to solve inhomogeneous differential equation. Therefore, it is used to find out the nonlinear relationship of the scattered signal and the target based on the total and incident signals. Eq. (1) can be rewritten in details using the Green function G0 (·)
$$ \mathrm{p}\left(\overrightarrow{\mathrm{r}}\right)={\mathrm{p}}^{\mathrm{inc}}\left(\overrightarrow{r}\right)+{\displaystyle \iint}\mathrm{T}\left(\overrightarrow{\mathrm{r}}\right)\mathrm{p}\left(\overrightarrow{{\mathrm{r}}^{\hbox{'}}}\right){\mathrm{G}}_0\left({k}_0,\left|\overrightarrow{r}-\overrightarrow{{\mathrm{r}}^{\hbox{'}}}\right|\right)\mathrm{d}\overrightarrow{{\mathrm{r}}^{\hbox{'}}} $$
(2)
Notice that when the background medium is homogenerous, G0 is the 0-th Hankel function of the first kind and described by
$$ \begin{array}{l}{\mathrm{G}}_0\left({k}_0,\left|\overrightarrow{\mathrm{r}}-\overrightarrow{{\mathrm{r}}^{\hbox{'}}}\right|\right)=\frac{- i}{4}{H}_0^{(1)}\left({k}_0\left|\overrightarrow{\mathrm{r}}-\overrightarrow{{\mathrm{r}}^{\hbox{'}}}\right|\right)=\\ {}\frac{- i}{4}\sqrt{\frac{2}{\pi {k}_0\left|\overrightarrow{\mathrm{r}}-\overrightarrow{{\mathrm{r}}^{\hbox{'}}}\right|}}{e}^{i\left({k}_0\left|\overrightarrow{\mathrm{r}}-\overrightarrow{{\mathrm{r}}^{\hbox{'}}}\right|-\pi /4\right)}.\end{array} $$
(3)
T (r) in Eq. (2) is the target function that needs to be estimated. It can be calculated as follows:
$$ \mathrm{T}\left(\mathrm{r}\right)=\left\{\begin{array}{l} k{(r)}^2\hbox{--} {k}_0^2={\omega}^2\left(\frac{1}{c^2}-\frac{1}{c_0^2}\right)\kern0.5em if\ r\le R\\ {}0\kern6.36em if\ r> R\ \end{array}\right. $$
(4)
Equation (4) indicates that the ideal target function depends on the frequency of the incident signal (ω = 2πf) and the sound speed difference of the background medium (c
0
) as well as the target medium (c). In order to calculate in details each pixel inside the region of interest, scattering wave equation (2) needs to be discreterized. The method of moment (MoM) which uses basic sinc function is used to solve this problem.
The total pressure field in the observed mesh area (N × N points) can be expressed by
$$ \overline{p} = \left(\overline{I}-\overline{C}. D\left(\overline{T}\right)\right){p}^{inc}, $$
(5)
where \( \overline{C} \) is the Green matrix showing the interactions among pixels, Ī is unit matrix, and D (·) returns a square diagonal matrix of the input vector. The scattered signal in the form of N
t
N
r
× 1 vector is described by
$$ {\overline{p}}^{sc}=\overline{B}. D\left(\overline{T}\right).\overline{p}, $$
(6)
where \( \overline{B} \) is the Green matrix showing the interaction of all pixels to the receiver. We have to determine two parameters \( \overline{p} \) and \( \overline{T} \) in Eqs. (5) and (6).
By rewritting these equations, we have [13]
$$ \varDelta {p}^{sc}=\overline{B}. D\left(\overline{p}\right).\varDelta \overline{T}=\overline{M}.\varDelta \overline{T}, $$
(7)
where \( \overline{M}=\overline{B}\cdot D\left(\overline{p}\right) \). For a transmitter and a receiver, we formulate a matrix \( \overline{M} \) and a scalar value Δp
sc. The target function \( \overline{T} \) has N
2 variables corresponding to the number of pixels in the region of interest. It can be estimated by:
$$ {\overline{T}}^n={\overline{T}}^{\left( n-1\right)}+\varDelta {\overline{T}}^{\left( n-1\right)}, $$
(8)
where n and n-1 are two consecutive discrete-time points. \( \varDelta \overline{T} \) is estimated by using Tikhonov’s regularization [14]:
$$ \varDelta \overline{T}= \arg \underset{\Delta \overline{T}}{ \min }{{\left\Vert \Delta \overline{p}\right.}^{sc}}_t-\overline{M_t}\Delta \overline{T}\left\Vert {}_2^2\right.+\gamma {\left\Vert \Delta \overline{T}\right\Vert}_2^2, $$
(9)
where \( \varDelta {\overline{p}}^{sc} \) is the difference between estimated and measured scattered signals whose size is (N
t
N
r
× 1). Besides, measurement results are assembled in a matrix form \( {\overline{M}}_t \) of (N
t
N
r
× N
2) elements and γ is the regularization factor that needs to be carefully chosen.
The DBIM procedure is presented in Algorithm 1.
Deterministic compressive sampling
Fundamentals of compressive sampling
Compressive Sampling (CS) or known as Compressed Sensing [15] allows for exactly recovering signal v ∈ ℝn from a small number of “random measurements” u ∈ ℝm (m < < n), which are defined as follows:
where Φ is a m × n matrix called sampling basis. The columns of Φ have entries equal to 1 at random positions; the entries at other positions are null. This structure means that measurements are done randomly.
The core problem of compressive sampling is that assuming v has sparse representation in an orthonormal basis Ψ, i.e.
in which, w, also known as sparse, only has s < m ≪ n non-zero coefficients. Note that s is the sparsity degree of the signal u and the number m of required measurements is normally higher than s. Compressive sampling theory shows that this sparse property allows accurate recovery w with overwhelming probability success [16]. In particular, sensing basis must have incoherent property to the model basis Ψ [17]. This property is guaranteed by the randomness of the non-zero components in Φ. Therefore, the problem can be written as follows:
$$ u = \Phi \Psi w= Aw, $$
(12)
where A ∈ R
mxn is full-rank, i.e. the m rows of A are linearly independent. By these settings, Eq. (12) is solved for w with w-sparse constraint. Once w is obtained, v can be calculated from (11). With A satisfying “Restricted Isometry Property (RIP)” condition which is first established by Candes et al. [16], the CS problem can be solved through l
0-minimization technique:
$$ \widehat{w}= \arg \underset{w\in {R}^n}{ \min}\left\Vert w\left\Vert {}_{l_0}\right.\right.\mathrm{subject}\ \mathrm{t}\mathrm{o}\ u= Aw, $$
(13)
where l
0 norm is defined as \( {\left\Vert w\right\Vert}_{l_0}:={ \max}_i\left|{w}_i\right| \). In this paper, instead of solving (13), the sparsest solution \( \widehat{v} \) of (12) can be found by solving the l
1-minimization problem [18]:
$$ \widehat{w}= \arg \underset{w\in {R}^n}{ \min}\left\Vert w\right.\left\Vert {}_{l_1}\right.\;\mathrm{subject}\ \mathrm{t}\mathrm{o}\ u= Aw, $$
(14)
in which, the l
1 norm is defined as \( \left\Vert w\right\Vert {l}_1:={\displaystyle {\sum}_{i=1}^n\left|{w}_i\right|.} \)
Equation (12) is introduced under the assumption that the exact form of the reconstructed signal is given. This is rarely the case in practice, because the measurements are often affected by noise. To reconstruct the signal in case of noisy measurements, we have:
where e represents the noise \( {\left\Vert e\right\Vert}_{l_2}\le \varepsilon \). l
1 problem in presence of noise can be expressed as follows [18]:
$$ \widehat{w}= \arg { \min}_{w\in {R}^n}{\left\Vert w\right\Vert}_{l_1}\mathrm{subject}\ \mathrm{t}\mathrm{o}\; u={\left\Vert u- Aw\right.}_{l_2}\le \varepsilon $$
(16)
Compressive sampling using deterministic filters
In this paper, we use a deterministic basis generated by a pseudo random sequence instead of pure random basis. The advantage of this method in comparison with CS is that it offers a simpler hardware design, because of the deterministic nature of the proposed measurement system. In fact, to facilitate the practical implementation, researchers focused on the measurement system design whose properties are not completely random or completely deterministic. The first work exploited this view is introduced in [11], followed by others [19, 20]. In the work [21], the authors offer some simple criteria to design a deterministic compressed sampling system that ensures the ability to successfully reconstruct sparse signals. Overall, the deterministic compressed sampling technique inherits some advantages over random one, such as more effective recovery time, clear structure, efficient storage and tighter recovery limit [20].
Because of simpler hardware implementation of a deterministic system compared to a random one’s, in this work, the matrix Φ is constructed using the deterministic method. To perform this, we consider a dynamic deterministic system with deterministic nonlinear characteristic; as its dynamic is very sensitive to initial conditions, its output has random-like property [22]. Note that the reconstruction accuracy of deterministic compressed sampling technique is ensured like the random one [23]. To construct Φ, we consider a Logistic map based dynamic structure that generates a very simple deterministic sequence which is transformed into a sequence that would have a Gaussian-like behavior [22, 23]. The following Logistic map sequence is used as follows
$$ q\left( n+1\right) = \rho q(n)\left(1- q(n)\right), $$
(17)
We then convert it by the Logit Transform to become Gaussian-like as
$$ {q}_G(n)= \ln \left[\frac{q(n)}{1- q(n)}\right] $$
(18)
Remarks: When the control parameter ρ in (17) equals to 4, the dynamics described by (17) is deterministic; this means that the dynamic of (17) is very sensitive in regard of the initial condition q (0). A very small change in q (0) would largely change q (n) in a relatively short time. More detailed construction of Φ from q
G
(n) can be found in [22]. Finally, the image reconstruction is performed using sparse approximation algorithms.
The deterministic compressive sampling DBIM
Measurement configuration
Assume that there are N
t
× N
r
transducers where N
t
is the number of transmitters and N
r
is the number of receivers (that gives N
t
× N
r
measurements). We need to reconstruct the target which is divided into N pixels vertically and horizontally, i.e. N
2 variables. There are three cases for the arrangement of transducers and transmitters that are required for consideration.
Firstly, if N
t
× N
r
is larger than N
2, the number of equations is greater than the number of variables. This is the over-determined problem and it is called a dense scattering domain. In this configuration, the amount of scattered information obtained around the target is large, from which we can calculate and easily reconstruct the target image. However, a large number of transmitters and receivers is required. This high number of devices makes the imaging system much more complex in terms of setting up the measurement equipment and of large required processing time consuming in target image reconstruction. Therefore, this configuration is less likely used in practice.
Secondly, if N
t
× N
r
is approximately or equals to N
2, the number of equations equals to the number of variables. This case is called a moderate scattering domain. For the high-resolution images in this paper, i.e. large N
2, we still need in this situation a relatively large amount of transducers to give enough measurement for correctly reconstructing images.
Thirdly, if N
t
× N
r
is smaller than N
2, i.e. the number of equations is smaller than the number of variables, this is the under-determined problem. When N
t
× N
r
is significantly smaller than N
2, we call this a sparse scattering domain, that is, transducers is sparse on the measurement system due to the small number of transducers and the amount of scattered information is less and sparser than the case of the dense and moderate scattering domains. Define the compressed ratio r as the ratio of N
t
× N
r
and N
2, i.e.
$$ r = \frac{N_t{N}_r}{N^2}. $$
(19)
It is clear that for the sparse scattering problem, r < 1. When r is small, the measurement system is simple and the collected data is of small size. Consequently, the system set-up and calculation are less complex. Henceforth, we need to focus on the case of the smallest value of r, with which, the imaging system is the simplest and it ensures the acceptable image reconstruction.
The implementation processes of the conventional method and the proposed method are shown in Fig. 2 and Fig. 3 respectively. In both figures, the Input represents the ideal target function and the Output represents the reconstructed target function. Nevertheless, in Fig. 2, the ideal target function of interest is reconstructed using Tikhonov regularizations while it is reconstructed using sparse-based approximation techniques such as l
1 optimization or greedy Algorithms in Fig. 3. Besides, the measurement configuration of linear transmitter-and-detector locations is used in conventional method and the measurement configuration of transmitters and detectors with linear transmitter locations and deterministic detector locations is used in the proposed method.
The signal accurate reconstruction of the CS is guaranteed if the sparse domain is not associated with the sampling domain [6, 7]. In [24], the authors compared the efficiency of sparse signal recovery methods, namely l
1-LSP, MOSEK, PDCO-CHOL, PDCO-LSQR, l
1-MAGIC and HOMOTOPY, and MOSEK. Their results show that l
1-LSP method outperforms others by its fast reconstruction and low computational complexity. Moreover, this also works effectively in large dense problems. With these observations, the DBIM method can be employed for not only detecting small scale targets (like breast cancer) but also for other larger-scale applications. Therefore, in this paper, we use l
1-regularized least-squares programs (l
1-LSP) [24] for reconstructing the image target.
l1-regularized least-squares and DCS-DBIM procedure
The l
1-LSP solves l
1-regularized least squares problems using the truncated Newton interior-point method and has the form as follows:
$$ min\left\Vert Aw- u\left\Vert {}^2\right.\right.+\zeta {\left\Vert w\right\Vert}_1, $$
(20)
where A is a data matrix in dense or sparse format with n columns and m rows; w is a vector of length n; and u is a vector of length m. In the DCS-DBIM approach, A, w, u are \( {\overline{M}}_t \), \( \varDelta \overline{T} \), and \( \varDelta {\overline{p}}_t^{sc} \) respectively. It is noted that the size of \( {\overline{M}}_t \) is dependent on the number of transmitters (N
t
), receivers (N
r
) and variables (N
2). In other words, it is changed in different scenarios. In the DCS-DBIM form, l
1 problem is expressed as
$$ \varDelta \overline{T}= \arg \underset{\Delta \overline{T}}{ \min }{{\left\Vert \Delta \overline{p}\right.}^{sc}}_t-\overline{M_t}\Delta \overline{T}\left\Vert {}_2^2\right.+\zeta {\left\Vert \Delta \overline{T}\right\Vert}_1, $$
(21)
where ζ is the regularization parameter in l
1 problem. The selection of regularization parameter ζ is crucial because it has a great influence on the image reconstruction quality. The large value of ζ will make the reconstructed image rough and the small value of ζ will result in high computational complexity. It is noted that the inverse solver matrix \( {\overline{M}}_t \) is changed at each iteration. In fact, in DBIM, the Green function is updated in each iteration; therefore, matrices B and C are changed, leading to the variation of \( {\overline{M}}_t \) at each iteration.
In this context, the regularization parameter ζ also changes at each iteration. It is chosen as a function of the forward error. For the computation process, the Rayleigh quotient is recursively iterated; the value σ0 of the first singular value of the inverse solver matrix \( \overline{M_t} \) is estimated; ζ is then selected as given in [25]. In the simulation scenario presented in this paper, the value of ζ for the first iteration is chosen as 10−15. The DCS-DBIM procedure is presented in Algorithm 2.