杜馨瑜 李永杰 尧德中

（电子科技大学神经信息教育部重点实验室，成都 610054）

引言

视觉是人类和高等动物赖以认识客观世界的主要感觉。在输入的全部感觉信息中，估计有70％以上与处理视觉信息有关［1］。一方面，现代信息技术的长足进步促进了视觉研究的深入；另一方面，对视觉基础研究的深入，又反过来对信息技术提供了新的研究思路。

借鉴视觉机制建立图像处理模型可以从以下两个方面考虑，一类是基于视觉单细胞神经元感受野电生理特性，如Marr等提出用双高斯滤波器（DOG）来解释人类视觉的底层行为，实现图像边缘提取［2］，邱芳土等对非经典感受野的三高斯模型改进用于图像增强［3］。Ghosh等引入δ函数来表示非经典感受野的高阶Laplace特性，用于图像增强，边缘检测和解释视觉心理［4］。另一类是在考虑视觉神经元系统的基础上，在视觉神经元的网络功能层次上进行图像处理建模，如脉冲耦合神经网络等。

Kayser等认为视觉系统高层次加工不仅仅是受单个细胞的影响［5，6］。因此，另一类模型在考虑视觉神经元系统的基础上，在视觉神经元的网络功能层次上进行图像处理建模。如，Munder等采用带有局部感受野的前向神经网络提取图像特征对行人进行识别［7］。Wohler等采用延时神经网络模拟感受野群的时空特性对序列图像进行分析［8］。Perez等采用前向神经网络模拟视觉系统多细胞感受野进行人脸识别和手写字体识别［9］。Sejnowski和Olshausen分别从独立成分分析和稀疏成分分析角度对视觉图像编码进行了探讨［10－11］。

采用小波或Gabor变换模拟视觉系统在以上两类模型中都有涉及。常用Gabor变换模拟视皮层单细胞神经元感受野的神经电生理特性，如朱舜山等利用Gabor变换模拟感受野，并用于图像的重建、压缩及解释稀疏编码［12］。Grigorescu等利用Gabor变换模拟感受野方位选择性进行边缘检测［13］。Gaudart等认为，Gabor变换并不能更好的协调实际视觉特性和模拟的激励响应，而应扩展单一的Gabor核到一般的小波核在视觉神经元网络层次上建模［14］，如Nuding等结合神经网络和小波分析，用非线性小波神经网络探讨视觉特性［15］。

马尔柯夫随机场（Markov random field，MRF）由于把物理学中分析现象空间与背景依赖性的概率理论，与图像处理中局部邻域像素关联性相联系而成为图像处理领域常用的模型之一。对MRF的详细描述见文献［16］。一些研究者更进一步采用多尺度MRF并在小波域进行建模。到目前为止，一些现有的小波域的马尔柯夫随机场金字塔模型，如文献［17－19］，仅仅是基于纯技术领域的图像处理模型，并没有与神经科学相联系。与之类似，也有许多模拟视觉系统的模型，但这些模型并未利用MRF［20－22］，或者这些模型并不是图像处理模型［23］。

本研究中，首先采用小波变换来模拟视皮层细胞感受野对输入图像进行稀疏表示；接者，采用多尺度MRF处理前一阶段的图像稀疏表示信息。通过配置不同的参数，模型可以完成不同的图像分割任务，如图像区域分割和边缘检测。

本研究的组织结构如下：在材料和方法部分，首先对视觉系统进行简要介绍，然后提出模型结构，并对模型中的自底向上和自顶向下两类信息流进行详细描述。在结果部分，一些实际应用中的生物医学图像分割实验展示了该模型的图像分割功能。最后，在讨论部分，对全文进行总结，并从神经科学的角度对模型进行初步的讨论。

1 材料和方法

1.1 材料

编程环境：Matlab；图像资料：采集到真实生物医学图像若干。

1.2 方法

1.2.1 视觉机制简介

视网膜是处理接收到的外界刺激信息（如图像）的第一个环节。然后，这些经视网膜预处理过的信息经由外侧膝状体传到基础视皮层（V1区）和高级视皮层进行更高层次的加工处理。视觉系统的特性可以概括为以下几点：（1）等级层次性：视觉信息是分级处理的，也就是说，视觉信息从视网膜到视皮层分级依次处理；（2）双向连接性：视觉神经元之间的连接绝大多数是双向的，即由前向连接和后向连接组成；（3）特征检测器：例如大多数的V1区神经元对具有朝向性的刺激很敏感；（4）学习机制：一般来说，视觉系统中的学习机制是自组织和非监督的；（5）选择注意机制：该特性可以使人类从复杂背景中快速发现感兴趣的目标。

1.2.2 从图像处理的角度描述小波域多尺度马尔柯夫随机场模型

许多图像分割的问题都可以看做是对图像的标记过程，也就是通过对原图像的像素或其他特征，分配一系列适当标签来完成特定的图像分割任务［16］。特别对小波域多尺度马尔柯夫随机场模型（A multiscale Markov random field model in wavelet domain，W3M）而言，用对输入图像（观测层）一层小波分解后系数（W）作为图像特征，用X代表各个隐含层的标签集，如图1（a）所示。从图像处理的角度，当给定观测层W时，图像分割问题可以转化为通过贝叶斯公式寻找关于X的最大后验概率（MAP）的问题，即

在多尺度模型中，不同尺度（金字塔的不同层）的标签集和相应的位置集分别表示为X={X（0），X（1），…，X（L）}和S={S（0），S（1），…，S（L）}，其中第0层表示金字塔的最底层，也就是最小尺度（最大分辨率）层；第L层表示金字塔的最顶层，也就是最大尺度（最小分辨率）层。从另一个方面来讲，所提出的W3M模型是一个四叉树和图的混合结构。其中四叉树表示自底向上的标记过程（如图1（b）所示），而图表示自顶向下的标记修正过程（如图1（c）所示）。从图1可以看到，图结构比四叉树稍微复杂。这是因为在图结构表示的自顶向下的过程中，3个父节点（如i，j，k）连接1个共有的子节点（如m），而在四叉树表示的自底向上的过程中，各个父节点没有共同的子节点。从总体而言，这两种结构充分利用上下文相关（邻域）信息，更真实地反应了实际视皮层神经元连接的复杂性。

图1 模型的基本结构。（a）观测层（特征提取层）W与标记层X对应关系简图；（b）自底向上过程的四叉树结构；（c）自顶向下过程的图结构（i，j，k为父节点，m为公共子节点）Fig.1 The basic structures of the model.（a）The sketch map of the observed layer W and the labeled layer X；（b）The quadtree in the bottom-up procedure；（c）The graph in the top-down procedure（The nodes i，j and k are fathers and the node m is the son）

本研究中W3M模型的四叉树与图的混合结构是建立在文献［24］所采用的结构基础上的，做了三方面的改进：一是采用小波变换以提取更为丰富的信息，而不是只采用灰度作为特征；二是W3M模型通过不同的参数设置完成不同的图像分割任务；三是赋予模型视觉神经科学方面的意义，从神经科学的角度说明W3M所表示的意义。

W3M模型的金字塔结构由L+1层组成，最上层定义为第L层，最下层定义为第0层。第n层的标签估计，根据参考文献［24］及贝叶斯公式，可写为

1.2.3 自底向上过程（图1（b）所示的四叉树结构）

式（2）右边的第一个似然函数项可以写为式（4）中，假设pw|x（0）（w|x（0））服从高斯分布。式（5）中，k为某一标签的值，M是标签的类别。θn∈［0，1］是标签转移概率，表示第n层某一节点经过与该节点同层的邻域节点相互竞争，向其共同的第n+1层父节点向上传递当前第n层该节点标签值的概率。这是尺度间多个子节点与一个父节点标签场向上的一阶Markov性的体现，表现为MRF。两种极端情况，θn=1和θn=0的意义分别为完全成功向上传递和完全失败传递。d－1（s）表示位置s处的子节点位置集合。式（5）的意义是：在第n+1层位置s处某一节点取标签值k的条件概率分布由该节点对应的所有子节点标签值的加权概率密度分布所决定，并且这些子节点的标签值可以是k，也可以是其它值。对于不同的图像分割任务，使用下式代替式（4）

图2 模型框架和计算流程。（a）模型框架（数字标号对应计算流程步骤号）；（b）计算流程Fig.2 The framework and the flowchart of the model.（a）The framework of the proposed W3M model（The number designations correspond to the steps in the computational flowchart.）（b）The computational flowchart

式中，λLL，λLH，λHL，λHH是加权系数，表示小波分解后低频、中频和高频不同分量的混合比例。这种混合可以理解为初级视皮层感受野对信号进行多种处理后在高级视皮层的整合［25］，调节不同的权重值可以完成不同的视觉任务。

1.2.4 自顶向下过程（图1（c）所示的图结构）

为了解释式（2）中的第2项，假设某一层的1个节点与3个父节点相连。αn∈［0，1］（n为当前金字塔层数）是该节点与它的某一父节点有相同标签的概率。这是尺度间多个父节点与一个子节点标签场向下的一阶Markov性的体现，表现为MRF。从另一方面，1－αn表示了该节点与它的某一父节点标签不一致的传递概率。因此，可以认为式（2）的第2项模拟了视觉系统自顶向下的反馈机制。把父节点标记为i，j，k，子节点标记为m，可以用下式来计算标签传递概率

考虑极端情况，αn=1代表了从父节点到子节点的标签完全确定性传递，各父节点的权重参考文献［24］的设置，依次取为3/7、2/7、2/7，表明假设各个父节点对该子节点的影响并不是等同的，其中一个较另外两个影响大；αn=0则表示在M类标签中完全随机传递，每一类标签都有相等的概率传递给子节点。显而易见，式（7）的意义在于W3M模型第n层位置s（对应节点m）的标签由该节点的3个父节点的标签以及模型参数αn所决定。δ为克罗内克函数。

根据以上各式，式（2）可以重新写为

在上述的过程中，对αn和θn初始值设为任意的（0，1）区间值，在自底向上和自顶向下过程中采用EM算法反复迭代。即：最初的E步骤，人为设定αn和θn值，根据底层观测场（小波分解系数）以及使用模糊聚类得到的初始底层标签场，按照上述的自底向上以及自顶向下过程执行一个循环，得到一个暂时确定的各层标签场；再进行M步骤，根据得到的各层标签场和底层观测场（小波分解系数）进行极大期望似然估计，对初始设定的αn和θn值进行修正。如此反复的进行E步骤和M步骤。因此最终的αn和θn值以及各层标签场由EM算法迭代终止时所决定，最终的分割结果由底层标签场决定。极大期望似然估计公式形式可以参见文献［24］的推导。由上所述，W3M模型的核心思想是：在金字塔最底层用小波提取输入图像特征并对这些特征进行相应的组合并初始标记；标签按四叉树结构和图结构分别进行自底向上和自顶向下的反复迭代传递；在每次迭代过程中，αn和θn被估计，并且根据式（8），图像各个位置的标签也同时被更新，最后在最底层输出图像分割结果。W3M模型的框架以及流程图分别如图2中（a）和（b）所示。

2 结果

2.1 实验

通过两组不同的实验，展示W3M模型所具有的图像区域分割和边缘检测能力，这些实验都是采用真实采集到的生物医学图像。

图3与图4是图像区域分割实验。图3（a）是待分割的显微图像。分割目标是区分出高亮荧光染色体以便于计数。这里，为了展示模型性能，采用原图像的灰度图作为待分割图像。图3（b）为采用模糊均值聚类的结果，可以看出有许多误分割。图3（c）为采用W3M模型的分割结果，可以看到所有高亮荧光染色体都得到正确的分割。这里，采用的是双正交小波（Matlab中表示为bior4.4），参数选择为：λLL=0.5、λLH=1.0、λHL=1.0、λHH=0。图4（a）是待分割的CT图像。分割的目标是骨皮质的髂骨和骶骨以及微钙化点。图4（b）是采用模糊均值聚类的结果，图4（c）是采用W3M模型分割的结果，可以看出在分割的细节方面，如微钙化点的分割，所提出的模型要好于模糊均值聚类。这里采用哈尔小波（Matlab中表示为Haar），参数选择为：λLL=1、λLH=1、λHL=1、λHH=0。

图3 显微图像分割结果。（a）原图像；（b）模糊聚类的结果；（c）W3M模型结果Fig.3 Comparison of image segmentation results on a micrograph image.（a）The input image；（b）The result of a fuzzy clustering method；（c）The result of the proposed model

图5为边缘检测实验的结果。目标是勾勒出白质边缘。上图是原始的一层磁共振扫描头部去脑壳后的图像。中间是用Canny算子边缘检测的结果，下图是W3M模型边缘检测结果。可以看出，所提出的模型在勾勒出的白质边缘上的视觉表现更为清晰。这里采用哈尔小波（Matlab中表示为Haar），参数选择为：λLL=0、λLH=1、λHL=1、λHH=1。

2.2 参数选择

图4 腹部CT图像分割结果。（a）原图像；（b）模糊聚类的结果；（c）W3M模型结果，包括骨皮质的髂骨和骶骨以及微钙化点Fig.4 Comparison of image segmentation results on an abdominal CT image.（a）The input image；（b）The result of a fuzzy clustering method；（c）The results of the proposed model are the cortical bone of iliac bone and sacrum，including the microcalcifications

图5 边缘检测结果比较（上图为头部MRI图像；中图为采用Canny算子结果；下图为采用W3M模型结果）Fig.5 Comparison of edge detection on a head MR image（Top：input image；Middle：result of Canny detector；Bottom：result of the proposed model）

通过大量的实验验证及相应的理论分析，主要选择两种小波基：双正交小波基bior4.4和哈尔小波基Harr。双正交小波bior4.4具有有限支撑、高阶消失钜，光滑性并且其对称性保证了对应的滤波器组具有线性相位，能更好地逼近原信号。而人类视觉对非线性相位的失真较为敏感，所以图像处理中的具有对称性的小波得到更广泛的应用。哈尔小波基具有对称性，并且其最短的支撑保证对图像空域的局部刻画能力强。因此关于边缘提取的实验选用Haar小波。在区域分割的实验中，由于受背景的干扰，过于细致的局部刻画能力反而会影响最终的效果，选择具有对称性，且支撑长度适中的双正交小波基bior4.4，就能达到较好的效果。

对于系数λLL、λLH、λHL、λHH（取值区间设定为［0，1］）的选择讨论如下。在图像区域分割实验中，由于分割结果以区域的形式存在，即结果要有对应于原图像平坦区域的部分，以及去除孤立干扰点，这要求我们倾向于保留原图像中能量相对集中的低频和中频成分进行分类，因此略去高频系数λHH，让水平分量系数和垂直分量系数相等（λLH=λHL=1.0），而低频系数λLL的大小取决于图像整体背景与分割目标在分割过程中的交互作用对最终分割效果的影响。低频系数λLL过小，使应该提取的部分目标误分为背景，若低频系数λLL过大，会使目标提取过度，把一部分背景误分为目标。按经验设置λLL=0.5。在图像边缘提取实验中，结果要对应于原图像的边缘区域，理所当然在参数设置中设置λLL=0.

3 讨论

本研究提出一个新的图像处理模型W3M以模拟人类某些视觉机制，如等级层次性、双向连接性、特征检测器和学习机制。把上述机制整合到一个统一的框架内，所提出的模型可以通过选择不同的参数来完成不同的图像分割任务。将模型与人类视觉机制的关系总结成表1，并作简要的讨论。

表1 视觉机制与W3M模型的关系Tab.1 The relationship between visiual mechanism and W3M

如引言中所述，视皮层细胞感受野具有局部性、方向性和带通滤波性，可以用小波变换来描述［14］。另外，视觉系统采用冗余压缩的方式来记录表示外部刺激，从而对感知到的信息提供了一种稀疏的表示［26］，而小波变换恰恰也是一种信号稀疏表示的工具［27－28］。本研究中，小波变换用作特征检测器对图像进行稀疏表示，不同的小波系数组合可以完成不同的图像分割任务，即式（6）所示，参数λ=［λLL、λHL、λLH、λHH］由不同的任务所决定。根据小波不同频段分解所代表的意义，对于图像区域分割，一般保留低频系数λLL以及中频系数λHL和λLH，设置高频系数λHH=0；对于边缘检测而言，一般设置低频系数λLL=0，保留其余系数。

以金字塔结构为基础的多尺度信息处理思想与视觉系统的等级层次性相对应。首先，视觉系统是由各层次的处理单元所组成，如视网膜、初级视皮层、高级视皮层等。并且各层之间都具有自底向上和自顶向下的双向连接。体现在所提出的W3M模型中，可以看到模型每一层都由代表层感受野的许多节点所组成，并且所有的节点都相互连接。且每一层节点所影响的范围（代表该层感受野）的大小随着金字塔结构自底向上逐层增大，这也与视觉系统中的感受野层次越高，影响范围越大相对应。

视知觉从哪里来是一个本源性问题。一种解释是直接来源于神经生理学的神经还原论［29］。该理论认为可以用作为特征检测器的滤波器组模拟感受野，对输入图像信息进行逐级提取，并抛弃冗余信息［30］。本模型采用小波作为特征检测器，并且在金字塔模型中用参数θn自底向上逐级竞争提取信息的思想与上述理论是相符合的。

另一种解释是格式塔理论的知觉组织原则。该理论强调整体知觉组织原则在视知觉中的重要性。MRF通过势团内的局部连接性（马尔柯夫性），以及实际计算中采用的等价吉布斯随机场全局势函数表现出的图像全局性，表现出视觉系统所具有的格式塔特性。MRF与格式塔心理学之间的关系也已由Zhu等建立的相关数学理论所描述［31］。本研究模型在不同的尺度下形成MRF，建立层与层之间因果关系，用参数αn表示视觉神经系统自顶向下的反馈连接特性，体现了具有更多全局信息的父节点对表现更多局部特征的子节点的控制，符合格式塔学派强调整体的观点。

另外，金字塔各层之间的侧向连接性提供了一种自学习机制［20］。正如参数θn，αn，一经初始设定，在本模型的算法中自动改变不需人为干预，体现了一种自学习特性。

4 结论

本研究提出的W3M模型在一定程度上模拟了视觉机制，通过不同的参数设置完成不同的图像分割任务，展示了在生物医学图像处理中的应用潜力。在后续研究中，将会用本模型对能够采集到的更多类型的真实生物医学图像进行测评；另外将采用更为先进的技术手段增加模型的运算速度。

［1］徐科.神经生物学纲要［M］.北京：科学出版社2005：209.

［2］Shapiro LG，Stockman GC.计算机视觉［M］.北京：机械工业出版社，2005：115.

［3］邱芳土，李朝义.同心圆感受野去抑制特性的数学模拟［J］.生物物理学报，1995，11（2）：214－220.

［4］Ghosh K，Sarkar S，Bhaumik K.Understanding image structure from a new multi-scale representation of higher order derivative filters［J］.Image Vision Comput，2007，25（8）：1228－1238.

［5］Kayser C，Einhauser W，Kongig P.Processing of complex stimuli and nature scenes in the visual cortex［J］.Current Opinion in Neurobiol，2004，14（4）：468－473.

［6］Kayser C，Einhauser W，Kongig P.Responses to Natural Secenes in Cat V1［J］.J Neurophysiol，2003，90（3）：1910－1920.

［7］Munder S，Gavrila DM.An experimental study on pedestrian classification［J］.IEEE Trans Pattern Anal Mach Intell，2006，28（11）：1863－1868.

［8］Wohler C，Anlauf JK.An adaptable time-delay eural-network algorithm for image sequence analysis［J］.IEEE Trans Neural Networ，1999，10（6）：1531－1536.

［9］Perez CA，Salinas CA，etc.Genetic design of biologically inspired receptive fields for neural pattern recognition［J］.IEEE Trans Syst Man Cy B，2003，33（2）：258－270.

［10］Bell AJ，Sejnowski TJ.The“independent components”of natural scenes are edge filters［J］.Vision Res，1997，37（23）：3327－3338.

［11］Olshausen BA，Field DJ.Sparse coding with an.overcomplete basis set：A strategy employed by V1?［J］.Vision Res，1997，37（23）：3311－3325.

［12］朱舜山，齐翔林，汪云九.基于视觉编码的图像处理研究［J］.生物物理学报，1996，12（2）：297－309.

［13］Grigorescu C，Petkov N，Westenberg MA.Contour detection based on nonclassical receptive field inhibition［J］.IEEE Trans Image Process，2003，12（7）：729－739.

［14］Gaudart L，Crebassa J，Petrakian JP.Wavelet transform in human visual channels［J］.Appl Opt，1993，32（23）：4119－4127.

［15］Nuding U，Zetzsche C.Learning the selectivity of V2 and V4 neurons using non-linear multi-layer wavelet networks［J］.Biosystems，2007，89（1－3）：273－279.

［16］Li SZ.Markov random field modeling in image analysis［M］.Berlin：Springer-Verlag，2009：1－20.

［17］NodaH，ShiraziM，KawaguchiE.MRF-based texture segmentation using wavelet decomposed images［J］.Pattern Recogn，2002，35（4）：771-782.

［18］刘国英，茅力非，王雷光，等.基于小波域分层Markov模型的纹理分割［J］.武汉大学学报（信息科学版），2009，34（5）：531－534.

［19］李旭超，朱善安，朱胜利.基于小波域层次 Markov模型的图像分割［J］.中国图象图形学报，2007，12（2）：308－314.

［20］George D，HawkinsJ.A hierarchicalBayesian modelof invariant pattern recognition in the visualcortex［C］//Prokhorov DV，eds.2005 IEEE International Joint Conference on Neural Networks（IJCNN’05）.Piscataway：IEEE Press，2005：1812－1817.

［21］Lee TS，Mumford D.Hierarchical bayesian inference in the visual cortex［J］.J Opt Soc Am，2003，20（7）：1434－1448.

［22］Serre T，Wolf L.Robust object recognition with cortex-like mechanisms［J］IEEE Trans Pattern Anal Mach Intell，2007，29（3）：411－426.

［23］Dean T.A computational model of the cerebral cortex［C］//Veloso M，Kambhampati S，eds.Proceedings of Twentieth National Conference on Artificial Intelligence（AAAI’05）.Menlo Park：AAAI Press，2005：938－943.

［24］Bouman C，Shapiro M.A Multiscale Random Field Model for Bayesian Image Segmentation［J］.IEEE Trans Image Process，1994，3（2）：162－177.

［25］尼克尔斯等著，杨雄里，等译.神经生物学——从神经元到脑（第4版）［M］.北京：科学出版社，2009：499－519.

［26］Attneave F.Some informational aspects of visual perception［J］.Psychol Revnt，1954，61：183-193.

［27］Bronstein A，Bronstein M.Sparse ICA for blind separation of transmitted and reflected images［J］.Int J Imag Syst Tech，2005，15（1）：84-91.

［28］Olshausen B，Field D.Sparse coding of sensory input［J］Curr Opin Neurobiol，2004，14（4）：481-487.

［29］Barlow HB.Single units and sensation：a neuron doctrine for perceptual psychology?［J］.Perception，1972，1（4）：371－394.

［30］Daugman J.Complete discrete 2-D Gabor transform by neural networks for image analysis and compression［J］.IEEE Trans Acoust，Speech，Signal Process，1988，36（7）：1169－1179.

［31］Zhu SC，Wu YN.From local features to global perception—A perspective of Gestalt psychology from Markov random field theory［J］.Neurocomputing，1999，26－27（6）：939－945.