Reconstruction of 3D ultrasound images based on Cyclic Regularized Savitzky–Golay filters

Ultrasonics 51 (2011) 136–147

Contents lists available at ScienceDirect

Ultrasonics journal homepage: www.elsevier.com/locate/ultras

Reconstruction of 3D ultrasound images based on Cyclic Regularized Savitzky–Golay filters Pollakrit Toonkum a, Nijasri C. Suwanwela b, Chedsada Chinrungrueng a,* a b

Department of Electrical Engineering, Faculty of Engineering, Chulalongkorn University, Bangkok 10330, Thailand Division of Neurology, Faculty of Medicine, Chulalongkorn University, Bangkok 10330, Thailand

a r t i c l e

i n f o

Article history: Received 7 May 2010 Received in revised form 9 July 2010 Accepted 9 July 2010 Available online 25 July 2010 Keywords: Savitzky–Golay filters Regularization Speckle reduction Reconstruction 3D ultrasound images

a b s t r a c t This paper presents a new three-dimensional (3D) ultrasound reconstruction algorithm for generation of 3D images from a series of two-dimensional (2D) B-scans acquired in the mechanical linear scanning framework. Unlike most existing 3D ultrasound reconstruction algorithms, which have been developed and evaluated in the freehand scanning framework, the new algorithm has been designed to capitalize the regularity pattern of the mechanical linear scanning, where all the B-scan slices are precisely parallel and evenly spaced. The new reconstruction algorithm, referred to as the Cyclic Regularized Savitzky– Golay (CRSG) filter, is a new variant of the Savitzky–Golay (SG) smoothing filter. The CRSG filter has been improved upon the original SG filter in two respects: First, the cyclic indicator function has been incorporated into the least square cost function to enable the CRSG filter to approximate nonuniformly spaced data of the unobserved image intensities contained in unfilled voxels and reduce speckle noise of the observed image intensities contained in filled voxels. Second, the regularization function has been augmented to the least squares cost function as a mechanism to balance between the degree of speckle reduction and the degree of detail preservation. The CRSG filter has been evaluated and compared with the Voxel Nearest-Neighbor (VNN) interpolation post-processed by the Adaptive Speckle Reduction (ASR) filter, the VNN interpolation post-processed by the Adaptive Weighted Median (AWM) filter, the Distance-Weighted (DW) interpolation, and the Adaptive Distance-Weighted (ADW) interpolation, on reconstructing a synthetic 3D spherical image and a clinical 3D carotid artery bifurcation in the mechanical linear scanning framework. This preliminary evaluation indicates that the CRSG filter is more effective in both speckle reduction and geometric reconstruction of 3D ultrasound images than the other methods. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Ultrasound imaging technique has been widely used in medical diagnosis because of its non-invasive, non-ionizing, real-time, portable and low-cost nature. Conventional two-dimensional (2D) ultrasound imaging is usually performed with hand-held probe which transmits ultrasonic pulses into the body and receives the echoes. These echoes are then used to create a grey-level image (B-scan) of a cross section of the body in the scan plane. However, the 2D nature of the B-scan image imposes some limitations to physician in visualizing, examining and assessing the progression of disease. Three-dimensional (3D) ultrasound imaging has been developed to overcome such limitations. One approach for 3D ultrasound reconstruction is to generate 3D ultrasound images from a sequence of 2D B-scans by attaching a position sensor to

* Corresponding author. Tel.: +66 2 218 6906; fax: +66 2 218 6912. E-mail addresses: [email protected] (P. Toonkum), fmednsu @md2.md.chula.ac.th (N.C. Suwanwela), [email protected] (C. Chinrungrueng). 0041-624X/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.ultras.2010.07.003

the freehand probe for labeling position and orientation of each B-scan image [1]. The acquired set of 2D B-scan images are then used to reconstruct a 3D image in a regular voxel (volume element) array by filling each pixel from the acquired 2D B-scans into its correct location in the volume, while the unobserved image intensities associated with unfilled voxels are to be estimated by interpolation or approximation from the intensities in the filled voxels [2]. The typical volume coordinate system is usually based on Cartesian volume (i.e., 3D rectangular grid), as depicted in Fig. 1 [3]. Although the freehand approach allows the user to manipulate the transducer and view the desired anatomical section, the approach requires careful and delicate end-user calibration to ensure the overall accuracy of an ultrasound position tracking system, and its image quality can be heavily degraded from problem of scanning gaps, particularly when imaging small structures at high resolution. Furthermore, the freehand scanning suffers from limitations and complications of the position sensors, i.e., acoustic sensor requires that an operator must maintain a clear line of sight between the probe and the position sensor’s fixed base station and that the corrections must be made for variation of the speed of sound in air due

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

Fig. 1. Typical volume coordinate system in 3D ultrasound reconstruction.

to changes in temperature and humidity, whereas electromagnetic field sensor requires short transmitter–receiver distance and no nearby metal, compelling the operator to be careful not to stray the sensors outside its operating region [4]. As an alternative to the freehand approach, dedicated 3D probes which incorporate an oscillating mechanism to sweep the scan plane is developed. There are various scanning geometries of mechanical 3D probes such as linear scanning, fan scanning and rotational scanning. In contrast to the freehand approach, the 3D scanning geometry of the mechanical 3D probes is known a priori and the transducer is automatically moved to acquire a sequence of 2D B-scan image slices, which helps solve the limitations and complications imposed on the freehand approach. In addition, the dedicated probe permits efficient computation of 3D reconstruction due to the pre-determined pattern of B-scan slices, thus resulting in short reconstruction times [5]. In addition to the pure mathematically approximation of the 3D ultrasound from B-scan slices, there is another problem in 3D ultrasound image reconstruction, which is related to the presence of random speckle noise within each image of 2D B-scans. The speckle noise is a random pattern occurred from interference of the acoustic wave within each image of 2D B-scans [6]. Such random patterns severely degrade the quality of ultrasound images, deviating the observed image intensities from their actual values, resulting in an inaccurate 3D reconstruction and interpretation. Therefore, to improve clinical accuracy, the reconstruction technique should be able to simultaneously approximate the missing data in unfilled voxels and perform speckle noise reduction in filled voxels assigned from sparse 2D B-scan ultrasound images, enabling the reconstruction technique to alleviate the geometric deformations of the human organs occurred during the data acquisition process. A number of algorithms for reconstructing 3D ultrasound images and reducing the speckle noise from sequences of 2D Bscan images have been reported and empirically evaluated in [7]. These algorithms can be grouped into three categories: the Voxel Nearest-Neighbor (VNN) interpolation, the Pixel Nearest-Neighbor (PNN) interpolation, and the Distance-Weighted (DW) interpolation. The most simple and straightforward algorithm to implement is the VNN interpolation. In the VNN interpolation, the 3D voxel array is traversed and each voxel is filled with the value of the nearest pixel contained in the set of B-scan images [8]. Although the method is simple and fast, it tends to produce artifacts on the 3D reconstruction, especially when the elevation scanning has a large sampling interval. For the PNN interpolation, each pixel in all 2D B-scan images is picked and its intensity is assigned to its nearest voxel. For a voxel

137

with multiple pixel assignments, its intensity is the average of all assigned pixel intensities. For a voxel with no assignment, its intensity is the average of all the initially filled voxels in its local neighborhood [9]. For the DW interpolation, all pixels in the local neighborhood centered at each voxel are weighted and averaged by the inverse distances between those pixels and the voxel center [10]. Due to the averaging tendency of the PNN and the DW interpolations, these reconstruction techniques are able to simultaneously approximate the missing data in unobserved voxels and perform speckle noise reduction in observed voxels assigned from sparse 2D B-scan images, providing good performance of 3D reconstruction with speckle and shadowing reduction. However, the employed averaging processes tend to eliminates the high frequencies and thereby smooth out the 3D image boundaries [11]. To preserve boundaries and enhance the contrast of the 3D reconstructed images, the Adaptive Distance-Weighted (ADW) interpolation algorithm has been recently proposed by [12]. The ADW interpolation method first determines for each voxel the mean and the variance of its neighboring pixels. If the mean and the variance associated with a voxel satisfied the homogeneity criterion, the voxel intensity is computed with an arithmetic mean filter. Otherwise, the intensity is interpolated by the DW method employing an adaptive weight determined by the mean and the variance of the neighboring pixels. Even though the ADW has demonstrated better capability in image interpolation, speckle suppression and edges preservation when comparing with the VNN and the DW method, it requires extremely high computation time due to the need to estimate the local mean and variance of each voxel. Nevertheless, all of the aforementioned 3D ultrasound image reconstruction algorithms have been developed and evaluated in the freehand scanning framework. In contrast to the freehand approach and due to the regularity property of the 2D B-scan images obtained from mechanical linear scanning, the acquired 2D B-scan images are precisely parallel to each other and accurately separated by predefined interval. Such regularity property of mechanical linear scanning pattern can be captured to permit efficient computation of the reconstruction algorithms. This paper aims to develop a new 3D ultrasound reconstruction algorithm in mechanical linear scanning framework for capitalizing its regularity property to improve the computational performance, and to evaluate the new algorithm by comparing with the most existing reconstruction algorithms. The new reconstruction algorithm, referred to as the Cyclic Regularized Savitzky–Golay (CRSG) filter, is a new variant of the Savitzky–Golay (SG) filter. Originally, the SG filter is the one-dimensional (1D) smoothing filter initially designed to render visible the relative widths and heights of spectral lines in noisy spectrometric data [13]. In addition to being extended to accept a 3D array of data as the input instead of a series of 1D data, the CRSG filter has been improved upon the original SG filter in two respects: First, the cyclic indicator function has been incorporated into the least square cost function to enable the CRSG filter to approximate nonuniformly spaced data of the unobserved image intensities contained in unfilled voxels and reduce speckle noise of the observed image intensities contained in filled voxels. Second, the regularization function has been augmented to the least squares cost function as a mechanism to balance between the degree of speckle reduction and the degree of detail preservation. In Section 2, we formulate the principle of the Cyclic Savitzky–Golay (CSG) filter used for reconstructing of 3D ultrasound images of the nonuniformly spaced B-scan data acquired from mechanical linear scanning. In Section 3, we develop the CRSG filter by incorporating the regularization function into the CSG filter. In Section 4, the reconstruction performances of the CRSG filter is evaluated and compared with those of the VNN interpolation post-processed by the Adaptive Speckle Reduction

138

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

(ASR) filter [14], the VNN interpolation post-processed by the Adaptive Weighted Median (AWM) filter [15], the DW interpolation, the ADW interpolation and the CSG reconstruction filter on both synthetic and clinical ultrasound data. Finally, the discussion and conclusions are presented in Section 5. 2. Formulation of the Cyclic Savitzky–Golay filters for 3D ultrasound reconstruction in the mechanical linear scanning framework 2.1. Volume coordinate system for the 3D mechanical linear ultrasound scanning framework In a mechanical linear scanning, a sequence of 2D B-scan images is obtained by mechanical movement of the transducer in a precise and predefined manner, corresponding to parallel cross sections of the human organ and perpendicular to the B-scan image plane. Fig. 2 depicts a volume coordinate system for a 3D voxel array, constructed from a sequence of 2D B-scan images acquired in mechanical linear scanning. The resolution within a 2D B-scan image, referred to as the in-plane resolution, is determined by the pulse bandwidth and transducer aperture, whereas the resolution in the direction perpendicular to the image slice, referred to as the elevation resolution, is determined by the transducer thickness and the elevation sampling interval [16]. In general, the in-plane resolution is much higher than the elevation resolution due to the transducer thickness and large elevation sampling interval. Planes of voxels corresponding to the acquired B-scan slices are shown as shaded voxels in Fig. 2. Notice that only one in every N scan planes in the 3D voxel array are filled with observed B-scan intensities. To reconstruct a 3D ultrasound image from a series of 2D B-scans obtained from such mechanical linear scanning approach, it is thus necessary that we interpolate or approximate the unobserved image intensities of void voxels which are sandwiched between intensity-filled voxels.

noise-free value ~f ðt i Þ for each ti, where ti is set to be t0 + Mi for some constant sample spacing M and i = . . ., 1, 0, 1, . . . . Notice that every data point to be estimated in the 1D SG filter needs to uniformly and equally spaced. For 3D ultrasound image reconstruction of data acquired from mechanical linear scanning, the relative distances between the center of neighborhood and the nearby filled voxel intensities are nonuniformly spaced along the elevation direction. To modify the original 1D SG filter for reconstructing 3D ultrasound images from B-scan data acquired from mechanical linear scanning, such 1D formulation has been improved in two respects. First, it has been extended to accept a 3D array of data as the filter input instead of a series of 1D data. Second, it has been modified by incorporating the indicator function to deal with nonuniformly elevationspaced data of the existing image intensities contained in neighborhood and to approximate the filter output at the fraction of sample spacing M other than at only the multiple of M. Consider a 3D voxel array t indexed by (i, j, k), where index i is defined over . . . , 1, 0, 1, . . . , index j over . . ., 1, 0, 1, . . . , and index k over . . . , 1, 0, 1, . . . . Let . . . , f1, f0, f1, . . . be a given sequence of observed 2D ultrasound B-scan images, where fl denotes a function representing the image intensity of the lth B-scan image, for l defined over . . ., 1, 0, 1, . . . . Define function f on voxel array t, representing the voxel intensities mapped from the observed sequence of B-scan intensity function fl. The intensity f at voxel tði; j; lN scan Þ is set to the value of fl(i, j), and the intensities f at other voxels to arbitrary value, such as 0. The filled voxels (i.e., voxels filled with the observed B-scan intensities fl) are portrayed by shaded voxels as illustrated in Fig. 2. In addition to function f, define function ~f over voxel array t representing the 3D ultrasound image which we aim to reconstruct from the observed B-scan intensities fl. Define voxel neighborhood Dfi;j;kg as a set of voxels centered at position (i, j, k) and be of the form:

Dfi;j;kg ¼ ftði þ m; j þ n; k þ oÞ;

m ¼ fM : Mg;

n ¼ fN : Ng; o ¼ fO : Ogg: 2.2. Formulation of the Cyclic Savitzky–Golay filters In this section, we formulate a new variant of the Savitzky–Golay (SG) filter, referred to as the Cyclic Savitzky–Golay (CSG) filter, for 3D ultrasound image reconstruction. According to the original 1D SG formulation, initially used to render visible the relative widths and heights of spectral lines in noisy spectrometric data [13], the filter performs a local 1D polynomial least squares fitting on a series of uniformly spaced data values f(ti), to approximate the

ð1Þ

Over neighborhood Dfi;j;kg , define function g{i, j, k} as a 3D polynomial of degree P in m, degree Q in n, and degree R in o; namely

g fi;j;kg ðm; n; oÞ ¼

Q X P X R X

afi;j;kg;fp;q;rg mp nq or :

ð2Þ

p¼0 q¼0 r¼0

To express Eq. (2) above in a vector format, we define index s denotes an integer ranging from 1 to (P + 1)(Q + 1)(R + 1), and let p,q, and r be indexing functions of the following forms:

pðsÞ ¼ bðs  1Þ=ððQ þ 1ÞðR þ 1ÞÞc;

ð3Þ

qðsÞ ¼ bðs  1Þ=ðR þ 1Þc mod ðQ þ 1Þ;

ð4Þ

rðsÞ ¼ ðs  1Þ mod ðR þ 1Þ;

ð5Þ

where bc denotes the floor function and mod the modulo function. According to these indexing functions, we arrange all a{i,j,k},{p,q,r}’s and mpnqor’s involved in (2) as vectors:

~ afi;j;kg ¼ ðafi;j;kg;fpðsÞ;qðsÞ;rðsÞg ;

s ¼ 1; . . . ; ðP þ 1ÞðQ þ 1ÞðR þ 1ÞÞT ;

ð6Þ

and

~ bfi;j;kg ðm; n; oÞ ¼ ðmpðsÞ nqðsÞ orðsÞ ;

s ¼ 1; . . . ;

T

ðP þ 1ÞðQ þ 1ÞðR þ 1ÞÞ :

ð7Þ

Such vectors ~ afi;j;kg and ~ bfi;j;kg thus allows us to correspondingly re-express polynomial g{i,j,k} defined in (2) compactly as: Fig. 2. Volume coordinate system in mechanical linear scanning.

afi;j;kg g: g fi;j;kg ðm; n; oÞ ¼ f~ bfi;j;kg ðm; n; oÞgT~

ð8Þ

139

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

Let ~f ði; j; kÞ denotes the output of the CSG reconstruction filter at voxel t(i, j, k). The CSG reconstruction filter computes the output ~f ði; j; kÞ by least squares fitting polynomial g {i,j,k} to the intensity f of all filled voxels contained in neighborhood Dfi;j;kg , and then setting output ~f ði; j; kÞ to g{i, j, k}(0, 0, 0). Since g{i,j,k}(0, 0, 0) is equal to a{i,j,k},{0,0,0} as defined in (2), computing output ~f ði; j; kÞ is equivalent to deriving coefficient a{i,j,k},{0,0,0} of polynomial g{i,j,k} over neighborhood Dfi;j;kg . That is,

~f ði; j; kÞ ¼ g fi;j;kg ð0; 0; 0Þ ¼ afi;j;kg;f0;0;0g :

ð9Þ

To determine solution a{i,j,k},{0,0,0} associated with the aforementioned polynomial least squares fitting problem, we minimize the following objective function:

 ð~afi;j;kg Þ ¼

X

Ifi;j;kg ðm; n; oÞ fefi;j;kg ðm; n; oÞg2 ;

ð10Þ

ðm;n;oÞ 2 Dfi;j;kg

Ifi;j;kg ðm; n; oÞ ¼

1; if 0;

/fi;j;kg ðoÞ ¼ 0;

otherwise;

t ¼ 1; . . . ;

T

ð2M þ 1Þð2N þ 1Þð2O þ 1ÞÞ :

ð19Þ

Let A be a (2M + 1)(2N + 1)(2O + 1)  (P + 1)(Q + 1)(R + 1) design matrix, of which entry (ts) is equal to

Ats ¼ fmðtÞgpðsÞ fnðtÞgqðsÞ foðtÞgrðsÞ :

ð20Þ

This design matrix allows us to further express ~ g fi;j;kg as

~ afi;j;kg : g fi;j;kg ¼ A~

ð21Þ

Using Eqs. (18), (19) and (21), we then rewrite objective function  in the matrix quadratic form as

ð~afi;j;kg Þ ¼

where e{i,j,k}(m, n, o) is defined as the difference between g{i,j,k}(m, n, o) and f(i + m,j + n,k + o). We have incorporated into objective function (10), indicator function Ifi;j;kg for assigning the weight of one to the squared error {e{i,j,k} (m, n, o)}2 associated with the filled voxels, and the weight of zero to those associated with the unfilled voxels. Mathematically, we choose to express above indicator function Ifi;j;kg as:



~ afi;j;kg ; bfi;j;kg ðmðtÞ; nðtÞ; oðtÞÞgT~ g fi;j;kg ¼ ðf~

 T   ~ efi;j;kg WrðkÞ ~ efi;j;kg ;

ð22Þ

where ~ efi;j;kg is defined to be ð~ f fi;j;kg  A~ afi;j;kg Þ, and Wr(k) to be a diagonal matrix of the form:

WrðkÞ ¼ diagfIrðkÞ ðoðtÞÞ;

t ¼ 1; 2; . . . ;

ð2M þ 1Þð2N þ 1Þð2O þ 1Þg:

ð23Þ

ð11Þ

Note that the sequence of arguments o(t) employed in weight matrix Wr(k) is independent of argument k of function r. The solution ~ afi;j;kg that minimizes objective function  is defined by

ð12Þ

~ f fi;j;kg g: afi;j;kg ¼ fAT WrðkÞ Ag1 fAT WrðkÞ~

where /{i,j,k} is an integer-valued function defined as:

/fi;j;kg ðoÞ ¼ jk þ oj mod Nscan ;

for all o 2 [O:O]. Based on the cyclic characteristics of function /{i,j,k}, it is proved in Appendix A that

Ii;j;k ðm; n; oÞ ¼ IrðkÞ ðoÞ;

ð13Þ

where rðkÞ ¼ k mod Nscan . This allows us re-express objective function (10) as

ð~afi;j;kg Þ ¼

X

IrðkÞ ðoÞ fefi;j;kg ðm; n; oÞg2 ;

As the output of the CSG reconstruction filter is set to be g{i,j,k}(0, 0, 0), which is equal to a{i,j,k},{0,0,0}, the definition of a{i,j,k},{0,0,0} in (24) allows us to symbolically express output g{i,j,k}(0, 0, 0) as a linear combination of the form:

g fi;j;kg ð0; 0; 0Þ ¼ afi;j;kg;f0;0;0g

ð14Þ ¼

ðm;n;oÞ2Dfi;j;kg

To minimize objective function ð~ afi;j;kg Þ in (14) for coefficient a{i,j,k},{0,0,0} in a straightforward manner is extremely time consuming as the number of voxels involved in objective function ð~ afi;j;kg Þ is usually large. To derive coefficient a{i,j,k},{0,0,0} efficiently, we choose to employ in the CSG reconstructing problem, the principle developed in the 1D SG filtering algorithm, which allows us to compute a{i,j,k},{0,0,0} by simply calculating a linear combination of observed image data f associated with filled voxels in neighborhood Dfi;j;kg . The derivation of such linear combination expression for computing coefficient a{i,j,k},{0,0,0} can be achieved by first expressing objective function {i,j,k} in (14) in a matrix quadratic form of variables a{i,j,k},{p,q,r}, and then solving for a symbolic expression of coefficient a{i,j,k},{0,0,0}. Similarly, by employing integer t ranging from 1 to (2M + 1)(2N + 1)(2O + 1), and indexing functions:

mðtÞ ¼ bðt  1Þ=ðð2N þ 1Þð2O þ 1ÞÞc  M; nðtÞ ¼ bðt  1Þ=ð2O þ 1Þc mod ð2N þ 1Þ  N;

ð15Þ ð16Þ

oðtÞ ¼ ðt  1Þ mod ð2O þ 1Þ  O;

ð17Þ

we arrange all g{i,j,k}(m, n, o)’s and all image data f(i + m, j + n, k + o)’s involved in objective function (14) correspondingly as vectors:

~ f fi;j;kg ¼ ðf ði þ mðtÞ; j þ nðtÞ; k þ oðtÞÞ; T

ð2M þ 1Þð2N þ 1Þð2O þ 1ÞÞ ; and

t ¼ 1; . . . ; ð18Þ

ð24Þ

ð2Mþ1Þð2Nþ1Þð2Oþ1Þ X

afi;j;kg;t f ði þ mðtÞ; j þ nðtÞ; k þ oðtÞÞ;

t¼1

ð25Þ where

n

o

afi;j;kg;t ¼ ðAT WrðkÞ AÞ1 ðAT WrðkÞ~et Þ : 1

ð26Þ

Notation ~ et denotes a unit vector of which element t is equal to one, and { }1 denotes the first element of vector. The expression in (25) suggests that for each voxel t(i, j, k), there exist a particular set of coefficients a{i,j,k},t which allows us to accomplish automatically the process of polynomial least squares fitting by simply calculating the linear combination of the observed image intensity f in voxel neighborhood Dfi;j;kg [17]. The fact that coefficients a{i,j,k},t depend only on design matrix A and weight matrix Wr(k), and that such matrix A and weight matrix Wr(k) are known in advance permit us to compute coefficients a{i,j,k},t prior to the filtering operation. Furthermore, as the definition of weight matrix Wr(k) associated with various voxels t(i, j, k) can be subsumed under just N scan distinct patterns due to its modulo characteristics, as proved in Appendix B, this together with the fact that design matrix A is constant for all voxel t(i, j, k), indicates that there exist only N scan distinct sets of a{i,j,k},t’s associated with the entire voxel array. This suggests that we need to pre-compute just only N scan distinct sets of coefficients a{i,j,k},t for the entire reconstructing operation, thus making the process highly computation efficient.

140

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

Consequently, the solution ~ afi;j;kg that minimizes the above objective function is equal to

3. Cyclic Regularized Savitzky–Golay filters Instead of minimizing the objective function {i,j,k} as defined in (14), the least squares fitting employed in the CRSG filter chooses to minimize the objective function of the form:

nfi;j;kg ð~ aÞ ¼ fi;j;kg ð~ aÞ þ kefi;j;kg ð~ aÞ;

ð27Þ

where k is a regularization parameter based on the noise level of image data f. The first function {i,j,k} is the original objective function of the CSG filter while the second function e{i,j,k} is a regularization function reflects the prior information of solution ~ afi;j;kg [18]. In our investigation, we choose to define the regularization function as

efi;j;kg ð~ aÞ ¼ Ifi;j;kg ðm; n; oÞfðD1 g fi;j;kg ðm; n; oÞÞ2 þ ðD2 g fi;j;kg ðm; n; oÞÞ2 þ ðD3 g fi;j;kg ðm; n; oÞÞ2 g;

~ d2fi;j;kg ¼ ðD2 g fi;j;kg ðmðtÞ; nðtÞ; oðtÞÞ;

ð29Þ

t ¼ 1; . . . ; T

ð2M þ 1Þð2N þ 1Þð2O þ 1ÞÞ ; ~ d3fi;j;kg ¼ ðD3 g fi;j;kg ðmðtÞ; nðtÞ; oðtÞÞ;

ð30Þ

t ¼ 1; . . . ; T

ð2M þ 1Þð2N þ 1Þð2O þ 1ÞÞ :

ð31Þ

Based on these vectors, the regularization function (28) can thus be expressed as

efi;j;kg ð~ aÞ ¼ WrðkÞ

  T  T  T ~ d1fi;j;kg ~ d1fi;j;kg þ ~ d2fi;j;kg þ ~ d3fi;j;kg : d2fi;j;kg ~ d3fi;j;kg ~ ð32Þ

Applying the steps similar to those used in deriving ~ g in (21), we can express ~ d1 ; ~ d2 and ~ d3 in matrix formats as

~ d1fi;j;kg ¼ B~ afi;j;kg ;

~ d2fi;j;kg ¼ C~ afi;j;kg ;

and ~ d3fi;j;kg ¼ D~ afi;j;kg ;

ð33Þ

where B is a (2M + 1)(2N + 1)(2O + 1)  (P + 1)(Q + 1)(R + 1) design matrix, of which entry ts is equal to

Bts ¼ pðsÞ mðtÞjpðsÞ1j nðtÞqðsÞ oðtÞrðsÞ ;

ð34Þ

where C is a (2M + 1)(2N + 1)(2O + 1)  (P + 1)(Q + 1)(R + 1) design matrix, of which entry ts is equal to

Cts ¼ qðsÞ mðtÞpðsÞ nðtÞjqðsÞ1j oðtÞrðsÞ ;

ð35Þ

and D is a (2M + 1)(2N + 1)(2O + 1)  (P + 1)(Q + 1)(R + 1) design matrix, of which entry ts is equal to

Dts ¼ rðsÞ mðtÞpðsÞ nðtÞqðsÞ oðtÞjrðsÞ1j :

g fi;j;kg ð0; 0; 0Þ ¼

ð2Mþ1Þð2Nþ1Þð2Oþ1Þ X

i,j,k(0, 0, 0)

of the

bfi;j;kg ðtÞ f ði þ mðtÞ; j þ nðtÞ; k

t¼1

þ oðtÞÞ;

ð39Þ

bfi;j;kg ðtÞ ¼ fðAT WrðkÞ A þ k ðBT WrðkÞ B þ CT WrðkÞ C þ DT WrðkÞ DÞÞ1 ðAT WrðkÞ~ et Þg1 :

ð40Þ

As coefficients b{i,j,k} can be computed prior to the reconstructing operation, calculating g{i,j,k} in (39) is equivalent to simply convolving coefficients b{i,j,k} with image data f in voxel neighborhood Dfi;j;kg . This CRSG filter thus possesses the same computation complexity as that of the CSG filter. 4. Performance evaluation

t ¼ 1; . . . ;

ð2M þ 1Þð2N þ 1Þð2O þ 1ÞÞT ;

Correspondingly, we can compute the output g CRSG reconstruction filter according to

ð38Þ

where

ð28Þ

where D1 is the partial derivative with respect to m, D2 is the partial derivative with respect to n and D3 is the partial derivative with respect to o. Such regularization function is designed to reduce the first order derivatives of polynomial g{i,j,k} to encounter the strength of the speckle noise [19]. Using indexing function m,n, and o as defined in (15)–(17), we arrange all D1 g{i,j,k}, D2 g{i,j,k} and D3 g{i,j,k} involved in the regularization function (28) as vectors

~ d1fi;j;kg ¼ ðD1 g fi;j;kg ðmðtÞ; nðtÞ; oðtÞÞ;

~ afi;j;kg ¼ ðAT WrðkÞ A þ k ðBT WrðkÞ B þ CT WrðkÞ C   f fi;j;kg : þ DT WrðkÞ DÞÞ1 AT WrðkÞ ~

ð36Þ

Using the definitions of ~ d1fi;j;kg , ~ d2fi;j;kg and ~ d3fi;j;kg as defined in (33), we can thus rewrite the objective function in (27) as

 aTfi;j;kg BT B ~ nfi;j;kg ð~ aÞ ¼ fi;j;kg ð~ aÞ þ k WrðkÞ ~ afi;j;kg þ ~ aTfi;j;kg CT C ~ afi;j;kg  afi;j;kg : þ~ aTfi;j;kg DT D~ ð37Þ

In this section, we evaluate the CRSG filter in reconstructing 3D ultrasound images of a synthetic sphere and of a clinical carotid artery bifurcation, and compare its results with those reconstructed from the commonly used VNN interpolation, the DW interpolation, and the recently developed ADW interpolation. In addition, we have also evaluated the CRSG filter compared with the CSG filter in order to investigate the degree of efficiency due to the regularization mechanism. Note that the CSG reconstruction filter is the special case of the CRSG filter where the regularization parameter k is set to 0. To enhance the VNN interpolation with the speckle noise reduction capability, we then post-process the 3D resultant images obtained from the VNN interpolation by speckle noise reduction filters. Two versions of such enhanced VNN interpolation are tested: the VNN interpolation post-processed by the Adaptive Speckle Reduction (ASR) filter [14], and the VNN interpolation post-processed by the Adaptive Weighted Median (AWM) filter [15]. Furthermore, similar to the CSG and the CRSG reconstruction algorithms, the VNN, the DW, and the ADW algorithms have been also implemented here by taking into account the regularity of the B-scan slices acquired from the mechanical linear scanning framework. Note that utilizing such regularity of the filled voxels allows all the VNN, the DW, and the ADW reconstruction algorithms to determine the exact locations of all the relevant filled voxels analytically, without searching the entire neighborhood for the filled voxels, thus speeding up the reconstruction process considerably. In the first evaluation, we aim to reconstruct a 3D sphere from a set of 32 parallel synthetic ultrasound cross sectional slices. We first mathematically synthesize 128 parallel ultrasound cross sectional slices, each with 128  128 pixels, by defining the tissue reflectivity of slice k as a disk of the form:

( T k ði; jÞ ¼

1

100; if fði  ic Þ2 þ ðj  jc Þ2 þ ðk  kc Þ2 g2 6 r; 10; otherwise;

ð41Þ

where i is defined over 1, . . . , 128, j over 1, . . . , 128, and k over 1, . . . , 128. The center coordinate ic, jc and kc of the sphere are set to be 64 and the sphere radius r is set to be 32. To simulate speckle noise, we multiply the tissue reflectivity at each pixel by a random value generated according to the Rayleigh probability distribution

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

of mean one [20]. We then uniformly select 32 slices from the 128 slices by choosing one out of every four slices, and use them to reconstruct the original 3D spherical image on a 3D regular grid of 128  128  128 voxels, The voxel width, the voxel height, and the voxel thickness are all set to be the same value as the pixel width of the cross sectional image. Note that parameter N scan in this evaluation is equal to 4. For each reconstruction algorithm, we measure its reconstruction capability of the synthetic ultrasound by the normalized mean square error (NMSEX) defined as

NMSEX ¼

P  2 ~ X ðf fi;j;kg  f fi;j;kg Þ ; P  2 X ðf fi;j;kg Þ

ð42Þ

where f fi;j;kg is the voxel intensity of the underlying image without speckle noise, ~f fi;j;kg is the reconstructed voxel intensity obtained from the reconstructing algorithm, and X is the domain of all voxels in the image volume. A smaller value of NMSEX reflects better capability in both interpolation of unfilled voxels and in reduction of the speckle noise in filled voxels. In this evaluation, we define the voxel neighborhood in all aforementioned reconstruction techniques to be of 9  9  9 voxels. For the ADW interpolation, the output g{i,j,k} at voxel t(i, j, k) is computed according to the formula:

g fi;j;kg ¼

P D W fi;j;kg ffi;j;kg P ; D W fi;j;kg

ð43Þ

where D denotes the domain of filled voxels falling into the predefined neighborhood region centered about t(i,j,k), f{i,j,k} denotes the intensity of the filled voxel in neighborhood D. The weighting factor W{i,j,k} is the relative weight for the filled voxel in neighborhood D, and is defined by

W fi;j;kg

8 0; > > > < ¼ 1; > > > : 1

dafi;j;kg

if if ; if

r2 l r2 l

6 H0 \ f fi;j;kg R ½l  r; l þ r;

6 H0 \ f fi;j;kg 2 ½l  r; l þ r;   r > H ; a ¼ b  r2  H 0 0 þ 1; l l 2

ð44Þ

141

where d{i,j,k} is the distance of the filled voxel in neighborhood D to its center, b is a parameter empirically set by the operator for adjusting the interpolation performance, l and r are the mean and standard deviation of all neighboring filled voxels, respectively. In this evaluation, we set the local homogeneity threshold H0 to be equal to 30.83 following the procedure recommended in [12] and empirically set b to 0.25. For our new CSG and CRSG filters, we define the polynomial gi,j,k to be of order 2 in all m, n and o and empirically set parameter k of the CRSG filter to be 0.1. For the speckle noise reduction filters which we use to post-process results of the VNN interpolation algorithm, the center-weight parameter and the a scaling constant of the AWM filter are respectively set to 20 and 1, while the ASR filter requires no parameter setting. Figs. 3 and 4 depict, respectively, exemplary cross sectional views and sagittal views of the synthetic 3D images obtained from various reconstruction algorithms. Comparing the images in Figs. 3 and 4 among various algorithms, it can be observed that the VNN interpolation post-processed by the ASR filter and by the AWM filter perform worst boundaries and appear to possess more reconstruction artifacts particularly on the circular disks, when compared with the images obtained from the other algorithms. When visually comparing the remaining resultant images, we see that the image obtained from the ADW, the CSG and the CRSG algorithms perform comparatively equivalent. The resultant images obtained from these algorithms appear visually to possess similar degree of reconstruction artifacts and image enhancement, while that obtained from the DW interpolation seem to suffer from over-smoothing at the circular disk boundaries. However, by further comparing the average NMSEX associated with all algorithms, we can clearly see that the CRSG reconstruction filter performs best with minimum average NMSEX, as summarized in Table 1. Here each average NMSEX in this table is derived from averaging an ensemble of the sphere corrupted by ten different realizations of Rayleigh speckle noise. Fig. 5 shows the 3D images of a spherical object obtained from each reconstruction algorithm. The images are rendered by using

Fig. 3. The exemplary cross sectional views of a spherical target reconstructed from (a) The VNN interpolation post-processed by the ASR filter. (b) The VNN interpolation post-processed by the AWM filter. (c) The DW interpolation. (d) The ADW interpolation. (e) The CSG reconstruction filter. (f) The CRSG reconstruction filter.

142

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

Fig. 4. The exemplary sagittal views of a spherical target reconstructed from (a) The VNN interpolation and post-processed by the ASR filter. (b) The VNN interpolation and post-processed by the AWM filter. (c) The DW interpolation. (d) The ADW interpolation. (e) The CSG reconstruction filter. (f) The CRSG reconstruction filter.

Table 1 The average NMSEX of the spherical object obtained from each reconstruction algorithm. Algorithms

NMSEX

VNN+ASR VNN+AWM DW ADW CSG CRSG

0.1570 0.1165 0.0942 0.0471 0.0465 0.0321

the thresholding and the ray-casting method [21], and we have set the intensity threshold level to be 55, which is equal to the average of the background intensity and the sphere intensity. The figure reveals that the sphere obtained from the CRSG reconstruction filter visually seems to be closet to the ideal sphere compared with the others while those obtained from the VNN interpolation seem to be far worst different from the ideal sphere. In comparison with the CSG reconstruction filter, the result of spherical surface obtained from the CRSG reconstruction filter appears much smoother due to the mechanism of regularization.

Fig. 5. 3D reconstruction of spherical object obtained from (a) The VNN interpolation post-processed by the ASR filter. (b) The VNN interpolation post-processed by the AWM filter. (c) The DW interpolation. (d) The ADW interpolation. (e) The CSG reconstruction filter. (f) The CRSG reconstruction filter.

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

143

Fig. 6. A series of slices of the carotid artery bifurcation acquired by the mechanical linear scanning system.

As a complimentary study to the first evaluation, we compare the performance of the aforementioned algorithms in reconstructing a clinical 3D ultrasound image of the carotid artery bifurcation. A series of 64 parallel B-scan images, each cropped to 90  90 pixels, is captured in a single sweep over the region of a carotid artery bifurcation by the mechanical linear scanning system on the ultrasound machine GE Logiq 9, as displayed in Fig. 6. Similarly to the first evaluation, we define all of the initial parameters in each reconstruction algorithm to have the same definition as aforementioned problem except in the ADW interpolation where the threshold value of local homogeneity is set to be 4.447. To evaluate the reconstruction performance of each algorithm in this clinical ultrasound images, we remove every even slices of all the 64 acquired slices and then employ various algorithms to reconstruct them back from the remaining 32 acquired slices. Fig. 7a–f depict the exemplary reconstructed slice of the ultrasound images of carotid artery bifurcation obtained from various reconstruction algorithms. The corresponding differences between the reconstructed results and the original acquired images are portrayed in Fig. 7g–l. These differences reflect the reconstruction efficiency; if the reconstructed result is closed to the acquired image, the difference should reveal only the noise which uncorrelated to the different tissues [22]. Inspecting the reconstructed results of all algorithms, we find that the one derived from the CRSG appears to visually possess less residual image structures than the other reconstruction algorithms. Fig. 8 illustrates the quantitative comparison of the voxel intensities along row 45 of Fig. 7. Fig. 8a and b depict the reconstruction profile images obtained from the VNN interpolation post-processed by the ASR filter and by the AWM algorithms, it can be observed that its resultant profiles suffer from the remaining of speckle noise due to the uncertainty of estimating the coefficients by using local statistics for classifying homogeneous and inhomogeneous regions. Comparing the remaining algorithms, the DW, the ADW, the CSG and the CRSG algorithms, as display respectively in Fig. 8c–f, can reduce the speckle noise well when compared with the VNN post-processed by the ASR filter and by the AWM filter. However, it can be observed that the DW algorithm smooth out the boundaries between vessels of the carotid artery bifurcation while the ADW, the CSG and the CRSG algorithms, as shown respectively in Fig. 9d–f, provide a better boundaries preservation. Although the reconstruction profiles derived from the ADW and the CSG algorithms seem to be able to preserve most of texture patterns, they still appear to remain a small effect of speckle noise especially in the homogeneous regions when compared with the CRSG algorithm.

Fig. 7. The exemplary reconstructed slices of carotid artery bifurcation, and the different images between its corresponding resultant image and the acquired image. (a) and (g) The VNN interpolation post-processed by the ASR filter. (b) and (h) The VNN interpolation post-processed by the AWM filter. (c) and (i) The DW interpolation. (d) and (j) The ADW interpolation. (e) and (k) The CSG reconstruction filter. (f) and (l) The CRSG reconstruction filter.

144

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

140

140 Observed profile VNN+ASR

100 80 60 40 20 0

Observed profile VNN+AWM

120

intensities along row 45

intensities along row 45

120

100 80 60 40 20

10

20

30

40 50 column

60

70

80

0

90

10

20

30

40 50 column

(a)

intensities along row 45

intensities along row 45

90

100 80 60 40

Observed profile ADW

120

20

100 80 60 40 20

10

20

30

40 50 column

60

70

80

0

90

10

20

30

40 50 column

(c)

60

70

80

90

(d)

140

140 Observed profile CSG

100 80 60 40 20

Observed profile CRSG

120

intensities along row 45

120

intensities along row 45

80

140 Observed profile DW

120

0

70

(b)

140

0

60

100 80 60 40 20

10

20

30

40 50 column

60

70

80

90

0

10

20

30

40 50 column

(e)

60

70

80

90

(f)

Fig. 8. The comparisons of intensities along row 45 of the exemplary reconstructed slices of carotid artery bifurcation in Fig. 7. (a) The VNN interpolation post-processed by the ASR filter. (b) The VNN interpolation post-processed by the AWM filter. (c) The DW interpolation. (d) The ADW interpolation. (e) The CSG reconstruction filter. (f) The CRSG reconstruction filter.

In additional to the visual assessment, we measure the performance of various reconstruction algorithms by computing the normalized mean square error (NMSEU) defined as

P NMSEU ¼

U ðffi;j;kg

P

 ~f fi;j;kg Þ2

U ðffi;j;kg Þ

2

;

ð45Þ

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147 Table 2 The average NMSEU of the images of the carotid artery bifurcation obtained from different reconstruction algorithms. Algorithms

NMSEU

VNN+ASR VNN+AWM DW ADW CSG CRSG

0.1041 0.0836 0.0618 0.0463 0.0545 0.0321

145

where U is the domain of all voxels associated with the even slices that have been removed for reconstruction. Note that the NMSEU is the normalized mean squared difference between the original and the reconstructed voxel intensities associated with all the even slices in the acquired B-scans. As indicated in Table 2, the image processed by the CRSG reconstruction algorithm consistently appears to has least NMSEU. To investigate the volume reconstruction of each algorithm, the aforementioned 64 parallel of clinical 2D B-scans are then used to reconstruct a 3D image on an array of 90  90  256 voxels. That is, we need to reconstruct 256 slices from the 64 originally

Fig. 9. Three-dimensional reconstruction of the carotid artery bifurcation are displayed in multiplanar reformatting. (a) The VNN interpolation post-processed by the ASR filter. (b) The VNN interpolation post-processed by the AWM filter. (c) The DW interpolation. (d) The ADW interpolation. (e) The CSG reconstruction filter. (f) The CRSG reconstruction filter.

146

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

acquired B-scan slices. Note that parameter N scan in this evaluation is equal to 4, indicating that there are every 3 new slices to be reconstructed between each pair of the originally acquired B-scan slices. Fig. 9 displays the 3D reconstruction of carotid artery bifurcation in multiplanar reformatting (MPR) by the various aforementioned algorithms. Fig. 9a and b depict, respectively, the reconstruction image obtained from the VNN interpolation postprocessed by the ASR filter and by the AWM filter. It can be observed that their resultant images suffer from the elevational artifacts which occur from interpolating the missing data by a constant defined by the nearest voxel intensity and then estimate the smoothed intensities from the redundant data in neighborhood. In comparison with the other reconstruction algorithms, the resultant image obtained from the DW algorithm as depicted in Fig. 9c tends to be overly blurred, while resultant images obtained from the ADW, the CSG, and the CRSG algorithms, as depicted in Fig. 9d–f, respectively, exhibit to better preserve the artery boundaries. However, when closely visually inspecting the results derived from the ADW and the CSG algorithms, it can be observed that there appear reconstruction artifacts between the acquired B-scan slices due to the lack of sufficient degree of smoothing along the elevation direction. Comparing the ADW and the CSG with the CRSG algorithm, the result derived from the CRSG algorithm seem to be better in compromising between the fine detail preservation and the degree of smoothness. Table 3 shows the computation time required by the various algorithms for reconstructing 90  90  256 voxels of the carotid artery bifurcation, when all of the reconstruction algorithms are implemented as programs written in MATLAB running on a personal computer with Pentium IV-2 GHz CPU. Table 3 indicates that the VNN interpolation post-processed by the ASR filter and by the AWM filter, and the ADW interpolation require much more time to reconstruct compared with the DW, the CSG and the CRSG algorithms, i.e. the DW, the CSG and the CRSG algorithms require only 40 s while the VNN interpolation post-processed by the ASR filter and by the AWM filter require 420 and 3086 s and the ADW interpolation requires 519 s. The extremely high computation time of the VNN interpolation post-processed by the ASR filter and by the AWM filter, and the ADW interpolation are due to the need to estimate the local mean and variance of each voxel during the reconstruction process. As for the DW, the CSG and the CRSG algorithms, their operations are formulated as a convolution sum. Even though the convolution coefficients required at various voxels may differ from one voxel to another voxel, the cyclic regularity among the locations of the observed B-scan voxels as acquired from mechanical linear scanning allows us to classify the sets of coefficients at various voxels used in the reconstruction into a finite number of sets, all of which can be pre-computed before the reconstruction operation. Consequently, all the DW, the CSG and the CRSG algorithms need to perform only the convolution sum during the reconstruction, thus resulting in less computation time when compared with the VNN interpolation post-processed by the ASR filter and by the AWM filter, and the ADW interpolation algorithms.

5. Conclusion The preliminary evaluation in Section 4 indicates that the CRSG filter is more effective in both speckle reduction and geometric reconstruction of 3D ultrasound images than most available reconstruction algorithms. This better performance is attributed to the following two factors: The first factor is that the CRSG filter reconstructs the 3D image estimate via the much flexible 3D polynomial least squares fitting on the observed B-scan image intensities. Such polynomial fitting enables the CRSG filter to better eliminate the speckle noise when compared with the VNN interpolation postprocessed by the ASR filter and the AWM filter, which attempt to interpolate the 3D image by a constant of nearest voxel intensity and then estimate the smoothed intensities from the redundant values in neighborhood. Furthermore, it also allows the CRSG filter to better avoid over smoothing image edges and fine details when compared with the DW interpolation algorithm, which attempts to interpolate the 3D image by averaging nearby weighted voxel intensities. Compared with the ADW interpolation algorithm, adaptive version of the DW interpolation algorithm, the CRSG reconstruction algorithm has demonstrated to possess better capability in image interpolation, speckle suppression, and computation time requirement than the ADW interpolation algorithm. The second factor allowing the CRSG filter to achieve more superior performance is that the CRSG filter employs the regularization function as a mechanism to balance between the degree of speckle reduction and the degree of detail preservation. Such technique helps us to ensure that the least square fitting of polynomial function computed by the CRSG reconstruction algorithm will be robust and conform to the strength of speckle noise in ultrasound images. In addition, since the computation of the CRSG filter is rather simple, i.e., involving mostly of a linear convolution operation, such new technique possesses a large potential in real-time ultrasound image reconstruction and in assisting automated segmentation. Acknowledgements This work was jointly supported by Thailand Research Fund under Grant Number RSA4580027; the Fund from the Cooperative Project between Department of Electrical Engineering and Private Sector for Research and Development; the 90th anniversary of Chulalongkorn University fund, Ratchadaphisek Somphot Endowment, Chulalongkorn University. The authors also wish to acknowledge Surot Phomwaranont, Leadership Product Manager of Ultrasound Division, GE Medical Systems (Thailand) Ltd., for many valuable suggestions and for providing the ultrasound images used in the evaluation. Appendix A. Partition of indicator function I fi;j;kg into N scan distinct groups Consider a 3D voxel array t indexed by (i, j, k), and voxel neighborhood Di;j;k defined by

Dfi;j;kg ¼ ftði þ m; j þ n; k þ oÞ; Table 3 The computation time of various algorithms in reconstructing 90  90  256 voxels of the carotid artery bifurcation. Algorithms

Time (s)

VNN+ASR VNN+AWM DW ADW CSG CRSG

420 3086 40 519 40 40

n ¼ fN : Ng;

m ¼ fM : Mg;

o ¼ fO : Ogg:

Let Ifi;j;kg denotes an indicator function defined over voxel neighborhood Dfi;j;kg , and be of the form:

 Ifi;j;kg ðm; n; oÞ ¼

1; if

/fi;j;kg ðoÞ ¼ 0;

0; otherwise;

ðA:2Þ

where /{i,j,k} is an integer-valued function defined as:

/fi;j;kg ðoÞ ¼ jk þ oj mod Nscan ;

ðA:3Þ

147

P. Toonkum et al. / Ultrasonics 51 (2011) 136–147

for all o 2 [O:O]. ~ be any integers and o be an integer in [O:O]. It can Let k and k be proved based on the cyclic characteristics of modulo function that

~ scan Þ þ oj mod Nscan : jk þ oj mod Nscan ¼ jðk þ kN

ðA:4Þ

As the expression on the left hand side of (A.4) defines function /{i,j,k} and that on the right hand side defines function /fi;j;kþkN , ~ scan g Eq. (A.4) indicates that the values of function /{i,j,k} and function /fi;j;kþkN are identical for any given o. It can thus be inferred that ~ scan g function /{i,j,k} is periodic with respect to parameter k and has the fundamental period of N scan . Due to the facts that indicator function Ifi;j;kg is completely defined by function /{i,j,k} and that function /{i,j,k} is periodic with k but totally independent of i and j, we conclude that

Ifi;j;kg ¼ Ifi0 ;j0 ;ðkþkN ~ scan Þg ;

ðA:5Þ

0 0 ~ Such equality thus allows us to clasfor all indexes i, i , j, j , k and k. sify indicator functions Ifi;j;kg associated with various {i, j, k} into N scan distinct categories of I0 ; . . . ; Ij ; . . . ; IðNscan 1Þ ; where

 Ij ðoÞ ¼

1; if 0;

jj þ oj mod Nscan ¼ 0;

otherwise;

ðA:6Þ

for j 2 f0; . . . ; ðN scan  1Þg. Note that we choose to express the argument of indicator function Ij as (o) instead of (m,n,o) since the value of Ij is totally independent of arguments m and n. The category Ij to which indicator function Ifi;j;kg belongs can be identified through an integer-valued function r of the form:

rðkÞ ¼ k mod Nscan ;

ðA:7Þ

where k is defined over . . ., 1, 0, 1, . . . . Such function r therefore allows us to establish that

Ii;j;k ðm; n; oÞ ¼ IrðkÞ ðoÞ;

ðA:8Þ

for all m 2 [M:M], n 2 [N:N], and o 2 [O:O]. Appendix B. Partition of weight matrix Wr(k) into N scan distinct groups Consider a diagonal matrix Wr(k) of the form:

WrðkÞ ¼ diag fIrðkÞ ðoðtÞÞ; t ¼ 1; 2; . . . ; ð2M þ 1Þð2N þ 1Þð2O þ 1Þg:

ðB:1Þ

~ define matrix W For any given integers k and k, ~ scan Þ accordrðkþkN ing to (B.1) as:

WrðkþkN ~ scan Þ ¼ diag fIrðkþkN ~ scan Þ ðoðtÞÞ;

t ¼ 1; 2; . . . ;

ð2M þ 1Þð2N þ 1Þð2O þ 1Þg:

ðB:2Þ

As indicator functions IrðkþkN ~ scan Þ and IrðkÞ are equal due to the cyclic characteristics of function r, we then substitute IrðkþkN in ~ scan Þ (B.2) with IrðkÞ , and thus obtain

WrðkþkN ¼ diag fIrðkÞ ðoðtÞÞ; ~ scan Þ

t ¼ 1; 2; . . . ;

ð2M þ 1Þð2N þ 1Þð2O þ 1Þg:

ðB:3Þ

Since the expression on the right hand of (B.3) is equal to the that of (B.1), we conclude that

WrðkþkN ¼ WrðkÞ : ~ scan Þ

ðB:4Þ

It can thus be inferred from this equation that the pattern of weight matrix Wr(k) is periodic with respect to parameter k and has the fundamental period of N scan . This fact therefore allows us to categorize weight matrix Wr(k) associated with various k into N scan distinct patterns of W0 ; . . . ; Wj ; . . . ; WðNscan 1Þ ; where

Wj ¼ diag fIj ðoðtÞÞ;

t ¼ 1; 2; . . . ;

ð2M þ 1Þð2N þ 1Þð2O þ 1Þg; for

ðB:5Þ

j 2 f0; . . . ; ðNscan  1Þg.

References [1] T.R. Nelson, D.H. Pretorius, Three-dimensional ultrasound imaging, Ultrasound Med. Biol. 24 (9) (1998) 1243–1270. [2] J.M. Sanches, J.S. Marques, A rayleigh reconstruction/interpolation algorithm for 3-d ultrasound, Pattern Recogn. Lett. 21 (10) (2000) 917–926. [3] R.S. José-Estérpar, M. Martín-Fernández, P.P. Caballero-Martínez, C. AlberolaLópez, J. Ruiz-Alzola, A theoretical framework to three-dimensional ultrasound reconstruction from irregularly sampled data, Ultrasound Med. Biol. 29 (2) (2003) 255–269. [4] R.J. Housden, A.H. Gee, G.M. Treece, R.W. Prager, Sensorless reconstruction of unconstrained freehand 3-d ultrasound data, Ultrasound Med. Biol. 33 (9) (2007) 408–419. [5] A. Fenster, D.B. Downey, 3-d ultrasound imaging: a review, IEEE Eng. Med. Biol. Mag. 15 (6) (1996) 41–51. [6] J.G. Abbott, F.L. Thurstone, Acoustic speckle: theory and experimental analysis, Ultrasound Imag. 1 (1979) 303–324. [7] R.N. Rohling, A.H. Gee, L. Berman, A comparison of freehand three-dimensional ultrasound reconstruction techniques, Med. Image Anal. 3 (4) (1999) 339–359. [8] R.W. Prager, R. Rohling, A. Gee, L. Berman, Rapid calibration for 3-d freehand ultrasound, Ultrasound Med. Biol. 24 (1998) 855–869. [9] H.A. McCann, J.C. Sharp, T.M. Kinter, C.N. McEwan, C. Barillot, J.F. Greenleaf, Multidimensional ultrasonic imaging for cardiology, Proc. IEEE 76 (9) (1998) 1063–1073. [10] C.D. Barry, C.P. Allott, N.W. John, P.M. Mellor, P.A. Arundel, D.S. Thomson, J.C. Waterton, Three-dimensional freehand ultrasound: image reconstruction and volume analysis, Ultrasound Med. Biol. 23 (1997) 1209–1224. [11] Q.H. Huang, Y.P. Zheng, Volume reconstruction of freehand three-dimensional ultrasound using median filters, Ultrasonics 48 (3) (2008) 182–192. [12] Q.H. Huang, M.H. Lu, Y.P. Zheng, Z.R. Chi, Speckle suppression and contrast enhancement in reconstruction of freehand 3D ultrasound images using an adaptive distance-weighted method, Appl. Acoust. 70 (1) (2009) 21–30. [13] A. Savitzky, M.J.E. Golay, Smoothing and differentiation of data by simplied least squares procedure, Anal. Chem. 36 (1964) 1627–1639. [14] J. Bamber, C. Daft, Adaptive filtering for reduction of speckle in ultrasound pulse-echo images, Ultrasonics 24 (1) (1986) 41–44. [15] W.M.T. Loupas, P. Allen, An adaptive weighted median filter for speckle suppression in medical ultrasound images, IEEE Trans. Circ. Syst. 36 (1) (1989) 129–135. [16] W. Huang, Y. Zheng, MMSE reconstruction for 3D freehand ultrasound imaging, Int. J. Biomed. Imag. 1 (3) (2008) 1–8. [17] C. Chinrungrueng, A. Suvichakorn, Fast edge-preserving noise reduction for ultrasound images, IEEE Trans. Nucl. Sci. 48 (3) (2001) 849–854. [18] H.K. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer Academic Publishers, 2000. [19] P. Toonkum, P. Boonvisut, C. Chinrungrueng, Real-time speckle reduction of ultrasound images based on regularized Savitzky–Golay filters, Proc. IEEE ICBBE. 1 (1) (2008) 2311–2314. [20] R.F. Wagner, S.W. Smith, J.M. Sandrik, H. Lopez, Statistics of speckle in ultrasound b-scans, IEEE Trans. Son. Ultrasonics 30 (3) (1983) 156–163. [21] J.D. Foley, A.V. Dam, S.K. Feiner, J.F. Hughes, Computer Graphics: Principle and Practice, Addison-Wesley Pubplishing company, Inc, 1995. [22] K. Krissian, C.F. Westin, R. Kikinis, K.G. Vosburgh, Oriented speckle reducing anisotropic diffusion, IEEE Trans. Image Proces. 16 (5) (2007) 1412–1424.