Speckle suppression and contrast enhancement in reconstruction of freehand 3D ultrasound images using an adaptive distance-weighted method

Applied Acoustics 70 (2009) 21–30

Contents lists available at ScienceDirect

Applied Acoustics journal homepage: www.elsevier.com/locate/apacoust

Speckle suppression and contrast enhancement in reconstruction of freehand 3D ultrasound images using an adaptive distance-weighted method Qinghua Huang a,c,*, Minhua Lu b, Yongping Zheng c, Zheru Chi d a

Department of Electronic Engineering, City University of Hong Kong, Kowloon, Hong Kong SAR, PR China Department of Biomedical Engineering, Shenzhen University, Shenzhen, PR China c Department of Health Technology and Informatics, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong SAR, PR China d Centre for Multimedia Signal Processing and Communications, Department of Electronic and Information Engineering, The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong SAR, PR China b

a r t i c l e

i n f o

Article history: Received 5 November 2007 Received in revised form 15 January 2008 Accepted 22 February 2008 Available online 15 April 2008 Keywords: 3D ultrasound Adaptive distance-weighted interpolation Volume reconstruction Speckle suppression Contrast enhancement

a b s t r a c t In three-dimensional (3D) freehand ultrasound (US), reconstructing a set of B-scans into a regular voxel array is the key procedure for consequent visualization and analysis. This paper presents a new adaptive interpolation algorithm for computing the voxel array to suppress speckle noises and enhance contrast. The local statistics of homogeneous regions including mean and variance were measured and the ratio of variance to mean was used as homogeneity criteria. For the computation of each voxel, the interpolation method was adaptively determined with respect to its local statistics. If the neighbouring pixels of a voxel satisfied the homogeneity criterion, its value was computed with an arithmetic mean filter. Otherwise, the voxel was probably locating in an inhomogeneous region and an adaptive distance-weighted (ADW) interpolation method was employed to compute its value. A resolution phantom and a subject’s forearm were reconstructed using the proposed algorithm and two other well-known methods – conventional distance-weighted (DW) and voxel nearest neighbourhood (VNN) interpolations. The comparison results demonstrated that the adaptive interpolation algorithm was able to suppress speckles, preserve edges and enhance contrast effectively for the volume reconstruction. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, three-dimensional (3D) ultrasound (US) has drawn great attention from both clinical and research fields, due to its advantages over traditional two-dimensional (2D) imaging, e.g., direct view of entire 3D anatomy, volume measurement and more accurate localization of anatomical parts [1–3]. Several techniques have been developed to obtain 3D US images. They were mainly grouped into two categories: using a dedicated 3D probe to acquire a small region of interest in real-time; and using a freehand scanning protocol with a spatial locator to sweep over a region of interest with arbitrary size and shape [4]. Dedicated 3D probes are relatively large and expensive in comparison with conventional 2D probes. Moreover, the field of view of such dedicated 3D probes was limited by the dimensions of piezoelectric elements in the probe housing, and their image resolution was not as good as their 2D counterparts [5]. In contrast, the tracked freehand scanning protocol is capable of overcoming the drawbacks of 3D probes, where the conventional 2D probe with a spatial locator at* Corresponding author. Tel.: +852 27887522; fax: +852 27844262. E-mail address: [email protected] (Q. Huang). 0003-682X/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.apacoust.2008.02.002

tached is moved by hand in an arbitrary manner and the acquisition of arbitrary volumes of interest can be accomplished. After the freehand scanning, a set of resulting 2D B-scans with high image resolution is reconstructed into a regular voxel array (volume) or a 3D surface using the corresponding positional information recorded by the spatial locator. Because of the low cost of imaging equipments and its capability of acquiring arbitrarily sized volumes, the freehand scanning for 3D US has been recognized as the cheapest and most flexible approach compared with other acquisition methods [3]. As the main purpose of 3D US is to provide an accurate 3D data set for qualitative and quantitative analyses, the volume reconstruction for irregularly spaced B-scans becomes the key procedure for tracked freehand systems. The accuracy of volume reconstruction critically depends on two stages, calibration and interpolation. Previous reports have addressed the issue of calibration in detail [5–7]. Although the quantitative analysis method without volume data set has been introduced [8], interpolation algorithms are commonly required to form a slice along any plane, a surface of a 3D object, or a 3D data set with regular voxel grids [9–13], and thus have been attracting more attentions.

22

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

Researchers in the field of 3D US have concerned about the issue of interpolation of scattered B-scan data in the past years [9,14– 19]. Because it was unavoidable that B-scan planes intersect and a voxel may be intersected by several B-scans, almost all 3D freehand systems had to deal with the problem of compounding [20]. In the 1990s, a number of attention was paid to the nearest neighbour (NN) approaches, including voxel nearest neighbour (VNN) and pixel nearest neighbour (PNN) interpolations (see a detailed survey conducted by Rohling et al. [9]). These two approaches are easy to implement. However, when the spacing of B-scans is large, reconstruction artefacts would be obvious using the VNN method and the gaps between B-scans could be highly smoothed if the PNN method was employed. Averaging as one of the important compounding operations has been considered as an optimal way for speckle reduction [17]. Barry et al. [14] proposed a distance-weighted (DW) interpolation method. Using the DW method, each voxel was computed by weighted average of its neighbouring pixels and the weights for the pixels were determined by their inverse distances to the voxel centre. The averaging operation in the DW method improved its performance in volume reconstruction with speckle and shadowing reduction. However, the averaging operation blurred image details, edges and contrast which originally occur in raw B-scans. Consequent efforts [15,18] have been made to improve the DW based interpolation for 3D freehand US. These studies demonstrated that the weights using higher power or the exponential of the distance between a neighbouring pixel and the voxel centre could produce 3D images with more texture patterns and enhanced contrast. Although improved image results with less smoothed details have been achieved, it remains a challenge for a certain averaging operation to simultaneously suppress speckles and preserve edges. Speckle suppression has long been investigated for medical ultrasonic images. Based on the local statistics of homogeneous speckle regions in B-mode ultrasonic images, many techniques [21–24] have been proposed to reduce speckles and preserve edges. Chen et al. [24] studied the speckle statistics and experimentally validated that the ratio of local variance to mean could be considered as a metrics to determine whether a local region was homogeneous containing only speckles [22,23]. Using the ratio, Huang and Zheng [19] proposed an adaptive squared distance-weighted (ASDW) interpolation method to both reduce speckle noises and preserve tissue edges. With the small ratio of local variance to mean, the ASDW was akin to a mean filter for speckle suppression; while it acted like a squared distanceweighted (SDW) to preserve more details in inhomogeneous regions if the ratio was relatively large. This method provided a good performance in speckle suppression, but failed to effectively enhance the contrast of small resolvable points. Another interpolation method made use of median filters instead of weighted mean filters for the volume reconstruction, and demonstrated good capability in edge preservation and noise reduction [25]. Unsatisfactorily, the reconstruction time is much longer than conventional methods, due to the great computation for data sorting. In this study, we proposed a new adaptive interpolation algorithm using local statistics of 2D images for the volume reconstruction of 3D freehand US. The local homogeneity of 2D US images, quantified by the ratio of the local variance to mean, was calculated and used to select reconstruction methods for the voxel to be constructed. The asymptotic value of the ratio of variance to mean measured at the relatively homogeneous regions was used as homogeneity threshold. If the ratio calculated for a specific voxel was below the threshold, it was considered locating in a more homogeneous region and would be computed with an arithmetic mean filter to reduce speckles. Otherwise, we used an adaptive distance-weighted (ADW) interpolation method. This could be classified as a new DW method, where the power of a pixel’s inverse

distance to the voxel centre was adaptively adjusted by the ratio of local variance/mean, thus effectively enhancing image contrast and preserving edges. The next section introduces the algorithm and experimental methods in detail. Section 3 demonstrates the in vitro and in vivo comparison results using the proposed method and two conventional interpolation methods, i.e., the DW and VNN. The discussion and conclusion for the adaptive interpolation method are presented in Section 4. 2. Methods 2.1. Acquisition system The freehand system was comprised of an electromagnetic spatial sensing device (miniBird, Ascension Technology Corporation, Burlington, VT, USA), a portable US scanner (SonoSite 180PLUS, SonoSite, Inc., Bothell, WA, USA), and a PC with a 3.0 GHz Pentium IV processor and 1 GB of RAM. A custom-designed software system programmed in VC++ was designed for data collection, volume reconstruction, visualization and data analysis. The receiver of the spatial sensing device was attached to a 7.5 MHz linear array probe. The video stream of the US scanner was digitized by a video capture card (NI-IMAQ PCI/PXI-1411, National Instruments Corporation, Austin, TX, USA) and then transferred to the PC. The spatial data, including the position and orientation of the spatial receiver, were simultaneously recorded by the freehand system through RS232 serial port in experiments. After a temporal calibration procedure [18], the B-scans with a pixel resolution of 0.08 mm were matched to the spatial data points at a rate of 20 frames per second. Spatial calibration using a cross-wire phantom [14,18] was performed to determine the location of the receiver with respect to the probe. 2.2. The measurement of speckle regions Loupas et al. [21] proposed the following signal-dependent noise model for speckle specifications in US images: pffiffi y ¼ s þ n s; ð1Þ where, y, s and n represent the observed US signal, noise-free signal, and noise, respectively. In a homogeneous region, it can be assumed that s is a constant. Therefore, the variance (r2) of the region has the following linear relationship to that of noise ðr2n Þ: r2 ¼ sr2n :

ð2Þ

Assume arithmetic mean l of the region is the expectation of s, Eq. (2) can be rewritten as r2 ¼ lr2n . If there are only speckle noises in a homogeneous region, Eq. (2) shows that the variance is proportional to the mean. This model implies that the relationship between local mean and variance can be used as a criterion to determine a homogeneous region. In this paper, we studied local statistics to find out the criterion of identifying a homogeneous region. The measurements on conventional 2D US images captured from a US resolution phantom (Model 44, CIRS Inc., USA) and a healthy male subject’s forearm were performed to determine the homogeneity of speckle regions. The local statistics including mean and variance at different locations in homogeneous speckle regions were measured. As shown in Fig. 1, the local mean and variance were calculated within typical speckle regions indicated by the rectangles. Karaman et al. [23] demonstrated that the measurement in a window size larger than 11  11 pixels would be good at approximating speckle statistics in homogeneous regions. In this study, we chose the window size of 15  15 pixels for the measurement of local statistics. Thirty-five measurements at different locations were conducted

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

23

Fig. 1. Two typical US images collected from the resolution phantom (a) and the forearm of the subject (b). The rectangles indicate homogeneous regions, where only speckles exist.

Fig. 2. Linear correlation of the speckle statistics for both of the US resolution phantom and the subject’s forearm. (a) The linear relationship between the local variance and mean for the phantom images, and (b) the linear relationship between the local variance and mean for the forearm images.

for both the phantom and forearm images with different brightness. Fig. 2 illustrates the measurement results. The asymptotic values of the ratio of variance/mean were 4.23 and 2.67, denoted as Hp and Hf in this paper, for the phantom and forearm images, respectively. They were used as the homogeneity thresholds for the volume reconstruction of the phantom and a part of the forearm accordingly in the current study. 2.3. Volume reconstruction A region of interest was scanned in a single sweep with a slow and steady motion in each experiment. Fig. 3 illustrates the outline of a sequence of B-scans in a typical sweep. Similar to the DW interpolation, our algorithm proceeded voxel by voxel. A spherical neighbourhood of radius R, centred about each voxel, was predefined. All pixels transformed into this neighbourhood and their

distances to the voxel centre were stored in a dynamic array in association with the current voxel. However, if there were a large number of raw B-scans or the voxel size was set as small as the pixel size, the memory of the PC would be quickly overwhelmed and the failure of reconstruction would result. To solve this problem, we used a voxel size relatively larger than the pixel size and reconstructed the whole voxel array slice by slice. Specifically, the voxel array was considered as a set of 2D slices along z-axis, i.e., the scanning direction in this paper. After Ns B-scans that could fully cover a slice in the voxel array were transformed into the volume coordinate system, the reconstruction of this slice began; and when the reconstruction was finished, the memory allocated for the slice was released for the reconstruction of the next slice. Thus, the memory of the PC was economically utilized though the computation time was prolonged. Ns should be chosen large enough to contain all spherical neighbourhoods of the voxels on any slices. In this study, Ns could be automatically determined by the radius of the spherical neighbourhood because the B-scan spacing was roughly even in a single sweep. The following interpolation algorithm was proposed to compute the voxel values. Once the data in a voxel’s dynamic array were ready, the local statistics of its neighbouring pixels were computed. If the ratio of variance/mean was larger than the homogeneity threshold, the voxel was regarded as locating in a structurally inhomogeneous region. Otherwise, it would be treated as locating in a homogeneous region. For a voxel locating in an inhomogeneous region, its neighbouring pixels were weighted by the inverse distances to its centre. These distances had a power adaptively determined by the ratio of local variance/mean. The weighted average of all neighbouring pixels was the interpolation output. In a previous study, Huang et al. [18] had revealed that higher power for the inverse distances in the DW method resulted in 3D images with more details and texture patterns. Based on those results, we suggested that the power of the distances was positively proportional to the ratio of local variance/mean in this paper. The proposed new method implied that the voxels in an inhomogeneous region would have a relatively large power of distances for inversely weighting their neighbouring pixels, which would enhance the contrast and well preserve small resolvable objects and edges. As the power of inverse distance for weighing pixels was determined with respect to the local statistics, this method was named as ADW interpolation in this study. For a voxel in homogeneous regions, an arithmetic mean filter was able to effectively reduce speckles [22,23]. Chen et al. [24] adopted a trimmed mean filter to achieve better speckle suppression results other than the regular mean filter. In their method, only the pixels with intensities within one standard deviation of the mean of the original distribution

24

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

Fig. 3. A typical single sweep of raw B-scans. Ns denotes the number of B-scans used for reconstructing one slice within the voxel array. Ns must be large enough to make the B-scans cover the range of all neighbourhoods for the voxels on the slice.

were averaged as the new output. In the current ADW method, the same trimmed arithmetic mean filter was implemented to the neighbouring pixels for speckle suppression in homogeneous regions. In summary, the ADW interpolation algorithm was formulated as follows: 8 > > > 0; < W k ¼ 1; > > > : 1a ; dk  2  2 r r > H0 ; a ¼ b   H0 þ 1; if l l

Pn k¼0 W k Iðpk Þ ; IðV j Þ ¼ P n k¼0 W k

r2 l

6 H0 \ I; ðpk Þ R ½l  r; l þ r;

r2 l

6 H0 \ I; ðpk Þ 2 ½l  r; l þ r;

r2 l

> H0 ;

ð3Þ where I(Vj) denotes the intensity of the jth voxel in the voxel array, n is the number of pixels falling into the predefined spherical region centred about Vj, I(pk) denotes the intensity of the kth neighbouring pixel of Vj, Wk is the relative weight for the kth pixel, dk is the distance of the kth pixel to the centre of the voxel Vj, l and r are the mean and standard deviation of all neighbouring pixels, respectively, H0 is the threshold value of local homogeneity, if r2 =l > H0 , the power a of the kth pixel’s distance to the centre of voxel Vj is determined with respect to the ratio of r2 =l, and b is a parameter empirically set by the operator for adjusting the interpolation performance. 2.4. In vitro experiment The US resolution phantom was scanned in the in vitro experiment. A set of 128 nearly parallel B-scans was captured in a single sweep. Each B-scan was cropped to contain 500  388 pixels. A voxel array with 150  116  120 voxels was reconstructed using the ADW, DW and VNN methods, respectively. For the ADW and DW methods, a spherical region with a radius of 0.4 mm was predefined for each voxel. The parameter b in the ADW method was set to be 0.25, 0.5, 0.75 and 1.0, and their results were compared. The homogeneity threshold H0 for the resolution phantom was set to 4.23 (Hp) as measured above in Section 2.2. Totally, 6 voxel arrays were reconstructed in the in vitro study. Four representative slices were chosen from the same location in the voxel arrays reconstructed using the ADW (b = 0.5), ADW (b = 1.0), DW and VNN methods. To compare the performance in speckle suppression and contrast enhancement among the three

methods, the intensities along row 30 on these slices were extracted. As shown in Fig. 4, row 30 passed through three round targets so that the contrast at the target boundaries and the homogeneous speckle regions reconstructed using the ADW, DW and VNN methods could be directly compared. In addition, the signal to noise-ratio (SNR) of the 3D homogeneous regions reconstructed using different methods was studied. A sub-volume with 30  30  30 voxels was selected at the same location but including no round target from each of the 6 voxel arrays, respectively, and their SNRs were computed for comparing the performances in reducing speckles among the three reconstruction methods. Larger SNR values indicated better performance in speckle reduction. Local contrast measure [26] was also performed to evaluate the contrast enhancement using the different methods. For a voxel V(i, j, k) at the coordinates (i, j, k) in a voxel array, the local contrast in its 2n + 1 neighbourhood can be expressed as Cði; j; kÞ ¼

MaxðIs ði; j; kÞÞ  MinðIs ði; j; kÞÞ ; MaxðIs ði; j; kÞÞ þ MinðIs ði; j; kÞÞ

ð4Þ

where Is(i, j, k) denotes a group of voxel intensities at the location (i, j, k) and the voxel’s neighbourhood. We used the averaged local contrast in a 3D inhomogeneous region to represent the performance of contrast enhancement. The method is described by P ði;j;kÞ2S Cði; j; kÞ ; ð5Þ CA ¼ N where S is a sub-volume cropped from a reconstructed volume and N is the voxel number within the sub-volume. In this phantom study, another sub-volume with 16  16  25 voxels mainly including a typical round target was selected at the same location from the 6 voxel arrays. The averaged local contrast measures were applied to the sub-volumes. It is obvious that sharper edges, more texture objects and details from B-scans should be preserved in a subvolume with a larger CA. 2.5. In vivo examination An in vivo examination was performed on the healthy male subject’s forearm. A dense set of 206 nearly parallel B-scans was captured and each B-scan was cropped to 500  430 pixels. The dimension of the voxel array was 240  212  145 voxels. Four voxel arrays using the ADW (b = 0.5), ADW (b = 1.0), DW and

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

25

Fig. 4. Four slices reconstructed for the resolution phantom using the: (a) DW, (b) VNN, (c) ADW (b = 0.5) and (d) ADW (b = 1.0) methods, with lines indicating the row 30 in the image.

VNN methods, respectively, were reconstructed. The radius of predefined neighbourhood for each voxel was 0.5 mm. The homogeneity threshold (H0) for this in vivo experiment was 2.67 (Hf) as measured above in Section 2.2. Similarly, four slices at the same location in the voxel arrays were extracted for comparison purpose. The intensities along column 90 on the slices were compared to evaluate the performance of speckle reduction and contrast enhancement of the proposed method. In addition, a quantitative evaluation test previously introduced in [9] was conducted in this in vivo study. The idea of the test was to evaluate the ability of a reconstruction technique in preserving true intensity values at the locations where a part of original data was removed. A good reconstruction method should be able to interpolate the removed grid points with values very close to the original data. A B-scan near the middle of the raw data set was selected for pixel removing. Different percentages of pixels were removed randomly from this selected B-scan. The remained data were used to reconstruct a voxel array with a voxel size equivalent to the pixel size. The average of the absolute differences between the interpolated grid points and the missing original pixel values was calculated for evaluating the reconstruction performance using the following equation: V¼

N 1 X jp  ri j; N i¼1 i

ð6Þ

where pi is the removed original pixel intensity, ri is the interpolated intensity at the location of pi and N is the number of removed pixels. A smaller V indicates a better performance of interpolation. Seven different data removing ratios were used in the current evaluation tests, i.e., 25%, 50%, 75%, 100%, 300%, 500% and 700%. The tests with the data removing ratios of 25%, 50%, 75% and 100% were performed by removing corresponding percentages of data from the selected B-scan n. For the 300% test, the data from the B-scan n ± 1 and B-scan n were removed. The 500% and 700% tests further removed all data from the B-scans n ± 2 and n ± 3, respectively. In this in vivo study, the evaluation test was performed to the 4 voxel

arrays and the interpolation errors for 10 randomly selected B-scans were averaged to compare the performance of the proposed method (b = 0.5) with that of the other two traditional methods, i.e., the DW and VNN. 3. Results 3.1. In vitro experiment Fig. 4 demonstrates a set of reconstructed slices of the US resolution phantom using the ADW, DW and VNN methods, respectively. It shows that the ADW method reduced most of the speckles in homogeneous regions and effectively preserved the edges of the round targets (Fig. 4c and d). In comparison, the DW method blurred the target edges (Fig. 4a) and the VNN seemed to overemphasize the original texture patterns of US images (Fig. 4b). Fig. 5 illustrates the intensities along row 30 on the slices in Fig. 4. The comparison between the ADW (b = 0.5 and 1.0) and DW methods was shown in Fig. 5a and b. It can be observed that the ADW method was able to enhance the sharp changes of intensity better than the conventional DW method, and offered desired improvement in speckle reduction in homogeneous regions. Meanwhile, compared with the VNN method, the currently proposed ADW method has significantly reduced speckle noises, as shown in Fig. 5c and d. Moreover, it can be seen that larger parameter b in the ADW method enhanced the sharp edges better. Table 1 presents the quantitative results for the SNR and the averaged local contrast measures. From the SNR results, the ADW method with four different values of the parameter b suppressed more speckles than the DW and VNN methods in homogeneous regions. Interestingly, different values of the parameter b resulted in no difference on the SNR measure. The possible explanation might be that for the homogeneous regions in the resolution phantom, the speckle patterns were very similar and the local ratio of variance/mean was mostly below the homogeneity threshold Hp. For the averaged local contrast measure, all the ADW measures with

26

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

Fig. 5. The comparison of intensities along row 30 of the four slices for the resolution phantom: (a) the ADW (b = 0.5) and DW, (b) the ADW (b = 1.0) and DW, (c) the ADW (b = 0.5) and VNN, and (d) the ADW (b = 1.0) and VNN.

Table 1 The quantitative results of the in vitro study for the measures of SNR in a homogeneous region and the averaged local contrast in a target region Metrics

DW

VNN

ADW b = 0.25

b = 0.5

b = 0.75

b = 1.25

SNR Contrast

15.51 0.1453

7.79 0.3056

16.45 0.1461

16.45 0.1519

16.45 0.1568

16.45 0.1608

different values of b produced higher contrast than the DW method in the target region, indicating a good performance in contrast enhancement, with larger value of b exhibiting better performance. Although the VNN method produced a much larger value of local contrast measure, it preserved all speckle noises at the same time with the lowest SNR value, hence showing the worst image quality. 3.2. In vivo examination Fig. 6 shows four typical slices reconstructed using the ADW (b = 0.5 and 1.0), DW and VNN methods, respectively. Fig. 6e–h illustrates partial image content enlarged from the original slices. Fig. 7 demonstrates the quantitative comparison between the

ADW and the other two methods by extracting the intensities along column 90. It can be observed that the DW method smoothed number of small resolvable objects as well as significant edges, while the proposed ADW method produced a good contrast for these small anatomical details and well preserved tissue edges, making the boundaries easier to be traced. However, as the DW method is traditionally regarded being good at speckle suppression [14], the ADW method did not present significant improvement in reducing speckle noises in low contrast regions, as shown in Fig. 7a and b. In comparison with the VNN method, the ADW method did much better in reducing speckle and other noises (see Fig. 7c and d). Although the image results using the VNN method seemed to be able to preserve most of texture patterns as shown in Fig. 6b and f, the noises were at the same time retained making the tissue boundaries difficult to be contoured. In addition, from Figs. 6 and 7, it can be indicated that larger value of b in the ADW method was better at contrast enhancement. The results for the evaluation test for reconstructing missing data points are summarized in Fig. 8. The VNN method generated the largest interpolation error though it showed the most texture patterns. Due to the greatest level of misalignments for relatively large gaps between B-scans, the VNN method generally has greater

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

27

Fig. 6. Four slices reconstructed using the DW, VNN and ADW (b = 0.5 and 1.0) methods for the subject’s forearm: (a) the DW, (b) VNN, (c) ADW (b = 0.5) and (d) ADW (b = 1.0) methods, with lines indicating the column 90. Images (e)–(h) illustrate the enlarged images of the rectangle indicated in image (a) from images (a)–(d), respectively.

28

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

Fig. 7. The comparison of intensities along column 90 of the four slices for the subject’s forearm: (a) the ADW (b = 0.5) and DW, (b) the ADW (b = 1.0) and DW, (c) the ADW (b = 0.5) and VNN, and (d) the ADW (b = 1.0) and VNN.

reconstruction error than the other methods [9,17]. The DW method produced improved interpolation result compared with the VNN. The ADW method (b = 0.5) further reduced the interpolation error, illustrating the best interpolation performance. This might result from the better ability of the ADW to effectively preserve edges and small objects in comparison with the DW method. 4. Discussion

Fig. 8. The comparison of the interpolation errors among the results of DW, VNN and ADW (b = 0.5) methods in the evaluation test for reconstructing missing data points.

In this paper, an adaptive DW method was proposed and its performance in speckle suppression and contrast enhancement was investigated for volume reconstruction of 3D freehand US. This method interpolated voxels using local statistics information. If a voxel was locating in a homogeneous region, the method made use of an arithmetic mean filter to reduce speckle noises. Otherwise, an adaptive power determination method was applied to the inverse pixel distances to the voxel centre in conventional weighted average interpolation method to better preserve edges and enhance local contrast. This method was successfully employed in current study for reconstructing a US resolution phantom and a subject’s forearm with a tracked 3D freehand system. The comparison results between this new method and the traditional DW and VNN methods showed that the proposed method was able

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

to effectively suppress speckle noises in homogeneous regions, and meanwhile preserve the contrast and edges in inhomogeneous regions and reduce interpolation errors. In order to find out the homogeneity of speckle regions in both the phantom and forearm images, the measurements of local statistics were performed in this study. The two thresholds (Hp and Hf) were critical for judging the homogeneity of the local neighbourhood of a voxel grid point. Generally, a neighbourhood was considered a homogeneous region if the ratio of variance/mean was below the threshold. Hence, a trimmed arithmetic mean filter was applied to sufficiently smooth speckles. As the speckle statistics may depend on the scanner specifications, the homogeneity threshold could be varied for different US scanners. The imaged material or tissue itself might also affect the speckle statistics; therefore the homogeneity threshold should be appropriately determined through enough tests. In addition, because the selection of homogeneous regions is manually conducted in this study, it is inconvenient for users and may result in incorrect selection. We will investigate a semi-automatic procedure that can globally analyse local statistics in an ultrasound image and classify the local regions into three groups with high, medium and low ratios of variance/mean, respectively. The regions belonging to the group of low ratio of variance/mean will be highlighted in the image. Then the user can make a decision which regions can be used as homogeneous regions for determining the threshold. A follow-up study could be the adaptive selection of the threshold calculated using the selected regions to optimize the reconstruction result. It is well known that the well-preserved edges, small point targets and the content of relatively low contrast lesions are crucial for clinical diagnosis. In the ADW method, the well-enhanced contrast in inhomogeneous regions made these image details more obvious. The adaptation strategy adopted in this study resulted in a proportional enhancement of contrast based on local statistics, thus a more accurate detection of small resolvable structures was realized. The parameter b could also have an effect on the interpolation performance. Theoretically, an extremely large value of b would lead to similar interpolation performance to the VNN method for inhomogeneous regions, since the pixel with the nearest distance to the voxel centre tended to carry the largest weight. According to the in vitro and in vivo studies, the parameter b was assigned to be from 0.25 to 1.0, and the different values worked well for contrast enhancement as demonstrated in the section of results. In addition, effective edge preservation was very important for the segmentation of anatomical structures, such as surface extraction [10–13], as the increased gradient at edge voxels could make an easier and more accurate tracing of tissue boundaries. Consequently, more accurate 3D measurements, such as volume estimation [3,10], can be realized for 3D objects. Furthermore, the improvement of edge preservation would be useful for the registration of multiple sweeps [4,27,28], due to that the well-preserved edges could be considered as landmarks and used for alignment of different sweeps. Moreover, the proposed adaptive method may also be potentially used to enhance the volume image obtained using the 3D ultrasound probe. However, as the local statistics of neighbouring pixels of each voxel should be measured prior to the computation of output value, as much as three times of the computation time of the DW was usually required using the ADW in above experiments. Fortunately, with the high performance of CPU in our system, the reconstructions were normally completed within 10 min. With further progress in computer technology, the time for volume reconstruction using the ADW method can be anticipated to be significantly reduced in the near future. In addition, speckle statistics of imaged human tissues should be studied prior to real examinations in this study. It would cause inconvenience for different clinical applications. Thus, one of our future tasks will be performed to investigate

29

the speckle statistics for different human tissues and provide a table which can be used for selecting appropriate homogeneity threshold for a specific examination. Moreover, the effect of the parameter b on different human tissues will be further quantitatively studied for optimizing the reconstruction results in various examinations. The proposed algorithm can also be applied to 3D volume data sets. One would ask why not perform the speckle suppression after the volume reconstruction has been completed. It is well known that the resolutions for voxels are normally much lower than those of original B-scans in freehand 3D ultrasound imaging. As aforementioned in this paper, the speckle suppression should fully make use of the statistics of local pixels. If we perform the suppression of speckles and preservation of small details when the 3D data have been formed, the proposed method has to deal with the data with lower resolutions and hence big loss of original information. Thus, the effect of post-processing of the 3D data may not be as good as that of direct processing of original 2D data. In summary, our recent work on the adaptive volume reconstruction of 3D US based on local homogeneity have been presented in this paper. The new method can adaptively select two reconstruction algorithms to interpolate a voxel array with suppressed speckle noises and enhanced contrast of edges and small details. According to the improved interpolation performance, we expect that this new method has potentials to be used in the reconstruction of 3D ultrasound images for both research and routine clinical practices. Acknowledgements This work was partially supported by The Hong Kong Polytechnic University (G-YE22) and the Research Grants Council of Hong Kong (PolyU 5245/03E). The authors thank the anonymous reviewers for their constructive comments in improving this paper. References [1] Nelson TR, Pretorius DH. Three-dimensional ultrasound imaging. Ultrasound Med Biol 1998;24:1243–70. [2] Fenster A, Downey DB, Cardinal HN. Three-dimensional ultrasound imaging. Phys Med Biol 2001;46(5):R67–99. [3] Gee AH, Prager RW, Treece GM, Berman L. Engineering a freehand 3D ultrasound system. Pattern Recogn Lett 2003;24(4–5):757–77. [4] Gee AH, Treece GM, Prager RW, Cash CJC, Berman L. Rapid registration for wide field of view freehand three-dimensional ultrasound. IEEE Trans Med Imag 2003;22(11):1344–57. [5] Mercier L, Langø T, Lindseth F, Collins LD. A review of calibration techniques for freehand 3-D ultrasound systems. Ultrasound Med Biol 2005;31(2):143–65. [6] Prager RW, Rohling RN, Gee AH, Berman L. Rapid calibration of 3-D freehand ultrasound. Ultrasound Med Biol 1998;24(6):855–69. [7] Treece GM, Gee AH, Prager RW, Cash CJC, Berman LH. High-definition freehand 3-D ultrasound. Ultrasound Med Biol 2003;29(4):529–46. [8] Huang QH, Zheng YP, Li R, Lu MH. 3-D measurement of body tissues based on ultrasound images with 3-D spatial information. Ultrasound Med Biol 2005;31(12):1607–15. [9] Rohling RN, Gee AH, Berman L. A comparison of freehand three-dimensional ultrasound reconstruction techniques. Med Image Anal 1999;3(4):339–59. [10] Treece GM, Prager RW, Gee AH, Berman L. Fast surface and volume estimation from non-parallel cross-sections, for freehand 3-D ultrasound. Med Image Anal 1999;3(2):141–73. [11] Treece GM, Prager RW, Gee AH, Berman L. Surface interpolation from sparse cross sections using region correspondence. IEEE Trans Med Imag 2000;19(11):1106–14. [12] Prager PW, Gee AH, Berman L. Stradx: real-time acquisition and visualization of freehand three-dimensional ultrasound. Med Image Anal 1999;3(2):129–40. [13] Zhang WY, Rohling RN, Pai DK. Surface extraction with a three-dimensional freehand ultrasound system. Ultrasound Med Biol 2004;30(11):1461–73. [14] Barry CD, Allot CP, John NW, Mellor PM, Arundel PA, Thomson DS, et al. Threedimensional freehand ultrasound: image reconstruction and volume analysis. Ultrasound Med Biol 1997;23(8):1209–24. [15] Meairs S, Beyer J, Hennerici M. Reconstruction and visualization of irregularly sampled three- and four-dimensional ultrasound data for cerebrovascular applications. Ultrasound Med Biol 2000;26(2):263–72.

30

Q. Huang et al. / Applied Acoustics 70 (2009) 21–30

[16] Sanches JM, Marques JS. A Rayleigh reconstruction/interpolation algorithm for 3D ultrasound. Pattern Recogn Lett 2000;21(10):917–26. [17] San José-Estépar R, Martı´n-Fernándes M, Caballero-Martı´nes PP, AlberolaLopes C, Ruiz-Alzola J. A theoretical framework to three-dimensional ultrasound reconstruction from irregularly sampled data. Ultrasound Med Biol 2003;29(2):255–69. [18] Huang QH, Zheng YP, Lu MH, Chi ZR. Development of a portable 3D ultrasound imaging system for musculoskeletal tissues. Ultrasonics 2005;43(3):153–63. [19] Huang QH, Zheng YP. An adaptive squared-distance-weighted interpolation for volume reconstruction in 3D freehand ultrasound. Ultrasonics 2006;44:E73–7. [20] Rohling RN, Gee AH, Berman L. Three-dimensional spatial compounding of ultrasound images. Med Image Anal 1997;1(3):177–93. [21] Loupas T, McDicken W, Allan P. An adaptive weighted median filter for speckle suppression in medical ultrasonic images. IEEE Trans Circuits Syst 1989;36(1):129–35.

[22] Koo JI, Park SB. Speckle reduction with edge preservation in medical ultrasonic images using a homogeneous region growing mean filter (HRGMF). Ultrason Imag 1991;13(3):211–37. [23] Karaman M, Kutay MA, Bozdagi G. An adaptive speckle suppression filter for medical ultrasonic imaging. IEEE Trans Med Imag 1995;14(2):283–92. [24] Chen Y, Yin RM, Flynn P, Broschat S. Aggressive region growing for speckle reduction in ultrasound images. Pattern Recogn Lett 2003;24(4–5):677–91. [25] Huang QH, Zheng YP. Volume reconstruction of freehand three-dimensional ultrasound using median filters. Ultrasonics 2008, in press. [26] Sun QL, Hossack JA, Tang JS, Acton ST. Speckle reducing anisotropic diffusion for 3D ultrasound images. Comput Med Imag Grap 2004;28:461–70. [27] Rohling RN, Gee AH, Berman L. Automatic registration of 3-D ultrasound images. Ultrasound Med Biol 1998;24:841–54. [28] Huang QH, Zheng YP. A new scanning approach for limb extremities using a water bag in freehand 3-D ultrasound. Ultrasound Med Biol 2005;31: 575–83.