Estimation of Thickness and Speed of Sound for Transverse Cortical Bone Imaging Using Phase Aberration Correction Methods: An In Silico and Ex Vivo Validation Study

Abstract

Delay-and-sum (DAS) beamforming of backscattered echoes is used for conventional ultrasound imaging. Although DAS beamforming is well suited for imaging in soft tissues, refraction, scattering, and absorption, porous mineralized tissues cause phase aberrations of reflected echoes and subsequent image degradation. The recently developed refraction corrected multi-focus technique uses subsequent focusing of waves at variable depths, the tracking of travel times of waves reflected from outer and inner cortical bone interfaces, the estimation of the shift needed to focus from one interface to another to determine cortical thickness (Ct.Th), and the speed of sound propagating in a radial bone direction (Ct.ν₁₁). The method was validated previously in silico and ex vivo on plate shaped samples. The aim of this study was to correct phase aberration caused by bone geometry (i.e., curvature and tilt with respect to the transducer array) and intracortical pores for the multi-focus approach. The phase aberration correction methods are based on time delay estimation via bone geometry differences to flat bone plates and via the autocorrelation and cross correlation of the reflected ultrasound waves from the endosteal bone interface. We evaluate the multi-focus approach by incorporating the phase aberration correction methods by numerical simulation and one experiment on a human tibia bone, and analyze the precision and accuracy of measuring Ct.Th and Ct.ν₁₁. Site-matched reference values of the cortical thickness of the human tibia bone were obtained from high-resolution peripheral computed tomography. The phase aberration correction methods resulted in a more precise (coefficient of variation of 5.7%) and accurate (root mean square error of 6.3%) estimation of Ct.Th, and a more precise (9.8%) and accurate (3.4%) Ct.ν₁₁ estimation, than without any phase aberration correction. The developed multi-focus method including phase aberration corrections provides local estimations of both cortical thickness and sound velocity and is proposed as a biomarker of cortical bone quality with high clinical potential for the prevention of osteoporotic fractures.

1. Introduction

The current standard method for bone strength assessment and fracture risk prediction is based on areal bone mineral density (aBMD) measured by dual-energy absorptiometry (DXA) [1]. Although aBMD is an important biomarker of bone quality, additional bone factors, including macro- and micro-structural bone parameters, as well as viscoelastic properties, are known to determine individual bone strength; therefore, to quantify these parameters for bone assessment, quantitative ultrasound (QUS) methods have been introduced as nonionizing alternatives. Early bone QUS technologies used dedicated hardware to measure acoustic properties, such as the speed of sound (SOS) and broadband ultrasound attenuation (BUA), at anatomical sites that contain mostly trabecular bone, such as the heel [2]. More recent QUS devices are aimed at imaging bone by using dedicated hardware electronics and ultrasound probes. An example of this by Lasaygues et al. developed ultrasonic image reconstruction methods to image the cortical diaphysis of long bones using quantitative ultrasonic tomography [3,4]. Another tomographic approach to image long bones is based on full-waveform inversion [5]. Additionally, Li et al. used Split-Step Fourier imaging to image bone fractures and to monitor bone healing [6]. Furthermore, a Born-based inversion method has been implemented on an ultrasonic wavefield imaging technique to reconstruct internal structures of long bones [7]. Limitations of these studies were that either the speed of sound or the thickness needed to be assumed a priori. Axial transmission devices can retrieve cortical parameters (i.e., porosity, thickness, and speed of sound), by measuring the propagating velocity of dispersive guided waves [8,9,10,11,12]; however, this technique is challenged by large soft tissue thickness, irregular bone shapes, and it does not provide direct image guidance.

A few recent technologies utilize sophisticated array-based pulse-echo imaging technology to estimate BMD in trabecular bones at major fracture sites (i.e., spine and proximal femur [13]), or to measure structural and material properties in the cortical bone (i.e., tibia and radius) [14,15]. Most medical ultrasound scanners on the market implement the standard delay-and-sum (DAS) beamforming method to reconstruct the brightness mode (hereinafter called B-mode) images. This technique uses a transducer array to transmit and receive focused ultrasound signals inside the body. Conventionally, the reconstruction of B-mode images using DAS is done by adding time specific delays to the individual ultrasound signals which are received at each element of the receiver array before summing all signals to create a beamformed received signal; therefore, the sensitivity of the beamformed signal can be maximized to a certain depth and direction. In medical ultrasound scanners, transmit and receive focusing is performed by assuming a constant speed of sound of soft tissue (1540 m/s) along the entire sound propagation path. This approximation provides satisfactory image quality for most soft tissues, because the true velocities only vary within 10% when compared with the assumed value [16]; however, this is not the case for mineralized tissues, such as cortical bones. The radial speed of sound in cortical bone is between the range of 2800 to 3500 m/s [17], which results in a substantial refraction at the soft tissue and cortical bone interface. In case of a wrong assumption on the constant sound velocity, the delay estimation, which is necessary to focus on a particular image location, is incorrect, subsequently leading to a phase-distorted DAS signal. As a result, in a conventional B-mode image reconstructed by medical ultrasound scanners, the internal bone structures appear blurred or cannot be reconstructed at all. Aside from radiofrequency echographic multispectrometry (REMS) technology [18], this conventional DAS beamforming is currently used for bone strength assessment and fracture risk prediction. There have been efforts to overcome this false assumption of a constant speed of sound in cortical bone. Renaud et al. [14] proposed the first in vivo image reconstruction of cortical bone using a conventional medical ultrasound scanner and seismic image reconstruction. This reconstruction method provides local estimations of Ct.Th and anisotropic sound velocity.

Consequently, the need for further methods brought about the multi-focus (MF) imaging technique, that was developed by our group to measure cortical thickness (Ct.Th) and the compressional sound velocity propagating in the radial bone direction (Ct.ν₁₁) [19]. Our method aims at imaging cortical bone at the central anteromedial tibia. This anatomical site is of clinical interest, as it is easy to access, and is composed of a thick and regular cortical bone shell. Alterations, such as reduced cortical thickness and the occurrence of large intracortical pores, have shown to be associated with reduced hip strength [20] and increased fracture risk [21]. The ultrasonic speed of sound in cortical bone has been proposed as a biomarker of bone quality since the 1970s [22], and is related to bone density and elastic constants, which are correlated to bone quality and fracture risk [23]. The MF method is based upon the consecutive focusing of ultrasound waves at varying depths, followed by the retrieval of focus locations and pulse-travel times of signals reflected from the periosteal (frontside) and endosteal (backside) cortical bone interfaces using conventional DAS beamforming (Figure 1). So far, the MF method has been validated on plates with constant thickness, positioned parallel to the probe array; however, it is important to extend this, as typical human cortical bones exhibit curvatures at their periosteal and endosteal interfaces. These curvatures introduce a distortion of the propagating wavefronts and the round-trip travel time, resulting in phase distorted beamformed signals. The objective of this research is to analyze the effect of the phase aberration caused by (a) bone curvature, (b) bone tilt with respect to the beam axis, and (c) material inhomogeneities due to the presence of cortical pores on the estimations of Ct.Th and Ct.ν₁₁. The phase aberration from bone surface curvature leads to a different round-trip travel time when compared with a flat plate bone model (surface time shift, ST), and was corrected using the concept of refraction compensation proposed by Yasuda et al. [24]. Bone tilt, with respect to the beam axis of the transducer array, shows orientation dependence of the received echoes compared with a flat bone interface. Additional time shifts caused by the orientation dependence of received echoes were corrected using autocorrelation analysis (ACF). Differences in the round travel time of the received echoes based on the interaction of ultrasound wave refractions with cortical pores were determined using cross correlation analysis (CC). We show the need to incorporate three phase-aberration correction (PAC) methods for non-plate shaped bone structures by means of numerical finite-difference time-domain (FDTD) simulation models, as well as measurements on a human tibia bone. Precision and accuracy values of estimated Ct.Th and Ct.ν₁₁ with and without corrections were compared.

2. Materials and Methods

An overview for all used abbreviations is summarized in

Table 1. List of abbreviations.

Abbreviation	Description
Ct.Th	Cortical thickness
Ct.ν11	Cortical compressional sound velocity propagating in the radial bone direction
νH2O	Speed of sound in water
dx	Lateral shift of center of mass of curved bone plate model relative to beam axis
r	Bone plate curvature radius
Ct.Po	Cortical porosity
Fz	Focus depth in z-direction
HF(Fz)	Amplitude of Hilbert-transformed envelope signal of beamformed frontside reflection at focus depth Fz
HB(Fz)	Amplitude of Hilbert-transformed envelope signal of beamformed backside reflection at focus depth Fz
FB	Front- and backside reflection
ΔTOF	Shift in time-of-flight between peak position of HF(Fz) and HB(Fz)
ΔFz	Shift in focus depth between peak position of HF(Fz) and HB(Fz)
Fz,B	Confocal focus depth position of backside reflection
θ	Semi-aperture angle of transmit and receive beams
keff	Correction factor keff for effective aperture keffθ
θcrit	Critical angle based on Snell’s law
Δθ	Difference of the semi-aperture angle to the critical angle
Txi, Rxi	Transmit or receive channel number
Rxref	Reference receive channel with maximum amplitude at envelope signal of pre-beamformed backside reflection
Vgb	Gated pre-beamformed backside reflection signals
VACF	Signal after using autocorrelation function (ACF)
|VACF|	Magnitude of the ACF signal
αACF	Inclination angle of the fitted ellipsoid on VACF to the major semi-axis
ΔtACF	Time shift correction based on αACF

.

2.1. Numerical Ultrasound Propagation Model

Ultrasound wave propagation in bone and water was simulated using a 2D finite-difference time-domain (FDTD) code (Simsonic, www.simsonic.fr, accessed on 10 March 2022) [25]. The simulation model considers elastic wave propagation including mode conversion, multiple scattering, frequency-independent absorption, refraction, and diffraction. A convergence study, as described in [19], provided stable results at grid sizes of 7 µm and time steps of 0.93 ns.

Table 2. Tissue material properties of bone and pores used for the numerical model. Mass density ρ, and c_ij (i.e., the coefficients of a transverse isotropic stiffness tensor were taken from [27] and the absorption value α was obtained from [26]) (adapted from [19] under the Creative Commons Attribution 4.0 license).

	Bone	Pores/Water
ρ [g/cm3]	1.93	1.00
c11 [GPa]	23.7	2.25
c22 [GPa]	23.7	2.25
c12 [GPa]	9.5	2.25
c66 [GPa]	6.6	0
ν11 [m/s]	3504	1500
α [dB/mm]	2.1	0.002

shows the material properties used for the models in this study. Material properties were used from an ex vivo study [26] and a previous acoustic microscopy study in a human femoral cortical bone [27]. All bone models were simulated as hollow cylinders immersed in water. The cylinders were defined by an outer curvature radius r and a wall thickness d. All bone models were placed 15 mm below a linear array with 64 transmitter and receiver elements (element and pitch sizes: 0.3 mm); therefore, the models assumed the sound propagation in the transverse image plane (i.e., perpendicular to the bone’s long axis, at the antero-medial midshaft of a tibia, where the outer bone surface is flat or slightly curved and the sound velocity of the tissue matrix can be assumed to be isotropic in the simulation plane). The transducer elements emitted broadband pulses with a center frequency of 4 MHz and a −6-dB bandwidth of 60%. Phase delays were applied to focus the transmit beam consecutively, at depths ranging from 13 mm to 40 mm, with an increment of 1 mm. The signals received by all elements were captured and downsampled to a sampling rate of 80 MHz for further processing. The sufficient aperture size of 64 was chosen based on a side study, which can be found in Appendix A.

2.1.1. Reference Bone Model: Flat Bone Plate

The reference model consisted of a 4 mm thick bone plate (Ct.Th^Ref = 4 mm) without pores. The material properties of the homogenous bone material results in a reference speed of sound of Ct.ν₁₁^Ref = 3504 m/s. The curvature radius of r = 10 m was used to simulate a flat bone plate (hereinafter simply called ‘flat bone plate’). The radius of 10 m was deemed sufficiently large to exhibit a negligible curvature within the simulation region (Figure 2a).

2.1.2. Bone Curvature

To investigate the effect of bone curvature, curved bone plate models were simulated and compared with the flat bone plate model (Figure 2b). Five curved bone models with radii of r = 60 mm, 50 mm, 40 mm, 30 mm, and 20 mm were simulated. The radius range was defined based on a previous study, in which human tibia midshaft bones of 55 postmenopausal women were measured by means of high-resolution peripheral computed tomography (HR-pQCT) [15]. In that study, the anteromedial tibia midshaft region had been chosen as the ultrasound measurement site due to the small amount of overlying soft tissue and the small curvature of the bone surface compared with other tibia regions. To estimate the curvature radius, circular fits were performed on the central anteromedial tibia region. Tibia bone curvature radii were found to be in the range between 12.6 mm and 68.8 mm with a mean radius of 30.3 mm. Three examples of the circular fits on the HR-pQCT scans are shown in Figure A2 of Appendix B, where the subjects with a minimum (Figure A2a), mean (Figure A2b), and maximum curvature radius (Figure A2c) were selected.

2.1.3. Bone Tilt

To study the effect of the angle of incidence, a bone surface tilt was incorporated by shifting the lateral position of the transducer array by dx (Figure 2c). The bone surface tilt was defined as the angle between the normal vector of the periosteal bone surface and the beam axis at their crossing point.

2.1.4. Material Inhomogeneity: Cortical Pores

To study the effect of material inhomogeneity, cortical pores were included in the curved bone plate models (Figure 2d). Previous ex vivo studies in human cortical bone reported cortical porosity (Ct.Po) and cortical pore diameter (Ct.Po.Dm) values between 2% and 22% and 7 and 95 µm, respectively [28,29,30]. Cortical pores were defined as circular pores with Ct.Po.Dm = 60 µm and varying pore densities, resulting in models with Ct.Po values ranging from 0% to 20% with an increment of 2%.

For the simulation models with cortical pores, transmission measurements were performed to calculate the reference speed of sound Ct.ν₁₁^Ref. An unfocused single-element transducer with a width of 0.3 mm emitted ultrasound waves with a center frequency of 4 MHz and a −6-dB bandwidth of 60%. The unfocused ultrasound wave traveled though the bone and the transmitted ultrasound wave was captured by a single element detector with a width of 0.3 mm, which was placed below the bone. The transducer and detector were placed at the beam axis of the reference MF simulation. In addition, a simulation was performed with the same configuration without the bone to measure the reference signal transmitted though water. The time-of-flight of the ultrasound wave transmitted through water TOF_H20 and bone TOF_bone was defined at the time of the maximum of the signal envelope. The Ct.ν₁₁^Ref of the bone models with pores were calculated using the following equation from [31].

$$Ct{.\mathit{\nu }}_{11}{}^{\mathit{Ref}} = \frac{\mathit{Ct}{.\mathit{Th}}_{\mathit{bone}}}{\frac{\mathit{Ct}{.\mathit{Th}}_{\mathit{bone}}}{{\mathit{\nu }}_{H_{2}O}} + \left({\mathit{TOF}}_{\mathit{bone}} - {\mathit{TOF}}_{H_{2}O}\right)},$$

(1)

with Ct.Th_bone = 4 mm.

2.2. Ex Vivo Measurement on a Human Tibia Bone

One left tibia bone from a human cadaver (female, age 85) was used for the ex vivo validation. The bone sample was received without the soft tissue and distal end (cut off at approximately 50%). The sample was collected by the institute of Anatomy, University of Lübeck, Germany, in accordance with the German law “Gesetz über das Leichen-, Bestattungs- und Friedhofswesen des Landes Schleswig-Holstein II Abschnitt, §9 Leichen-öffnung, anatomisch”, from 2 April 2005. A 30 mm disk was cut from the tibia midshaft using a band saw (EXACT GmbH, Remscheid, Germany). A HR-pQCT scan was performed (XtremeCT II, Scano Medical AG, Bassersdorf, Switzerland) with a total scan length of 10.2 mm in the axial direction and an isotropic voxel size of 60.7 µm. Cortical thickness at the anteromedial tibia section was extracted using a custom protocol adapted from Iori et al. [32] and used as reference value. Cortical porosity was calculated from the HR-pQCT scan using the algorithm proposed by Burghardt et al. [33]. A site-matched multi-focus measurement was performed using a medical scanner SonixTouch equipped with a 3D linear array transducer 4DL14-5/38 (consisting of a 1D 128 element array) and a SonixDAQ single-channel data acquisition system (Ultrasonix, Richmond, BC, Canada). The SonixDAQ allows the pre-beamformed single-channel radio frequency (RF) data acquisition of all channels without any signal processing. The sample was immersed in water and the transducer array was positioned perpendicular to the bone’s long axis. The multi-focus measurement sequence consisted of a series of conventional B-mode imaging sequences with 128 lateral scan positions. At each scan position, sound waves were focused on the radial bone where the direction was into the tibia sample using a 64-element transmit aperture. Subsequent B-mode images were acquired using 17 gradually increasing focus depths (starting from 14 mm with a step size of 2 mm). The transducer elements were excited with a “+−” signal at a system transmit frequency of 4 MHz to optimize the penetration depth. Single-channel RF data were captured with all 128 array elements at a sampling rate of 40 MHz with a 12-bit resolution.

2.3. Signal Processing

2.3.1. Reference Bone Model: Flat Bone Plate

Details of the multi-focus signal processing steps have been described previously [19]. From the delay and beamformed (DAS) Hilbert-transformed envelope signal, the amplitudes [H_F(F_z) and H_B(F_z)] and pulse travel times [TOF_F(F_z) and TOF_B(F_z)] of the signals reflected from the front- and backsides of the plate were tracked for each beam focus position F_z. The time-of-flight difference between front- and backside reflections was defined as ΔTOF = TOF_B(F_z) − TOF_F(F_z). Spline interpolation was used to estimate H_F(F_z) and H_B(F_z) at an F_z increment of 0.1 mm. The interpolated data, and the front- and backside focus positions F_F and F_B, respectively, were retrieved from the peak positions of H_F(F_z) and H_B(F_z), and ΔF_z (i.e., the shift needed to focus either on the front- or backside of the plate, and to estimate the time delay between front- and backside reflections ΔTOF). Ct.Th and Ct.ν₁₁ were estimated using Equation (3) in [19] with sound velocity in water ν_H20:

$$\mathit{Ct}.\mathit{Th} = \frac{\mathsf{\Delta }{F}_z}{0.5·\frac{\mathrm{Ct}{.\mathit{\nu }}_{11}}{{\mathit{\nu }}_{H_{2}O}}·\left(1 - \frac{\mathit{Ct}{.\mathit{\nu }}_{11}{}^2}{{\mathit{\nu }}_{H_{2}O}{}^2}\right)·\left(1 - \cos \left(k_{\mathit{eff}}\theta \right)\right)-\frac{\mathit{Ct}{.\mathit{\nu }}_{11}}{{\mathit{\nu }}_{H_{2}O}}},$$

(2)

where θ is the semi-aperture angle of the transmitting and receiving beams, and k_eff is an effective aperture contributing to the beam focusing on the backside of the plate. The effective aperture accounts for the increased conversion of compressional waves into shear waves with increasing inclination angles and the absence of compressional wave transmission into the bone tissue for inclination angles larger than the critical angle θ_crit [19]:

$${\mathit{\theta }}_{\mathit{crit}} = {\sin }^{-1}\left(\frac{{\mathit{\nu }}_{H_{2}O}}{\mathit{Ct}{.\mathit{\nu }}_{11}}\right).$$

(3)

In contrast to our previous study [19], we have used an aperture size of 64 elements and adjusted the factor to estimate the effective aperture k_eff from 0.1 to 0.122:

$$k_{\mathit{eff}} = \left(\begin{array}{cc}1& \mathit{if} \theta < {\theta }_{crit} - 10°\\ 0.122·\mathsf{\Delta }\theta & \mathit{if} \theta > {\theta }_{crit} - 10°\end{array}\right).$$

(4)

More details on the estimation for k_eff can be found in Appendix A.

2.3.2. Phase Aberration Correction

Phase aberrations caused by bone curvature, bone tilt and material inhomogeneities are corrected for signals reflected from the backside cortical bone interface. Three phase-aberration correction (PAC) methods are used: (1) The curved bone surface geometry results in different round-trip travel times compared with the flat bone model. A time-shift correction based on the periosteal bone surface geometry (hereinafter called ‘surface time correction’ ST), was used to correct for the additional ultrasound wave propagation paths in the water due to the bone curvature. The ST correction used the concept of refraction compensation proposed by Yasuda et al. [24]. Further details are summarized in Figure A3a in Appendix C. (2) For tilted bone models, the reflected wavefront exhibits a tilt with respect to the beam axis (Figure 2c). To correct the phase aberration caused by surface inclination, an autocorrelation function (ACF) analysis on the reflected backside echoes was performed (Figure A3b–d in Appendix C). (3) Local variations of the sound velocity caused by material inhomogeneities lead to small fluctuations of the transit time measured at individual receiver elements and subsequently to a distortion of the summed signal; therefore, the following method was used to estimate the backside focus depth. The arrival times for all receiver elements was estimated using a cross-correlation (CC) method. The receiver channel that measured the highest signal amplitude was used as the reference signal. The inter-element arrival times exhibit either a concave, flat, or a convex shape, depending on the distance of the beam focus relative to the backside bone interface. A second-order polynomial was fitted to the inter-element arrival times, and the confocal focus depth was determined by finding the zero-crossing point of the second order fit coefficients (Figure A3e,f in Appendix C). This zero-crossing point was used to determine ΔF_z, and to estimate Ct.Th^MF and Ct.ν₁₁^MF using Equation (2).

2.4. Statistics

For each model, the retrieved Ct.Th^MF and Ct.ν₁₁^MF values were compared with the reference Ct.Th^Ref = 4 mm and Ct.ν₁₁^Ref. Simulation models without cortical pores had the reference speed of sound of Ct.ν₁₁^Ref = 3504 m/s. The bone models with cortical pores Ct.ν₁₁^Ref were extracted from the transmission measurements. Pearson linear regression analysis was performed to compare the parameters obtained using the multi-focus method with reference values. For all models with a 64-element aperture, the relative error (RE), precision, and accuracy values for each PAC method were determined and compared with the values without any PAC. Precision was defined as the coefficient of the variation of the difference between the predicted Ct.Th^MF, Ct.ν₁₁^MF and the reference values for Ct.Th^Ref, Ct.ν₁₁^Ref. Accuracy was determined by means of the root mean square error (RMSE) compared with the reference values. All analyses were performed using MATLAB R2019b, including the Signal Processing, Curve Fitting, and Statistics Toolboxes (The Mathworks, Natick, MA, USA).

3. Results

3.1. Numerical Simulations

A total of 22 bone models were simulated (

Table A2. Summary of simulation models. ‘r’ and ‘dx’ in the model abbreviations represent the curvature radius of the bone model and the lateral shift of the transmit and receive arrays relative to the beam axis, respectively. ‘Po’ represents the porosity value when pores were simulated.

Effect of.	Model Abbreviation	Curvature Radius r (mm)	Lateral Shift to Beam Axis dx (mm)	Bone Surface Tilt (°)	Porosity [%]
	flat plate	10,000	0	0	0
curvature	r60dx0	60	0	0	0
	r50dx0	50	0	0	0
	r40dx0	40	0	0	0
	r30dx0	30	0	0	0
	r20dx0	20	0	0	0
curvature	r40dx1.11	40	1.11	1.4	0
and	r40dx2.11	40	2.11	3.1	0
tilt	r40dx3.11	40	3.11	4.5	0
	r40dx4.11	40	4.11	5.9	0
	r40dx5.11	40	5.11	7.4	0
	r40dx6.11	40	6.11	8.9	0
curvature	r40dx0Po2	40	0	0	2
and	r40dx0Po4	40	0	0	4
porosity	r40dx0Po6	40	0	0	6
	r40dx0Po8	40	0	0	8
	r40dx0Po10	40	0	0	10
	r40dx0Po12	40	0	0	12
	r40dx0Po14	40	0	0	14
	r40dx0Po16	40	0	0	16
	r40dx0Po18	40	0	0	18
	r40dx0Po20	40	0	0	20

in Appendix D). The reference sound velocities Ct.ν₁₁^Ref of the porous models, as determined by transmission simulations, are summarized in

Table A3. Results of Ct.Th^Ref and Ct.ν₁₁^Ref of transmission simulation.

Model	Ct.ν11Ref [m/s]
r40dx0Po2	3428.6
r40dx0Po4	3321.8
r40dx0Po6	3189.4
r40dx0Po8	3127.0
r40dx0Po10	3038.0
r40dx0Po12	2953.8
r40dx0Po14	2848.7
r40dx0Po16	2774.6
r40dx0Po18	2704.2
r40dx0Po20	2681.6

in Appendix D. The estimated Ct.Th^MF and Ct.ν₁₁^MF values for all models and the relative errors are summarized in

Table A4. Results of Ct.Th^MF, Ct.ν₁₁^MF and relative errors (RE) for each PAC method.

Model	Correction	Ct.ThMF [mm]	RECt.Th [%]	Ct.ν11MF [m/s]	RECt.ν11 [%]
flat plate	No	4.00	0.01	3490	0.41
(reference)	ST	4.00	0.01	3490	0.41
	ST + ACF	4.00	0.01	3490	0.41
	ST + ACF + CC	4.00	0.01	3490	0.41
r60dx0	No	3.77	5.77	3330	5.83
	ST	3.99	0.32	3470	0.98
	ST + ACF	3.99	0.32	3470	0.98
	ST + ACF + CC	4.00	0.01	3490	0.41
r50dx0	No	3.72	6.91	3255	7.11
	ST	3.98	0.48	3475	0.83
	ST + ACF	3.98	0.48	3475	0.83
	ST + ACF + CC	4.00	0.01	3495	0.26
r40dx0	No	3.65	8.70	3180	9.25
	ST	3.96	1.07	3460	1.26
	ST + ACF	3.96	1.07	3460	1.26
	ST + ACF + CC	4.02	0.51	3510	0.16
r30dx0	No	3.57	10.80	3120	10.97
	ST	3.94	1.49	3440	1.83
	ST + ACF	3.93	1.65	3445	1.69
	ST + ACF + CC	4.06	1.48	3545	1.16
r20dx0	No	3.37	15.66	2955	15.67
	ST	3.82	4.49	3340	4.69
	ST + ACF	3.85	3.78	3360	4.12
	ST + ACF + CC	3.85	3.78	3360	4.12
r40dx1.11	No	3.67	8.26	3210	8.40
	ST	3.96	0.91	3455	1.41
	ST + ACF	3.96	1.07	3460	1.26
	ST + ACF + CC	4.03	0.67	3505	0.02
r40dx2.11	No	3.69	7.66	3235	7.68
	ST	3.93	1.65	3445	1.69
	ST + ACF	3.96	0.91	3455	1.41
	ST + ACF + CC	3.96	0.91	3455	1.41
r40dx3.11	No	3.67	7.50	3205	7.83
	ST	3.79	5.23	3325	5.12
	ST + ACF	3.95	1.19	3465	1.12
	ST + ACF + CC	3.95	1.19	3465	1.12
r40dx4.11	No	3.64	8.99	3175	9.40
	ST	3.58	10.45	3140	10.39
	ST + ACF	3.93	1.86	3445	1.69
	ST + ACF + CC	3.93	1.86	3445	1.69
r40dx5.11	No	3.61	9.84	3150	10.11
	ST	3.47	13.28	3040	13.25
	ST + ACF	3.80	4.91	3335	4.83
	ST + ACF + CC	3.80	4.91	3335	4.83
r40dx6.11	No	3.47	13.37	3030	13.53
	ST	3.29	17.74	2895	17.39
	ST + ACF	3.54	11.48	3110	11.25
	ST + ACF + CC	3.54	11.48	3110	11.25
r40dx0Po2	No	3.60	10.02	3065	10.60
	ST	3.88	3.08	3310	3.46
	ST + ACF	3.90	2.45	3330	2.88
	ST + ACF + CC	3.93	1.83	3350	2.29
r40dx0Po4	No	4.39	9.67	3640	9.58
	ST	3.82	4.57	3165	4.72
	ST + ACF	3.81	4.47	3170	4.57
	ST + ACF + CC	3.76	6.08	3135	5.62
r40dx0Po6	No	3.57	10.79	2830	11.27
	ST	3.85	3.79	3050	4.37
	ST + ACF	3.85	3.87	3055	4.21
	ST + ACF + CC	3.92	1.92	3120	2.18
r40dx0Po8	No	3.54	11.49	2805	10.30
	ST	3.79	5.27	3015	3.58
	ST + ACF	3.83	4.18	3010	3.74
	ST + ACF + CC	3.89	2.87	3055	2.30
r40dx0Po10	No	5.16	29.12	3950	30.02
	ST	3.98	0.48	3030	0.26
	ST + ACF	3.95	1.24	3015	0.76
	ST + ACF + CC	3.98	0.48	3030	0.26
r40dx0Po12	No	5.58	39.69	4090	38.47
	ST	3.86	3.46	2835	4.02
	ST + ACF	3.87	3.23	2830	4.19
	ST + ACF + CC	3.92	1.94	2870	2.84
r40dx0Po14	No	4.78	19.54	3440	20.76
	ST	3.78	5.40	2760	3.11
	ST + ACF	3.78	5.40	2760	3.11
	ST + ACF + CC	3.78	5.40	2760	3.11
r40dx0Po16	No	4.82	20.55	3870	39.48
	ST	3.79	5.36	2635	5.03
	ST + ACF	3.79	5.14	2630	5.21
	ST + ACF + CC	3.90	2.58	2695	2.87
r40dx0Po18	No	3.82	4.52	2585	4.41
	ST	3.98	0.56	2680	0.89
	ST + ACF	3.97	0.79	2685	0.71
	ST + ACF + CC	3.83	4.21	2640	2.37
r40dx0Po20	No	5.65	41.24	3750	39.84
	ST	2.22	44.45	2675	0.25
	ST + ACF	3.33	16.66	2730	1.80
	ST + ACF + CC	3.03	24.14	2670	0.43

in Appendix D. Without PAC, all deviations from the ideal flat plate geometry led to deteriorations of precision and accuracy. In most situations, PAC improved both the precision and accuracy (

Table 3. Precision of Ct.Th^MF and Ct.ν₁₁^MF after each PAC method.

	Model	No PAC	ST	ST + ACF	ST + ACF + CC
Ct.ThMF	Curved bone plate	4.3%	1.7%	1.4%	2.0%
	Curved tilt bone plate	2.3%	7.3%	4.3% (4.2%) *	4.3% (1.1%) *
	Material inhomogeneity	18.5%	1.4%	4.7%	7.2% (1.9%) **
Ct.ν11MF	Curved bone plate	4.3%	1.6%	1.4%	2.0%
	Curved tilt bone plate	2.3%	7.1%	4.2% (2.5%) *	4.3% (0.8%) *
	Material inhomogeneity	15.9%	7.9%	7.9%	8.2% (7.9%) **

* Exclusion of bone models with tilt angles over 7°. ** Exclusion of bone model with porosity 20%.

and

Table 4. Accuracy of Ct.Th^MF and Ct.ν₁₁^MF after each PAC method.

	Model	No PAC	ST	ST + ACF	ST + ACF + CC
Ct.ThMF	Curved bone plate	10.2%	2.2%	1.9%	1.8%
	Curved tilt bone plate	9.6%	10.3%	5.2% (1.3%) *	5.2% (1.2%) *
	Material inhomogeneity	23.2%	14.6%	6.3%	8.3% (3.5%) **
Ct.ν11MF	Curved bone plate	10.4%	2.4%	2.1%	1.9%
	Curved tilt bone plate	9.8%	10.1%	5.1% (1.4%) *	5.1% (1.2%) *
	Material inhomogeneity	25.3%	3.4%	3.4%	2.8% (3.0%) **

* Exclusion of bone models with tilt angles over 7°. ** Exclusion of bone model with porosity 20%.

), which will be described in more detail in the following sections.

3.1.1. Effect of Bone Curvature

All three PAC methods showed improvements of Ct.Th^MF and Ct.ν₁₁^MF estimations with respect to precision and accuracy (Table 3 and Table 4). Although the ST correction alone showed the strongest improvement, the combination of all three PAC only yielded slight further improvements.

3.1.2. Effect of Bone Tilt

To correct for the bone tilt, using ST correction was not sufficient, and it even degraded accuracy and precision values (Table 3 and Table 4). The wavefront inclination caused by the tilted surface was effectively corrected using the ACF; however, for bone models with tilt angles above 7°, the CC correction method failed, because no zero-crossing point for the estimation of confocal focus depth could be retrieved (Figure A4c in Appendix D). After excluding these models, precision values for Ct.Th^MF and Ct.ν₁₁^MF were 1.1% and 0.8%, respectively, and accuracy values were 1.2% for both the Ct.Th^MF and Ct.ν₁₁^MF estimations.

3.1.3. Effect of Material Inhomogeneities

The presence of pores strongly degraded accuracy and precision values without PAC. The ST correction strongly improved precision and accuracy. Additional ACF correction had no effect in the evaluated simulations, because all porous bone models were modeled without a tilt. The CC further improved precision and accuracy values for Ct.Th^MF and Ct.ν₁₁^MF. For the bone model with the highest porosity value of 20%, all PAC methods did not result in a precise and accurate estimation of Ct.Th^MF and Ct.ν₁₁^MF (Table A4 and Figure A5 in Appendix D).

3.1.4. Overall Effect of PAC

Table 5. Precision and accuracy of Ct.Th^MF and Ct.ν₁₁^MF for each PAC method for all 64-element aperture models. Reference thickness for all models is Ct.Th^Ref = 4 mm. Models without cortical pores have Ct.ν₁₁^Ref = 3504 m/s. Reference values Ct.ν₁₁^Ref for models including cortical pores can be found in Table A3 in Appendix D.

	Correction	Precision	Accuracy
Ct.ThMF	No	17.4%	17.1%
	ST	10.3%	11.2%
	ST + ACF	4.1%	5.2%
	ST + ACF + CC	5.7% (2.1%) *	6.3% (2.6%) *
Ct.ν11MF	No	11.6%	18.5%
	ST	9.5%	5.9%
	ST + ACF	9.5%	3.7%
	ST + ACF + CC	9.8% (9.6%) *	3.4% (2.3%) *

* Exclusion of bone models with tilt angles over 7° and/or porosity above 20%.

shows the precision and accuracy values for all 22 simulation models. Note that precision values are defined as the coefficient of variation of the difference between the estimated and reference value of the flat bone plate model and accuracy is defined as RMSE as a percentage. That means the smaller the precision and accuracy value, the more precise and accurate the parameter estimation is with respect to the reference value. Overall, the combination of the three PAC methods results in an improved precision and accuracy estimation Ct.Th^MF and Ct.ν₁₁^MF. Precision values for Ct.Th^MF and Ct.ν₁₁^MF were 5.7% and 9.8%, respectively. Accuracy values for Ct.Th^MF and Ct.ν₁₁^MF were 6.3% and 3.4%, respectively.

3.2. Ex-Vivo Multi-Focus Measurement

Reference cortical thickness and porosity values of the human tibia bone at the central anteromedial part were found to be Ct.Th^Ref = (2.65 ± 0.61) mm and 15.1%, respectively, using HR-pQCT. For the MF measurement, Ct.Th^MF and Ct.ν₁₁^MF were determined in a manually selected region of interest (Figure 3a green ROI). The maximum amplitudes over all focus depths with and without PAC methods (ST + ACF + CC) are shown in Figure 3a,b. The endosteal surface of the human tibia sample is more blurred in the maximum projection image without PAC methods (Figure 3a). Without PAC, Ct.Th^MF and Ct.ν₁₁^MF could be retrieved at 14 lateral scan positions, whereas with PAC, cortical parameter estimations were achieved at 29 scan positions. The mean and standard deviation of Ct.Th^MF without PAC was (2.39 ± 0.25) mm, which was significantly different from the reference value. In contrast, the estimation of Ct.Th^MF with PAC of (2.71 ± 0.22) mm was not significantly different from the reference value. The estimated cortical speed without and with PAC were (2870 ± 95) m/s and (2857 ± 52) m/s, respectively.

4. Discussion

In this study, we have extended the estimations of thickness and speed of sound in cortical bone in a transverse plane using the multi-focus approach to realistic bone geometries. For this, several phase aberration corrections were proposed. The effects of bone curvature, surface inclination relative to the beam axis, and the presence of intracortical pores’ parameter estimations were analyzed.

4.1. Numerical Simulation

4.1.1. Effect of Bone Curvature

For curved bone models positioned parallel to the probe array (without bone tilt), the ST correction was sufficient and corrected the additional geometrical time shifts for curved bone interfaces compared to a flat bone plate. For the correction, it was assumed that ultrasound waves propagate in a straight direction, as described in ray theory [34].

4.1.2. Effect of Bone Tilt

Additional phase aberration corrections on the backside echoes were necessary for the curved models with bone tilt relative to the beam axis, to correct the orientation dependence of the reflected wavefront. Here, autocorrelation function was used on the backside echoes to estimate the inclination of the backside echoes; therefore, an ellipsoid was fitted on the magnitude of the backside signals after autocorrelation analysis.

4.1.3. Effect of Material Inhomogeneities

Cortical pores result in scattering and subsequent diffusion of the ultrasound waves. This causes local fluctuations of the arrival time of the received backside echoes compared with the reference flat bone plate model. The backside confocal depth, which is required for the simultaneous estimation of both thickness and sound velocity, has been estimated in our previous work by detecting the peak position of the DAS beamformed backside echoes with respect to the focusing depth [19]. Phase aberration induced by cortical pores cause a decrease in the intensity of the beamformed signal. With increasing porosity, the confocal peak arising from the backside reflection becomes less sharp and the peak position is harder to detect; therefore, we have developed another method to extract the confocal backside position by analyzing the curvature of the backside echo wavefront prior to the summation of all receive channels at each focus depth. The curvature of the wavefront was extracted by analyzing the cross correlation of the backside echoes relative to the backside echo with the highest signal amplitude. The change of the wavefront curvature from a convex shape (negative curvature) to a concave shape (positive curvature) was used to extract the focus position. The zero-crossing point was calculated using a linear fit of the retrieved curvature values over the focus depth. Incorporating cross correlation analysis prior to the summation of the beamformed signals improved the accuracy of the estimation of the backside confocal position, as well as precision and accuracy in simulations including pores (Table 5). Moreover, this method improved the backside signal detection rate and the accuracy of the estimation of cortical thickness in the ex-vivo measurement.

4.1.4. Combination of Phase Aberration Methods

For the transition to in-vivo applications of the multi-focus method, the combination of all three PAC methods is necessary, because all the investigated deviations from an ideal flat homogenous plate are present in real cortical bone. Overall, the combination of the three PAC showed a strong improvement of precision and accuracy values for cortical thickness and speed of sound estimations than when compared to the values without PAC.

4.2. Ex Vivo Measurement

The endosteal surface of the human tibia sample was tracked with and without PAC methods; however, more endosteal surface locations were retrieved when PAC was used. The cortical thickness measured by ultrasound was consistent with the reference value measured by HR-pQCT. The cortical sound velocity of (2857 ± 53) m/s was in the range of the cortical speed of sound values typically found in human cortical bone [17]. Our previous study showed a dependency of cortical speed of sound on cortical porosity (Po) Ct.ν₁₁^fit = 0.39·Po² − 51.4·Po + 3485 [m/s], Figure 6a in [19]). By inserting the reference cortical porosity value of 15.1% obtained from HRpQCT into this equation a speed of sound value of Ct.ν₁₁^fit = 2804 m/s was determined for the human tibia sample. The MF-based estimation was in the range of the expected speed of sound value; however, this observation needs to be confirmed in a larger sample size. In conclusion, the ex vivo measurement on a human tibia sample suggests the ability to measure cortical thickness and speed of sound using the MF approach by incorporating PAC methods.

4.3. Transition to In Vivo Application

Cortical bone has been proposed as significant predictor of a bone’s mechanical strength because mechanical force given to a bone is carried primally by cortical bone [35]. Clinical studies showed an improvement of fracture prediction by measuring cortical thickness [36,37,38]. HR-pQCT is the most precise modality to measure cortical thickness at the tibia with a precision of 1.6% [38]. Our study showed a thickness precision estimation of 5.7%. We expect that the clinical precision of the MF approach could be larger than for controlled simulations; however, HR-pQCT uses ionizing radiation and is extremely expensive compared with ultrasound imaging. Wydra et al. [39] proposed a similar refraction measurement method and reported precision values for Ct.Th of 8.5% for measurements on porous plate-shaped skull bone phantoms. In contrast, our study considered bone curvature and bone tilt with a better precision value of 5.67%, which can be attributed to the PAC methods, the use of a higher frequency (4 vs. 2.25 MHz), and the consideration of an effective aperture [19].

In addition to cortical thickness, ultrasonic wave-speed in cortical bone has been proposed as a biomarker for bone quality [10,40,41,42]. Bidirectional axial transmission techniques use a probe with several ultrasonic transmitters and receivers to measure waves traveling in the longitudinal direction of long bones. An in vivo study by Minonzio et al. used a bidirectional axial transmission technique (BDAT) to estimate the cortical thickness and porosity, and they reported those parameters as suitable biomarkers for fracture discrimination in postmenopausal women [43]. The QUS device Bindex^® calculates the apparent cortical thickness at the distal radius and tibia using BDAT and reported the correlation with BMD (r ≥ 0.71, p < 0.001, 0.20 < R² < 0.55) [44]. Talmant et al. [41] showed that the velocity of the first arriving signal (νFAS) is a significant biomarker for fracture discrimination and to predict fracture risk in vivo. Inter-operator precision (repeated measurements by different operators) for FAS velocities were reported at ~7%, respectively. In our study we report the precision value for different simulation models (precision of radial cortical speed of sound was 9.8%), which have been simulated only once. Compared with axial transmission techniques, the multi-focus measurement estimates cortical thickness and speed of sound within the imagined plane and provides image guidance.

Another approach to measure Ct.Th and Ct.ν₁₁ using corrected refraction was proposed by Renaud et al. [14] using a single-element excitation, full-array waveform capture, and an adapted Kirchhoff migration developed by seismologists to image the earth subsurface. The method was validated in vivo on two young healthy subjects. No precision or accuracy values were reported. In two separate studies, Karjalainen et al. [11,45] proposed the estimation of an apparent Ct.Th from TOF between periosteal and endosteal bone interface at the tibia using a constant predefined speed of sound in cortical bone of 3565 m/s. This approach fails to capture the microstructural changes in porous bone structures and changes in Ct.ν₁₁. In contrast, our method estimates Ct.Th and Ct.ν₁₁ independently; however, in this study, a very simple pore structure was assumed. Further studies should therefore target bone models with more realistic pore diameter distributions.

Recently, Iori et al. proposed a cortical backscatter model to retrieve the intracortical pore size distribution non-invasively in the tibia midshaft [46]. These findings were further supported by another study on the same set of bones, which suggested that cortical thinning and backscatter parameters describing the presence and accumulation of large cortical pores in the tibia provide similar or better predictions of proximal femur stiffness and ultimate force than aBMD [20]. The cortical backscatter (CortBS) method has been applied for the first time in vivo by Armbrecht et al. [15] on postmenopausal women with low bone mineral density. The study reported a better discrimination performance for vertebral and non-vertebral fragility factures using cortical backscatter parameters (0.69 ≤ AUC ≤ 0.73) compared with DXA based aBMD (0.54 ≤ AUC ≤ 0.55). As the CortBS and multi-focus measurement modalities can be implemented in the same device, future in-vivo studies should be performed to evaluate if such a multiparametric assessment of macro- and microstructural (i.e., Ct.Th and intracortical pore size, respectively) and viscoelastic (i.e., Ct.ν₁₁ and attenuation coefficient α(f)) cortical bone properties can improve the discrimination and risk prediction performance for distinct types of fragility fractures. The combined estimation of Ct.Th, Ct.ν₁₁ and pore size distribution using nonionizing and noninvasive technique may have a high clinical potential to prevent osteoporotic fractures.

4.4. Limitations

Several limitations of the proposed PAC methods were observed in this study. The methods fail for bone inclination angles larger than 7° with respect to the beam axis, as well as for the bones with high porosity values (20% or more). For bone models with tilt angles larger than 7°, most backside echoes were not captured by the receiver array resulting in a much smaller DAS beamformed signal and the transition from convex to concave shape of the backside signal wavefront disappeared. Subsequently, the zero-crossing point could not be retrieved (Figure A4c in Appendix D). As the bone surface inclination in the imaging plane can be reliably reconstructed, the application of the PAC methods can be easily restricted to locations, in which the surface inclination is within ±7°. Second, the simulation study was restricted to one scan position for one multi-focus measurement, while the ex vivo measurement performed the multi-focus measurement at 128 scan positions along the lateral distance; therefore, future in silico studies should simulate multi-focus measurements with more scan positions along the lateral distance and include simulation models with real bone curvature, tilt, and porosity. Moreover, compound imaging with beam steering [15,46] should be used to ensure that the bone area of interest is probed with sufficiently small beam inclinations. Third, for high porosity values, large amounts of scattering of ultrasound waves resulted in a strong attenuation and distortion of the backside signal, yielding an imprecise estimation of the confocal backside position (Figure A5c in Appendix D).

Another limitation is the use of simplified bone models. For in vivo transition, the effect of heterogeneous cortical pores and heterogeneous backside surface on the phase aberration should be investigated. Cortical pores lead to increased scattering, and therefore, increased phase aberration, which could be corrected with cross correlation analysis. Furthermore, the effect of changes regarding the speed of sound in soft tissue should be considered in the future, based on realistic simulation models. Conventional image reconstruction assumes an invariant speed of sound of 1540 m/s. Although the higher and variable velocity in bone was considered, soft tissue velocities can also vary by up to 10% between subjects depending on the relative distribution of skin, fat, and connective tissue along the bone length [16]. This leads to additional wave distortion, defocusing of bone regions, and misalignments of beamformed signals. Anderson et al. [47] showed on a tissue-mimicking phantom that a speed of sound error up to ±8% degrades the lateral resolution of the image by up to a factor of three. The mismatch between the assumed and actual speed of sound could be compensated for by evaluating the focus quality using the coherence factor proposed by Hasegawa et al. [48] or by using the minimum average sum of absolute differences between all pre-beamformed radio frequency channel data proposed by [49]. Renaud et al. proposed an autofocused method to estimate the optimal speed of sound of the overlaying soft tissue [50]. Additional phase aberration corrections by tissue structure may improve lateral resolution, signal quality and the accuracy and precision of the measured time-of-flight through the bone for in vivo transition measurements.

Another limitation of this study is the non-repeated measurement design. Only one simulation was performed for each simulation model and the ex vivo measurement was performed once; therefore, for the reproducibility and precision of the multi-focus method for realistic bone simulations, ex vivo and in vivo measurements should be investigated in the future, by repeating the measurements by repositioning of the transducer between each measurement.

5. Conclusions

This study demonstrates the assessment of cortical thickness and speed of sound in the radial direction using refraction- and phase-aberration corrected MF imaging. Conventional DAS beamforming was improved using phase aberration correction methods to account for bone curvature, bone tilt, and bone material homogeneities from cortical pores. The method was developed and validated using in silico simplified bone models with and without pores, and one ex vivo measurement was performed on a human tibia cadaver. For a reliable in vivo estimation of cortical thickness and speed of sound values, the real bone structures and soft tissue velocity inhomogeneity must be considered. The derived parameters showed an improvement in precision and accuracy using phase aberration corrections and demonstrated good agreement with reference values.

6. Patent

K.R. has the patent “CortBS: Ultrasonic method for determining pore dimensions in cortical bone” pending.