Redefining Ultrasound Compounding: Computational Sonography

11/05/2018 ∙ by Rüdiger Göbl, et al. ∙ 2

Freehand three-dimensional ultrasound (3D-US) has gained considerable interest in research, but even today suffers from its high inter-operator variability in clinical practice. The high variability mainly arises from tracking inaccuracies as well as the directionality of the ultrasound data, being neglected in most of today's reconstruction methods. By providing a novel paradigm for the acquisition and reconstruction of tracked freehand 3D ultrasound, this work presents the concept of Computational Sonography (CS) to model the directionality of ultrasound information. CS preserves the directionality of the acquired data, and allows for its exploitation by computational algorithms. In this regard, we propose a set of mathematical models to represent 3D-US data, inspired by the physics of ultrasound imaging. We compare different models of Computational Sonography to classical scalar compounding for freehand acquisitions, providing both an improved preservation of US directionality as well as improved image quality in 3D. The novel concept is evaluated for a set of phantom datasets, as well as for in-vivo acquisitions of muscoloskeletal and vascular applications.



There are no comments yet.


page 1

page 7

page 9

page 10

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In recent years, 3D ultrasound (US) imaging has gained interest [1, 2] both due to the greater availability of 3D ultrasound systems (tracked, mechanical driven or matrix arrays)[3, 4]

as well as to the improvements in computational power allowing for faster and online processing of the acquisitions. To date, the most common way to obtain a 3D US volume remains the interpolation of 2D images into a 3D volume 

[5]. We focus here on tracked 3D freehand ultrasound, a technique that yields 3D US volumes by first acquiring 2D images along with their position and orientation and then fusing them into a volumetric regular grid through a process known as 3D reconstruction or compounding [5]. Up to now, conventional compounding techniques (e.g. forward, backward or functional interpolation) have neglected an essential property of ultrasound images, namely, its direction-dependent nature. It is known that imaging a point from different directions may lead to different intensities. This direction-dependency is caused by the different attenuation paths followed by the reflected signals, as well as by the dependency of the reflected intensity on the tissue’s impedance and orientations as defined by Snell’s law. This directional behavior in part explains the information loss resultant from current 3D reconstruction methods[6] and is the main reason for restricted acquisition protocols (e.g. straight probe motion) recommended in practice. It is to be noted that such protocols modify the physician’s common practice based on image feedback and interactive motion of the probe.

In this work, we introduce Computational Sonography (CS) as a novel paradigm for the 3D reconstruction and representation of tracked freehand 3D US that opens new possibilities for visualizing, interacting and processing volumetric data. In more detail, we introduce two major changes to current reconstruction paradigms (cf.Fig. 5). First, we consider and preserve the directional dependency of US in the 3D reconstruction. And second, we enable each element in the 3D volume to store a complex model instead of just a single scalar value. We argue that with such complex models it becomes possible to retain the directional information previously discarded and make it vailable for further use. We want to point out that the “Computational” term in CS, is an analogy to the concept of Computational Photography [7], and stands for the fact that once such a more complex volumetric representation has been stored, we deploy computational algorithms in order to extract and display the relevant information on demand.

Figure 4: Illustration of the directionality dependency in US imaging using images of a phantom depicted in (a). (b) When the probe is perpendicular to the wires these appear clearly in the image. (c) When the probe is tilted the wires may not be visible.

The goal of this paper is to demonstrate how the concept of Computational Sonography can be translated into new direction-preserving compounding methods. In practice, we present two methods that reconstruct the amount of reflected signal expected from every direction at each point in a volume. An important design specification is the model to represent these values, where desired properties are the directionality-preservation but also compactness, and fixed-length per voxel. We study two representation models complying with the requirements: tensors, as initially explored in a conference paper 

[8], and a novel spherical representation. In summary, our contributions to 3D sonography are:

  • A reconsideration of tacit assumptions in the process of fusing multiple US views regarding i) the directionality of US, and ii) the reconstruction of scalar volumes.

  • The introduction of two models to represent multi-directional data, namely tensor-CS and spherical-CS.

  • The analysis of two direction-preserving compounding algorithms and their comparison to classical compounding.

  • The demonstration of key capabilities of CS.

Figure 5: Computational Sonography concept.

Finally, to show the practical value, we quantitatively evaluate the fidelity of the CS reconstructions (vs. information loss). The qualitative and quantitative evaluations are performed for both phantom and in-vivo data in comparison to state-of-the-art compounding methods.

This work follows an initial conference publication [8] but includes the following extensions: in addition to the tensor model we propose a more expressive spherical representation, which encodes the directionality of ultrasound-information with minimal information-loss. We also further analyse in-depth and quantitatively compare the two models against conventional compounding methods for both phantom and in-vivo acquisitions. In particular, we have extended the evaluation of the fidelity from 15 in-vivo and in-vitro cases to 40 cases.

1.1 Related Work

Only few methods have considered the directionality information in the processing 3D US volumes. Efforts in our group include a method to solve quadratic optimization problems in freehand 3D US in a compounding-free fashion [9], as well as the design of directional dependent filters to vessels in 3D US images [10]. For 3D multiple-view US reconstruction, Schulte zu Berge et al.[11] use the orientation information to cluster 2D slices, compound one volume per cluster and fuse the clusters. Different to CS, this method reconstructs only a scalar field, the directional information is lost for further processing.

Towards reconstructing more complex models per voxel, Klein et al.[12]

cluster the samples in a voxel according to their viewing angles but combine them into a volume of mixture models of Nakagami distributions, resulting in a vector containing the distribution parameters per voxel. However, the directionality information is also lost after reconstruction. Recently, Lee

et al.[13] and Papadacii et al.[14] proposed tensorial-based reconstructions (respectively Elastic Tensor Imaging (ETI) and Backscatter Tensor Imaging (BTI)) of US images from multiple view angles. Based on shear wave imaging, ETI measures the anisotropy of tissue elasticity and represents it with a tensor. In BTI, the local spatial coherence of the US signals observed for different orientations and focal regions was used to reconstruct 3D tensor volumes. Although these methods also recover tensor fields, their focus is on recovering the fiber orientation of the underlying tissues. Additionally, they both rely on a specialized hardware setup and a fixed scanning protocol. In contrast, our approach focuses on information preservation of the full acquisition and can be seamlessly integrated into the clinical workflow as well as applied directly to clinical 3D freehand US acquisitions. Finally, beyond tensor fields, we have also proposed a new representation for 3D US, Spherical-CS.

Representing US directional dependent information is similar to the problem of describing diffusion patterns in Diffusion MRI. In this field, both tensors (Diffusion Tensor Imaging-DTI[15]) and spherical representations (Q-ball[16]) have been studied. However, apart from the representation model, the methods to fit and visualize tensors and q-balls are specialized to the physics of MRI acquisition and therefore cannot be applied in a straightforward manner to 3D US data.

In conclusion and to the best of our knowledge, this is the first orientation preserving reconstruction method for 3D ultrasound. The proposed CS-models bear some similarity to those used in diffusion MRI, but their computation from ultrasound data is new.

2 Computational Sonography Concept

We investigate the problem of building a direction-dependent 3D US volumetric reconstruction from a set of tracked 2D US images consisting of a set of tracked US samples111US images are composed of rays, each formed by samples resulting from measuring the reflected signal at regular times.. We denote the set of acquired samples and perform the reconstruction on a 3D volumetric grid in a lattice , with . We can then decompose the reconstruction problem in two tasks:

  1. Finding the relevant subset of samples contributing to the reconstruction of each voxel .

  2. Fitting the chosen samples to the orientation-dependent model stored in each voxel .

While we follow the strategy of [17] for the first stage (cf.subsec. 3.3), we focus hereafter on the reconstruction task. The problem can be formalized as finding a reconstruction function mapping from the set of measured samples to a model for every voxel in ,


Let each US sample be characterized by a tuple , where is the measured reflected signal, and and stand respectively for the sample’s position and acquisition direction. In conventional compounding methods, maps the samples to a scalar value . Accordingly, the mean reconstruction can be formalized as


Other common scalar reconstructions include the median, k-Nearest Neighbor or Inverse Distance Weighted compounding methods.

With computational sonography, we propose two new reconstruction functions , which fuse samples while preserving ultrasound’s directionality; i.e. they capture the fact that the reconstruction may show different intensity values for different directions. The two proposed methods, one tensor-based and another spherical-based , differ in the fusion as well in the final model used for compact, but direction preserving, representation of the reconstruction, as detailed in sec. 3.

The result of the two CS reconstruction procedures is a compact representation of the full acquisition, which preserves the directional nature of US for visualization, as well as other possible common post-processing tasks.

3 Methodology

In order to design a direction-preserving reconstruction function , we propose to store a representation model instead of a scalar value in each element of a target lattice. We present two reconstruction functions and representation models. The first, initially introduced in [8], is a tensor-based approach, Tensor Computational Sonography (Tensor-CS). The second, is a novel approach relying on the concept of spherical grids (Spherical-CS).

3.1 Tensor-Models

For Tensor-CS, we choose to fill the output volumetric grid with tensors, that is we store a tensor model in each voxel. The reconstruction function of Tensor-CS is implicitly defined as:


where the tensor is assumed to be symmetric. The direction preservation is achieved by means of the constraint , which states that tensor should compactly encode the intensity value observed from the given direction . In practice, a least squares approach is employed to find the values of the tensor , where we minimize


As shall be symmetric, at least six sampling points determine the tensor coefficients. In practice, we use typically between 12 and 500 samples, as this showed to provide best results in our experiments. We solve the least squares problem using Cholesky decompositions implemented on a GPU.

The main advantage of the tensor model is that it allows for a direction dependent modeling of the ultrasound signal while requiring only six real values to be stored per voxel, making Tensor-CS very suitable for parallel processing. The tensor model is also suited for an interactive orientation-dependent visualization (see subsec. 3.4). Finally, tensor fields can be further post-processed by means of tensor-based regularization and visualization methods, such as those used in DTI [18, 19, 20]. On the other hand, the tensors’ compactness comes at the price of resolution loss and some simplifications. For instance, the tensor modeling assumes that signals are point symmetric. While the symmetry is irrelevant for the most common freehand acquisitions, it becomes a problem for acquisitions from all orientations, and in particular, in the presence of acoustic shadowing. This phenomenon is illustrated in Fig. a and Fig. b.

Figure 10: Symmetry and two dimensional Computational Sonography. We simulate two distributions of measured samples, shown as gray dots. (a) The first distribution is point-symmetric. (b) In the second distribution the symmetry is broken by a “shadow". In (c), only measurements from a restricted angle are considered. The reconstruction with Tensor-CS (red) fits nicely to noisy point-symmetric data (a), but has problems with dissymmetries, as they could result from shadowing (b). These artefacts are somewhat mitigated when acquisitions are taken only from a restricted angle as in (c). The reconstruction obtained with Spherical-CS (blue) accurately approximates the signal in the three cases, as shown by the blue lines from (a) to (c).

3.2 Spherical Models

To overcome some of the limitations of Tensor-CS regarding the resolution and enforced symmetry, we propose a new representation for 3D US data. More precisely, for each voxel we model the ultrasound signal reflected from all directions as a continuous real function defined on the unitary sphere . This model allows a unique intensity value for every direction , or equivalently, . As a result, the Spherical-CS representation consists of one spherical function per element of the lattice , that is, . Recall that the samples relevant to voxel are defined by , i.e. by their intensity , position and imaging direction . We then try to approximate the spherical function from the samples’ tuples . Although one could store every sample relevant to each voxel to approximate with the highest fidelity, the resulting representation would be very inefficient, due to a varying number of samples. Instead, we choose as domain a spherical grid, computing a discrete approximation of . The domain is defined as a partitioning the sphere into segments , such that . Then, we approximate a scalar value per cell. Finally, we collect the values from the different cells to store a model in each voxel. Thereby, the sought reconstruction function is:


In practice, each entry of the vector is computed only from the samples falling inside the corresponding cell . First we determine the subset of samples associated to each segment . Then, the value assigned to the cell is obtained from the samples intensity values through one interpolation schemes of


where stands for the number of samples in segment . They are derived from classical compounding methods, but applied only to the associated sample subset

This discretization scheme permits adjusting the desired resolution of the representation. In this way, we achieve a compact, yet accurate representation of the acquisition allowing for subsequent processing or extraction of information according to the needs of the user. The proposed models are also well capable of representing arbitrary functions including non-point-symmetric data (cf.Fig. 10).

A challenge with the spherical representation arises from the discretization of the sphere, which is itself an open subject of research [21]. We are interested in a discretization that: i) leads to a regular partition of the area on the sphere, and ii) enables computational efficient mappings from a directions to the respective segment and backwards, as these mappings will be needed for the reconstruction as well as for later access, processing or visualization.

Next, we analyse and discuss three possible discretization approaches, as well as their mappings. To this end, we first define a grid on the unit sphere as a set of points , where each point is described by its spherical coordinates, namely the azimuthal and elevation angles, i.e. . Then, the Voronoi cells of the grid points inherently define a partition of the sphere.

3.2.1 Latitude-longitude grids

Given a sphere parametrized with spherical coordinates, the latitude-longitude discretization is simply obtained through an equidistant sampling of the azimuth and elevation angles. As a result, it is trivial to map direction vectors to their respective point on the spherical grid by converting Euclidean coordinates to spherical coordinates and visce-versa. More precisely, the coordinates of the points on a latitude-longitude grid with a resolution of radians are given in spherical coordinates by


As it is apparent in Fig. a, the inconvenience of latitude-longitude grids is that the density of points varies significantly across the surface of the sphere, which results in cells of variable area, and is undesired for computational sonography.

3.2.2 Icosahedral geodesic grids

Geodesic grids are built by iteratively subdividing a polyhedron. It is common to use a icosahedron as a starting point. At every iteration, each polygon of the current polyhedron is further divided in cells, and projected onto the sphere. This subdivision is continued until the desired resolution is reached. Following this discretization, the obtained partitions have near-equal area. For a detailed description we refer to the work of Thuburn [22]. Fig. b shows a geodesic grid based on an icosahedron subdivision. Similar to the latitude-longitude grids, the spherical coordinates of the grid-points are known by construction and the mapping to their represented directions is trivial. However, no constant-time mapping of arbitrary direction vectors to the corresponding grid-point (or directional segment) exists for these grids [23]. Since measurements can originate from arbitrary directions, icosahedral geodesic grids are not practical from a computational point of view.

3.2.3 Fibonacci geodesic grids

The third type of discretization we studied here are spherical Fibonacci grids, which are generated by drawing points along a spiral on the surface of the sphere. The construction rules describe each grid-point in terms of the golden ratio [24], as follows:


Fig. c shows a spherical Fibonacci grid with 42 grid-points. Mapping a grid-point to a direction vector is equivalent to converting spherical to Euclidean coordinates. Most importantly, however, the inverse mapping recently proposed by Keinert et al.[23] provides a mapping from any direction vector to the nearest grid-point or directional segment with constant run-time. We refer the reader to the work for the derivation of the inverse mapping. This property makes Fibonacci geodesic grids the most suited for computational sonography.

Figure 15: Voronoi Diagrams (a)-(c) showing the extent of each cell in the partition of the sphere. (a) Latitude-longitude grid with 30 degree resolution. (b) Icosahedral spherical grid with 42 grid-points. (c) Spherical Fibonacci grid with 42 grid-points. (d) Exemplary tensor evaluated in all directions.

3.3 Sample Selection

As discussed in sec. 2, the compounding process is composed of two stages. So far we focused on the reconstruction step and assumed the subset of measurements relevant to fit one voxel model is known. This subset is determined using a sample selection strategy initially described in [8], and inspired by [17]. For completeness, we describe the process next.

The task is to select all samples which contribute to the reconstruction of a specific voxel, centered at coordinates . Practically, we select a sample if the voxel center lies inside the region of influence of the ultrasound sample. This region is defined as an ellipsoid centered around the sample’s position and with one of its axis aligned to the sample’s direction . Formally, if:


where represent the half of the length of the ellipsoid’s axes in each direction.

Figure 16: Sample Selection Process. Three rays pass through the voxel centered at . A sample is considered as relevant to if it falls within the region of action of a sample defined by an ellipsoid (cf. (12)). Only one sample per ray is taken into account.

Fig. 16 shows the sample selection process exemplarily. Based on (12), multiple samples of each scan line could be selected depending on the proximity of a ray to the target position . As the main interest for the reconstruction is to collect complementary information, only the nearest sample of each ray is retained in .

3.4 Computational Aspects of CS

After the samples have been selected, and a model (either tensor or discretized spherical function) has been computed per voxel, the 3D reconstruction process is finished, and the result is a 3D lattice of models. In this section we discuss the “computational" aspect of CS, where we will make use of algorithms to extract the relevant information from the lattice both on-demand and according to the application. Within the scope of this paper we restrict the algorithms to the extraction of the information needed for visualization.

3.4.1 Tensor visualization

From the Tensor-CS reconstruction, it is possible to extract directional images, i.e. images of the reflected ultrasound signal from a given direction. The directional intensity for a voxel from direction , with

, can be estimated from the tensors with the second order approximation:


Thereby, (13) allows us to interactively visualize the image as if taken from an arbitrary viewpoint, even if not previously imaged. Such interactive direction-dependent visualization can serve as a basis for realistic simulations for training.

Secondly, from the Tensor-CS we can also retrieve a variety of 3D scalar representations suitable for simple 3D renderings. The first option is to compute the absolute value of the tensor traces and use them as the intensities. In this way the intensity at voxel becomes

. The tensor trace amounts to display the accumulated signal intensity from the three major components (eigenvectors) of the tensor. A second possibility is to extract the intensity from the tensor’s dominant orientation corresponding to the largest eigenvalue of the tensor


3.4.2 Visualization of discretized spherical functions

In the case of the Spherical-CS, and in particular for the geodesic grids, we studied the visualization of two quantities: the mean and the maximum of the intensity values stored in the different cells of the geodesic grids.


These values are suitable for showing respectively, the most representative intensity in each voxel and the intensity from the direction with highest response. Examples for these visualizations can be seen in subsec. 4.4. These choices are motivated to maximize the speed of the visualization task. It is worth noting though that the representation is open to more complex reconstruction approaches.

4 Experimental Validation

To evaluate the different reconstruction capabilities of Computational Sonography we present next a series of experiments on both phantom and in-vivo datasets (subsec. 4.1). First, we perform a quantitative evaluation of the CS representation quality in terms of information preservation and in comparison to classical compounding (subsec. 4.2). We then investigate the robustness of the representation within a voxel to increasing variations of the intensity values, or to an increasing coverage of the acquisition directions (subsec. 4.3). We also visually illustrate the direct use of the CS representation to create scalar 3D volumes (subsec. 4.4). Finally, we present exemplary uses of the additional information present in CS volumes in subsec. 4.5. The evaluation is based on a comparison of the following spatial compounding methods:

  • Mean compounding.

  • Tensor-CS.

  • Spherical-CS with latitude-longitude grid.

  • Spherical-CS with spherical Fibonacci grid.

For all the Spherical-CS reconstructions we used the mean of sample intensities in each cell (6). For the Tensor-CS one constraint not enforced so far is the expected positive-definiteness of the tensors (), arising from the the fact that intensities are by nature positive.

4.1 In-vitro and In-vivo Data Acquisition

TypeSeq. No. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Total
Vessel phantom 125 217 556 633 725 1166 149 269 100 561 592 1000 6093
Carotid 160 50 292 336 437 479 224 379 25 197 209 187 204 3179
Femoral 89 377 312 78 209 188 306 81 31 249 91 232 204 90 187 2724
Table 1: Dataset description including number of frames per sequence. The dataset consists of three types of image sequences including a vessel phantom as well as carotid and femoral in-vivo sequences. For every type we recorded several sequences, each with a varying coverage of the orientations.

We acquired a number of phantom and in-vivo datasets using a 3D-Ultrasound system consisting of an Ultrasonix SonixMDP machine with a linear transducer and a NDI Polaris Vicra optical tracking system. Every sweep is composed of several B-mode images acquired in a sequence while smoothly moving the probe. All acquired 3D-US sweeps differ in the number of images and in the span of directions covered. We followed a protocol where we gradually varied from sweeps taken for a single direction (as commonly done in today’s 3D US freehand protocols) to sweeps densely covering the sphere.

Figure 17: Latex vessel phantom used for static acquisitions.

In practice, we acquired 12 sweeps of a static latex vessel phantom (cf. Fig. 17) and 28 in-vivo sweeps, including 13 vascular scans (carotid arteries) and 15 scans in an osteological setting (femoral bones). The number and size of the different acquisition sequences are detailed in Table 1 while the acquisition parameters are summarized in Table 2.

Parameter Value
US probe L14-5/38 Linear
Frequency 6.6MHz
Depth carotid 4.0cm
Depth femoral 5.5cm
Depth vessel phantom 4.0cm
Sample selection (12)
Reconstruction resolution 0.5mm spacing
Table 2: Acquisition parameters

4.2 Representation Quality

Using all in-vitro and in-vivo datasets above we evaluate the representation capabilities of CS. First, we evaluate the representation quality in terms of how well a reconstructed volume preserves the information from the original ultrasound images, which is standard evaluation criterion for 3D reconstruction methods [5]. In particular, we focus on quantifying the preservation of the directional intensity, i.e. the intensity observed from different directions. To this end, and in contrast to the evaluation of classical compounding methods based on restricted scanning protocols (i.e. parallel images), our datasets contain information from various probe directions.

To quantify the information preservation, we reproject the volumes at the positions and directions from where the original US images were acquired. For a volume reconstructed from a sequence including samples , we measure the Mean Squared Error (MSE) between the input and the reprojected intensity values of every sample as follows:


where is the intensity of sample and is the intensity computed from the reprojection of the compounded volume at the position and direction of the original sample . The error is averaged over all input samples. For all experiments, 3D volumes were reconstructed with a 0.5mm spatial resolution. For the latitude-longitude grid a resolution of was used, whereas for the Fibonacci-geodesic reconstruction the grid was built with cells. For scalar volumes, the reprojection simply corresponds to interpolation. In the case of tensor-CS, we reproject using (13), whereas for spherical-CS we use the value of intensity stored in the associated grid cell. Examples of the reprojections for each method are shown in image Fig. 30. It can be seen that the 2D projections from computational-sonography volumes better preserve the information from the original B-mode images as CS avoids averaging information from different directions. Although the tensor-CS recovers appropriate intensity values for large portions of the image, the optimization fails in voxels where the information exceeds the low-complexity of the second order model or the symmetry. This is reflected as white regions in the images.

Figure 30: Reprojections. 2D slices obtained by projecting the volumes reconstructed with the three different methods in comparison to the original B-mode image.
Figure 31: Violin plot of observed representation errors in a carotid acquisition (Seq. 5, 437 frames).
Mean Tensor-CS Lat-Long-CS Geodesic-CS
Carotid 0.024  0.034 0.024  0.092 0.021  0.031 0.014 0.028
Femur 0.028  0.039 0.038  0.157 0.022  0.034 0.015 0.030
In-vitro 0.025  0.031 0.061  0.242 0.020  0.029 0.008 0.024
Overall 0.026  0.036 0.040  0.170 0.021  0.032 0.013 0.028
Table 3: Mean and SD of Squared Representation Errors. For the in-vivo and in-vitro acquisitions, comparing representations of B-mode intensities, for voxels for which a tensor could be estimated.

We report in Table 3

the mean and standard deviation of the representation error in (

16) for all volumes computed from the in-vivo and in-vitro acquisitions. For a fair comparison with the tensor-CS, we considered only the representation error of voxels for which a tensor could be estimated. The results indicate that representing the B-mode data by tensors leads to errors comparable to mean compounding, which is in accordance to the results in our prior work [8]. Spherical-CS based on Latitude-longitude grids does not substantially reduce the representation error, whereas geodesic-based spherical-CS results in the lowest average errors for all datasets. Fig. 31 shows violin plots of the representation errors found in a carotid acquisition (Seq. 5, 437 frames, ROI dir.var.  0.0257  0.0061) and emphasizes the similarity of error distributions between the presented representations.

4.3 Robustness to Intensity Variability and Spherical Coverage

In the following experiments, we investigate the robustness of the compared methods to two potential sources of inconsistency within the samples of a voxel: intensity variability and directional coverage.

To evaluate the robustness against the intensity-span

, we calculate the representation error for an increasing variance of intensity values associated to a voxel. Given intensities normalized to

, we computed the variance of the intensities for each voxel , and used a histogram to group voxels with similar variance. Figures (a)a and (b)b show the representation error for the different quantized levels of intensity-variations. It can be clearly seen that all Computational Sonography methods results in significantly lower errors, especially for larger values of intensity variance. Moreover for Spherical-CS the error is stable.

Figure 36: Robustness against intensity variability and directional coverage. (a) & (b) Mean square representation error for increasing intensity-spans within each voxel, for the in-vivo sweeps in (a) and the vessel phantom sweeps in (b). (c) & (d) Mean square representation error for an increasing directional coverage (angle-span) within each voxel. (c) shows the errors for the in-vivo acquisitions, (d) for the phantom sweeps.

We also analyze robustness to the spherical-span, that is, how the variance in the span of directions covered by the acquisition affects the representation error. We repeat the procedure above, but quantizing the spherical variance , as defined by Mardia and Jupp [25], where is the so called mean resultant length, which depends on the clustering of the vectors around the the mean direction. Note that is in

, corresponding to acquisitions going from a single direction to uniformly distributed over the sphere. We then estimate the representation error for the different levels of spherical variance. Figures

(c)c and (d)d shows how an increasing variance of directions affects the tensor representation the most, which is explained by the fact that the symmetry assumption is more likely broken in voxels compounding samples with high spherical variance. For the mean-compounding and spherical-CS the increase of spherical spans does not seem to have a noticeable impact on the reconstruction error.

The two graphs provide an experimental evidence to the fact that classical compounding approaches result in reconstruction errors when applied to unconstrained free-hand acquisitions. The reason for this is not directly the increased spherical variance but rather the ultrasound physics causing the same structure to have large intensity variations in different views. With regard to the two different CS representations, the tensor-CS suffers from a limited 2nd order model complexity, which reflects in a poor performance under large intensity variations On the other side, the strongest advantage of spherical-CS comes for acquisitions with heterogeneous ultrasound information (i.e. strong intensity variations) or high spherical variance. In fact, the error of spherical-CS reconstructions decreases for larger spherical variances, better achieving the goal of preserving the directional information, although at the cost of increased model complexity. These results further confirm computational sonography as a viable alternative to traditional reconstruction methods, allowing for high quality 3D free-hand reconstructions from arbitrary scanning trajectories instead of restricted scanning protocols.

4.4 Qualitative comparison

Exemplary results of the different reconstruction methods are shown in Fig. 62. They include axial and lateral slices for each of the sequence types: femoral, carotid and vessel phantom. For tensor-CS, we show the absolute value of the tensor trace. For spherical-CS, we use two forms of visualization: first, the intensities reflecting the mean of the cell values excluding empty cells, and second, the color-coded direction of the cell with the highest value. The color coding maps each vector in the volume to a unique color (cf.Fig. 37)222In practice, we first estimate the per-acquisition main direction by means of averaging the direction vectors in the volume. The color-coding then shows deviations from the main direction with increasing saturation, the hue is determined by the angle on the plane perpendicular to the main direction. Because the map considers 3D direction vectors, similar colors in different views stand for similar maximum directions.

Figure 37: Direction color-map w.r.t. the per-dataset main direction (black) and the direction determining the orientation of the map (red).
Figure 62: 3D Reconstructions from a femoral (Seq. 6, 166 frames, ROI dir.var.  0.0635  0.0165) and a carotid (Seq. 5, 437 frames, ROI dir.var.  0.0257  0.0061) in-vivo acquisitions as well as from a vessel phantom (Seq. 6, 1166 frames, ROI dir.var.  0.0421  0.0151). From left to right we compare the mean reconstruction, the absolute value of the tensor trace in tensor-CS, the mean of all grid-cells in a geodesic grid and in the last column we show the estimated surface normals in color coding (scaled with the geodesic mean). The red arrows show the main direction the color coding is based upon.

The images reveal first that Tensor-CS has an unstable behavior in unconstrained large scale evaluations. Given its reduced representation complexity it is prone to fail when the data from different directions is of higher order than our 2nd order model (the intensities vary significantly across different directions in a complex pattern). With the current implementation this effect is translated into a failure of the tensor fitting, shown as white voxels in the image. This effect seems to be located in boundary areas and locations where numerical instabilities are evoked by low intensities measured from a high number of directions, as can be seen in the phantom views. Despite this behaviour, it is apparent that whenever the tensor calculation succeeds the details of soft tissues are increased. Moreover, the tissue interfaces, e.g. in the two views of the carotid sequence, are best defined in Tensor-CS. In the case of Spherical-CS, the cell-mean visualization filters artefacts coming from directions with secondary reflections. This can be clearly seen for the femur lateral view: while the mean compounding features reflections behind the femoral interface, these reflections do not appear in the spherical-CS mean reconstruction. Also, a considerable blur reduction is observed when comparing spherical-CS to the mean compounding, in particular for the two views of the carotid and the lateral vessel sequence. The color visualization of the maximum Spherical-CS cells provides information about the interfaces’ normals, as can be seen in the femoral lateral view but more clearly as the normal of the carotid changes in the carotid lateral view or in the vessel phantom.

4.5 Free view point evaluation

Figure 65: Tensor-CS based Free-view visualization on a carotid scan (Seq. 5, 437 frames). The red arrows indicate the scanline directions of the reprojections. (a) shows a reprojection in original image direction, (b) the intensities derived form the tensor model in a different direction.

The directional information preserved in the CS models can be used to estimate intensities from arbitrary directions. This is achieved by evaluating (13) for tensor-CS with a chosen direction. Fig. 65 shows this for a carotid scan (Seq. 5), where the left is extracted from the original view-point and the image resulting from another direction on the right. It can be seen that the contrast between the interfaces and the tissue intensities is changed between the two directions.

5 Discussion

The qualitative results in combination with the quantitative evaluation show that geodesic spherical CS provides the best representation for arbitrary US-sweeps, while complying with the requirements of data fidelity and constant storage. The mappings from samples’s directions to the sphere are slightly more complex than for Tensor-CS, but this is compensated by the fact that as opposed to Tensor-CS, Spherical-CS does not rely on an optimization procedure to reconstruct the cell-values for a given voxel. As this optimization is based on all samples in a voxel, it can only be performed in a two-step reconstruction as we have presented it. The spherical-CS on the other hand could be implemented in an iterative reconstruction process, thus allowing for run-times comparable to state of the art scalar reconstructions. In addition, the current implementation of Tensor-CS has shown an unstable behaviour for general acquisitions with wide spherical variation. We believe that part of the instability could be solved by using more sophisticated optimization algorithms resulting in additional computational cost. However, the inherent limitations of the 2nd order symmetric model would still remain.

The concept of Computational Sonography opens new visualization paths for 3D US, relying on CS ability to store and then compute values on demand. In addition to the visualizations in Fig. 62 and according to the application, one may want to display the most meaningful cell for each voxel, to ensure faithfulness to the original measures. Online interactive visualization methods can be used to visualize the true directional-dependency of the ultrasound volume even after reconstruction and as such can serve as a basis for realistic simulations for training. Finally, we may also use the directional information to guide online acquisitions towards a viewpoint where the structure is most visible or to create surfaced visualization of interfaces.

Beyond visualization, Computational Sonography aims at compactly and reliably encoding all the information available during a free-hand acquisition in order to make it available for later processing. For instance, the CS concept is also useful for segmentation where the direction of the maximum cell gives hints about the presence of surfaces. Also for registration, CS allows us enables to compute the better similarity functions between a 2D slice and the reconstructed 3D volume, when accounting for the directional dependency of the images.

6 Conclusion

We showed the performance of several methods for the representation of the directional dependent information found in freehand US, using the conventional mean compounding as a baseline. Two storage models were proposed. A tensor-CS model, which we initially proposed in [8] and a novel approach based on spherical grids. We illustrated the adaption capabilities of the two representations, showing that the information can be extracted on demand according to the application. Finally, an exhaustive quantitative evaluation of the representation error resulted in a favorable performance of the spherical grids over the tensors and mean compounding methods. These results serve as evidence to our claim that with our new representation freehand acquisitions do not any longer need to rely on linear protocols.

Future research includes exploring other applications that require more complex algorithms to extract the relevant information, comparable to compressed sensing approaches. Additionally, one could devise methods using non-local information, for example to extract surfaces of tissue interfaces.


  • [1] L. J. Brattain, N. V. Vasilyev, and R. D. Howe, Enabling 3D Ultrasound Procedure Guidance through Enhanced Visualization.   Berlin, Heidelberg: Springer Berlin Heidelberg, 2012, pp. 115–124.
  • [2] H. Neshat, D. W. Cool, K. Barker, L. Gardi, N. Kakani, and A. Fenster, “A 3d ultrasound scanning system for image guided liver interventions,” Medical Physics, vol. 40, no. 11, 2013.
  • [3] A. Fenster, J. Bax, H. Neshat, D. Cool, N. Kakani, and C. Romagnoli, “3d ultrasound imaging in image-guided intervention,” in 2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Aug 2014, pp. 6151–6154.
  • [4] A. Vegas and M. Meineri, “Three-dimensional transesophageal echocardiography is a major advance for intraoperative clinical management of patients undergoing cardiac surgery: A core review,” Anesthesia and Analgesia, vol. 110, no. 6, pp. 1548–73, june 2010.
  • [5] O. V. Solberg, F. Lindseth, H. Torp, R. E. Blake, and T. A. N. Hernes, “Freehand 3d ultrasound reconstruction algorithms–a review,” Ultrasound in Medicine and Biology, vol. 33, no. 7, pp. 991 – 1009, 2007.
  • [6] R. Prager, A. Gee, G. Treece, and L. Berman, “Freehand 3d ultrasound without voxels: volume measurement and visualisation using the stradx system,” Ultrasonics, vol. 40, no. 1–8, pp. 109 – 115, 2002.
  • [7] R. Ng, M. Levoy, M. Brédif, G. Duval, M. Horowitz, and P. Hanrahan, “Light field photography with a hand-held plenoptic camera,” Computer Science Technical Report CSTR, vol. 2, no. 11, 2005.
  • [8] C. Hennersperger, M. Baust, D. Mateus, and N. Navab, “Computational sonography,” Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015, pp. 459–466, 2015.
  • [9] C. Hennersperger, D. Mateus, M. Baust, and N. Navab, “A quadratic energy minimization framework for signal loss estimation from arbitrarily sampled ultrasound data,” in Medical Image Computing and Computer-Assisted Intervention MICCAI 2014, ser. Lecture Notes in Computer Science, P. Golland, N. Hata, C. Barillot, J. Hornegger, and R. Howe, Eds.   Springer International Publishing, 2014, vol. 8674, pp. 373–380.
  • [10] C. Hennersperger, M. Baust, P. Waelkens, A. Karamalis, S.-A. Ahmadi, and N. Navab, “Multi-scale tubular structure detection in ultrasound imaging,” Medical Imaging, IEEE Transactions on, vol. 34, no. 1, pp. 13–26, Jan 2015.
  • [11] C. S. z. Berge, A. Kapoor, and N. Navab, Orientation-Driven Ultrasound Compounding Using Uncertainty Information.   Cham: Springer International Publishing, 2014, pp. 236–245.
  • [12] T. Klein, M. Hansson, and N. Navab, “Modeling of multi-view 3d freehand radio frequency ultrasound,” in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2012.   Springer, 2012, pp. 422–429.
  • [13] W.-N. Lee, M. Pernot et al., “Ultrasound elastic tensor imaging: comparison with mr diffusion tensor imaging in the myocardium,” Physics in medicine and biology, vol. 57, no. 16, p. 5075, 2012.
  • [14] C. Papadacci, M. Tanter, M. Pernot, and M. Fink, “Ultrasound backscatter tensor imaging (BTI): analysis of the spatial coherence of ultrasonic speckle in anisotropic soft tissues,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, vol. 61, no. 6, pp. 986–996, June 2014.
  • [15] D. Le Bihan, J.-F. Mangin, C. Poupon, C. A. Clark, S. Pappata, N. Molko, and H. Chabriat, “Diffusion tensor imaging: concepts and applications,” Journal of magnetic resonance imaging, vol. 13, no. 4, pp. 534–546, 2001.
  • [16] D. S. Tuch, “Q-ball imaging,” Magnetic resonance in medicine, vol. 52, no. 6, pp. 1358–1372, 2004.
  • [17] C. Hennersperger, A. Karamalis, and N. Navab, Vascular 3D+T Freehand Ultrasound Using Correlation of Doppler and Pulse-Oximetry Data.   Cham: Springer International Publishing, 2014, pp. 68–77.
  • [18] V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, “Fast and simple calculus on tensors in the log-euclidean framework,” in International Conference on Medical Image Computing and Computer-Assisted Intervention.   Springer, 2005, pp. 115–122.
  • [19] O. Christiansen, T.-M. Lee, J. Lie, U. Sinha, and T. F. Chan, “Total variation regularization of matrix-valued images,” International Journal of Biomedical Imaging, vol. 2007, 2007.
  • [20] G. Rosman, Y. Wang, X.-C. Tai, R. Kimmel, and A. M. Bruckstein, “Fast regularization of matrix-valued images,” in

    Efficient Algorithms for Global Optimization Methods in Computer Vision

    .   Springer, 2014, pp. 19–43.
  • [21] A. Staniforth and J. Thuburn, “Horizontal grids for global weather and climate prediction models: a review,” Quarterly Journal of the Royal Meteorological Society, vol. 138, no. 662, pp. 1–26, 2012.
  • [22] J. Thuburn, “A pv-based shallow-water model on a hexagonal icosahedral grid,” Monthly Weather Review, vol. 125, no. 9, pp. 2328–2347, 1997.
  • [23] B. Keinert, M. Innmann, M. Sänger, and M. Stamminger, “Spherical fibonacci mapping,” ACM Transactions on Graphics (TOG), vol. 34, no. 6, p. 193, 2015.
  • [24] R. Swinbank and R. James Purser, “Fibonacci grids: A novel approach to global modelling,” Quarterly Journal of the Royal Meteorological Society, vol. 132, no. 619, pp. 1769–1793, 2006.
  • [25] K. V. Mardia and P. E. Jupp, Directional Statistics

    , ser. Wiley series in probability and statistics.   John Wiley & Sons, Inc., 2000.