CN103150573A - 基于多分辨率分形特征的神经树突棘图像分类方法 - Google Patents

Abstract

本发明公开了一种基于多分辨率分形特征的神经树突棘图像分类方法，包括以下步骤：（1）对神经树突棘图像进行特征提取获取神经树突棘图像的多分辨率分形特征；（2）采用线性判别分析（LDA）基于神经树突棘图像的多分辨率分形特征进行分类。该方法分类精度高、分类结果稳定。

Description

基于多分辨率分形特征的神经树突棘图像分类方法

技术领域

本发明属于智能图像分析领域，尤其是神经树突棘图像的高可靠性分类方法，具体涉及一种基于多分辨率分形特征的神经树突棘图像分类方法。

背景技术

近年来，细胞显微镜成像技术得到了很快发展。脑神经元成像对于理解神经元形态和神经元移动机制至关重要，尤其在对很多疾病的研究中，探查神经元对复杂感知和运动过程的动态响应是一个重要手段，脑神经元成像在这些研究中体现出了其无法替代的重要地位。毫无疑问，显微镜成像已成为突破各种神经疾病研究的重要工具。

尽管在显微镜成像技术上已经取得了很多突破，而与之相对应的图像计算分析方法的发展却相对滞后。例如，对于治疗神经损伤退化和紊乱疾病，如阿兹海默病和帕金森病，分析神经突生长是一个很有前景的方向，然而，人工对神经突增长进行分析是非常耗费人力和时间的，海量的图像使得人工分析不现实。因此，自动神经元图像分析技术的发展至关重要，而分类已成为自动神经元图像分析的一个总要组成部分。如何有效地区分各种类型的神经元细胞这个课题近年来吸引了大量的来自不同领域的研究人员，其中包括很多来自计算机视觉和人工智能领域的研究人员。

视网膜神经节细胞承担着把视觉信号传递到大脑的任务，获取和分析视网膜神经细胞的不同形态和动态响应机制对于理解它们的功能而言至关重要。研究人员在对猴子和猫的视网膜神经细胞图像的分析表明，分形方法对于图像中神经元结构的描述和区分具有很好的效果。有研究结果表明，计算大鼠视网膜神经元图像中的全局分形维数能够提供对细胞进行分类和功能判别的依据。本发明因此而来。

发明内容

本发明目的在于提供一种基于多分辨率分形特征的神经树突棘图像分类方法，解决了现有技术中图像分类效果差、分类精度差等问题。

为了解决现有技术中的这些问题，本发明提供的技术方案是：

一种基于多分辨率分形特征的神经树突棘图像分类方法，其特征在于所述方法包括以下步骤：

（1）对神经树突棘图像进行特征提取获取神经树突棘图像的多分辨率分形特征；

（2）采用线性判别分析（LDA）基于神经树突棘图像的多分辨率分形特征进行分类：通过已分类的神经树突棘图像的多分辨率分形特征进行最大似然估计（Maximum likelihood estimation）得到每一类图像特征的高斯分布参数，包括：先验概率、均值和协方差矩阵；当新的神经树突棘图像加入时，通过计算新的神经树突棘图像的图像特征和各个类的高斯分布的马氏距离（Mahalanobis distance）予以分类，新的神经树突棘图像将划分到与其马氏距离最近的类别当中。

优选的，所述方法步骤（1）中特征提取的步骤包括：

A)读取神经树突棘图像并将其转换为256色灰度图像作为初始分辨率的灰度图像；

B)使用OTSU算法确定二值化阈值T，把灰度图像转换为二值黑白图像（Binary image）；使用细化算法从灰度图像得到骨架图像（Skeleton image）；

C)分形维数可以作为度量图像纹理粗糙度的一种度量指标，盒子维或“盒计数”算法是通过划分图像形成网格，统计出网格中包含的盒子数。

本发明采用差分盒子维法，具体说明如下：

将M×M大小的图像分割成S×S的子块（M/2≥s）1,s为整数），令R=s/m。将图像作为三维空间中的曲面，x,y表示平面位置，z轴表示灰度值。xy平面被分割成许多ss的网格。在每个网格上，是一列s×s×s的盒子。设图像灰度在第(i,j)网格中的最小值和最大值分别落在第k和第l个盒子中，则：

nr(i,j)=l-k+1是覆盖第(i,j)网格中的图像所需的盒子数，而覆盖整个图像所需的盒子数为Nr而

N_r=∑n_r(i,j)；

分形维数为

$D=\lim \left(\frac{\log \left(N_r\right)}{\log (1/r)}\right);$

针对不同的r，计算Nr，应用最小二乘法，即可求得分形维数D。

D)通过插值算法，由当前分辨率缩放系数r从256色灰度图像得到下一分辨率的灰度图像，继续进行步骤B）～D)，依次循环获得所有分辨率下的二值图像和骨架图像的分形维数FD；

E）连接该幅神经树突棘图像在所有分辨率下的二值图像和骨架图像的分形维数构成图像特征描述。

优选的，所述方法使用OTSU算法确定二值化阈值T的方法是：

Otsu算法的基本思想是用某一假定的灰度值t将图像的灰度分成两组，当两组的类间方差最大时，此灰度值t就是图像二值化的最佳阈值。

设图像有M个灰度值，取值范围在0～M-1，在此范围内选取灰度值t，将图像分成两组G₀和G₁，G₀包含的像素的灰度值在0～t，G₁的灰度值在t+1～M-1，用N表示图像像素总数，n_i表示灰度值为i的像素的个数。

已知：每一个灰度值i出现的概率为p_i＝n_i/N；假设G₀和G₁两组像素的个数在整体图像中所占百分比为W₀,W₁，两组平均灰度值为u₀,u₁，可得: $W_0=\underset{i=1}{\overset{t}{\mathrm{\Σ}}}{p}_i$ $W_1=\underset{i=t+1}{\overset{M-1}{\mathrm{\Σ}}}{p}_i=1-W_0$

$u_0=\underset{i=1}{\overset{t}{\mathrm{\Σ}}}{\mathrm{ip}}_i$ $u_1=\underset{i=t+1}{\overset{M-1}{\mathrm{\Σ}}}{\mathrm{ip}}_i$

图像的总平均灰度为：u＝W₀×u₀+W₁×u₁

类间方差为：g(t)＝W₀(u₀-u)²+W₁(u₁-u)²＝W₀W₁(u₀-u₁)²

最佳阈值为：T=argmax(g(t))，即使得类间方差g(t)最大时所对应的t值。

优选的，所述方法细化算法的方法是：

细化算法是采用迭代的手段不断剥除目标的边界点，保留边界上满足骨架特性的点，直到没有多余边界可剥除，剩下的就是骨架。本发明选择中轴变换细化算法。

中轴变换实际上是在不影响原图拓扑性的基础上，通过抽取对表达原图形状最关键的点，使得原图中宽度大于1个像素的线条变成单像素的过程。

中轴变换通常利用二值形态操作判断像素点的八邻域的情况。在文献2，中轴变化算法充分考虑了上下左右4个方向，通过两个步骤、8个条件来判断是否删除中心像素点，其算法过程描述如下：

(1)分析图像任意3×3区域，若同时满足以下条件，消除边界点：

a.周围目标像素数在2个到6个之间(包括2和6)；

b.中心像素周围的8个像素中同时存在背景像素和目标像素(保证当前处于边界位置)，背景像素有且只有1个连通，目标像素也有且只有1个连通；

c.中心像素的上邻、左邻、下邻像素至少有一个为背景像素；

d.中心像素的左邻、下邻、右邻像素至少有一个为背景像素。

(2)在将所有满足步骤1中4个条件的边界点去除后，再进行进一步的删除，分析图像任意3×3区域，若同时满足以下条件，消除边界点：

a.周围目标像素数在2个到6个之间(包括2和6)；

b.中心像素周围的8个像素中同时存在背景像素和目标像素(保证当前处于边界位置)，背景像素有且只有1个连通，目标像素也有且只有1个连通：

c.中心像素的上邻、左邻、右邻像素至少有一个为背景像素：

d.中心像素的上邻、下邻、右邻像素至少有一个为背景像素。

经过反复迭代处理直至结果不再发生变化(达到收敛)位置，中轴变换结束。从上述的限制条件可以看出，迭代过程充分考虑了消除线段端点、过度深入区域内部以及割断骨架等多种可能引起误处理的情况。

优选的，所述方法中插值算法的步骤是

本发明采用双三次插值算法生成多分辨率图像。双三次插值又称立方卷积插值。该算法利用待采样点周围16个点（4×4邻域）的灰度值作三次插值，不仅考虑到4个直接相邻点的灰度影响，而且考虑到各邻点间灰度值变化率的影响。

当计算一个新的像素点位置（i+u,j+v）亮度时，双三次插值公式如下:

f(i+u，j+v)＝ABC

其中，A、B，C均为矩阵，形式为：

A＝[S(1+u)S(u)S(1-u)S(2-u)]；

$B=\left[\begin{array}{cccc}f(i-1,j-2)& f(i,j-2)& f(i+1,j-2)& f(i+2,j-2)\\ f(i-1,j-1)& f(i,j-1)& f(i+1,j-1)& f(i+2,j-1)\\ f(i-1,j)& f(i,j)& f(i+1,j)& f(i+2,j)\\ f(i-1,j+1)& f(i,j+1)& f(i+1,j+1)& f(i+2,j+1)\end{array}\right];$

C＝[S(1+v)S(v)S(1-v)S(2-v)]^T。

其中，f(i,j)为原图像在(i，j)位置处的亮度值；S(w)为插值核，其函数为：

$S\left(w\right)=\left\{\begin{array}{c}1-2{\left|w\right|}^2+{\left|w\right|}^3,\left|w\right|<1\\ 4-8\left|w\right|+5{\left|w\right|}^2-{\left|w\right|}^3,1\&le;\left|w\right|<2\\ 0,\left|w\right|\&GreaterEqual;2\end{array}\right..$

优选的，所述方法中通过最大似然估计获得先验概率、均值和协方差矩阵的方法是：

使用高斯模型描述图像特征数据，混合模型被定义为M个高斯密度函数的线性组合： $P\left(x\right)=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_i{N}_i(x;{\mathrm{\&mu;}}_i,{\mathrm{\&Sigma;}}_i)$

其中，N_i(x;π_i,∑_i)是均值为μ_i，协方差为∑_i的高斯分布，π_i是混合参数，是第i个高斯分布的权重，表征先验概率，且：

且0≤π_i≤1

N_i(x；π_i,∑_i)的概率密度函数为：

$N_i\left(x\right)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^{d/2}{\left|{\mathrm{\&Sigma;}}_i\right|}^{1/2}}\exp \{-\frac{1}{2}{(x-{\mathrm{\&mu;}}_i)}^T{\mathrm{\&Sigma;}}_i^{-1}(x-{\mathrm{\&mu;}}_i)\}$

将高斯混合密度函数中所有待定的参数记为θ，则似然函数为：

$P\left(X\right|\mathrm{\&theta;})=\underset{i=1}{\overset{N}{\mathrm{\&Pi;}}}P\left(x_i\right|\mathrm{\&theta;})\&DoubleRightArrow;{\mathrm{\&theta;}}^{*}=\underset{\mathrm{\&theta;}}{\arg \max }P\left(X\right|\mathrm{\&theta;})$

对于每一个观察到的x，用Z(z|1)来表示每一个高斯分布，因z是未知变量，即无法知道x属于哪个分布。因此，使用EM（Expectation-Maximization）算法得到分布参数：

首先计算似然函数：

$P(X,Z|\mathrm{\&theta;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p(x_i,z_i|\mathrm{\&theta;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left({\mathrm{\&pi;}}_{i}N(x_i;{\mathrm{\&mu;}}_{z_i},{\mathrm{\&Sigma;}}_{z_i})\right)$

关于条件分布p(Z|X,θ^old)的期望：

Q(θ，θ^old)＝E[log p(X,Z|θ)|X,θ^old]

其次，是将似然函数最大化以获得新的参数值，使用θ^new更新θ^old，使QYθ|θ^oldY最大化，即θ^new＝arg max QYθ|θ^oldY。

重复上两个步骤直到迭代停止条件达到，即充分小时算法结束，此时得到的参数集合θ为最终确定的高斯分布参数。

优选的，所述方法中马氏距离的计算方法是：

设经线性判别和最大似然估计后得到的某一数据类别G（m维）有均值向量为μ＝(μ₁,μ₂,K,μ_m)^t,协方差矩阵为∑＝(σ_ij)，则样品X＝(x₁,x₂,K,x_m)^t与类G的马氏距离定义为：

d²(X，G)＝(X-μ)′∑^-1(X-μ)

当m=1时，

$d^2(x,G)=\frac{{(x-\mathrm{\&mu;})}^{\&prime;}(x-\mathrm{\&mu;})}{{\mathrm{\&sigma;}}^2}=\frac{{(x-\mathrm{\&mu;})}^2}{{\mathrm{\&sigma;}}^2}$

本发明基于多分辨率分形特征的神经树突棘图像分析和分类方法属于智能图像分析领域，尤其是神经树突棘图像的分析研究，解决了细胞神经元图像的特征描述和分类问题。本发明采用的线性判别分析（LinearDiscriminant Analysis），简称判别分析，是统计学上的一种分析方法，用于在已知的分类之下遇到有新的样本时，选定一个判别标准，以判定如何将新样本放置于哪一个类别之中。

本发明依据视网膜神经细胞原型和变异型的树状形态以及轴突投影模式对其进行分类。发明人在深入研究细胞神经元图像的分形特性基础上，提出了基于多分辨率的分形图像特征描述方法，并利用获取的神经元描述特征对神经元进行分类，获得了较好的分类效果。目前，尚未出现基于多分辨率分形的神经元图像特征提取方法。

相对于现有技术中的方案，本发明的优点是：

本发明技术方案提供了一种基于多分辨率分形特征的神经树突棘图像分类方法，包括以下步骤：（1）对神经树突棘图像进行特征提取获取神经树突棘图像的多分辨率分形特征；（2）采用线性判别分析（LDA）基于神经树突棘图像的多分辨率分形特征进行分类。该方法分类精度高、分类结果稳定。

附图说明

下面结合附图及实施例对本发明作进一步描述：

图1为本发明对神经树突棘图像二值化结果，左图为原图像，右图为二值化结果；

图2为本发明对神经树突棘图像提取骨架结果；左图为原图像，右图为提取的图像骨架；

图3为本发明方法具体实施流程图；

图4为本发明基于多分辨率分形的图像特征提取后的处理结果图。

图5为采用多分辨率的分形特征对测试图像库的分类精度比较结果。

图6为本发明（FLD）与其他分类器的分类稳定度比较结果。

具体实施方式

以下结合具体实施例对上述方案做进一步说明。应理解，这些实施例是用于说明本发明而不限于限制本发明的范围。实施例中采用的实施条件可以根据具体厂家的条件做进一步调整，未注明的实施条件通常为常规实验中的条件。

实施例

本实施例基于多分辨率分形特征的神经树突棘图像分类方法，包括以下步骤：

（1）对神经树突棘图像进行特征提取获取神经树突棘图像的多分辨率分形特征；

（2）采用线性判别分析（LDA）基于神经树突棘图像的多分辨率分形特征进行分类：通过已分类的神经树突棘图像的多分辨率分形特征进行最大似然估计（Maximum likelihood estimation）得到每一类图像特征的高斯分布参数，包括：先验概率、均值和协方差矩阵；当新的神经树突棘图像加入时，通过计算新的神经树突棘图像的图像特征和各个类的高斯分布的马氏距离（Mahalanobis distance）予以分类，新的神经树突棘图像将划分到与其马氏距离最近的类别当中。

本发明具体实施流程图可见图3。

如图3所示，本实施例的具体实施细节为：

（1）将彩色图像转换为256色灰度图像，使用如下公式完成转换：

GRAY=0.2989×R+0.5870×G+0.1140×B。

（2）使用OTSU算法对256色灰度图像进行二值化：

算法的基本思想是用某一假定的灰度值t将图像的灰度分成两组，当两组的类间方差最大时，此灰度值t就是图像二值化的最佳阈值。

设图像有M个灰度值，取值范围在0～M-1，在此范围内选取灰度值t，将图像分成两组G₀和G₁，G₀包含的像素的灰度值在0～t，G₁的灰度值在t+1～M-1，用N表示图像像素总数，n_i表示灰度值为i的像素的个数。

已知：每一个灰度值i出现的概率为p_i＝n_i/N；假设G₀和G₁两组像素的个数在整体图像中所占百分比为W₀,W₁，两组平均灰度值为u₀,u₁，可得:

$W_0=\underset{i=1}{\overset{t}{\mathrm{\&Sigma;}}}{p}_i$ $W_1=\underset{i=t+1}{\overset{M-1}{\mathrm{\&Sigma;}}}{p}_i=1-W_0$

$u_0=\underset{i=0}{\overset{t}{\mathrm{\&Sigma;}}}{\mathrm{ip}}_i$ $u_1=\underset{i=t+1}{\overset{M-1}{\mathrm{\&Sigma;}}}{\mathrm{ip}}_i$

图像的总平均灰度为：u＝W₀×u₀+W₁×u₁

类间方差为：g(t)＝W₀(u₀-u)²+W₁(u₁-u)²＝W₀W₁(u₀-u₁)²

最佳阈值为：T=argmax(g(t))，即使得类间方差g(t)最大时所对应的t值。图像二值化结果可见图1。

（3）使用细化算法对256色灰度图像进行骨架图像提取：

细化算法是采用迭代的手段不

的边界点，保留边界上满足骨架特性的点，直到没有多余边界可剥除，剩下的就是骨架。本发明选择中轴变换细化算法。

中轴变换实际上是在不影响原图拓扑性的基础上，通过抽取对表达原图形状最关键的点，使得原图中宽度大于1个像素的线条变成单像素的过程。

中轴变换通常利用二值形态操作判断像素点的八邻域的情况。在文献2，中轴变化算法充分考虑了上下左右4个方向，通过两个步骤、8个条件来判断是否删除中心像素点，其算法过程描述如下：

1)分析图像任意3 3区域，若同时满足以下条件，消除边界点：

a.周围目标像素数在2个到6个之间(包括2和6)；

b.中心像素周围的8个像素中同时存在背景像素和目标像素(保证当前处于边界位置)，背景像素有且只有1个连通，目标像素也有且只有1个连通；

c.中心像素的上邻、左邻、下邻像素至少有一个为背景像素；

d.中心像素的左邻、下邻、右邻像素至少有一个为背景像素。

2)在将所有满足步骤1中4个条件的边界点去除后，再进行进一步的删除，分析图像任意3×3区域，若同时满足以下条件，消除边界点：

a.周围目标像素数在2个到6个之间(包括2和6)；

b.中心像素周围的8个像素中同时存在背景像素和目标像素(保证当前处于边界位置)，背景像素有且只有1个连通，目标像素也有且只有1个连通：

c.中心像素的上邻、左邻、右邻像素至少有一个为背景像素：

d.中心像素的上邻、下邻、右邻像素至少有一个为背景像素。

经过反复迭代处理直至结果不再发生变化(达到收敛)位置，中轴变换结束。从上述的限制条件可以看出，迭代过程充分考虑了消除线段端点、过度深入区域内部以及割断骨架等多种可能引起误处理的情况。

图像细化实例可见图2。

（4）盒计数算法的具体说明：

分形维数可以作为度量图像纹理粗糙度的一种度量指标，盒子维或“盒计数”算法是通过划分图像形成网格，统计出网格中包含的盒子数。

本实施例采用差分盒子维法，具体说明如下：

将M×M大小的图像分割成S×S的子块（M/2≥s）1,s为整数），令R=s/m。将图像作为三维空间中的曲面，x,y表示平面位置，z轴表示灰度值。xy平面被分割成许多S×S的网格。在每个网格上，是一列S×S×S的盒子。设图像灰度在第(i,j)网格中的最小值和最大值分别落在第k和第l个盒子中，则：

n_r(i,j)=l-k+1是覆盖第(i,j)网格中的图像所需的盒子数，而覆盖整个图像所需的盒子数为N_r而

N_r=∑n_r(i,j)

分形维数为

$D=\lim \left(\frac{\log \left(N_r\right)}{\log (1/r)}\right)$

针对不同的r，计算Nr，应用最小二乘法，即可求得分形维数D。

在实际应用中，对原始图像使用长度为r的盒子对图像进行网格分割，r的取值分别为4,16和64，计算每个不同r大小时图像中目标（二值化和细化后白色部分）所覆盖的盒子数目N(r)。分形维数可使用最小二乘法由公式（1）得到。

（5）由插值算法得到多分辨率图像

本实施例采用双三次插值算法生成多分辨率图像。双三次插值又称立方卷积插值。该算法利用待采样点周围16个点（4×4邻域）的灰度值作三次插值，不仅考虑到4个直接相邻点的灰度影响，而且考虑到各邻点间灰度值变化率的影响。

当计算一个新的像素点位置（i+u,j+v）亮度时，双三次插值公式如下:

f(i+u，j+v)＝ABC

其中，A、B，C均为矩阵，形式为：

A＝[S(1+u)S(u)S(1-u)S(2-u)]；

$B=\left[\begin{array}{cccc}f(i-1,j-2)& f(i,j-2)& f(i+1,j-2)& f(i+2,j-2)\\ f(i-1,j-1)& f(i,j-1)& f(i+1,j-1)& f(i+2,j-1)\\ f(i-1,j)& f(i,j)& f(i+1,j)& f(i+2,j)\\ f(i-1,j+1)& f(i,j+1)& f(i+1,j+1)& f(i+2,j+1)\end{array}\right];$

C＝[S(1+v)S(v)S(1-v)S(2-v)]^T。

其中，f(i,j)为原图像在(i,j)位置处的亮度值；S(w)为插值核，其函数为：

$S\left(w\right)=\left\{\begin{array}{c}1-2{\left|w\right|}^2+{\left|w\right|}^3,\left|w\right|<1\\ 4-8\left|w\right|+5{\left|w\right|}^2-{\left|w\right|}^3,1\&le;\left|w\right|<2\\ 0,\left|w\right|\&GreaterEqual;2\end{array}\right..$

重复步骤2-5，直到完成所有分辨率下图像的分形维数特征提取。

连接一幅图像在所有分辨率下的二值、骨架分形维数构成图像特征描述，由线性判别分析进行分类。

Fisher线性判别(HD)的基本原理如下，对于线性判别函数可以将d维矢量视作特征空间X_d中的以a₁，a₂，…,a_d为分量的一个矢量，则a^Tx表示矢量x在以a为方向的轴上投影的‖a‖倍。希望所求的a使投影后同类样本相距较近，即同类样本密集；不同类样本相距较远。FLD是求解满足类间离散度和总类内离散度之比最大的投影方向，然后在一维空间中确定判决规则。设给定两类d维训练样本各类样本均值矢量m_i和总的样本均值矢量m分别为：

$m_i=\frac{\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ij}}}{N_i},i=\mathrm{1,2}$

$m=\frac{\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N_i}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ij}}}{\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}{N}_i}$

设各类类内离散度矩阵S_i（i=1,2）,总的类内离散度阵S_w=S₁+S₂。类间离散度阵S_b为

S_b＝(m₁-m₂)(m₁-m₂)^T

定义类问离散度与类内离散度之比为Fisher准则函数

$J_F\left(a\right)=\frac{a^T{S}_{b}a}{a^T{S}_{w}a}$

并使其最大。将标量J_F对矢量a求导并令其为零矢量，然后利用二次型关于矢量求导公式可得：

$a^{*}=S_w^{-1}(m_1-m_2)$

此时的a*可使Fisher准则函数取最大值，即最佳投影方向。

本发明使用判别分析中常用的马氏距离分类方法：

设经线性判别和最大似然估计后得到的某一数据类别G（m维）有均值向量为μ＝(μ₁,μ₂,K,μ_m)^t,协方差矩阵为∑＝(σ_ij)，则样品X＝(x₁,x₂,K,x_m)^t与类G的马氏距离定义为：

d²(X，G)＝(X-μ)′∑^-1(X-μ)

当m=1时，

$d^2(x,G)=\frac{{(x-\mathrm{\&mu;})}^{\&prime;}(x-\mathrm{\&mu;})}{{\mathrm{\&sigma;}}^2}=\frac{{(x-\mathrm{\&mu;})}^2}{{\mathrm{\&sigma;}}^2}$

样本与哪一类别的马氏距离最近，将被划分为此类。

本实施例具体应用时，得到的图像特征如图1所示，图1的左列为各分辨率256色灰度神经元图像，中列为OTSU阈值化处理的二值图像，右列为经细化算法得到的骨架图像。将其应用于一个具有2个类别的视网膜神经元图像测试库的参数设定如下：

初始图像为1024×1024大小的RGB彩色图像。采用双三次方插值算法得到多分辨率图像，得到的新图像分辨率为：1024×1024、768×768、512×512和256×256；对于每个分辨率的图像，采用使用OTSU算法自动确定图像二值化阈值；同时，采用细化算法得到骨架图像。各分辨图像及其二值化、骨架提取结果实例可见图4。

应用盒计数方法对于二值化图像和骨架图像得到当前分辨率下的分形维数；每幅图像各分辨率下的分形维数连接形成对图像的特征描述用于分类。

其中分类器采用Fisher线性判别分析方法。图6为分类精度的比较，本发明采用的线性判别分析方法（FLD）较其他分类器取得了较好的分类精度结果。其他比较的分类器为：k近邻（kNN）、逻辑回归（Logistic）、多层感知器（MLP）和支持向量机（SVM）。

图6中可以看到使用线性判别分析对神经树突棘图像进行分类，分类准确率为85%，而使用kNN、逻辑回归、多层感知器和支持向量机的分类精度分别为63%、76%、65%和80%，使用线性判别分析方法对神经树突棘图像分类取得了较好的分类精度。

同时，为验证本发明的分类稳定性，对各分类器的分类稳定程度进行了比较，结果如图5所示。可以看出，本发明采用的方法得到了较好的稳定性。图中矩形框代表了各分类器在100次分类实验中的稳定程度，矩形框越小代表分类器分类越稳定，可以看到，线性判别分析分类器在实验中获得了较好的分类稳定性。同时，矩形框的在纵坐标的位置越高，代表平均分类精度越高，可以看到，线性判别分析在保证稳定分类的前提下能同时保证分类精度。

上述实例只为说明本发明的技术构思及特点，其目的在于让熟悉此项技术的人是能够了解本发明的内容并据以实施，并不能以此限制本发明的保护范围。凡根据本发明精神实质所做的等效变换或修饰，都应涵盖在本发明的保护范围之内。

Claims (7)

1.一种基于多分辨率分形特征的神经树突棘图像分类方法，其特征在于所述方法包括以下步骤：

（1）对神经树突棘图像进行特征提取获取神经树突棘图像的多分辨率分形特征；

（2）采用线性判别分析（LDA）基于神经树突棘图像的多分辨率分形特征进行分类：通过已分类的神经树突棘图像的多分辨率分形特征进行最大似然估计（Maximum likelihood estimation）得到每一类图像特征的高斯分布参数，包括先验概率、均值和协方差矩阵；当新的神经树突棘图像加入时，通过计算新的神经树突棘图像的图像特征和各个类的高斯分布的马氏距离（Mahalanobis distance）予以分类，新的神经树突棘图像将划分到与其马氏距离最近的类别当中。

2.根据权利要求1所述的神经树突棘图像分类方法，其特征在于所述方法步骤（1）中特征提取的步骤包括：

A)读取神经树突棘图像并将其转换为256色灰度图像作为初始分辨率的灰度图像；

B)使用OTSU算法确定二值化阈值T，把灰度图像转换为二值黑白图像（Binary image）；使用细化算法从灰度图像得到骨架图像（Skeleton image）；

C)使用盒计数算法分别计算两种图像的分形维数：

将M×M区域大小的图像分割成S×S的子块，其中M/2≥s，1,s均为整数，将图像作为三维空间中的曲面，x,y表示平面位置，z轴表示灰度值；xy平面被分割成许多S×S的网格，在每个网格上，是一列S×S×S的盒子；设图像灰度在第(i,j)网格中的最小值和最大值分别落在第k和第l个盒子中，则n_r(i,j)=l-k+1是覆盖第(i,j)网格中的图像所需的盒子数，而覆盖整个图像所需的盒子数N_r为：

N_r=∑n_r(i,j)；

针对不同的r，计算N_r，应用最小二乘法，即可求得分形维数D，所述分形维数D为：

$D=\lim \left(\frac{\log \left(N_r\right)}{\log (1/r)}\right)---\left(1\right);$

D)通过插值算法，由当前分辨率缩放系数r从256色灰度图像得到下一分辨率的灰度图像，继续进行步骤B）～D)，依次循环获得所有分辨率下的二值图像和骨架图像的分形维数；

E）连接该幅神经树突棘图像在所有分辨率下的二值图像和骨架图像的分形维数构成图像特征描述。

3.根据权利要求2所述的神经树突棘图像分类方法，其特征在于所述方法使用OTSU算法确定二值化阈值T的方法是：

设图像有M个灰度值，取值范围在0～M-1，在此范围内选取灰度值t，将图像分成两组G₀和G₁，G₀包含的像素的灰度值在0～t，G₁的灰度值在t+1～M-1，用N表示图像像素总数，n_i表示灰度值为i的像素的个数；假设已知每一个灰度值i出现的概率为p_i＝n_i/N；G₀和G₁两组像素的个数在整体图像中所占百分比为W₀,W₁，两组平均灰度值为u₀,u₁，可得：

$W_0=\underset{i=0}{\overset{t}{\mathrm{\&Sigma;}}}{p}_i;$ $W_1=\underset{i=t+1}{\overset{M-1}{\mathrm{\&Sigma;}}}{p}_i=1-W_0;$ $u_0=\underset{i=0}{\overset{t}{\mathrm{\&Sigma;}}}{\mathrm{ip}}_i;$ $u_1=\underset{i=t+1}{\overset{M-1}{\mathrm{\&Sigma;}}}{\mathrm{ip}}_i;$

则图像的总平均灰度为u＝W₀×u₀+W₁×u₁；

类间方差为g(t)＝W₀(u₀-u)²+W₁(u₁-u)²＝W₀W₁(u₀-u₁)²；

最佳阈值为T=argmax(g(t))，即使得类间方差g(t)最大时所对应的t值。

4.根据权利要求2所述的神经树突棘图像分类方法，其特征在于所述方法细化算法的方法是选用中轴变换细化算法，其算法过程包括：

(1)分析图像任意3×3区域，若同时满足以下条件，消除边界点：

a.周围目标像素数在2个到6个之间；

b.中心像素周围的8个像素中同时存在背景像素和目标像素，保证当前处于边界位置，背景像素有且只有1个连通，目标像素也有且只有1个连通；

c.中心像素的上邻、左邻、下邻像素至少有一个为背景像素；

d.中心像素的左邻、下邻、右邻像素至少有一个为背景像素；

(2)在将所有满足步骤1中4个条件的边界点去除后，再进行进一步的删除，分析图像任意3×3区域，若同时满足以下条件，消除边界点：

a.周围目标像素数在2个到6个之间；

b.中心像素周围的8个像素中同时存在背景像素和目标像素，保证当前处于边界位置，背景像素有且只有1个连通，目标像素也有且只有1个连通；

c.中心像素的上邻、左邻、右邻像素至少有一个为背景像素；

d.中心像素的上邻、下邻、右邻像素至少有一个为背景像素；

经过上述步骤反复迭代处理直至结果不再发生变化即达到收敛位置，中轴变换结束。

5.根据权利要求2所述的神经树突棘图像分类方法，其特征在于所述方法中插值算法选用双三次插值算法生成多分辨率图像，通过利用待采样点周围16个点（4×4邻域）的灰度值作三次插值，当计算一个新的像素点位置（i+u,j+v）亮度时，双三次插值公式为：

f(i+u,j+v)=ABC；

其中，A、B，C均为矩阵，形式为：

A=[S(1+u)S(u)S(1-u)S(2-u)]；

$B=\left[\begin{array}{cccc}f(i-1,j-2)& f(i,j-2)& f(i+1,j-2)& f(i+2,j-2)\\ f(i-1,j-1)& f(i,j-1)& f(i+1,j-1)& f(i+2,j-1)\\ f(i-1,j)& f(i,j)& f(i+1,j)& f(i+2,j)\\ f(i-1,j+1)& f(i,j+1)& f(i+1,j+1)& f(i+2,j+1)\end{array}\right];$

C＝[S(1+v)S(v)S(1-v)S(2-v)]^T；

其中，f(i,j)为原图像在(i，j)位置处的亮度值；S(w)为插值核，其函数为：

$S\left(w\right)=\left\{\begin{array}{c}1-2{\left|w\right|}^2+{\left|w\right|}^3,\left|w\right|<1\\ 4-8\left|w\right|+5{\left|w\right|}^2-{\left|w\right|}^3,1\&le;\left|w\right|<2\\ 0,\left|w\right|\&GreaterEqual;2\end{array}\right..$

6.根据权利要求2所述的神经树突棘图像分类方法，其特征在于所述方法中通过最大似然估计获得均值、协方差矩阵，先验概率的方法是使用高斯模型描述图像特征数据，混合模型P(x)被定义为M个高斯密度函数的线性组合；其中：

$P\left(x\right)=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_i{N}_i(x;{\mathrm{\&mu;}}_i,{\mathrm{\&Sigma;}}_i);$

其中，N_i(x;π_i,∑_i)是均值为μ_i，协方差为∑_i的高斯分布，π_i是混合参数，是第i个高斯分布的权重，表征先验概率，且：

N_i(x;π_i,∑_i)的概率密度函数为：

$N_i\left(x\right)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^{d/2}{\left|{\mathrm{\&Sigma;}}_i\right|}^{1/2}}\exp \{-\frac{1}{2}{(x-{\mathrm{\&mu;}}_i)}^T{\mathrm{\&Sigma;}}_i^{-1}(x-{\mathrm{\&mu;}}_i)\};$

将高斯混合密度函数中所有待定的参数记为θ，则似然函数为：

$P\left(X\right|\mathrm{\&theta;})=\underset{i=1}{\overset{N}{\mathrm{\&Pi;}}}P\left(x_i\right|\mathrm{\&theta;})\&DoubleRightArrow;{\mathrm{\&theta;}}^{*}=\underset{\mathrm{\&theta;}}{\arg \max }P\left(X\right|\mathrm{\&theta;});$

对于每一个观察到的x，用Z(z|1)来表示每一个高斯分布；使用EM（Expectation-Maximization）算法得到分布参数：

首先计算似然函数：

$P(X,Z|\mathrm{\&theta;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log p(x_i,z_i|\mathrm{\&theta;})=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log \left({\mathrm{\&pi;}}_{i}N(x_i;{\mathrm{\&mu;}}_{z_i},{\mathrm{\&Sigma;}}_{z_i})\right);$

关于条件分布p(Z|X,θ^old)的期望：

Q(θ，θ^old)＝E[log p(X,Z|θ)|X,θ^old]；

其次，是将似然函数最大化以获得新的参数值，使用θ^new更新θ^old，使QYθ|θ^oldY最大化，即θ^new＝arg max QYθ|θ^oldY；

重复上述两个步骤直到迭代停止条件达到，即充分小时算法结束，此时得到的参数集合θ为最终确定的高斯分布参数。

7.根据权利要求2所述的神经树突棘图像分类方法，其特征在于所述方法中马氏距离的计算方法是：

设经最大似然估计后得到的某一数据类别G（m维）有均值向量为μ＝(μ₁,μ₂,K,μ_m)^t,协方差矩阵为∑＝(σ_ij)，则样品X＝(x₁,x₂,K,x_m)^t与类G的马氏距离定义为：

d²(X，G)＝(X-μ)′∑^-1(X-μ)；

当m=1时，

$d^2(x,G)=\frac{{(x-\mathrm{\&mu;})}^{\&prime;}(x-\mathrm{\&mu;})}{{\mathrm{\&sigma;}}^2}=\frac{{(x-\mathrm{\&mu;})}^2}{{\mathrm{\&sigma;}}^2}.$

Images