### System

The freehand 3D ultrasound system consists of three modules: a conventional 2D ultrasound scanner (DC-7, Mindray Medical International Ltd., Shenzhen, China) used to acquire ultrasound image, an electromagnetic spatial sensing device (Aurora
[23], NDI, Ontario, Cannada) acquiring the position and orientation of ultrasound images and surgical instruments, and a workstation with custom-designed software used for data collection, volume reconstruction, and visualization
[1]. Figure
1 illustrate this system framework. The spatial information (position (*x*,*y*,*z*) and orientation (*R*
_{
x
},*R*
_{
y
},*R*
_{
z
})) of ultrasonic probe and surgical instrument embeding sensor is recorded by Aurora system connects to the workstation through its USB port so that it can be acquired by a custom-designed software by use of Aurora system API. Besides, the real-time ultrasound video is captured by a video capture card (RGB-133, VTimage Inc., Shenzhen, China) installed in the workstation.

### Data collection

Before volume reconstruction, data collection is a very important step which influences accuracy and efficiency of reconstruction. Data collection consists of two steps: calibration between spatial data and 2D B-scan images, and selecting ROI (region of interest). In this paper, the set of B-scan image *I*
_{
i
} and its position *T*
_{
i
} are collected by the data collection method of
[1].

### Volume reconstruction

A 3D volume data is reconstructed from the collected data, including the 2D B-scan images and its spatial information. In this study, the algorithm for volume reconstruction is composed of two stages: bin-filling and regression.

#### Bin-filling stage

The bin-filling stage is to map the pixel in 2D B-scans into the voxel in 3D volume data based on its corresponding positional information. In this freehand 3D ultrasound system, the mapping of the coordinate system from the 2D B-scans to 3D volume is named after forward mapping and is defined as

\begin{array}{l}{V}_{r}=M\times {V}_{p}\end{array}

(1)

where *V*
_{
p
} is the physical position, *M* is the forward transformation matrix, and *V*
_{
r
} is the resulting voxel location in the reconstructed 3D volume data. The forward transformation matrix *M* must be decomposed and implemented to find an matrix for this transformation, which is discussed in detail in
[1].

#### Regression stage

Since the scanning of freehand-style is subjective, the collected B-scan images are usually irregular and highly sparse. Therefore, there are some gaps in the reconstructed volume after the bin-filling stage, as addressed in Figure
2(a). The goal of the regression stage is to make the nonparametric estimation
[20] for the whole volume data from the previous sparse volume data.

The sparse volume data after the bin-filling stage are given by

{Y}_{i}=r({X}_{i})+{\epsilon}_{i},\phantom{\rule{1em}{0ex}}i=1,\cdots ,P

(2)

where *r*(·) is the regression function, *X*
_{
i
} = (*X*
_{
i 1},*X*
_{
i 2},*X*
_{
i 3})^{T} is the three-dimension coordinate of data, *ε*
_{
i
}s are the independent and identically distributed zero mean noise values.

Specifically, if *X* is near the sample at *X*
_{
i
}, we can approximate it with a *N*-term Taylor series

\begin{array}{ll}\phantom{\rule{5pt}{0ex}}r({X}_{i})& \approx r(X)+{r}^{\prime}(X)({X}_{i}-X)+\frac{1}{2!}{r}^{\mathrm{\prime \prime}}(X){({X}_{i}-X)}^{2}\\ \phantom{\rule{1em}{0ex}}+\cdots +\frac{1}{N!}{r}^{N}(X){({X}_{i}-X)}^{N}\\ ={\beta}_{0}+{\beta}_{1}({X}_{i}-X)+{\beta}_{2}{({X}_{i}-X)}^{2}+\cdots +{\beta}_{N}{({X}_{i}-X)}^{N}\end{array}

(3)

A least-squares formulation capturing this idea is to solve the folowing optimization problem:

\begin{array}{ll}\phantom{\rule{5pt}{0ex}}\underset{\{{\beta}_{n}\}}{\mathit{\text{min}}}\sum _{i=1}^{P}& \left[{Y}_{i}-{\beta}_{0}-{\beta}_{1}\left({X}_{i}-X\right)-{\beta}_{2}{\left({X}_{i}-X\right)}^{2}\right.\\ \phantom{\rule{14.0pt}{0ex}}{\left(\right)close="]">-\cdots -{\beta}_{N}{\left({X}_{i}-X\right)}^{N}}^{}2& \frac{1}{h}K\left(\frac{{X}_{i}-X}{h}\right)\end{array}\n

(4)

where *K*(·) is the kernel function which penalizes distance away from the local position where the approximation is centered, and the smoothing parameter *h* (bandwidth) controls the strength of this penalty. In particular, the function *K* is a symmetric function which attains its maximum at zero, satisfying

\underset{{R}^{1}}{\int}\mathit{\text{tK}}(t)\mathit{\text{dt}}=0,\phantom{\rule{1em}{0ex}}\underset{{R}^{1}}{\int}{t}^{2}K(t)\mathit{\text{dt}}=c

(5)

where *c* is a constant value. The choice of the particular form of the function *K* is usually Gaussian, exponential, or other forms, which comply with the above constraints. Because the choice of the kernel has little impact on the accuracy of estimation. Therefore, the Gaussian kernel, being computational complexity, is selected in this paper.

For the estimation problem based on Least Square Method upon showed in Equation 4, the order *N* effect the accuaracy and complexity of local approximation of the volume data. Therefore, it must be appropriately chosen. In the nonparametric statistics literature, locally constant, linear, and quadratic approximations (corresponding to *N*=0,1,2) have been considered most widely
[24–27].

The Kernel function *K* is now a function of 3 variables. Given a nonsingular positive definite 3 × 3 bandwidth matrix *H*, which is defined

{K}_{H}(X)=\frac{1}{{|H|}^{\frac{1}{2}}}K\left({H}^{-\frac{1}{2}}X\right).

(6)

Often, one scales each covariate to have the same mean and variance, then use the kernel

{h}^{-3}K\left(\frac{||X||}{h}\right)

(7)

where *K* is any one-dimensional kernel. therefore, there is a single bandwidth parameter *h*. At a target value *X* = (*X*
_{1},*X*
_{2},*X*
_{3})^{T}, the local sum of squares is given by

\sum _{i=1}^{n}{w}_{i}(X){\left({Y}_{i}-{\beta}_{0}-\sum _{j=1}^{3}{\beta}_{j}\left({X}_{\mathit{\text{ij}}}-{X}_{j}\right)\right)}^{2}

(8)

where

{w}_{i}(X)=K\left(\frac{||{X}_{i}-X||}{h}\right).

(9)

The estimator is

\hat{r}(X)=\hat{{\beta}_{0}}

(10)

where
\hat{\beta}={\left(\hat{{\beta}_{0}},\cdots ,\hat{{\beta}_{3}}\right)}^{T} is the value of *β* = (*β*
_{0},⋯,*β*
_{3})^{T} that minimizes the weighted sums of squares. The solution
\hat{\beta} is

\hat{\beta}={\left({X}_{x}^{T}{W}_{x}{X}_{x}\right)}^{-1}{X}_{x}^{T}{W}_{x}Y

(11)

where

\mathit{\text{Xx}}=\left(\begin{array}{cccc}1& {X}_{11}-{X}_{1}& {X}_{12}-{X}_{2}& {X}_{13}-{X}_{3}\\ 1& {X}_{21}-{X}_{1}& {X}_{22}-{X}_{2}& {X}_{23}-{X}_{3}\\ 1& {X}_{31}-{X}_{1}& {X}_{32}-{X}_{2}& {X}_{33}-{X}_{3}\end{array}\right)

(12)

and *W*
_{
x
} is the diagonal matrix whoses (*i*,*i*) element is *w*
_{
i
}(*X*). Therefore,

\hat{r}(X)=\hat{{\beta}_{0}}={e}_{1}^{T}{\left({X}_{x}^{T}{W}_{x}{X}_{x}\right)}^{-1}{X}_{x}^{T}{W}_{x}Y

(13)

where *e*
_{1} is a column vector with all elements equal to zero but the first element equal to one,
\hat{r}(X) is the final values of reconstructed volume data.