The optimization of acoustic fields for ablative therapies of tumours in the upper abdomen

Consider a scatterer defined by a closed smooth surface S. Let the medium surrounding the scatterer be defined as the exterior domain V_ext. Let the boundary condition at infinity satisfy Sommerfeld's acoustic radiation condition. The propagation of time-harmonic acoustic waves in a homogeneous isotropic inviscid medium is described by the Helmholtz equation

$$\begin{equation} \nabla ^2 p( {\vec r} ) + {\bm k}^2 p( {\vec r} ) = 0\quad \quad \forall \vec r \in V_{{\rm ext}} ,\end{equation} \tag{ 1 }$$

where p is the acoustic pressure, k = ω/c is the complex acoustic wave number where ω is the angular frequency, c the complex sound speed in the external medium and $\vec r$ is the position vector. A source term may be included in the right-hand side of (1) (incident pressure wave), but it is convenient to split the total pressure $p( {\vec r} )$ as the sum of the incident pressure $p_i ( {\vec r} )$ (i.e. in absence of the scatterer) and the pressure scattered by the surface S, $\forall \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} \in V_{{\rm ext}}$.

For scattering problems, the integral representation of the solution to (1) is given by the Helmholtz integral equation (Courant and Hilbert 1962, p 258):

$$\begin{equation} \fl\mathop \smallint \limits_S^{} \left[ {p( {\vec r_q } )\frac{{\partial G( {\vec r{\rm |}\vec r_q } )}}{{\partial n_q }} - \frac{{\partial p( {\vec r_q } )}}{{\partial n_q }}G( {\vec r{\rm |}\vec r_q } )} \right]{\rm d}S = \frac{1}{2}p( {\vec r} ) - p_i ( {\vec r} )\quad \quad \forall \vec r \in S,\end{equation} \tag{ 2 }$$

$$\begin{equation} \fl\mathop \smallint \limits_S^{} \left[ {p( {\vec r_q } )\frac{{\partial G( {\vec r{\rm |}\vec r_q } )}}{{\partial n_q }} - \frac{{\partial p( {\vec r_q } )}}{{\partial n_q }}G( {\vec r{\rm |}\vec r_q } )} \right]{\rm d}S = p( {\vec r} ) - p_i ( {\vec r} )\quad \quad \forall \vec r \in V_{{\rm ext}} ,\end{equation} \tag{ 3 }$$

where the Green's function in a three-dimensional space is given by

$$\begin{equation} G( {\vec r{\rm |}\vec r_q } ) = \frac{1}{{4\pi \|\vec r - \vec r_q \|}}{\rm e}^{ - {\rm i}{\bm k}\|\vec r - \vec r_q \|} \end{equation} \tag{ 4 }$$

and $p( {\vec r_q } )$ is the pressure on the surface S at point $\vec r_q$ and $\vec n_q$ is the outward normal.

It is well known that the problem described by the discretization of integral equation (2) suffers from non-uniqueness at frequencies of excitation approaching an eigenvalue of one of the (fictitious) modes of the cavity inside the scatterer. In cases for which the dimensions of the scatterer are large compared with the wavelength in the propagating medium (i.e. in the exterior domain V_ext), such as in transthoracic HIFU applications, this is likely to occur. The matrix formed by discretizing equation (2) then becomes close to singular. When the wave number is purely real, the method which appears to offer the best compromise in terms of application of the Helmholtz integral equation to exterior acoustic problems involving scatterers of arbitrary shape remains the Burton–Miller formulation (Burton and Miller 1971), which solves for a linear combination of equation (2) and its derivative with respect to the outward normal vector on S at $\vec r$. However, for a wave number whose imaginary part is non-zero, solving the discretized form of equation (2) may be sufficient, as the fictitious modes of the cavity inside the scatterer are likely to be damped out (Colton and Kress 1983, p 81 and p 84). For the complex sound speed, frequency of excitation and scatterer dimensions considered in this paper, it will be demonstrated in section 3 that a discretization of the Helmholtz integral equation is indeed sufficient and that the coupling term involving the normal derivative formulation need not be considered. The methodology described below for solving the inverse problem may nevertheless be generalized to the Burton–Miller formulation.

By discretizing equation (2), for all position vectors $\vec r$ on S corresponding to each node on the mesh of the surface, a linear system of equations may be generated. Assuming generic properties of 1500 m s⁻¹ for the speed of sound in the surrounding medium and 1000 kg m⁻³ for its density, the ratio of characteristic acoustic impedances of the medium and of bone is around 0.2, using the speed of sound and density values for ribs as reported by Aubry et al (2008). As part of this work, it was therefore deemed not unreasonable to assume the ribs to be perfect reflectors. For a perfectly rigid scatterer, the normal derivative terms in equation (2) are zero.

$$\begin{equation} [ H ]\{ {p_{{\rm surf}} } \} = - \{ {p_i } \},\end{equation} \tag{ 5 }$$

where [H] is the boundary element matrix and the acoustic pressures on the surface S have been relabelled {p_surf}.

Let the incident field be a linear combination of plane circular pistons rigidly vibrating in an infinite baffle. The incident field on S is a linear combination of the source velocities of each piston. Equation (5) may therefore be re-written as

$$\begin{equation} [ H ]\{ {p_{{\rm surf}} } \} = - [ \beta ]\{ U \},\end{equation} \tag{ 6 }$$

where the elements of [β] may be obtained analytically in the far-field or by solving the Rayleigh integral in the near-field. [β] is of dimension M × N, where M is the number of nodes on the surface S after its discretization and N is the number of plane pistons. {U} is the vector of source velocities and is of dimension N × 1.

Consider equation (3), where the position vectors are located in the exterior volume V_ext. Again, the term on the left-hand side involving the normal derivative of the acoustic pressure on S is zero for a perfectly rigid surface. When evaluated numerically, the integral may be expressed as a weighted sum of the pressures on S. Additionally, the incident pressure at any given field location is a linear combination of the source velocities. Hence, for a specified number of locations in the exterior volume, we have

$$\begin{equation} [ Q ]\{ {p_{{\rm surf}} } \} = \{ {p_{{\rm ext}} } \} - [ \gamma ]\{ U \},\end{equation} \tag{ 7 }$$

where [Q] is a matrix of weighting coefficients and the acoustic pressures in the exterior volume V_ext have been relabelled {p_ext}. The coefficients in [γ] may be obtained from the Rayleigh integrals relating each plane piston on the multi-element array to the each location in V_ext.

The vector of surface pressures may be eliminated by combining equations (6) and (7):

$$\begin{equation} ( {[ \gamma ] - [ Q ][ H ]^{ - 1} [ \beta ]} )\{ U \} = \{ {p_{{\rm ext}} } \}.\end{equation} \tag{ 8 }$$

Equation (6) may be re-written as follows:

$$\begin{equation} - [ H ]^{ - 1} [ \beta ]\{ U \} = \{ {p_{{\rm surf}} } \}\end{equation} \tag{ 9 }$$

or

$$\begin{equation} [ A ]\{ U \} = \{ {p_{{\rm surf}} } \}\end{equation} \tag{ 10 }$$

where [A] = −[H]⁻¹[β]{U}.

Equation (8) may be re-written as follows:

$$\begin{equation} [ C ]\{ U \} = \{ {p_{{\rm ext}} } \},\end{equation} \tag{ 11 }$$

where [C] = [γ] − [Q][H]⁻¹[β].

All quantities in equation (11) will generally be complex. When investigating an inverse problem involving complex quantities, it is often convenient to reformulate in terms of purely real variables by rewriting equation (11) as follows:

$$\begin{equation} \left[ {\begin{array}{@{}cc@{}} {{\rm Re}[ C ]} & { - {\rm Im}[ C ]} \\ {{\rm Im}[ C ]} & {{\rm Re}[ C ]} \\ \end{array}} \right]\left\{ {\begin{array}{@{}c@{}} {{\rm Re}\{ U \}} \\ {{\rm Im}\{ U \}} \\ \end{array}} \right\} = \left\{ {\begin{array}{@{}c@{}} {{\rm Re}\{ {p_{{\rm ext}} } \}} \\ {{\rm Im}\{ {p_{{\rm ext}} } \}} \\ \end{array}} \right\}\end{equation} \tag{ 12 }$$

or

$$\begin{equation} [ {\hat C} ]\{ {\hat U} \} = \{ {\hat p_{{\rm ext}} } \}.\end{equation} \tag{ 13 }$$

We wish to obtain a set of real and imaginary parts of source velocities which best fit a prescribed field pressure distribution in a least-squares sense such that

• the acoustic pressure magnitudes on S do not exceed a threshold defined by $p_{{\rm surf}_{{\rm max}} }$,
• the source velocity magnitudes do not exceed the upper bound of the dynamic range of each element U_max.

This may be expressed as a least-squares minimization problem with nonlinear constraints:

$$\begin{equation} \mathop {\min }\limits_{\hat U} \frac{1}{2}\|[ {\hat C} ]\{ {\hat U} \} - \{ {\hat p_{{\rm ext}} } \}\|_2^2 \quad {\rm such}\;{\rm that}\quad \left\{ {\begin{array}{@{}cc@{}} {| {\{ {p_{{\rm surf}} } \}} | \le p_{{\rm surf} _{_{{\rm max}} } } } \\ {| U | \le U_{{\rm max}} } \\ \end{array}} \right.\!\!\!\!.\end{equation} \tag{ 14 }$$

Within the context of this work, $p_{{\rm surf}_{{\rm max}} }$ was chosen as 50% of the maximum value of the pressure magnitude on the surface of the ribs resulting from all elements vibrating in phase and with unit amplitude velocity. This 50% value is not necessarily of clinical significance and was arbitrarily chosen to illustrate the technique. It is understood that in clinical applications, a damage threshold related to a dose quantity would have to be experimentally established. The prescribed exterior field pressure distribution {p_ext} was chosen as the incident field (i.e. without the ribs) resulting from the above excitation condition along a 51 × 51 point 2.5 cm × 2.5 cm square grid located in the y–z plane centred at the geometric focus of the array. U_max was chosen as 1 m s⁻¹. In this particular case, the ribs did not cause significant aberrations in the focal plane; hence, the locations at which {p_ext} were evaluated were limited to a grid in the y–z plane. Were significant aberrations to occur in the spherical focusing case, such as splitting at the focus, one could choose a set of points for {p_ext} located over a 3D region such as a cuboid volume.

The HIFU transducer modelled as part of the underlying work was a spherically shaped bowl populated with N = 256 plane circular elements mounted onto its surface. The elements were each of a = 3 mm radius. A radius of curvature of D = 18 cm was used. The outer diameter of the HIFU transducer was chosen as 16 cm. The elements were pseudo-randomly spatially distributed on the surface of the array. The arrangement is identical to that described by Gélat et al (2011). The excitation frequency was assumed to be 1 MHz. A frontal view of the array is shown in figure 1.

Zoom In Zoom Out Reset image size

The Helmholtz integral equation boundary element approach was implemented on a dedicated computer cluster, parallelized over 100 cores and solved iteratively employing a Generalized Minimal Residual (GMRES) scheme in a methodology analogous to Ochmann et al (2003) and Shen and Liu (2007). In all GMRES calculations, the initial trial pressure vector was chosen as the incident field at nodal positions on the surface of the scatterer. 40 iterations were sufficient for the residual norm to fall below 1% of its initial value. The approach used here is essentially the same as that described by Gélat et al (2011) but differs as follows.

• The coupling coefficient in the Burton–Miller (normal derivative) formulation was set to zero.
• The speed of sound in the exterior domain was specified as a complex number.
• Operations were vectorized to solve equation (5) for multiple right-hand sides simultaneously, thus enabling the coefficients of [A] to be obtained more efficiently.

The routines employed to generate the boundary element matrices were obtained from the PAFEC (Program for Automatic Finite Element Calculations) VibroAcoustics software with permission from PACSYS Ltd.

The NAG® Numerical Library e04us routine was used to carry out the constrained minimization. It is designed to minimize an arbitrary smooth sum of squares function subject to constraints, which may include simple bounds on the variables, linear constraints and smooth nonlinear constraints. It employs a sequential quadratic programming method described by Bonnans et al (2006, p 490). The partial derivatives of the constraints and of the cost function with respect to the minimization variables were supplied. Details of the cost function and Jacobian matrices may be found in the appendix.

Initial values of the minimization variables were all specified as $1/\sqrt 2$ m s⁻¹. The problem was scaled by a factor of 10⁻⁶ to avoid numerical instability. In addition to this, the constraints involving the source velocities were scaled by a factor of 10⁶ so that they were of the same order of magnitude as the surface pressure magnitude constraints.

The constrained optimization results based on the methodology described in section 2.2 were compared with the binarized apodization method. The technique employed is similar to that described by Quesson et al (2010). Ray tracing from the geometric focus of the array onto the surface of each element of the HIFU array was carried out. Each element was discretized into 276 points positioned along a regular Cartesian grid within the plane of each element. Rays traced from each point to the focus which came within 0.25 mm of any node on the rib mesh were discarded. This distance corresponds to half the maximum quadratic patch dimension on the rib mesh. Elements with over 50% of their surface shadowed were deactivated. Other elements were assigned a 1 m s⁻¹ normal velocity and zero phase.

The DORT method consists of the construction of the wave fronts that are invariable under a time-reversal process (Prada 2002). For an array comprising N transducers functioning in both transmitting and receiving mode, the time-reversal operator (TRO) matrix can be obtained from the matrix of the inter-element transfer functions [K(ω)], for a time-invariant, linear system. The TRO is then defined as [K(ω)]*[K(ω)], where * denotes the conjugate transpose. The TRO is Hermitian and with positive eigenvalues and its eigenvectors are invariant wave fronts of the time-reversal process. The diagonalization of the TRO is equivalent to the SVD of the array response matrix [K(ω)], from a mathematical point of view. This can be written as follows (Prada 2002).

$$\begin{equation} [ {K( \omega )} ] = [ {V( \omega )} ][ {{\rm \Lambda }( \omega )} ][ {W( \omega )} ]^* ,\end{equation} \tag{ 15 }$$

where

• [V(ω)] and [W(ω)] are unitary matrices and the columns of [W(ω)] are the eigenvectors of the TRO.
• [Λ(ω)] is a real diagonal matrix of the singular values.

In practice, the array response matrix [K(ω)] is measured by emitting a pulse on each array element successively and measuring the corresponding echoes from the scatterers on the N transducers. Cochard et al (2009) employ a time window to select the back-scattered echo from the ribs, and carry out a Fourier transform of the selected signal at each frequency within the transducer bandwidth. Under continuous wave excitation conditions, [K(ω)] can be obtained from equation (8) by considering that the total pressure $\forall \vec r \in V_{{\rm ext}}$ is the sum of the incident pressure and the scattered pressure.

$$\begin{equation} ( {[ Q ][ H ]^{ - 1} [ \beta ]} )\{ U \} = \{ {p_{{\rm scattered}} } \},\end{equation} \tag{ 16 }$$

where {p_scattered} is the vector of scattered pressures at positions on the surface of the HIFU array elements. 276 locations on the surface of each element were considered. These were positioned along a regular Cartesian grid on the surface of each element. {p_scattered} was then evaluated N = 256 times where all rows of {U} were set to zero except for the row corresponding to the source being excited, which was successively set to unity. The scattered pressure was then averaged over the surface of each element, resulting in the N × N matrix of inter-element transfer functions [K(ω)].

After carrying out an SVD on the array response matrix, the singular values of the diagonal matrix [Λ(ω)] were plotted, and a threshold was determined to separate the singular vectors into two categories: those associated with the higher singular values which focus on the ribs, and those associated with the lower singular values which do not send energy to the scatterer. This procedure is described by Cochard et al (2009).

A focusing vector U_focus was then defined. This was chosen as unit amplitude and zero phase velocities for all elements on the HIFU array. Since the intercostal spacing is large compared with the wavelength in the propagating medium at 1 MHz, this surface velocity distribution is already successful at focusing (Gélat et al 2011). Song et al (2007) describe how the emissions focusing on the scatterer may be removed from those that focus on the target using an orthogonal projection in the frequency domain. By projecting the focusing vector orthogonally onto the first set of eigenvectors associated with the higher singular values, we have

$$\begin{equation} \{ {U_{{\rm projected}} } \} = \{ {U_{{\rm focus}} } \} - \mathop \sum \limits_{i = 1}^{i_{{\rm max}} } ( {\{ {W_i } \}^* \{ {U_{{\rm focus}} } \}} )\{ {W_i } \},\end{equation} \tag{ 17 }$$

where

• {W_i} is the ith column of [W],
• i_max is the number of eigenvectors associated with the higher singular values.

A forward calculation was then carried out using the values of {U_projected} after which the acoustic pressures on the surface of the ribs were obtained as were the field pressures at selected locations around the focal region.

A validation was carried out on simple scatterers to justify setting the coupling coefficient in the Burton–Miller formulation to zero. The results displayed here are for a 1 MHz unit amplitude plane wave scattered by a perfectly rigid sphere and a perfectly rigid cylinder with hemi-spherical end-caps. The spherical scatterer was assumed to be of 0.5 cm radius. The cylinder was 22 cm in length and 1 cm in diameter and centred about the x-axis. The boundary element calculations were compared against the analytical solution provided by Morse and Ingard (1968, p 419 and p 401) for the cases of a plane wave scattered by a rigid sphere and by an infinite rigid cylinder. The scatterers were meshed using isoparametric eight-noded quadratic patches ensuring at least three elements per wavelength for a wavespeed of 1500 + 4.405i m s⁻¹, corresponding to an attenuation coefficient of 0.123 Np cm⁻¹ at 1 MHz. The density of the medium was assumed to be 1000 kg m⁻³. The incident fields were plane waves of unit amplitude travelling in the positive x direction for the sphere and in the positive z direction for the cylinder. Comparisons between the analytical solutions and BE calculations in the shadow zone are displayed in figures 2 and 3, for the sphere and the cylinder, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2 shows agreement of peak acoustic pressures within 1.3% of the analytical solution against the Helmholtz integral equation implementation of the BE formulation on a spherical scatterer in a dissipative medium. Figure 3 shows agreement within 0.2% of the analytical solution for an infinite cylinder compared with the BE implementation of the Helmholtz integral equation for a finite cylinder, at 0.2 cm from the surface of the cylinder along the z-axis in the shadow zone. For the finite cylinder, oscillations about the analytical solution are observed. These are due to constructive and destructive interference from the waves diffracted by the end-caps of the cylinder. Figures 2 and 3 display no artefacts of numerical instability associated with the non-uniqueness of the solution of the discretized version of equation (2), demonstrating that any fictitious internal modes of the scatterers are indeed being damped out. It was verified that, by setting the imaginary part of the sound speed to zero, artefacts of the non-uniqueness were clearly visible.

After completing this validation, the incident field for the array shown in figure 1 was computed using a linear superposition of plane circular piston sources, assumed to be rigidly vibrating in an infinite baffle. This approach is similar to that used by Gavrilov and Hand (2000). The main axis of the HIFU array was assumed to be the Cartesian z-axis, the global origin being positioned at the geometric focus of the array. The array was positioned with respect to the ribs to simulate an intercostal treatment between ribs 10 and 11 on the right side, approximately 3 cm deep into rib cage. This configuration was in part chosen to assess the feasibility of the technique, although it is relevant to clinical applications where shallow tumours may occur. These may be more problematic to treat than more deep-seated ones because the higher energy density at the skin increases the risk of skin burn. This arrangement is illustrated in figure 4.

Zoom In Zoom Out Reset image size

The acoustic field positions of interest are located at distances from the individual array sources much greater than the radius of each individual element (a = 3 mm). The analytical expression described by Kinsler et al (1982, p 179) was therefore used for the far-field approximation. The incident pressure field (i.e. unobstructed by ribs) in the y–z plane is shown in figure 5 for all elements vibrating at an amplitude of 1 m s⁻¹ and with zero phase, where the focal pressure is 4.4 MPa.

Zoom In Zoom Out Reset image size

A simulation solving equation (5) for multiple right-hand sides was launched so that obtaining the product [H]⁻¹[β] would not necessitate 256 separate runs. Due to RAM limitations, it was however not possible to solve for 256 right-hand sides simultaneously, so four sets of jobs with 64 right-hand sides were run. Each of these jobs took about one day to complete.

From a knowledge of [H]⁻¹[β], the acoustic pressure on the surface of the ribs may be obtained using equation (6) for any given distribution of source velocities. In order to formulate the cost function and the constraints, the case of all elements vibrating in phase and with unit velocity amplitude (spherical focusing) was investigated in the presence of ribs. The acoustic pressure magnitudes at selected locations in the y–z plane are shown in figure 6. The focal pressure is 3.8 MPa.

Zoom In Zoom Out Reset image size

To save on post-processing times, acoustic pressures within 8 mm of the surface of the ribs were not calculated. This would necessitate using a higher integration order in equation (3) to ensure accuracy, which would lead to longer run times.

The resulting acoustic pressure magnitude on the surface of the ribs is shown in figure 7.

Zoom In Zoom Out Reset image size

Figures 5 and 6 show that inserting the ribs between the array and the focus causes the acoustic pressure magnitude at the focus to drop from 4.4 to 3.8 MPa. Figure 7 shows that the maximum pressure amplitude on the surface of the ribs is 1.8 MPa for the spherical focusing case, with all elements of the array vibrating at an amplitude of 1 m s⁻¹ and with uniform phase.

To fulfil the requirement of a 50% reduction in the rib surface pressure magnitude compared to the spherical focusing case, the value $p_{{\rm surf} _{_{\max }} }$ in equation (14) was chosen as 0.9 MPa. U_max was selected to be 1 m s⁻¹. The vector of pressures in the exterior volume featured in the cost function (i.e. the 'desired' field pressure distribution) was generated from pressure field values in absence of ribs at 2601 equally spaced locations in the y–z plane such that −2.5 cm ≤ y ≤ 2.5 cm and −2.5 cm ≤ z ≤ 2.5 cm, the focus of the array being at the global origin.

The surface mesh of the ribs contains 250 370 nodes. It was impractical to impose a constraint at all these locations on the surface of the scatterer as this would result in numerical instability. Surface locations associated with pressure magnitudes below 40% of the maximum pressure in figure 7 (i.e. 40% of 1.8 MPa) were not included in the constraints. This resulted in only 2635 constraints associated with the pressure magnitude on the ribs being required, along with the 256 constraints for the magnitude of the element velocities.

With the above input conditions, a solution was obtained for which both sets of constraints were satisfied. The real and imaginary parts of the velocities were subsequently rescaled and the surface and field pressures calculated. Figures 8(a) and (b) respectively display the magnitudes and phases of the element velocities of the array resulting from the constrained minimization. It was verified that the velocity magnitudes did not exceed 1 m s⁻¹.

Zoom In Zoom Out Reset image size

Figure 9 shows the acoustic pressure magnitude in the y–z plane resulting from the source velocity distribution displayed in figures 8(a) and (b).

Zoom In Zoom Out Reset image size

Figure 10 shows the acoustic pressure magnitude on the surface of the ribs, resulting from the source velocity distribution displayed in figures 8(a) and (b), which was verified not to exceed 0.9 MPa at all 250 370 nodal positions.

Zoom In Zoom Out Reset image size

The constrained minimization algorithm is unsuccessful at rendering a focal peak pressure of 4.4 MPa as specified by the cost function in equation (14). This result is not unexpected since constraints have been imposed on the element velocity magnitudes and the presence of ribs will inevitably prevent some of the energy from being deposited at the focus. Nevertheless, a 50% drop in the maximum acoustic pressure on the surface of the ribs was achieved with the acoustic pressure magnitude at the focus reduced by only 4% in figure 9 compared with the spherical focusing case in figure 7.

Currently, the technique is computationally expensive. About four days of pre-processing are required to obtain the required matrices used in the inverse problem. Run times associated with the actual constrained minimization are not long (approximately 5 min on a 16 GB RAM, two-core Dell® Precision™ workstation with 3 GHz Intel® Xeon™ processors). It is acknowledged that this method is currently unsuitable for real time clinical applications and is more likely to be of relevance when incorporated into a treatment planning software and associated treatment protocol.

The ray tracing procedure described in section 2.3 led to the deactivation of 28% of the elements of the HIFU array. The velocity distributions are displayed in figures 11(a) and (b).

Zoom In Zoom Out Reset image size

Figure 12 shows the acoustic pressure magnitude in the y–z plane resulting from the source velocity distribution displayed in figures 11(a) and (b).

Zoom In Zoom Out Reset image size

Figure 13 shows the acoustic pressure magnitude on the surface of the ribs, resulting from the source velocity distribution displayed in figures 11(a) and (b). This results in a maximum pressure magnitude on the surface of the ribs of 1.0 MPa.

Zoom In Zoom Out Reset image size

The binarized apodization approach based on ray tracing provides a 41% reduction in the peak acoustic pressure on the surface of the ribs compared with the reference case shown in figure 7. This, however, results in a 3.2 MPa acoustic pressure magnitude at the focus, which corresponds to an 18% drop compared with the spherical focusing case. These results suggest that binarized apodization may not offer the flexibility of the constrained minimization approach in terms of maximizing focal pressures whilst keeping rib surface pressure below a specified threshold. By comparing figures 8(a) and 11(a), it is worth noting that the ray tracing approach identifies a larger number of 'problematic' sources than the least-squares minimization, thus highlighting further limitations of this implementation of binarized apodization. In reality, the behaviour of the sources and their interaction with the ribs is more complex than that predicted by ray acoustics and it is likely that this specific array–rib geometry results in the constructive interference observed on the surface of rib 11 in figure 13.

The array response matrix [K(ω)] described in section 2.4 was formulated using equation (16). An SVD was carried out and the singular values of the diagonal matrix [Λ(ω)] were obtained. The results are plotted in figure 14.

Zoom In Zoom Out Reset image size

The threshold i_max which separated the singular vectors into those associated with the higher singular values which focus on the ribs, and those associated with the lower singular values which do not send energy to the scatterer, was set to 25. This maximized the ratio of maximum pressure at the focus to maximum pressure on the ribs. The element velocities shown in figures 15(a) and (b) were obtained using the spherical focusing case as the focusing vector and by projecting the focusing vector orthogonally onto the first set of eigenvectors associated with the higher singular values according to equation (17). These were rescaled so that the maximum velocity amplitude was 1 m s⁻¹, since this criterion defined the upper limit of each element's dynamic range.

Zoom In Zoom Out Reset image size

Figure 16 shows the acoustic pressure magnitude in the y–z plane corresponding to the focusing vectors shown in figures 15(a) and (b).

Zoom In Zoom Out Reset image size

Figure 17 shows the acoustic pressure magnitude on the surface of the ribs, resulting from the source velocity distribution displayed in figures 15(a) and (b), which results in a maximum pressure magnitude of 1.1 MPa on the surface of the ribs.

Zoom In Zoom Out Reset image size

The DORT method offers a 37% reduction in the peak acoustic pressure on the surface of the ribs compared with the reference case in figure 7 for spherical focusing. Nevertheless, this implementation results in a 2.7 MPa acoustic pressure magnitude at the focus, which is a 29% drop compared to the spherical focusing case and 26% lower than the constrained optimization case. These results suggest that, although the DORT method does not require any a priori knowledge of the rib topology to avoid energy being deposited on the surface scatterer, the ratio of focal pressure magnitude to maximum pressure magnitude on the surface of the ribs is inferior to that resulting from the constrained minimization approach in this particular example.

Axial scans of the acoustic pressure magnitudes centred at the global origin resulting from spherical focusing, the constrained minimization, binarized apodization and from the DORT method are displayed in figures 18, 19 and 20, along the x-, y- and z-axes, respectively. These scans were normalized against the spherical focusing case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For each focusing method, the maximum values of the acoustic pressures on the surface of the ribs relative to the maximum field pressures are displayed in table 1. The values of the maximum pressures at the focus relative to the spherical focusing case are also displayed in table 1.

The ribs considered as part of these calculations were assumed to be perfectly rigid. As well as being highly reflecting, it is also well known that human ribs are highly attenuating. Although a full visco-elastic model would be necessary to tackle this problem, a surface impedance condition in the $\partial p( {\vec r} )/\partial n$ term in equation (2), $\forall \vec r \in S$ could provide a means of simulating attenuation at the surface of the ribs and consequently estimate the acoustic dose-rate (Duck 2009) deposited. Efforts are currently being made to incorporate this in the BE routines. A limitation of this approach is that potential effects of mode conversion due to the tissue–bone interface would be unaccounted for. Also, the method cannot currently monitor heat deposition on the surface of the ribs. An optimization scheme which would ensure that the thermal dose (Duck 2009) on the surface of the ribs was kept below a given threshold would be considerably more involved, as this would necessitate interaction with a heat transfer model.

The surrounding medium was assumed to be homogeneous with an attenuation coefficient of 0.123 Np cm⁻¹ at 1 MHz. It should be noted that acoustic pressures predicted as part of this work may be considerably underestimated, as the water path between the array and the abdomen which would normally be present during treatment, and in which attenuation effects will be much smaller, was unaccounted for. Additionally, the upper bound of 1 m s⁻¹ on the element velocities was not intended to be a realistic value but was instead chosen for demonstrative purposes. Overall, this method could be generalized to include a combination of exterior and interior boundary element formulations in order to investigate the propagation of ultrasound in multiple regions. This would necessitate meshing the closed surfaces associated with each region (e.g. skin, fat, tissue, bone).

The effects of nonlinear propagation were not considered in this paper. It is however well known that HIFU fields can result in highly nonlinear behaviour in the focal region (Wu et al 2004) leading to distortion of the acoustic waveform and transfer of energy from the fundamental frequency towards higher order harmonics. Nevertheless, from the point of view of treatment planning, it has been shown that correcting for aberrations introduced by the presence of bone need not require nonlinear propagation models and that linear approaches generally suffice (Pinton et al 2011). Furthermore, nonlinear behaviour is likely to be mainly confined to the central focal lobe. Hence, although a more rigorous treatment of the focal region may be required, useful information can be obtained from linear models when investigating energy depositions at the surface of the ribs. A full-wave nonlinear 3D model for propagation of ultrasound in tissue based on a finite volume approach is currently being developed alongside the boundary element method described in this paper. This model allows for shock wave formation and for sharp discontinuities in the media interfaces (e.g. tissue/bone).

This work was carried out as part of a wider project whose goal is the design of a prototype transcostal HIFU system for the treatment of liver and pancreatic tumours. This will feature a multi-element array with electronic steering capabilities together with a system for ultrasonic ablation/cavitation monitoring. The technique presented in this paper is expected to help validate the treatment planning software and inform the treatment protocol. In vivo applications will be investigated and are expected to rely upon treatment plans generated using virtual propagation models based on anatomical data that has been obtained from CT or MRI scans of the patient. Treatment monitoring will be performed using an ultrasonic probe either integrated into the HIFU transducer, or placed in a separate position (e.g. to scan from a subcostal approach). Registration of the ribs from CT images to the patient will be carried out as follows. First, a diagnostic ultrasound probe, tracked by an optical 3D localizer, will be used to acquire a 3D ultrasound image of the rib surface (the 3D pose of each of a series of B-mode ultrasound images is calculated using the 3D tracking data, and the images combined to form a 3D image of the region of interest). The CT bone surface may then be registered directly and automatically to the ultrasound image using the method described by Penney et al (2006). Alternative methods could also be used since the CT scans are pre-segmented to produce a geometric representation of the bone surface as part of treatment planning and acoustic modelling tasks. Since the HIFU transducer is also tracked in 3D by the localizer, and registering the ribs using this procedure allows the ribs to be localized with respect to the 3D physical reference co-ordinate system of the localizer, the correct (planned) physical position and orientation of the HIFU transducer relative to the ribs may be determined.

Initially, effects of respiratory movement and organ deformation will be ignored and breath-hold or gating will be used. A subsequent version of the treatment planning software will feature motion compensation capabilities based on methods developed by Rijkhorst et al (2011).

The cost function F is a function of 2N real variables (elements of the vector $\{ \hat U\}$), which returns a real scalar. It is defined as follows:

$$\begin{equation} F( {\{ {\hat U} \}} ) = \|[ {\hat C} ]\{ {\hat U} \} - \{ {\hat p} \}\|^2 \end{equation} \tag{ A.1 }$$

or

$$\begin{equation} F( {\hat U_1 ,\hat U_2 , \ldots ,\hat U_{2N} } ) = \mathop \sum \limits_{i = 1}^{2M'} \left( {\mathop \sum \limits_{j = 1}^{2N} \hat C_{i,j} \hat U_j - \hat p_i } \right)^2 .\end{equation} \tag{ A.2 }$$

The constraints, denoted as c, are nonlinear and are required to be input in the format

$$\begin{equation} c_i \le 0.\end{equation} \tag{ A.3 }$$

Equation (8) is rewritten as follows:

$$\begin{equation} \left[ {\begin{array}{@{}cc@{}} {{\rm Re}( {[ A ]} )} & { - {\rm Im}( {[ A ]} )} \\ {{\rm Im}( {[ A ]} )} & {{\rm Re}( {[ A ]} )} \\ \end{array}} \right]\left\{ {\begin{array}{@{}c@{}} {{\rm Re}( {\{ U \}} )} \\ {{\rm Im}( {\{ U \}} )} \\ \end{array}} \right\} = \left\{ {\begin{array}{@{}c@{}} {{\rm Re}( {\{ {p_{{\rm surf}} } \}} )} \\ {{\rm Im}( {\{ {p_{{\rm surf}} } \}} )} \\ \end{array}} \right\}\end{equation} \tag{ A.4 }$$

or

$$\begin{equation} [ {\hat A} ]\{ {\hat U} \} = \{ {\hat p_{{\rm surf}} } \}.\end{equation} \tag{ A.5 }$$

The ith row and the jth columns of Re([A]) and Im([A]) are defined as A_{R i,j} and A_{I i,j}, respectively. Similarly, let the jth row of Re({U}) and Im({U}) be defined as U_{R i,j} and U_{I i,j}, respectively.

For all M nodes on the surface S, we require the acoustic pressure magnitude to be less than a specified threshold $p_{{\rm surf}_{{\rm max}} }$. There are therefore M constraints associated with this condition. These can be expressed as follows.

For i ≤ M

$$\begin{equation} \fl c_i = \left( {\mathop \sum \limits_{j = 1}^N A_{Ri,j} U_{Rj} - \mathop \sum \limits_{j = 1}^N A_{Ii,j} U_{Ij} } \right)^2 + \left( {\mathop \sum \limits_{j = 1}^N A_{Ii,j} U_{Rj} + \mathop \sum \limits_{j = 1}^N A_{Ri,j} U_{Ij} } \right)^2 - p_{{\rm surf}_{{\rm max}} }^2 .\end{equation} \tag{ A.6 }$$

In practice, numerical instability may arise if M is too large. In some cases, a subset of the nodes on S may have to be chosen.

In addition, the source velocity magnitudes cannot be greater than the maximum amplitude U_max allowed by the dynamic range. This is equivalent to adding the N further constraints.

For M < i ≤ M + N

$$\begin{equation} c_i = U_{Ri - M}^2 + U_{Ii - M}^2 - U_{{\rm max}}^2 .\end{equation} \tag{ A.7 }$$

The gradient of the cost function with respect to its input variables may be expressed as follows.

For 1 ≤ k ≤ 2N,

$$\begin{equation} \frac{{\partial F}}{{\partial \hat U_k }} = 2\mathop \sum \limits_{i = 1}^{2M'} \hat C_{i,k} \mathop \sum \limits_{j = 1}^{2N} ( {\hat C_{i,j} \hat U_j - \hat p_i } ).\end{equation} \tag{ A.8 }$$

The gradient of the constraints with respect to the input variables forms a 2N by M + N Jacobian matrix (2N variables and M + N constraints). Its elements are as follows.

For k ≤ N and i ≤ M,

$$\begin{equation} \fl\frac{{\partial c_i }}{{\partial \hat U_k }} = 2A_{Ri,k} \left( {\mathop \sum \limits_{j = 1}^N A_{Ri,j} U_{Rj} - \mathop \sum \limits_{j = 1}^N A_{Ii,j} U_{Ij} } \right) + 2A_{Ii,k} \left( {\mathop \sum \limits_{j = 1}^N A_{Ii,j} U_{Rj} + \mathop \sum \limits_{j = 1}^N A_{Ri,j} U_{Ij} } \right).\end{equation} \tag{ A.9 }$$

For k > N and i ≤ M,

$$\begin{eqnarray} &&\fl\frac{{\partial c_i }}{{\partial \hat U_k }} = - 2A_{Ii,k - N} \left( {\mathop \sum \limits_{j = 1}^N A_{Ri,j} U_{Rj} - \mathop \sum \limits_{j = 1}^N A_{Ii,j} U_{Ij} } \right) + 2A_{Ri,k - N} \left( {\mathop \sum \limits_{j = 1}^N A_{Ii,j} U_{Rj} + \mathop \sum \limits_{j = 1}^N A_{Ri,j} U_{Ij} } \right). \nonumber\\ \end{eqnarray} \tag{ A.10 }$$

For k ≤ N and i > M,

$$\begin{equation} \frac{{\partial c_i }}{{\partial \hat U_k }} = 2U_{Rk} \quad {\rm if}\;i = k + M,\end{equation} \tag{ A.11 }$$

$$\begin{equation} \frac{{\partial c_i }}{{\partial \hat U_k }} = 0\quad {\rm if}\;i \ne k + M.\end{equation} \tag{ A.12 }$$

For k > N and i > M,

$$\begin{equation} \frac{{\partial c_i }}{{\partial \hat U_k }} = 2U_{Ik - M} \quad {\rm if}\;i = k + M - N,\end{equation} \tag{ A.13 }$$

$$\begin{equation} \frac{{\partial c_i }}{{\partial \hat U_k }} = 0\quad {\rm if}\;i \ne k + M - N.\end{equation} \tag{ A.14 }$$