The effects of image homogenisation on simulated transcranial ultrasound propagation

Transcranial transmission of ultrasound is increasingly used in a variety of clinical and research applications, including high intensity ablation, opening the blood brain barrier, and neural stimulation. Numerical simulations of ultrasound propagation in the head are used to enable effective transcranial focusing and predict intracranial fields. Such simulations require maps of the acoustic properties of the head, which can be derived from clinical CT images. However, the spatial resolution of these images is typically coarser than the scale of heterogeneities within the skull bone, which are known to exert a major influence on ultrasound propagation. In the present work, the impact of image related homogenisation on transcranial transmission from a single element transducer is examined using a dataset of co-registered clinical resolution CT and micro-CT images of skull sections. Reference acoustic property maps are derived from micro-CT images of cortical bone tissue. The influence of imaging resolution is examined by progressively downsampling the segmented acoustic property maps, and through comparison with maps derived from co-registered clinical CT images. The influence of different methods of segmenting the bone volume from the clinical CT images, and for resampling the clinical and micro-CT data are also examined. Image related homogenisation is demonstrated to have a substantial effect on the transcranial transmission of ultrasound, resulting in underestimations of simulated transmission loss and time-of flight. Effects on time-of flight are due to the loss of the internal scattering microstructure of the skull, while changes in transmitted ultrasound amplitude are due to both loss of microstructure and other smoothing effects. Inflating the simulated attenuation coefficient of the skull layer reduces the error in transmitted pressure amplitude to around 40%, however this is unable to correct fully for errors in time of flight and the pressure distribution of the transmitted field.

Large scale numerical simulations of transcranial ultrasound propagation are an important tool in achieving effective transcranial ultrasonic focusing. Time-reversal simulations, which involve modelling the propagation of ultrasound from the location of the desired focus inside a model of the head, through the skull to a virtual transducer surface, can be used to derive drive signals for large multi-element arrays (Aubry et al 2003). Forward propagation simulations modelling the propagation of ultrasound into the skull from a known ultrasound source are also of use for predicting intracranial fields and anticipating treatment impact , Lee et al 2015, Top et al 2016.
In order to simulate the propagation of ultrasound within the head, accurate maps of the acoustic properties of the various tissues, which can be derived from medical images, are necessary (Aubry et al 2003). X-ray CT images are typically used due to their excellent bone contrast, and the potential for accurate measurement of density from CT Hounsfield units (HU) (Kyriakou et al 2014). Heterogeneous models derived from clinical CT images have previously been used for successful transcranial time-reversal focusing in a number of studies (Marquet et al 2009, Chauvet et al 2013, using a range of different methods for deriving acoustic property maps (Aubry et al 2003, Pichardo et al 2011. However, the use of these images to produce acoustic property maps is subject to some limitations. Notably, the spatial resolution of clinical CT images is limited to an approximate scale of 0.5 mm, which leads to a loss of the fine internal microstructure of the skull and smoothing of interfaces between bone and soft-tissue due to partial volume effects (Pinton et al 2011). The trabecular portion of the skull was identified by Fry and Barger (1978) as a major influence on transcranial propagation, while Pinton et al (2011) confirmed that scattering from the trabecular microstructure is the dominant source of attenuation in the skull. In their method for deriving acoustic properties from clinical CT images, Aubry et al (2003) attempt to account for this by artificially increasing the acoustic absorption value assigned to the skull bone to account for the scattering that occurs in the trabecular layer, based on the attenuation measurements made by Fry and Barger (1978). However, this can only account for losses due to scattering, and not changes in the shape of propagating ultrasonic fields.
In the present paper, simulations of ultrasound propagation through the skull are performed to examine the impact of smoothing and loss of microstructure detail due to the limits of medical image resolution. Acoustic property maps are derived from a dataset comprising CT images of skull bone obtained at clinical resolutions (herein referred to as clinical CT images or clinical CT data) and corresponding micro-CT data of sections of skull bone co-registered to the original clinical CT images. These simulations are used to examine the impact of progressive downsampling of high resolution acoustic property maps based on micro-CT data, which mimics the effects of image related homogenisation. Simulations through maps based on micro-CT data are compared with simulations through heterogeneous acoustic property maps derived from co-registered clinical CT data. Different methods for segmenting the CT data and resampling it onto simulation grids at the resolution of the micro-CT images are compared. The degree to which assigning whole skull attenuation to the skull layer en-bloc is able to correct for the loss of microstructure in clinical CT data is also examined. Simulations through homogeneous medium maps derived from the same clinical CT data are carried out to examine the difference from heterogeneous maps relative to the simulations through micro-CT data. The impact of elastic wave propagation through the microstructure of skull bone is examined using simulations in 2D.

Clinical-CT and micro-CT datasets
Two human cadavers were obtained from the Wake Forest School of Medicine Whole Body Donation Program (Lillie et al 2015). Clinical CT scans of each cadaveric head were obtained with a slice thickness of 0.625 mm and in-plane resolution of 0.488-0.625 mm. Following acquisition of the clinical CT images, small samples of bone (four from skull A, three from skull B) were extracted for high-resolution imaging in a micro-CT scanner and reconstructed at an isotropic resolution of 50 µm. Micro-CT data was then manually aligned and co-registered with the clinical CT datasets using GEOMAGIC STUDIO (v 12.1.0 Geomagic, Research Triangle Park, NC, USA), described in Lillie et al (2015). This alignment was used to create an affine transform matrix for each micro-CT image, that allowed them to be aligned and registered with the corresponding clinical CT datasets (Lillie et al 2015).
Following visual inspection, two of the datasets were rejected from the study due to their geometry being unsuitable for the planned transmission simulations. This left a total of two micro-CT images for skull A, and three for skull B. These comprised a section of frontal and parietal bone for each skull, in addition to a portion of temporal bone from skull B. 3D renderings of the individual bone samples are shown in figure 1, alongside sections through the raw micro-CT data, demonstrating the three-layered structure and internal heterogeneity of the skull. Here, it can be seen that the section of temporal bone has a different structure than the samples of frontal and parietal bone, comprising mainly dense cortical bone with very little trabecular microstructure. The temporal bone sample serves as a comparison to the other bone fragments, and is used to determine the extent to which any effects observed are due to the loss of bone microstructure, rather than other aspects of the smoothing and homogenisation process.
Micro-CT and clinical CT datasets were aligned inside the Matlab environment. Due to the need to interpolate the clinical CT data onto spatial grids of the same spatial sampling as the micro-CT data for simulations, the affine transform matrix was inverted and used to rotate the clinical CT data. Interpolation of the clinical CT onto the micro-CT resolution grids was carried out using both nearest neighbour and linear interpolation, to determine which method results in more accurate transcranial simulations.

Assignation of acoustic medium properties
Acoustic property maps were derived from the micro-CT datasets by segmenting the 3D image volume into cortical bone and soft-tissue components based on the CT intensity. Segmentations of all micro-CT datasets were carried out using the Seg3D segmentation toolbox (CIBC 2016) and exported directly into the Matlab environment. A 40% threshold intensity was used to create the segmentation, followed by the application of Seg3D's connected component algorithm to remove spurious components external to the bone volume. A separate, solid segmentation delineating the bone segment as a whole was also created to allow soft-tissue within the bone layer (bone marrow) to be separated from the rest of the volume. The acoustic properties of cortical bone were assigned to the initial segmentation of the bone, soft-tissue within the bone volume was assigned the acoustic properties of bone marrow and the rest of the domain was assigned the properties of water. To ensure Figure 1. 3D renderings of the micro-CT skull sections produced in the Seg3D segmentation toolbox, alongside sections through the raw micro-CT data, demonstrating the internal structure of the skull bone. that acoustic property maps derived from micro-CT data could be compared to those derived from clinical CT data using the algorithm described below, the absorption value assigned to cortical bone was based on the measurements by Pinton et al (2011), and the same sound speed and density were assigned to water and bone marrow. The acoustic property values used are shown in table 1.
Acoustic properties were assigned to clinical CT data using the algorithm developed by Aubry et al (2003), which uses HU values to interpolate between acoustic sound speed and density values for cortical bone and softtissue. In accordance with previous usage of this algorithm by Deffieux and Konofagou (2010), the first step is segmentation of clinical CT data into skull and non-skull segments using a threshold segmentation of 50% of the peak intensity in the image, referred to herein as the 'tight' segmentation. Segmentations of the clinical CT data were carried out in the Matlab environment separately for linear and nearest neighbour interpolations following rotation and interpolation of the clinical CT images to align them with the grids of micro-CT data. To determine the influence of accounting for smoothing around the edges of acoustic property maps, a separate 'loose' segmentation was created by dilating the tight segmentation by 1 mm in all directions. Within the loose segmentation volume, intensity values of the clinical CT data were then rescaled to vary between 0 and 1, and used to interpolate between the sound speed and density values used for cortical bone and marrow/water in the micro-CT data. This ensures that the basis for assigning acoustic property values remained consistent between clinical and micro-CT datasets. The tight segmentation was then applied to the resulting acoustic property maps, to determine the impact of enforcing a sharp medium boundary at the edges of the skull. The range of different acoustic property maps produced from both clinical and micro-CT data are shown in figure 2.
Simulations through acoustic property maps derived from the clinical CT data were carried out both with and without attenuation, modelled within the simulations as absorption, to determine how effectively it accounted for the reduced scattering due to loss of microstructure detail at clinical CT resolutions. As in Deffieux and Konofagou (2010), attenuation was assigned homogeneously to the bone volume using the corresponding skull segmentation. The whole skull attenuation value assigned was that measured by Pinton et al (2011), who also provide the absorption value for cortical bone.
In addition, a series of homogeneous acoustic property maps were created based on the average values of sound speed and density across the heterogeneous maps derived from clinical CT data. The homogeneous properties were derived separately for each skull segment, interpolation methods, and segmentation boundary. The average value and standard deviation of the sound speed and density of the bone layer across all the homogeneous simulations were 2490 ± 167 ms −1 and 1725 ± 127 kgm −3 , respectively.
For all simulations where absorption was modelled (including attenuation modelled as absorption), the power law exponent was set to y = 2. As all simulations used either CW sources, or tonebursts with a narrow frequency spectrum, this does not influence the absorption of ultrasound. However, when absorption is included, it is accompanied by physical dispersion (required for the wave to obey causality), which causes a frequency dependent change in sound speed determined by the power law exponent, shown in figure 3 Cox 2010, 2014). Using y = 2 results in no dispersion, which allows direct comparison between lossless and absorbing simulations through clinical CT data.

Homogenisation of property maps
The maps of acoustic properties derived from micro-CT data were progressively homogenised in order to mimic the effects of imaging resolution. Homogenisation was applied directly to the maps of acoustic properties rather than the original micro-CT images, and consisted of spatially downsampling the maps by averaging medium properties across neighbouring groups of grid points in order to imitate the voxels of images obtained at a lower spatial resolution. These downsampled medium maps were then resampled back onto the simulation grid at the micro-CT resolution using either linear interpolation or nearest neighbour interpolation, resulting in both smooth and pixelated property maps. The initial acoustic property maps, with a spatial resolution of 50 µm, were homogenised to effective voxel sizes of up to 1 mm in 50 µm increments. Figure 4 shows examples of the resulting medium maps at different effective resolutions alongside those derived from clinical CT. Table 1. Acoustic property values used in the derivation of property maps from CT images. Sound speed and density values for softtissue and cortical bone are taken from the values for water used in Deffieux and Konofagou (2010). Absorption values for cortical bone, and attenuation values for whole bone are taken from Pinton et al (2011). The absorption value for bone marrow is taken from the IT'IS database (Hasgall et al 2015).

Overview
Simulations of ultrasound propagation were carried out to evaluate the impact of image related homogenisation of acoustic property maps derived from micro-CT data. Numerical simulation was carried out using the k-Wave toolbox. k-Wave uses a k-space corrected pseudospectral time domain (PSTD) to solve the coupled first-order conservation equations for the propagation of ultrasound in heterogeneous media, including power law absorption, and also includes a PSTD scheme for modelling elastic wave propagation (Treeby et al 2012.
The results for micro-CT data homogenised at different spatial sampling levels were compared with simulations through both heterogeneous and homogeneous maps derived from clinical CT data. Simulations through the highest resolution, non-homogenised micro-CT data, where acoustic properties were assigned based on cortical bone segmentations, were considered reference simulations and errors were computed relative to them. Previous numerical studies using tissue properties derived from micro-CT images have shown good agreement with experimental measurements made with the same samples (Luo et al 1999, Bossy et al 2005, Nagatani et al 2006, Bossy et al 2007, Moilanen et al 2007, Machado et al 2011. Simulation through maps based on clinical CT data when attenuation was not simulated were compared with simulations through non-absorbing micro-CT data, while simulations through attenuating clinical CT maps were compared with simulations through absorbing micro-CT data.

Figure 2.
A comparison of the acoustic property maps generated from non-homogenised micro-CT dataset, and the clinical CT data. The clinical CT data has been resampled to the same resolution as the micro-CT data using both linear and nearest neighbour interpolation. The maps based on clinical CT data are also constrained using both 'tight' and 'loose' volume segmentations, which were also used to define the homogeneous attenuation maps for the clinical CT data. In keeping with the main uses of numerical simulation of transcranial ultrasound propagation, two testing schemes were used. The first was based on simulations designed to examine the intracranial ultrasound field produced by applied ultrasound, designated the 'forward' problem. The second was based on time-reversal simulations used to derive drive signals for transducer elements of large transcranial arrays, designated the 'reverse' problem. The layout for both forward and reversal simulations is shown in figure 5. This demonstrates the position of the transducers relative to the bone and the source and sensor distributions placed on the internal side of the bone map.
The same simulated transducer element was used in both forward and reversal simulations (in transmit and receive modes, respectively) based on individual elements of large hemispherical transcranial arrays like that presented by Clement et al (2000). The element was modelled as a single element focused bowl transducer with a diameter of 8 mm and a radius of curvature of 10 cm, operating at a frequency of 1 MHz, which lies at the upper end of ultrasound frequencies feasible for transcranial transmission (Deffieux and Konofagou 2010). Forward simulations were carried out using ten-cycle toneburst (TB) sources, while reversal simulations used both continuous wave (CW) and ten-cycle TB ultrasound sources.
2D elastic simulations were carried out in the MATLAB environment with CUDA hardware acceleration (Treeby et al 2014) on a Dell PowerEdge R730 compute server with 2 × 6-core Xeon E5-2620 2.4 GHz CPUs, 64 GB of 1866 MHz memory, and an Nvidia Titan X GPU with 3072 CUDA cores and 12 GB of memory. 3D fluid simulations were run using the MPI version of k-Wave on the IT4I Salomon supercomputing cluster (Jaros et al 2016). Each 3D simulation was carried out on a compute node with an 8-core Intel Xeon E5-4627v2 3.3 GHz CPU, and 256 GB of RAM. The grid size and simulation time for each bone sample is shown in table 2. All simulations used a grid spacing of 50 μm giving approximately 30 points per wavelength for 1 MHz in water (Robertson et al 2017a). The timestep was set to 0.65 ns, giving 1538.5 PPP and a CFL number of 0.04.

Forward propagation simulation
The forward simulation modelled the propagation of ultrasound into the skull from the single element transducer. The use of a 50 µm spatial discretisation step prevented extending the simulation domain to include the theoretical focus of the transducer, as the grid size would become too large to execute the simulations. Therefore the simulated field was evaluated by recording the time-varying pressure field over a 2D transverse plane with a size of 11 × 11 mm, located 2 cm from the centre of the virtual transducer (see figure 5(a)).
For each grid point across the recording plane, the TB data signal envelope was computed and used to calculate the amplitude and time-of-flight of the TB peak across the plane. An L2 error norm in the temporal peak amplitude at each grid point, as well as the average error in time-of-flight, were computed across the recording field for each simulation relative to the corresponding reference simulation.

Time-reversal simulation
For the reversal simulation, the source term was designed to imitate an outwardly propagating spherical wave originating from the theoretical focus of the transducer element. The source was modelled as a curved surface spanning the entire simulated domain, with a radius of curvature of 8 cm. The center of the curve was placed at 2 cm from the center of the transducer element on the internal side of the bone map, as shown in figure 5(b). After the simulation was completed, the time-varying pressure recorded across the virtual transducer was integrated into a single signal. In practical focusing simulations, this signal would normally be time-reversed, and used to derive drive signals for a transcranial sonication. This could involve using the raw time-reversed signal as an input, or analysing it spectrally to calculate phase and amplitude drive parameters. The integrated signal was analysed accordingly. For CW simulations, it was decomposed spectrally to give phase and amplitude values, while for TB data the signal envelope was computed, and the errors in the amplitude and time-of-flight of the envelope peak were calculated. An L2 error was also calculated for the entire TB signal, which is more representative of the error that would be introduced in full time-reversal of the recorded signal. It should be noted that, to allow analysis of phase values for the CW simulations, the results for the progressive homogenisation were unwrapped to determine the trend in phase. For plotting, phase values for the individual simulations based on clinical CT data were shifted by increments of 2π to move them onto the same plotting axis.

Elastic simulations
Additional testing was carried out in 2D to examine the impact of shear mode propagation on transmission through acoustic property maps derived from the micro-CT data. Previous simulation studies through cortical bone have demonstrated that shear waves do not play a significant role when the ultrasound waves are close to normal incidence (Robertson et al 2017b).The reversal simulation setup was converted into 2D for testing, and simulations were carried out using TB sources only.
Simulations were carried out using the elastic PSTD code in the k-Wave toolbox , and simulations were conducted both with and without elastic properties assigned to bone layer to determine the impact of shear mode propagation and the homogenisation of elastic properties. When simulated, a shear sound speed  of 1500 m s −1 was assigned to the cortical bone segment based on the measurements made by White et al (2006), and an absorption value of 19.5 dB cm −1 for ultrasound absorption at 1 MHz, based on the power law derived by . Shear sound speed and absorption were set to zero throughout the rest of the domain. Homogenisation of the acoustic property maps was carried out using linear interpolation, as shown in figure 4. Elastic simulations were not carried out using the clinical CT data, as a reliable method for converting clinical CT images of the skull into maps of elastic acoustic properties is not available.

Forward propagation
Progressive homogenisation of acoustic property maps derived from micro-CT data leads to an overestimation in the amplitude of transmitted ultrasound, an underestimation of time-of-flight, and a loss of fine spatial detail in the intracranial field. Figure 6 summarises the error metrics computed for the forward simulations, averaging the results across frontal and parietal bone samples. These are compared with the results of simulations through temporal bone, which demonstrated reduced errors relative to the bone samples with more internal microstructure, especially in terms of the effect on the time-of-flight. The errors shown are calculated for nonabsorbing simulations only, as absorbing data showed identical trends only with a slight reduction in error. When downsampling acoustic property maps derived from micro-CT data, the use of a nearest neighbour algorithm for resampling onto the high-resolution grid consistently results in reduced error in amplitude, including for simulations through temporal bone. The use of a nearest neighbour algorithm for interpolation of clinical CT data also results in reduced error when the property maps are loosely segmented, but not when tightly segmented. Errors in toneburst amplitude and time-of-flight resulting from simulations through clinical CT data are summarised in tables 3 and 4, comparing average errors in simulations through frontal and parietal bone with simulations through temporal bone. Tight segmentation of property maps gives reduced error when attenuation is not simulated, and the reduced error for attenuating, loosely segmented property maps is likely due to the larger segmentation volume over which attenuation is introduced. These results also appear to demonstrate that homogeneous models derived from clinical CT data tend to perform approximately the same as heterogeneous models derived from clinical CT data, at least in terms of the error metrics computed here.
Results for clinical CT maps that result from linear interpolation onto the simulation grid are shown alongside the homogenised micro-CT data in figure 6 for comparison with the micro-CT data. Errors calculated for simulations through heterogeneous property map derived from clinical CT data are plotted on the axis in a position corresponding to the resolution of the CT images, while results for the homogeneous medium maps are plotted separately. Simulations through maps derived from clinical CT data show similar errors in average timeof-flight to the micro-CT data homogenised at the same resolution, but the error in amplitude over the recording plane is much larger. Accounting for attenuation leads to a large reduction transmitted amplitude error, and a small reduction in time-of-flight error.
Examples of the peak pressure maps across the intracranial scan planes are shown in figure 7 for simulations through the parietal bone sample from skull B. They show fields resulting from propagation through acoustic property maps derived from the micro-CT data at various levels of downsampling, and through heterogeneous and homogeneous maps derived from clinical CT data (non-absorbing maps are shown). This demonstrates the increase in transmitted amplitude with homogenisation of the micro-CT maps, and the improved performance obtained from using a nearest neighbour algorithm. In terms of the fields resulting from transmission through clinical CT data, while tightly segmenting the acoustic property maps and using a nearest-neighbour algorithm did restore some finer spatial detail to the fields, it is not apparent that the new field shape more accurately represents the reference simulations.

Time-reversal
Generally, the results for the time-reversal simulations follow the same trends as those observed for the forward simulation tests. Homogenisation of micro-CT leads to an overestimation of the transmitted amplitude for both CW and TB simulations. It leads to an underestimation of simulated time-of-flight for TB data, which manifests in the CW data as an increasing error in phase. Figure 8 summarises these results, showing the average of each error metric across frontal and parietal bone segments, with results for the temporal bone plotted separately to illustrate any difference. As with forward simulations, the error in amplitude is lower for temporal bone, where it shows a steady increase with homogenisation, rather than the rapid increase between 0.05 mm and 0.5 mm observed for frontal and parietal bone. The difference is more pronounced for errors in time-of-flight and phase, to the extent that the simulations through temporal bone show an increase in time-of-flight, rather than the reduction observed for frontal and parietal bone.
Average errors for TB simulations resulting from maps derived from clinical CT data for frontal and parietal samples are shown in tables 5 and 6. For the heterogeneous maps without attenuation, the errors in time-offlight, phase (not shown), and amplitude are similar to those obtained from micro-CT data homogenised at approximately the same spatial resolution. As with forward simulations, the use of nearest neighbour interpolation and tightly segmenting acoustic property maps leads to reduced error, while adding attenuation leads to a large reduction in amplitude error. For TB simulations, introducing attenuation leads to a slightly worse error in time-of-flight, but almost no change in phase for single frequency, CW simulations. This could be due to frequency dependent absorption leading to a change in the pulse shape for TB simulations, or due to greater attenuation of longer ultrasound paths through the bone leading to a shift in the signal recorded at the transducer surface. As with the forward simulation setup, propagation through homogenous models was at least as erroneous as simulations through heterogeneous models. For reversal simulations, it should be noted that the relationship between the error in the signal recorded at the virtual transducer and the resulting impact on the intracranial fields when using these simulations to focus is not straightforward.
If it is necessary to reconstruct the exact time varying pressure signal, then errors in the amplitude and timeof-flight assigned to the transducer element will lead to corresponding errors in the amplitude and timing of the pressure field at the target. However, if the error in the amplitude and phase (or time-of-flight) is constant across each element, the relative phase and amplitude will be conserved. This means that the intracranial field should be effectively spatially focused, albeit potentially with an incorrect amplitude. In this respect, it is not the absolute error in amplitude and phase/time-of-flight that is important, but the variation in this error between different transducer elements. Figure 9 shows the standard deviation in the error in both amplitude and time-of-flight for reversal simulations using a TB source, for all bone samples (including temporal bone). It demonstrates how the standard deviation in both amplitude and time-of-flight error increases with increasing homogenisation. Amplitude error shows a standard deviation of over 60% for homogenisation on the scale of clinical CT images, while time-of flight varies by over 0.3 μs. The variation in the error resulting from simulations through linearly interpolated clinical CT data is also shown. Heterogeneous property maps show similar variance in error to the homogenised micro-CT data, while homogeneous property maps show greater variance in the non-attenuating case.

Elastic simulations
Results for simulations through 2D property maps of the skull, with and without shear mode propagation modelled, are shown in figure 10. Error in the signal across the transducer surface is computed relative to simulations through non-homogenised micro-CT data without elastic medium properties assigned. For the non-homogenised micro-CT data, the introduction of shear acoustic properties leads to an average reduction in the amplitude of 33% for frontal and parietal bone, and an average reduction in time-of-flight of 0.45 μs. However, the difference in amplitude error between elastic and non-elastic simulations does not increase with homogenisation, and the difference in time-of-flight error reduces with increasing homogenisation.
For temporal bone simulations, the difference between elastic and non-elastic simulations is much smaller, consistent with previous studies (Robertson et al 2017b). Simulations through the non-homogenised map result in an error in amplitude of 1.2%, and an error in time-of-flight of 0.02 μs when introducing shear mode propagation. This indicates that the influence of shear mode properties is mainly the result of wave propagation in the  trabecular bone layer. Overall, these results indicate that accounting for elastic wave propagation in the bone tissue can have a significant effect when modelling the propagation of ultrasound across the internal microstructure of the skull. However, the inclusion of elastic wave propagation does not appear to modify the change in transmission through the skull that occurs when acoustic property maps are progressively homogenised.

Discussion
In the present paper, it has been demonstrated that homogenisation of acoustic property maps derived from micro-CT data leads to an underestimation of the transmission loss when simulating the propagation of ultrasound through skull bone. It also introduces errors in time-of flight or phase on the scale of the acoustic period of 1 MHz ultrasound. It was observed that the use of a nearest neighbour algorithm for the downsampling of micro-CT property maps results in reduced error compared to linear downsampling. However, this trend only holds for simulations through heterogeneous maps derived from clinical CT data when loosely segmented, but not when tightly segmented. At the same time, simulations through tightly segmented clinical CT data consistently give reduced error before attenuation is considered. Together these results highlight the importance of defining the edges of the bone layer as a discontinuous, rather than smoothly varying, interface. This is likely because, when acoustic properties vary over a finite distance (rather than forming sharp, discontinuous interfaces), the resulting reflection and transmission coefficients demonstrate frequency dependent behaviour determined by the spatial rate of change of acoustic properties (Cieszko et al 2012). Ultrasound with wavelengths close to the spatial extent of the smoothly varying interface (or shorter) experience reduced reflection and increased transmission, consistent with the results shown here.
To determine the extent to which the impact of image related homogenisation is caused specifically by loss of trabecular microstructure, comparative simulations were carried out through a sample of temporal bone. Simulations through temporal bone demonstrated a large reduction in the error in simulated time-of-flight and phase. Amplitude error increased more slowly with increasing homogenisation, but was still significant. This suggests that the effects of homogenisation on simulated time-of-flight are due to the loss of heterogeneous internal microstructure not found in temporal bone, possibly due to changes in the mean pathlength through the  skull layer. Likewise, it suggests that the error in simulated amplitude is due to a combination of loss of internal structure and other effects of homogenisation such as smoothing. To place these results in context, they can be compared with a recent paper by Grimal et al (2014), who examined the impact of changes in the porosity of cortical bone in the femoral neck on simulated time-of-flight and amplitude. They observed a reduction in simulated time-of-flight, and an increase in transmitted amplitude when the internal heterogeneity of the cortical bone layer was removed, consistent with the results here. They also observed that the changes in the transmitted amplitude that occur due to the loss of the internal scattering structure can be corrected for by accounting for the porosity and average pore size in the average acoustic properties assigned to the cortical bone layer.
Simulations were also carried out through both heterogeneous and homogeneous acoustic property maps derived from clinical CT images co-registered to the micro-CT data. For the time-reversal simulation setup, the errors resulting from simulation through heterogeneous medium property maps derived from clinical CT data are on the same scale as the errors resulting from micro-CT data homogenised at the same resolution. For the forward propagation simulations, the error across the intracranial field is consistently larger for simulations through clinical CT data. Modelling attenuation (inflated to account for scattering losses) leads to a large reduction in the error in transmitted pressure amplitude across all simulations, indicating that the attenuation value measured by Pinton et al (2011) is effectively able to account for a portion of the impact of image related homogenisation. For simulations using a toneburst source, the addition of attenuation actually leads to a change in the error in time-of-flight, increasing it for reversal simulations, and reducing it for forward simulations. The reasons for the change in time-of-flight and why forward and reversal simulations demonstrate opposite trends, is unclear.
With respect to the homogeneous acoustic property maps, the error relative to the reference simulations is generally similar to that obtained using heterogeneous medium maps, although the variation in these errors (shown in figure 9) seems to be higher. In contrast, homogeneous acoustic property maps have previously been identified as giving reduced performance when using simulations of ultrasound propagation in time-reversal focusing (Jones and Hynynen 2016). This suggests that the improved performance of heterogeneous maps is not due to their ability to account for small scale acoustic property variations inside the skull layer, but instead due to the fact that a heterogeneous model will account for gross variations in acoustic properties across the skull surface. It is likely that the homogeneous models used in the current study demonstrate better performance because they were based on average acoustic property values across a small volume of bone tissue. When an entire skull map is assigned a single sound speed and density value, it fails to account for the variation in average property values beneath each transducer element. The fact that homogeneous models ensure a discontinuous interface between bone and soft-tissue may also contribute to their performance (Robertson et al 2017b).
A limitation of the present study is the limited exploration of the impact of elastic wave propagation, which was shown to have an impact on simulation through micro-CT data. Previous studies have indicated that shear mode propagation has a negligible impact on propagation through acoustic property maps derived from clinical CT images when applied at normal incidence (White et al 2006). While the results shown in figure 10 have demonstrated that elastic mode propagation does not affect the impact of image related homogenisation, it may limit the usefulness of the comparison with simulations through clinical-CT data. Further work, including comparison of micro-CT simulations with practical experiments would be of use to validate the use of micro-CT data as a gold standard in the simulation of transcranial ultrasound propagation. An additional limitation is the use of a single frequency to examine the impact of image related homogenisation. It is expected that the impact of image homogenisation will vary depending on the ultrasound wavelength, as the impact of the trabecular bone is known to vary with the frequency of the applied ultrasound (Fry and Barger 1978). In the present paper, 1 MHz was chosen because it lies at the upper end of the range of frequencies effective for transcranial focusing (Deffieux and Konofagou 2010). While it is unlikely that quantitative results presented here can be directly extrapolated to different ultrasound frequencies, it is possible that by considering the scale of the homogenisation voxel relative to the wavelength of the applied ultrasound that general trends can be inferred.

Conclusion
Simulations of the propagation of ultrasound through acoustic property maps derived from micro-CT data were used to examine the impact of image related homogenisation effects. Homogenisation results in a significant reduction in transmission loss due to smoothing of the interface of bone and soft-tissue, and loss of internal microstructure. The loss of internal microstructure also results in a change in simulated time of flight, likely due to a change in the mean pathlength through the bone. These results were compared with simulations through co-registered clinical CT images, which demonstrated the importance of correcting for the loss of scattering by assigning an inflated attenuation coefficient. Simulations through both clinical and micro-CT data highlighted the impact of enforcing a sharp boundary between bone and soft tissue.
The results presented here suggest that homogenisation of acoustic property maps on the scale of clinical CT images may have a significant impact on the simulation of ultrasound through skull bone. While the impact on the effectiveness of simulation-based time-reversal focusing is less clear, given the history of successful focusing using clinical CT data, these results have clear relevance to current efforts to focus ultrasound onto the brain for novel technologies.