2. 中国科学院大学 计算机与控制学院, 北京 100049
2. School of Computer and Control Engineering, University of Chinese Academy of Sciences, Beijing 100049, China
B型心脏超声图像分割是指根据像素和先验知识将图像分割成不同区域并得到壁、腔、瓣、流四种结构的过程[1]。由于超声图像具有高噪声、多伪影、成像缺陷以及边界模糊等特点,在自然图像分割中效果较好的分割算法无法直接用于心脏超声图像。目前,超声图像的主流分割方法有基于模型的方法、基于阈值的方法、基于图论的方法以及基于聚类的方法等。
其中,基于模型的方法如基于卷积神经网络(Convolutional Neural Network,CNN)的模型[2]和蛇形模型[3],虽然效果较好,但是需要大量样本进行训练;另外,基于CNN的方法目前难以融合人们已经构建出来的先验知识,而蛇形模型针对不同的目标区域需要训练不同的模型。基于阈值的方法如大津算法,阈值的计算时间较长,难以满足实时要求。近年来,相关研究者在此基础上作了大量改进,然而由于超声图像具有高噪声、多伪影、成像缺陷以及边界模糊等特征[4],基于阈值的方法的分割效果仍然无法满足实际需求。基于图论的方法如图割法等需要人为指定前景和背景的种子点;改进的分割方法如GrowCut算法[5]可以根据图像灰度直方图自动生成初始种子模板,但仍需用户的二次交互,在对大量图像进行分割时工作量较大。基于聚类方法是对图像像素进行聚类,在图像分割上获得了广泛应用。如文献[6]提出了改进的K均值聚类的图像分割方法,但这种方法仅适用于彩色图像; 文献[7]实现了K-means算法用于分割灰度图像,当聚类的数目比较少时得到的效果比较好,但聚类中心数较多时效果会不如归一化分割; 文献[8]将自适应K-means用于X-ray图像中分割乳腺,取得了满意的结果,并将算法集成到了计算辅助决策系统。 但由于超声图像的分辨率较X-ray图像的分辨率要低,因此超声图像分割难度更大,且分割算法常需要与去噪、边缘增强等技术结合。 此外,图像处理技术中的各向异性扩散也可用于提高图像分割精度。如文献[9]提出了基于各向异性扩散的活动轮廓模型并用于分割心脏核磁共振图像,实验表明可克服噪声和伪影的干扰;文献[10]针对心脏核磁共振图像,使用各向异性滤波对图像作预处理,使用模糊C-均值聚类提高了分割灰度重叠、目标不连续和边界模糊时的分割效果。但在超声图像上结合各向异性扩散和K-means聚类进行图像分割的方法仍有待进一步研究。
B型心脏超声图像不同的像素区域对应B型心脏的壁、腔、瓣结构,像素灰度值相对大的区域对应壁或瓣结构,灰度值相对低的区域对应腔结构,且区域内部具有相对稳定的属性,像素无复杂的颜色变化。本文基于超声图像的这些先验知识,借鉴在其他图像上取得了较好分割效果的聚类思想,同时以峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)作为图像质量评价标准,结合各向异性扩散[11]和一维K-means算法[12],提出了一种新的B型心脏超声图像分割算法,称为k-gray算法。该算法不需要大量数据进行训练,即使对单张图像也可以进行分割。k-gray算法主要分为两个阶段:首先采用各向异性扩散作去噪、平滑区域内部像素、增强边缘等预处理;然后使用K-means聚类,将图像按照区域或像素值分成几个互不重叠的区域。同时本文在数学上证明了一维K-means可使得图像的PSNR值达到最大值,并且在真实心脏超声图像上与大津算法等进行了对比实验,结果显示k-gray取得了更高的PSNR值和更好的视觉分割效果。
1 基本原理 1.1 各向异性扩散对于一幅B型心脏超声图像,由于超声图像采集设备的采样点有限且采样必须满足Shannon-Nyquist采样定理,B型心脏超声图像存在模糊现象;而且由于成像设备本身的热效应等因素,B型心脏超声图像还存在大量的噪声[13]。同时观察B型心脏超声图像,还可以发现图像会有明暗相间的现象。如图 1(a)是原始B型心脏超声图像的A4C切面,采用Canny算子进行边缘检测的结果如图 1(b)所示。在分割前需要进行预处理去除噪声、平滑区域内部像素。各向异性扩散又叫Perona-Malik扩散,该方法在去除图像噪声的同时能保持图像的重要内容和一些细节信息;实验也证明各向异性扩散方法可以有效地去除噪声,同时增强边缘。
原始图像可以用方程表示为:
$\mathit{\boldsymbol{I}} = {\mathit{\boldsymbol{I}}_0}(x,y)$ | (1) |
在定义时间参数t和扩散系数r(x,y,t)的情况下,图像的扩散方程为:
${\mathit{\boldsymbol{I}}_t} = {\rm{div}}\left( {r\left( {x,y,t} \right) \cdot \Delta \mathit{\boldsymbol{I}}} \right) = r\left( {x,y,t} \right) \cdot \mathit{\boldsymbol{I}} + \nabla r \cdot \nabla \mathit{\boldsymbol{I}}$ | (2) |
初值条件为:
$\mathit{\boldsymbol{I}}\left( {x,y,t = 0} \right) = {\mathit{\boldsymbol{I}}_0}\left( {x,y} \right)$ | (3) |
其中:It为二维函数簇I(x,y,t)对参数t的偏导数;∇为梯度算子;Δ为拉普拉斯算子。
当扩散系数r(x,y,t)为常数时,扩散方程为各向同性扩散;当扩散系数r(x,y,t)变化时,扩散方程为各向异性扩散。在图像处理中,设扩散系数为图像亮度梯度的函数,即:
$\mathit{\boldsymbol{I}}(x,y,t) = g\left( {\nabla \mathit{\boldsymbol{I}}\left( {x,y,t} \right)} \right)$ | (4) |
取:
$g\left( {\nabla \mathit{\boldsymbol{I}}} \right) = C/1 + {\left( {\nabla \mathit{\boldsymbol{I}}/H} \right)^{1 + \alpha }}$ | (5) |
$\left\| {\nabla \mathit{\boldsymbol{I}}} \right\| = \sqrt {\mathit{\boldsymbol{I}}_x^2 + \mathit{\boldsymbol{I}}_y^2} $ | (6) |
其中:α、C和H是常数,且C∈[0,1],H表示梯度阈值,决定了像素是属于噪声还是真正的边缘[14]。扩散系数的取值范围为:r(x,y,t)∈[0,1]。图像的梯度Ix和Iy的计算方法可以采用Sobel算子与图像卷积的方式。
1.2 一维K-means像素聚类假设超声图像已转换为灰度图,且像素值的取值区间为[0, 255]。首先,将图像的二维坐标空间上的像素值拉成一维数组,对所有像素值按照从0到255非降排序; 然后,采用自适应的一维K-means算法对像素值进行聚类,指定聚类的最少中心数和最多中心数,根据贝叶斯信息量决定最佳聚类中心数,并在给定聚类中心数后采用动态规划策略计算出最佳聚类中心;最后输出类中心,并将图像上的每个像素值赋值为其类中心像素值大小,到此整个图像像素聚类步骤完成。 对一个一维数组进行K-means聚类的过程使用动态规划的方法来求解。
假设n个已按非降序排序的元素x1,x2,…,xn聚类为k个类时取得最小类内平方和,其类中心为μ1,μ2,…,μk,最小类内距离平方和记为D(n,k),则有:
${\mu _1} < {\mu _2} < \cdots < {\mu _k}$ | (7) |
$D\left( {n,k} \right) = \sum\limits_{i = 1}^k {\mathop \sum \limits_{{x_j} \in {c_i}} {{\left( {{x_j} - {\mu _i}} \right)}^2}} {\rm{.}}$ | (8) |
若属于第k类的元素的最小下标为m,则前m-1个元素的最小类内平方和必在聚类中心数为k-1时取得,且为D(m-1,k-1) ,否则与聚类为k时取得最小类内平方和相矛盾。于是,一维的K-means聚类问题可以使用动态规划方法来求解。动态规划方程为:
$\begin{array}{l} D\left( {n,k} \right) = \\ \mathop {\min }\limits_{k \le m \le n} \mathop {\min }\limits_{k \le m \le n} \left\{ {D\left( {m - 1k - 1} \right) + d\left( {{x_m},{x_{m + 1}}, \ldots ,{x_n}} \right)} \right\},\\ 1 \le n \le k,1 \le m \le k{\rm{.}} \end{array}$ | (9) |
其中:
$d\left( {{x_m},{x_{m + 1}}, \ldots ,{x_n}} \right) = \mathop \sum \limits_{j = m}^n {\left( {{x_j} - {\mu _k}} \right)^2}$ | (10) |
B型心脏超声图像的像素与自然图像相比具有一些不同的特点:1) 超声图像的像素通道数虽然是三通道,但是通道之间的信息关联度非常大,将三通道的超声图像转为灰度图像基本丢失的信息很少。2) 心脏超声图像的像素值具有一定的含义,像素值大的区域一般为器官组织区,呈现为图像上的亮区;而像素值小的区域多为心脏的腔室或血管通道内,表现为图像上的暗区;且同一区域内像素值相差不大,区域间像素值相差明显。3) 超声图像分割对于区域内部的细节要求不高,但是对于边界信息要求较高。
本文在各向异性预处理后使用一维K-means算法对超声图像的像素值进行聚类,从而将图像化为由几个像素中心值确定的灰阶图像,它具有几个明显的作用:1) 简化了图像上大部分细节,使得图像区域性更强;2) 将像素值分成几个灰阶后,方便后续处理,如对比度增强;3) 将图像化为几个灰阶图像后,可以从图像上看出B型超声图像的各个区域的亮度层次信息。由于经过一维K-means像素聚类后,图像的像素仅有k个灰度值,因此把此超声图像分割算法称为k-gray算法,算法流程描述如图 2所示。
图 2所示流程描述的k-gray算法主要包含两个步骤:第1步是各向异性扩散处理,直到满足条件时结束;第2步是一维K-means像素聚类,将像素聚成k个类,然后根据类别修改图像像素。此两个步骤中,一维K-means主要是对输入的向量进行聚类,得到k个类中心,k-gray的后续阶段将根据一维K-means的结果重新修改图像像素值。一维K-means算法起到两个作用:在各向异性扩散的基础上固定聚类数k时使得信噪比最大;一维K-means算法相当于对像素作一个量化,平滑图像区域内部的细节信息。
最佳聚类中心数k的确定方法采用的是贝叶斯信息准则(Baysian Information Criterion,BIC)度量,它是基于贝叶斯后验概率最大化的一种模型选择标准,它选择使得BIC具有最大值的k作为聚类中心数。
$BIC = 2\ln \left( L \right) - k\ln \left( n \right)$ | (11) |
其中:L表示图像像素拉成一维向量后的似然概率函数;k表示聚类的中心数;n表示像素点个数。
B型心脏超声图像可以看成是带有加性噪声退化了的图像,在根据先验知识对图像进行复原过程中,图像与其原始未退化图像之间的差异性越来越小,同时迭代处理前后两幅图像之间差异性也在减小,否则算法将不收敛,从而可以采用处理前后两幅图像之间差异性的指标来衡量图像质量变化,比如可以采用均方误差(Mean Square Error,MSE);进一步考虑整幅图像像素值大小,可以采用峰值信噪比(PSNR)来计算同一图像在处理前后的质量变化;且考虑到峰值信噪比较均方误差准确,较其他图像客观质量评价指标如结构相似性(Structural SIMilarity,SSIM)等简单、快速,因此本文使用峰值信噪比作为图像客观质量评价标准。均方误差和峰值信噪比的计算方法如下:
$MSE = \frac{1}{{mn}}\mathop \sum \limits_{i = 0}^{m - 1} \mathop \sum \limits_{j = 0}^{n - 1} {\left[ {{I_{i,j}} - {G_{i,j}}} \right]^2}$ | (12) |
$PSNR = 10\lg \left( {MAX_I^2/MSE} \right)$ | (13) |
其中:Ii,j表示原始图像在像素位置为(i,j)处的像素值,Gi,j表示迭代一次后图像在像素位置为(i,j)处的像素值; MAXI表示图像I的像素最大值;I、G都是灰度图像。
可以证明,一维K-means算法得到的结果就是使得峰值信噪比值最大的结果。
证明 各向异性处理得到的图像经过一维K-means聚类后,中心像素数变为k个像素,每个像素被修改为离它最近的类中心像素值,k的确定是依据模型选择标准之一的BIC。同时根据先验知识,将所有像素值小于10的像素置0,像素值大于240的像素值置为255,于是有:MAXI=255。
通过将二维图形化为一维向量后,有:
$\begin{array}{l} x = \\ \left( {{I_{0,0}},{I_{0,1}}, \ldots ,{I_{0,n - 1}},{I_{1,0}},{I_{1,1}}, \ldots ,{I_{1,n - 1}}, \ldots ,{I_{m - 1,0}},{I_{m - 1,1}}, \ldots ,{I_{m - 1,n - 1,}}} \right) \end{array}$ | (14) |
$\begin{array}{l} y = \\ \left( {{G_{0,0}},{G_{0,1}}, \ldots ,{G_{0,n - 1}},{G_{1,0}},{G_{1,1}}, \ldots ,{G_{1,n - 1}}, \ldots ,{G_{m - 1,0}},{G_{m - 1,1}}, \ldots ,{G_{m - 1,n - 1}}} \right) \end{array}$ | (15) |
y具有k个不同值的分量,记为μ1,μ2,…,μk。于是:
$\begin{array}{l} MSE = \frac{1}{{mn}}\mathop \sum \limits_{i = 0}^{m - 1} \mathop \sum \limits_{j = 0}^{n - 1} {\left[ {{I_{i,j}} - {G_{i,j}}} \right]^2} = \\ \frac{1}{L}\mathop \sum \limits_{l = 0}^{L - 1} {\left[ {{x_l} - {y_l}} \right]^2} = \mathop \sum \limits_{i = 1}^k \mathop \sum \limits_{{x_j} \in {c_{\rm{i}}}} {\left( {{x_j} - {\mu _i}} \right)^2} = \\ D\left( {n,k} \right) \end{array}$ | (16) |
式(16) 最后一项D(n,k)正是在给定k时一维K-means的目标函数式。因此,当采用动态规划策略求得D(n,k)的最小值时,均方误差也取得最小值,从而根据式(13) 可得峰值信噪比的最大值。
2.2 实现细节在对图像实现各向异性扩散预处理时,采取迭代计算的方式得到处理后的图像,各向异性扩散的迭代方程为:
$\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{I}}_{t + 1}} = {\mathit{\boldsymbol{I}}_t} + }\\ {\frac{\lambda }{4}\left[ {\sum\limits_{p \in \left\{ {N,S} \right\}} {{g_1}} \left( {\nabla {I_{s,p}}} \right) \cdot \nabla {I_{s,p}} + \sum\limits_{p \in E,W} {{g_2}} \left( {\nabla {I_{s,p}}} \right) \cdot \nabla {I_{s,p}}} \right]} \end{array}$ | (17) |
$\nabla {I_{s,p}} = {I_t}\left( p \right) - {I_t}\left( s \right);p \in \left\{ {{\rm{N,S,E,W}}} \right\}$ | (18) |
其中:系数1/4表示求4个方向的平均值;λ表示扩散速率;函数g1(∇I)和g2(∇I)表达式形式一样,差别在于常数H不同,分别对应垂直方向和水平方向上的扩散函数。B型心脏超声图像的成像区域为扇形区域,坐标系为笛卡尔坐标系,如果以扇形中心点为原点,化为极坐标并旋转90°后,可以看到像素在水平方向上的带状分布,如图 3(a)所示,此时图像在纹理特征上与地震图像有明显相似之处,地震图像的纹理结构如图 3(b)[14]所示。
因此,根据B型心脏超声图像在水平和垂直方向上的图像亮度梯度不同,上述的垂直和水平方向上的扩散函数分别采用g1(∇I)和g2(∇I)表示。
式(5) 中的参数H的确定方法有两种:1) 使用文献[15]中的knee算法计算每个方向上的图像梯度差的绝对值,然后取90%;2) 使用下文中将提到的一维K-means算法得到图像的灰度中心,计算相邻两个灰度中心差的绝对值和,然后取平均值。
λ通常取值为1。迭代次数t会严重影响各向异性扩散的效果,迭代次数过少会使图像的噪声没有被平滑,迭代次数过多又会使图像的边缘被平滑,因此本文算法每隔10次迭代计算一次峰值信噪比,在峰值信噪比取最大值时终止各向异性扩散的迭代。
在采用一维K-means算法对超声图像进行像素聚类时,有两个实现的细节:1) 将图像目标区域按比例缩放到一定大小,保证不丢失准确率的同时像素数不至于过大,比如大小为60×80,减少K-means聚类时间。2) 根据先验知识,经过各向异性扩散后,对于像素聚类中心值小于threshold的像素值,可以认为是心脏的腔结构,直接将其像素置为0;对于像素聚类中心值大于255-threshold的像素值,可以认为是心脏的壁或瓣结构,直接将其像素置为255。阈值threshold因成像设备不同而异,需事先统计得出。
3 实验结果实验算法通过C++实现,采用OpenCV进行图像的存取。输入一张B型心脏超声图像简图,添加噪声和一些像素上的干扰,然后使用迭代算法进行分割,得到的效果如图 4所示。通过人为添加的噪声来模拟B型心脏超声图像的噪声,可视化k-gray算法的两个处理步骤即各向异性扩散和一维K-means聚类的结果,可以发现k-gray算法可以有效地去除人为添加的噪声,且对k-gray算法得到的结果经过图像增强后与原图的视觉差异较小。
输入一张真实的B型心脏超声图像,使用各向异性扩散进行预处理,用K-means算法对像素进行聚类,得到的效果如图 5所示。从图中可以发现,在真实B型心脏超声图像上,k-gray算法得到的结果经过二值化后,其壁、瓣、腔结构的边界较原图 1(a)更明显。
此外,从数据库中随机抽取4张不同背景、不同亮度的真实B型心脏超声图像,分别使用人工分割、大津(Otsu)算法、区域增长算法、图割法及k-gray算法得到的效果如图 6所示。图 6中:(a)是人工分割的图像;(b)是采用OpenCV库中单阈值的大津算法分割的结果;(c)是将所有像素值大于100的像素点作为种子点的区域增长算法分割的结果;(d)是图割法分割的结果;(e)是k-gray算法分割的结果。大津算法、区域增长算法和图割法分割得到的图像均为二值图像。
从图 6中可以看出,单阈值的大津算法阈值和区域增长算法既可能出现过分割也可能出现欠分割的情形,适应性较差;图割法受背景影响较大,且分割过程中需人工指定前景和背景;k-gray算法能较好地适应不同背景、亮度的图像,且k-gray算法在二尖瓣及壁结构处出现的像素断裂区域很少,而其他三种算法在二尖瓣处出现了像素断裂及在壁结构处出现较大区域的像素不连续现象。
表 1是四种分割算法在含20张B型心脏超声图像数据库中随机抽取4张,计算峰值信噪比的统计平均值的结果。可以看出k-gray算法得到的PSNR值较大津算法提高11.5%,较区域增长算法提高11.94%,较图割算法更稳定。
从图 2的流程中和整个分割过程中可以看到,k-gray算法对输入的心脏超声图像没有任何额外的要求,分割的区域可以是任意形状、大小,参数的确定也无需事先经过大量训练;从图 5~6及表 1的结果对比中可以看到,k-gray算法分割得到的心脏超声图像质量比单阈值的大津算法、区域增长算法及图割法都要好。
4 结语B型心脏超声图像分割就是要分割出壁、腔、瓣或流结构。在对现有成果研究分析的基础上,结合B型心脏超声图像先验知识,提出了一种基于像素聚类的B型超声图像分割算法。采用峰值信噪比作为图像质量评价指标,算法在理论上可以保证峰值信噪比达到最大值,且在实际中分割得到的壁、腔、瓣结构较其他无需训练的分割算法更准确,可用于分割B型心脏超声图像中任意目标。本文算法对图像全局区域像素进行聚类,由于B型心脏超声图像各个区域性质有差异,如何将图像分成各个子区域,在子区域内使用一维K-means进行像素聚类有待进一步研究;此外,在图像质量评价指标中,有待进一步研究如何与能区别对待不同像素点的评价指标相结合。
[1] | 中华医学会麻醉学分会.围手术期经食管超声心动图监测操作的专家共识[EB/OL].(2015-02-11)[2016-05-27]. http://www.csaol.cn/a/xuehuigongzuo/linchuangzhinan/2015/0907/2555.html. ( Anesthesiology Branch of Chinese Medical Association. The expert consensus on the operation of the perioperative period of esophageal echocardiography[EB/OL]. (2015-02-11)[2016-05-27]. http://www.csaol.cn/a/xuehuigongzuo/linchuangzhinan/2015/0907/2555.html. ) |
[2] | CARNEIRO G, NASCIMENTO J C, FREITAS A. The segmentation of the left ventricle of the heart from ultrasound data using deep learning architectures and derivative-based search methods[J]. IEEE Transactions on Image Processing, 2012, 21 (3) : 968-982. doi: 10.1109/TIP.2011.2169273 |
[3] | MARSOUSI M, ALIREZAIE J, AHMADIAN A, et al. Segmenting echocardiography images using B-spline snake and active ellipse model[C]//Proceedings of the 2010 IEEE Annual International Conference on Engineering in Medicine and Biology Society. Piscataway, NJ:IEEE, 2010:3125-3128. |
[4] | MAZAHERI S, SULAIMAN P S B, WIRZA R, et al. Echocardiography image segmentation:a survey[C]//Proceedings of the 2013 International Conference on Advanced Computer Science Applications and Technologies. Piscataway, NJ:IEEE, 2013:327-332. |
[5] | 兰红, 闵乐泉. 多阈值优化交互式分割算法及其在医学图像中的应用[J]. 计算机应用, 2013, 33 (5) : 1435-1438. ( LAN H, MIN L Q. Interactive segmentation algorithm optimized by multi-threshold with application in medical images[J]. Journal of Computer Applications, 2013, 33 (5) : 1435-1438. doi: 10.3724/SP.J.1087.2013.01435 ) |
[6] | 杨明川, 吕学斌, 周群彪. 不完全K-means聚类与分类优化结合的图像分割算法[J]. 计算机应用, 2012, 32 (1) : 248-251. ( YANG M C, LYU X B, ZHOU Q B. Image segmentation algorithm based on incomplete K-means clustering and category optimization[J]. Journal of Computer Applications, 2012, 32 (1) : 248-251. ) |
[7] | TATIRAJU S, MEHTA A. Image segmentation using k-means clustering, EM and normalized cuts[EB/OL].[2016-05-20].http://www.ics.uci.edu/~dramanan/teaching/ics273a_winter08/projects/avim_report.pdf. |
[8] | PATEL B C, SINHA G R. An adaptive k-means clustering algorithm for breast image segmentation[J]. International Journal of Computer Applications, 2010, 10 (4) : 35-38. doi: 10.5120/ijca |
[9] | 祖克举, 周昌雄, 张尤赛. 基于各向异性扩散活动轮廓模型的左心室MRI分割[J]. 计算机测量与控制, 2007, 15 (3) : 339-341. ( ZU K J, ZHOU C X, ZHANG Y S. Left ventricle MRI segmentation based on anisotropic diffusion snake model[J]. Computer Measurement & Control, 2007, 15 (3) : 339-341. ) |
[10] | 曾文权, 何拥军, 崔晓坤. 基于各向异性滤波和空间FCM的MRI图像分割方法[J]. 计算机应用研究, 2014, 31 (1) : 316-320. ( ZENG W Q, HE Y J, CUI X K. MRI image segmentation method based on anisotropic diffusion and spatial FCM[J]. Application Research of Computers, 2014, 31 (1) : 316-320. ) |
[11] | PERONA P, MALIK J. Scale-space and edge detection using anisotropic diffusion[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1990, 12 (7) : 629-639. doi: 10.1109/34.56205 |
[12] | WANG H, SONG M. Ckmeans.1d.dp:optimal k-means clustering in one dimension by dynamic programming[J]. The R Journal, 2011, 3 (2) : 29-33. |
[13] | BUADES A, COLL B, MOREL J M. A review of image denoising algorithms, with a new one[J]. Multiscale Modeling & Simulation, 2005, 4 (2) : 490-530. |
[14] | TSIOTSIOS C, PETROU M. On the choice of the parameters for anisotropic diffusion in image processing[J]. Pattern Recognition, 2013, 46 (5) : 1369-1381. doi: 10.1016/j.patcog.2012.11.012 |
[15] | PETROU M, PETROU C. Image processing:the fundamentals[M]. Hoboken, NJ: John Wiley and Sons, 2010 : 549 -550. |