时域与频域的图像增强及Matlab实现

中 国 矿 业 大 学

本科生毕业设计

姓 名：学 院：专 业：设计题目：专 题：指导教师： 陶士昌 学 号： 08083576 计算机科学与技术 电子信息科学与技术

时域与频域的图像增强及Matlab实现 图像处理 梁志贞 职 称： 副教授

2012年 6月 徐州

中国矿业大学毕业设计任务书

学院 计算机科学与技术 专业年级 信科08-3 学生姓名 陶士昌

任务下达日期：2012年 1 月 10 日

毕业设计日期： 2012年 1 月 4 日 至 2012 年 6 月10日

毕业设计题目： 时域与频域的图像增强及Matlab实现

毕业设计专题题目：图像处理

毕业设计主要内容和要求：

毕业设计（论文）的目的是对毕业生所学的专业基础知识和研究能力、自学能力以及各种综合能力的检验，要进一步巩固和加强学生基本知识的掌握和基本技能的训练，加强对学生的多学科理论、培养刻苦钻研、勇于探索的精神。

图像增强（image enhancement）是一种按特定的需要突出一幅图像的某些信息，同时，削弱或去除某些不需要的信息，从而有目的地强调图像整体或局部特征，加强图像判读和识别效果的处理方法。其主要目的是使处理后的图像对某些特定的应用比原来的图像更加有效。

常用的图像增强技术主要包括直方图修改处理、图像平滑化处理、图像锐化处理等。在实际应用中，常常是几种方法联合处理，以便达到预期的增强效果。

从纯技术上讲，图像增强技术基本上可分为两大类：频域增强法和空域增强法。空域增强法是直接对图像中的像素进行处理，基本上是以灰度映射变换为基础的。频域增强法的基础是卷积定理，关键在于图像空域与频域的转换类型

本课题首先要介绍了图像增强方面发展的状况和常用图像增强的基本理论,然后根据所学知识对已有的算法用Matlab编程实现。

院长签字： 指导教师签字：

中国矿业大学毕业设计指导教师评阅书

指导教师评语（①基础理论及基本技能的掌握；②独立解决实际问题的能力；③研究内容的理论依据和技术方法；④取得的主要成果及创新点；⑤工作态度及工作量；⑥总体评价及建议成绩；⑦存在问题；⑧是否同意答辩等）：

成 绩： 指导教师签字： 年 月

日中国矿业大学毕业设计评阅教师评阅书

评阅教师评语（①选题的意义；②基础理论及基本技能的掌握；③综合运用所学知识解决

实际问题的能力；③工作量的大小；④取得的主要成果及创新点；⑤写作的规范程度；⑥总体评价及建议成绩；⑦存在问题；⑧是否同意答辩等）：

成 绩： 评阅教师签字： 年 月日

中国矿业大学毕业设计答辩及综合成绩 答 辩 情 况 回 答 问 题 提 出 问 题 有一有原正 基本 没有 般性则性确 正确 回答 错误 错误 答辩委员会评语及建议成绩： 答辩委员会主任签字： 年 月 日 学院领导小组综合评定成绩： 学院领导小组负责人： 年 月 日

摘 要

在图像的采集、处理、存储、显示或传输的过程中，由于受到多种因素的影响，如光电系统失真、噪声干扰、曝光不足或过量、相对运动、传输误码等，往往使图像与原始景物之间或图像与原始图像之间产生某种差异，从而引起图像的降质或退化。降质或退化的图像通常模糊不清，使人观察起来不满意，或者使机器从中提取的信息减少甚至造成错误。因此必须对图像进行增强。图像增强的方法基本上可分为空间域方法和频率域方法两大类。空间域方法是在原图像上（空间域）直接进行数据运算，对像素的灰度值进行处理。频率域方法是在图像的变换域上进行处理，增强感兴趣部分的频率分量，然后进行反变换得到增强后的图像。

论文在介绍图像增强原理的基础上，利用MATLAB工具进行了算法的设计和实现。实验证明在质量较差的图像中，选择不同的算法对图像增强有不同的效果。

关键词：图像增强；时域；频域；MATLAB

ABSTRACT

In the process of collection, processing, storage, display or transmission of the image, the images have some differences with the original scene or the original image, because of many factors such as PV systems distortion, noise, underexposure, overexposure, relative motion, transmission error, etc. These results are called image degradation or degeneration. The degraded or degenerated image is usually fuzzy. It makes us feel unsatisfied or the machine fail to extract full information and even results in errors. So, the image enhancement is needed. The image enhancement methods are basically divided into two categories: spatial domain and frequency domain. Spatial domain methods directly deal with the original image data and the operations are performed on the pixel gray value. The frequency domain methods are explored on the transformation domain of the image processing to enhance the interesting part of the frequency components and the enhanced image is obtained by performing the inverse transformation.

Based on the principle of the image enhancement, this paper introduces some algorithms and implements them by the MATLAB tools. The experiment proves that different methods will give different results when processing low quality image.

Keywords: image enhancement; spatial domain; frequency domain; MATLAB; image processing

目 录

1绪论 ................................................................................................................................................................... 1 1.1课题研究的背景和意义 ............................................................................................................................... 1

1.2 图像增强的国内外研究情况 .............................................................................................................. 2

1.2.1 图像增强的国外发展情况 ...................................................................................................... 2 1.2.2图像增强技术国内发展状况 ............................................................................................... 3

2 时域图像增强方法 .......................................................................................................................................... 4

2.1 引言 ...................................................................................................................................................... 4 2.2 时域增强的定义和步骤 ...................................................................................................................... 4 2.3 灰度变换 .............................................................................................................................................. 4

2.3.1 灰度扩展（缩减） .................................................................................................................. 4 2.3.2 分段线性变换 .......................................................................................................................... 5 2.3.3 非线性变换 .............................................................................................................................. 5 2.3.4 MATLAB函数及示例 ............................................................................................................... 6 2.4 直方图修正 .......................................................................................................................................... 6

2.4.1 图像的直方图 .......................................................................................................................... 6 2.4.2 直方图均衡化 .......................................................................................................................... 8 2.4.3 MATLAB函数及实现 ................................................................................................................. 9 2.5 空域滤波增强 ..................................................................................................................................... 11

2.5.1 空域滤波基本原理和分类 ..................................................................................................... 11 2.5.2 平滑滤波器 ............................................................................................................................ 12 2.5.3 锐化滤波器 ............................................................................................................................ 16

3 频域图像增强方法 ........................................................................................................................................ 20

3.1 二维离散傅里叶变换（DFT）简介 .................................................................................................. 20 3.2频域增强定义和步骤 ......................................................................................................................... 21 3.3 低通滤波 ............................................................................................................................................ 21

3.3.1 理想低通滤波器 .................................................................................................................... 21 3.3.2 Butterworth 低通滤波器 .................................................................................................... 22 3.3.3 指数低通滤波器 .................................................................................................................. 22 3.3.4 梯形低通滤波器 .................................................................................................................... 22 3.3.5 MATLAB算法及其实现 ........................................................................................................... 22 3.4高通滤波 ............................................................................................................................................. 26

3.4.1 理想高通滤波器 .................................................................................................................... 26 3.4.2 Butterworth高通滤波器 ..................................................................................................... 26 3.4.3 指数高通滤波器 .................................................................................................................... 27 3.4.4 MATLAB算法及其实现 ........................................................................................................... 27

4 结论与对未来展望 ........................................................................................................................................ 31 参考文献 ............................................................................................................................................................ 33 英文原文 ............................................................................................................................................................ 34 中文译文 ............................................................................................................................................................ 44 致 谢 ................................................................................................................................................................ 52

中国矿业大学2012届本科生毕业论文 第1页

1绪论

1.1课题研究的背景和意义

随着人类社会的进步，科学技术的发展，人们对信息处理和信息交流的要求越来越高。图像信息具有直观、形象、易懂和信息量大等特点，它是在人们日常生活、生产中接触最多的信息种类之一。近年来，随着信息社会数字化的发展，数字图像处理（digital image processing）无论是在理论研究方面还是在实际应用方面都取得了长足的发展。尤其是计算机技术的应用、遥感技术的发展、数字通信的兴起、互联网的普及、数字处理芯片性能的提高以及应用数学理论与方法的更新，对数字信号处理起了关键性的推动作用；而数字图像处理技术的发展又有力地促进和加速了上述各项技术的进步[1]。

人类传递信息的主要媒介是语言和图像。有研究表明，大约有70%的信息是通过人眼获得图像的图像信息，所以图像信息是十分重要的信息传递媒体和方式。在当代科学研究、军事技术、航天、气象、医学、工农业生产等领域中，人们越来越多地利用图像信息来认识和判断事物解决实际问题。例如，人们利用人造卫星所拍摄的地面照片来分析地球资源、气象态势和污染情况，利用“神州”宇宙飞船所拍摄的月球表面照片来分析月球的地形；在医学上，通过计算机断层图像，医生可以观察和诊断人体内部是否有病变组织；在公安侦破案件中，采用指纹图像提取和比对来识别罪犯；在军事上，目标的自动识别和跟踪都有赖于高速图像处理。

图像传递系统包括图像采集、图像压缩、图像编码、图像存储、图像通信、图像显示这六个部分。在实际应用中每个部分都有可能导致图像品质变差，使图像传递的信息无法被正常读取和识别。例如，在采集图像过程中由于光照环境或物体表面反光等原因造成图像整体光照不均，或是图像采集系统在采集过程中由于机械设备的缘故无法避免的加入采集噪声，或是图像显示设备的局限性造成图像显示层次感降低或颜色减少等等，往往使图像与原始景物或者原始图像之间产生某种差异，常将这种差异称之为降质或退化。降质或退化的图像通常模糊不清，使人观察起来不满意，或者使机器从中提取的信息减少甚至造成错误。因此，必须对图像进行改善。

改善的方法主要包括图像增强和图像复原。图像增强是从人的主观出发，可以不考虑图像降质的原因，只将图像中感兴趣的部分加以处理或突出有用的图像特征，故改善后的图像并不一定要去逼近原图像，如增加图像的对比度、提取图像中目标物轮廓、衰减各类噪声、均衡图像灰度等。图像复原技术是从客观出发，针对图像降质的具体原因，设法补偿降质因素，从而使改善后的图像尽可能地逼近原始图像。显然，图像复原的主要目的是提高图像的逼真度，这里我不做讨论和研究。

图像增强处理的方法基本上可分为时域（空间域）方法和频率域方法两大类。空间域方法是在原图像上（时域）进行数据运算，对像素的灰度值进行处理。如果是对图像作逐点运算，称为点运算，如果是在像素点邻域内进行运算，称为局部运算或邻域运算。频率域方法是在图像的变换域上进行处理，增强感兴趣部分的频率分量，然后进行反变换，得到增强后的图像。

中国矿业大学2012届本科生毕业论文 第2页

1.2 图像增强的国内外研究情况

1.2.1 图像增强的国外发展情况

在1921年年底提出了一种基于光学还原的新技术。在这一时期由于引入了一种用编码图像纸带去调制光束达到调节底片感光程度的方法，使灰度等级从5个灰度级增加到15个灰度等级，这种方法明显改善了图像复原的效果。到20世纪60年代早期第一台可以执行数字图像处理任务的大型计算机制造出来了，这标志着利用计算机技术处理数字图像时代的到来。1964年，研究人员在美国喷气推进实验室(JPL)里使用计算机以及其它硬件设备，采用几何校正、灰度变换、去噪声、傅里叶变换以及二维线性滤波等增强方法对航天探测器―徘徊者7号‖发回的几千张月球照片进行处理，同时他们也考虑太阳位置和月球环境的影响，最终成功地绘制出了月球表面地图。随后他们又对1965年―徘徊者8号‖发回地球的几万张照片进行了较为复杂的数字图像处理，使图像质量进一步提高。这些成绩不仅引起世界许多有关方面的注意而且JPL本身也更加重视对数字图像处理地研究和设备的改进，并专门成立了图像处理实验室IPL。在IPL里成功的对后来探测飞船发回的几十万张照片进行了更为复杂的图像处理，最终获得了月球的地形图、彩色图以及全景镶嵌图。从此数字图像增强技术走进了航空航天领域。

20世纪60年代末和20世纪70年代初有学者开始将图像增强技术用于医学图像、地球遥感监测和天文学等领域。X射线是最早用于成像的电磁辐射源之一，在1895年X射线由伦琴发现。20世纪70年代Godfrey N. Hounsfield先生和Allan M. Cormack教授共同发明计算机轴向断层技术：一个检测器围绕病人，并用X射线源绕着物体旋转。X射线穿过身体并由位于对面环中的相应检测器收集起来。其原理是用感知的数据去重建切片图像。当物体沿垂直于检测器的方向运动时就产生一系列的切片，这些切片组成了物体内部的再现图像。到了20世纪80年代以后，各种硬件的发展使得人们不仅能够处理二维图像，而且开始处理三维图像。许多能获得三维图像的设备和分析处理三维图像的系统已经研制成功了，图像处理技术得到了广泛的应用。进入20世纪90年代，图像增强技术已经逐步涉及人类生活和社会发展的各个方面。计算机程序用于增强对比度或将亮度编码为彩色，以便解释X射线和用于工业、医学及生物科学等领域的其他图像。地理学用相同或相似的技术从航空和卫星图像中研究污染模式。在考古学领域中使用图像处理方法已成功地复原模糊图片。在物理学和相关领域中计算机技术能增强高能等离子和电子显微镜等领域的实验图片。直方图均衡处理是图像增强技术常用的方法之一。1997年Kim 提出如果要将图像增强技术运用到数码相机等电子产品中，那么算法一定要保持图像的亮度特性。在文章中Kim提出了保持亮度特性的直方图均衡算法(BBHE)。Kim的改进算法提出后，引起了许多学者的关注。在1999年Wan等人提出二维子图直方图均衡算法(DSIHE)。接着Chen和Ramli提出最小均方误差双直方图均衡算法(MMBEBHE)。为了保持图像亮度特性，许多学者转而研究局部增强处理技术，提出了许多新的算法：递归均值分层均衡处理(RMSHE)、递归子图均衡算法(RSIHE)、动态直方图均衡算法(DHE)、保持亮度特性动态直方图均衡算法(BPDHE)、多层直方图均衡算法(MHE)、亮度保持簇直方图均衡处理(BPWCHE)等等[2]。

中国矿业大学2012届本科生毕业论文 第3页

1.2.2图像增强技术国内发展状况

在借鉴国外相对成熟理论体系和技术应用体系的条件下，国内的增强技术和应用也有了很大的发展。总体来说，图像增强技术的发展大致经历了初创期、发展期、普及期和应用期4个阶段。初创期开始于20世纪60年代，当时的图像采用像素型光栅进行扫描显示，大多采用中、大型机对其进行处理。在这一时期由于图像存储成本高，处理设备造价高，因而其应用面很窄。20世纪70年代进入了发展期，开始大量采用中、大型机进行处理，图像处理也逐渐改用光栅扫描显示方式，特别是出现了CT和卫星遥感图像，对图像增强处理提出了一个更高的要求。到了20世纪80年代，图像增强技术进入普及期，此时的计算机已经能够承担起图形图像处理的任务。20世纪90年代进入了应用期，人们运用数字图像增强技术处理和分析遥感图像，以有效地进行资源和矿藏的勘探、调查、农业和城市的土地规划、作物估产、气象预报、灾害及军事目标的监视等。在生物医学工程方面，运用图像增强技术对X射线图像、超声图像和生物切片显微图像等进行处理，提高图像的清晰度和分辨率。在工业和工程方面，主要应用于无损探伤、质量检测和过程自动控制等方面。在公共安全方面，人像、指纹及其他痕迹的处理和识别，以及交通监控、事故分析等都在不同程度上使用了图像增强技术。图像增强是图像处理的重要组成部分，传统的图像增强方法对于改善图像质量发挥了极其重要的作用。随着对图像技术研究的不断深入和发展，新的图像增强方法不断出现。例如一些学者将模糊映射理论引入到图像增强算法中，提出了包括模糊松弛、模糊熵、模糊类等增强算法来解决增强算法中映射函数选择问题，并且随着交互式图像增强技术的应用，可以主观控制图像增强效果。同时利用直方图均衡技术的图像增强也有许多新的进展：例如提出了多层直方图结合亮度保持的均衡算法、动态分层直方图均衡算法。这些算法通过分割图像，然后在子层图像内做均衡处理，较好地解决了直方图均衡过程中的对比度过拉伸问题，并且可以控制子层灰度映射范围，增强效果较好[3]。

中国矿业大学2012届本科生毕业论文 第4页

2 时域图像增强方法

2.1 引言

空域变换增强根据对图像的每次处理是对单个像素进行的或是对小的子图像(模板)进行的，可分为2组：基于像素(点)的和基于模板的。在基于像素的处理(也叫点处理)中，增强过程对每个像素的处理与其他像素无关；而模板处理则是指每次处理操作都是基于图像中的某个小区域进行的。

2.2 时域增强的定义和步骤

时域方法是在原图像上（空间域）直接进行数据运算，对像素的灰度值进行处理。用公式描述如下：G?x,y??F?x,y?H?x,y?，其中是F?x,y?原图像；H?x,y?为空间转换函数；G?x,y?表示进行处理后的图像。

2.3 灰度变换

灰度变换，或灰度级修正是在空间域对图像进行增强的一种简单而有效的方法，可以根据对图像不同的要求而采用不同的修正方法。灰度级修正属于（像素）点运算，它不改变原图像中像素点的位置，只改变像素点的灰度值，并且逐点进行，和周围其他像素无关。设输入图像为f?x,y?，经过变换后输出图像为g?x,y?，修正函数或变换函数为T???，则有

g(x,y)?T?f?x,y?? （2.1）

通过不同的映射变换，达到不同的灰度修正的效果。

一般成像系统只具有一定的亮度响应范围，亮度的最大值与最小值之比称为亮度对比度。由于成像系统的限制，常会出现对比度不足（或过大）的弊病，使人眼观察图像时视觉效果很差。因此需要调整对比度，即灰度变换。有三种常用的灰度变换：线性、分段线性及非线性变换。 2.3.1 灰度扩展（缩减）

这是一种最简单的灰度变换，假定原图像f?x,y?的灰度范围为?a,b?，希望变换后图像g?x,y?的灰度范围扩展至?c,d?，则线性变换可表示为

d?c??gx,y? ? f ? x , y ? ? a ? ? c f ( x , y ) ? ? a , b ? （2.2）

b?a

此式可用图2.1(a) 所示。

中国矿业大学2012届本科生毕业论文 第5页

2.3.2 分段线性变换

为了突出感兴趣的目标或灰度区间，相对抑制那些不感兴趣的灰度区域，可采用分段线性变换。常用的是3段线性变换法，其数学表达式如下 c?f?x,y?0?f?x,y??a?a ??d?c? x? ?xy? g ? , ? ? ? f ? x , y ? ? a ? ? c a f , y b （2.3） ?b?a? ?G?d?f?x,y??b??db?f?x,y??F?F?b?

式（2.3）对灰度区间?o,a?和?b,M?进行压缩，对灰度区间?a,b?进行扩展，如图2.1（b）。

通过调整折现拐点的位置及控制分段线性的斜率，可对任一灰度区间进行扩展或压缩。这种变换适用于在黑色或白色附近有噪声干扰的情况[4]。 g ? x ? ,yg?x,y?

G G d d c c o f(x,y)f(x,y) ab o F aF b

(a)简单线性变换 （b）分段线性变换 g?x,y?

G

c

o f(x,y)F （c） 非线性对苏变换

图2.1 灰度范围的变换 2.3.3 非线性变换

当用某些非线性变换函数，如用对数函数、指数函数等作为式（2.1）的变换函数时，

可实现图像灰度的非线性变换。例如，对数变换的一般式为

[(x,y)?1] （2.4） g(x,y)?c?klnf

这里的k、c是用于调整曲线位置和形状的参数，其图形如图2.1（c）所示。它使f?x,y? 的灰度范围得到扩展，高灰度范围得到压缩，以使图像的灰度分布与人的视觉特性相匹配。

中国矿业大学2012届本科生毕业论文 第6页

指数变换的一般式为

g(x,y)?bc?f?x,y??a??1 （2.5） 式中a、b、c 3个参数也是用来调整曲线的位置和形状，它的效果与对数变换相反，对图像的高灰度区给予较大的扩展。 2.3.4 MATLAB函数及示例

在MATLAB中，函数imadjust(I, [low_in high_in],[low_out high_out]) 来实现图像的灰度值变换。low_in 和high_in 参数分别用来指定输入图像需要映射的灰度范围，low_out和high_out指定输出图像所在的灰度范围。代码如下[5]

I=imread('cameraman.tif'); J=imadjust(I,[0 0.2],[0.5 1]); imshow(I);

figure,imshow(J) 结果如图2.2所示

(a) 原图 (b) 灰度变换后图

图2.2 灰度变换示例

2.4 直方图修正

2.4.1 图像的直方图

直方图是图像中像素灰度分布统计特性的一种图形表示方式，可定义为 H?k??nk 或 P?k??nk/N （归一化） （2.6） 式中k为灰度级；nk为第k级灰度的像素数；N为一幅图像的总像素数；H?k?和P?k?分别为该图像的直方图和归一化直方图（频数）。图2.3是图片cameraman.tif及其直方图统计

和数字图像类似，模拟图像（连续图像）的直方图在概念上也表示图像灰度的分布情况类似于概率密度。

中国矿业大学2012届本科生毕业论文 第7页

设连续图像为f?x,y?，D表示某一灰度，A?D? 表示图像中灰度大于D的面积（称之为面积函数），则该图像的直方图H?D?可以用它的面积函数的微分来定义，即

limA?D??D??A?D?dA?D? （2.7） H?D????D?0dDdD

(a) cameraman原图 (b) 直方图统计

图2.3 cameraman图像及直方图

(a) 位置一 (b) 位置一直方图

(c) 位置二 (d) 位置二直方图

图2.4 不同位置的图像及直方图

中国矿业大学2012届本科生毕业论文 第18页

-1 -2 -1 -1 0 1 0 0 0 -2 0 2 1 2 1 -1 0 1 (a)对水平边缘响应大 （b）对垂直边缘响应大

图2.13 Sobel算子基本方法

式（2.19）和式（2.20）分别对应图2.13所示的两个滤波模板。为了简化运算，也可以用来 x ? S y 替代式（2.21）的计算，从而得到锐化后的图像。由此可知，Sobelg? S算子不像普通梯度算子那样用两个像素的差值，而用两列或者两行加权和的差值，这就有一下两个优点：

① 由于引入了平均因素，因而对图像中的随机噪声有一定的平滑作用。 ② 由于它是相隔两行或两列的差分，故边缘两侧的元素得到了增强，边缘显得粗而亮。

（a）原图 （b）sobel算子滤波

（c）prewitt算子滤波 （d）log算子滤波

图2.14三种边缘增强算子效果

在MATLAB中常用空域微分算子sobel算子、prewitt算子、高斯-拉普拉斯算子等来实现非线性锐化滤波。下面的例子来显示几种边缘增强算子在图像增强中的不同效果。 I=imread('eight.tif');

figure,imshow(I);

中国矿业大学2012届本科生毕业论文 第19页

h1=fspecial('sobel'); K1=filter2(h1,I); figure,imshow(K1); h2=fspecial('prewitt'); K2=filter2(h2,I); figure,imshow(K2); h3=fspecial('log'); K3=filter2(h3,I); figure,imshow(K3);

中国矿业大学2012届本科生毕业论文 第20页

3 频域图像增强方法

3.1 二维离散傅里叶变换（DFT）简介

在连续信号的处理中，傅里叶变换为人们深入理解和分析信号的特性提供了一种强有力的手段。为了能进行定量数值计算，可以通过抽样使原来连续分布的信号变成离散信号。并且由抽样定理可以知道，当满足一定条件时，就可以从有限的抽样精确地恢复出原连续信号，由此大大地降低了分析处理的复杂程度。对于离散信号，与之相对应的是离散傅里叶变换[11]。

离散傅里叶变换建立了离散时域（空间域）与离散频域之间的联系。当信号或图像处理在空域上直接处理时，计算量会随着取样点数的增加而极具增加。采用DFT方法可减少计算量：将输入的数字信号首先进行DFT，把时域（空域）中的卷积或相关运算简化为在频域上的相乘处理，然后进行DFT反变换，恢复为时域（空域）的信号。这样，计算量减少，提高了处理速度。此外，DFT还存在快速算法，可以进一步减少计算量，提高处理速度。

定义二维离散信号{f(x,y)|x?0,1,?,M?1;y?0,1,?,N?1}的DFT的变换为

?uxvy?M?1N?1?j2????1 F ? , y ? e ? M N ? （3.1）? ? f ?xu ,vMNx?0y?0

?? ?uxvy?M?1N?1j2????1 f ? x , y ? ? F (u ,v ? M N ? （3.2） )eMNu?0v?0

式中，u,x??0,1,?,M?1?;v,y??0,1,?,N?1?。

在DFT变换对中，F(u,v)为离散信号f(x,y)的频谱，一般情况下是复数，设F?u,v?为其幅度谱，??u,v?为其相位谱，F(u,v)还可用式（3.3）表达，即

F(u,v)?F(u,v)exp?j??u,v???R(u,v)?jI(u,v) （3.3）

式中R(u,v)和I(u,v)分别为二维复频谱F(u,v)函数的实部和虚部[12]。

F(u,v)的幅度函数为

F?u,v??R2?u,v??I2?u,v? （3.4）

??

除了二维离散傅里叶变换，还有DCT变换、沃尔什变换、哈达玛变换、正交的K-L变换等。都是时域到频域的变换方法。受篇幅问题，这里不再做介绍了。

?I?u,v?? ? ? u , ? ? （3.5） ?arctanv

?R?u,v??F(u,v)的相位函数为

中国矿业大学2012届本科生毕业论文 第21页

3.2频域增强定义和步骤

频域增强是利用图像变换方法将原来的图像空间中的图像以某种形式转换到其它空间中，然后利用该空间的特有性质方便地进行图像处理，最后再转换回原来的图像空间中，从而得到处理后的图像。 频域增强的主要步骤是：

(1) 选择变换方法，将输入图像变换到频域空间；

(2) 在频域空间中，根据处理目的设计一个转移函数并进行处理； (3) 将所得结果用反变换得到图像增强。 频域增强的方法有两个关键点：

(1) 将图像从空间域转换到频域所需的变换及将图像从频域空间转换回空间域所需的变换；

(2) 在频域空间对图像进行增强的操作。

根据适用范围的不同，常用的频域增强方法有低通滤波和高通滤波。

3.3 低通滤波

图像在传递过程中，由于噪声主要集中在高频部分，为去除噪声改善图像质量，滤波器采用低通滤波器H?u,v?来抑制高频成分，通过低频成分，然后再进行逆傅立叶变换获 得滤波图像，就可达到平滑图像的目的。在傅里叶变换域中，变换系数能反映某些图像的特征，如频谱的直流分量对应于图像的平均亮度，噪声对应于频率较高的区域，图像实体位于频率较低的区域等，因此频域常被用于图像增强。在图像增强中构造低通滤波器，使低频分量能够顺利通过，高频分量有效地阻止，即可滤除该领域内噪声。由卷积定理，低通滤波器数学表达式为：G?u,v??F?u,v?H?u,v?式中，F?u,v?为含有噪声的原图像的傅里叶变换域；H?u,v?为传递函数；G?u,v?为经低通滤波后输出图像的傅里叶变换。假定噪声和信号成分在频率上可分离，且噪声表现为高频成分。H?u,v?滤波滤去了高频成分，而低频信息基本无损失地通过。

选择合适的传递函数H?u,v?对频域低通滤波关系重大。常用频率域低通滤波器H?u,v?有四种[13]：理想低通滤波器、Butterworth 低通滤波器、指数低通滤波器、梯形低通滤波器。

3.3.1 理想低通滤波器

设傅立叶平面上理想低通滤波器离开原点的截止频率为D0，则理想低通滤波器的传递函数为：

D?u,v??D0?1 H ?u , v ? ? ? （3.6）

??0Du,v?D0?

式中，D?u,v??u2?v2

??1/2表示点?u,v?到原点的距离，D0表示截止频率点到原点的距离。

中国矿业大学2012届本科生毕业论文 第22页

3.3.2 Butterworth 低通滤波器

n 阶Butterworth 滤波器的传递函数为[14]：

1??Hu,v? 2n?D?u,v?? （3.7） 1????D0?

Butterworth低通滤波器是一种物理上可以实现的滤波器。它的特性是连续性衰减，而不像理想滤波器那样陡峭变化。 3.3.3 指数低通滤波器

指数低通滤波器是图像处理中常用的另一种平滑滤波器。它的传递函数为：

nD?u,v?? （3.8） D0, v ? H u e

??3.3.4 梯形低通滤波器

梯形低通滤波器是理想低通滤波器和完全平滑滤波器的折中。它的传递函数为：

?1D?u,v??D0 ??D?u,v??D1??Hu,v?D? D 0 ? ? u ,v ? ? D 1 （3.9）

D?D01? ?0D?u,v??D0?

3.3.5 MATLAB算法及其实现

（1）理想低通滤波器[15]

图3.1理想低通滤波器特性曲线

中国矿业大学2012届本科生毕业论文 第23页

理想是指小于D0的频率可以完全不受影响地通过滤波器，而大于D0的频率则完全通不过，如公式（3.1）所示。

给出H的1个剖面图（设D对原点对称），这里理想是指小于D0的频率可以完全不受影响地通过滤波器，而大于D0的频率则完全通不过。因此D0也叫截断频率。尽管理想低通滤波器在数学上定义得很清楚，在计算机模拟中也可实现，但在截断频率处直上直下的理想低通滤波器是不能用实际的电子器件实现的。

如果使用这些“非物理”的理想滤波器，其输出图像会变得模糊和有“振铃（ring）”现象出现。我们可借助卷积定理解释如下。

为简便，考虑1-D的情况。对1个理想低通滤波器，其h(x)的一般形式可由求式（3.1）的傅立叶反变换得到，其曲线可见图3.2（a）。现设f(x)是1副只有1个亮像素的简单图像，见图3.2（b）。这个亮点可看作是1个脉冲的近似。在这种情况下，f(x)和h(x)的卷积实际上是把h(x)复制到f(x)中亮点的位置。比较图3.2（b）和图3.2（c）可明显看出卷积使原来清晰的点被模糊函数模糊了。对更为复杂的原始图，如我们认为其中每个灰度值不为零的点都可以看作是1个其值正比于该点灰度值的1个亮点，则上述结论仍可成立。

图3.2空间模糊示意图

由图3.2还可以看出h(x,y)在2-D图像平面上将显示出一系列同心圆环。如对1个理想低通滤波器的H(u,v)求反变换，则可知道h(x,y)中同心圆环的半径是反比于D0的值的。所以如果D0较小，就会使h(x,y)产生数量较少但较宽的同心圆环，并使g(x,y)模糊得比较厉害。当增加D0时，就会使h(x,y)产生数量较多但较窄的同心圆环[16]，并使g(x,y)模糊得比较少。如果D0超出F(u,v)的定义域，则h(x,y)在其对应的空间区域为1，h(x,y)与f(x,y)的卷积仍是f(x,y)，这相当于没有滤波。

I1=imread('saturn.png');

I2=imnoise(I1,'salt & pepper'); imshow(I1);

figure,imshow(I2); s=fftshift(fft2(I2));

[M,N]=size(s); %分别返回s的行数到M中，列数到N中 n1=floor(M/2); %对M/2进行取整

n2=floor(N/2); %对N/2进行取整

中国矿业大学2012届本科生毕业论文 第24页

%ILPF滤波，d0=5，15，30 （程序中以d0=5为例）

d0=5; %初始化d0

for i=1:M for j=1:N

d=sqrt((i-n1)^2+(j-n2)^2); %点（i,j）到傅立叶变换中心的距离 if d<=d0 %点（i,j）在通带内的情况 h=1; %通带变换函数

else %点（i,j）在阻带内的情况 h=0; %阻带变换函数 end

s(i,j)=h*s(i,j); %ILPF滤波后的频域表示

end end

s=ifftshift(s); %对s进行反FFT移动

%对s进行二维反离散的Fourier变换后，取复数的实部转化为无符号8位整数 s=uint8(real(ifft2(s)));

figure; %创建图形图像对象

imshow(s); %显示ILPF滤波后的图像 title('ILPF滤波（d0=5)') 结果如图3.3所示

（a）saturn原图 （b）加入椒盐噪声 （c）理想低通滤波后

图3.3理想低通滤波效果

（2） 巴特沃斯低通滤波器

由于Butterworth低通滤波器在抑制噪声的同时，图像边缘模糊程度大大减小，且没有振铃效应。基于以上特点，用Butterworth低通滤波器将低频分量和高频分量分离，低频分量进行均衡后，再将两部分融合，实现图像的增强[17]。

一个阶为n，截断频率为D0的巴特沃斯低通滤波器的转移函数为式（3.2）。

中国矿业大学2012届本科生毕业论文 第25页

阶为1的巴特沃斯低通滤波器剖面示意图见图3.4。由图可见低通巴特沃斯滤波器在高低频率间的过渡比较光滑，所以用巴特沃斯滤波器得到的输出图其振铃效应不明显。 一般情况下，常取是H最大值降到某个百分比的频率为截止频率。在上面式中，当D(u,v)=D0时，H(u,v)=0.5(即降到50%)。另一个常用的截止频率值是使H降到最大值的1/21/2时的频率。

图3.4巴特沃斯低通滤波器剖面示意图

程序：

I1=imread('coins.png');

I2=imnoise(I1,'salt & pepper'); f=double(I2);

g=fft2(f); %傅里叶变换

g=fftshift(g); % 转换数据区域

[N1,N2]=size(g); n=2; d0=50;

n1=fix(N1/2); n2=fix(N2/2); for i=1:N1 for j=1:N2

d=sqrt((i-n1)^2+(j-n2)^2); h=1/(1+0.414*(d/d0)^(2*n)); result(i,j)=h*g(i,j); end

end

result=ifftshift(result); X2=ifft2(result); X3=uint8(real(X2)); subplot(1,3,1),imshow(I1)

中国矿业大学2012届本科生毕业论文 第26页

subplot(1,3,2),imshow(I2) subplot(1,3,3),imshow(X3) 运行结果如图3.5

（a）coins原图 （b）加入椒盐噪声 （c）butterworth低通滤波

图3.5巴特沃斯低通滤波效果

3.4高通滤波

图像中的细节部分与其频率的高频分量相对应，所以高通滤波可以对图像进行锐化处理。高通滤波器与低通滤波器的作用相反，它使高频分量顺利通过，而消弱低频。 图像的边缘、细节主要位于高频部分，而图像的模糊是由于高频成分比较弱产生的。采用高通滤波器可以对图像进行锐化处理，是为了消除模糊，突出边缘[18]。因此采用高通滤波器让高频成分通过，使低频成分削弱，再经逆傅立叶变换得到边缘锐化的图像。常用的高通滤波器有：理想高通滤波器、Butterworth 高通滤波器、指数高通滤波器。 3.4.1 理想高通滤波器

二维理想高通滤波器的传递函数为：

D(u,v)?D0?0(u H , v ) ? ? 1 D ( u ,v ) （3.10） ? D0?

3.4.2 Butterworth高通滤波器

n 阶巴特沃斯高通滤波器的传递函数定义如下：

1 H ( u , v ) ? （3.11） 2n?D0? 1?(2?1)??D(u,v)??

中国矿业大学2012届本科生毕业论文 第27页

3.4.3 指数高通滤波器

指数高通滤波器的传递函数为：

H(u,v)?exp{?0.347[D0/D(u,v)]n} （3.12）

3.4.4 MATLAB算法及其实现

（1）理想高通

理想高通的程序和理想低通相似，只是程序中h的取值正好相反。 I=imread('coins.png'); figure(1),imshow(I); title('原图像'); s=fftshift(fft2(I)); figure(2);

imshow(abs(s),[]);

title('图像傅里叶变换所得频谱'); figure(3);

imshow(log(abs(s)),[]);

title('图像傅里叶变换取对数所得频谱'); [a,b]=size(s); a0=round(a/2); b0=round(b/2); d=10;

p=0.2;q=0.5; for i=1:a

for j=1:b

distance=sqrt((i-a0)^2+(j-b0)^2); if distance<=d h=0; else h=1; end;

s(i,j)=(p+q*h)*s(i,j); end; end;

s=uint8(real(ifft2(ifftshift(s)))); figure(4); imshow(s);

title('高通滤波所得图像'); figure(5); imshow(s+I);

中国矿业大学2012届本科生毕业论文 第28页

title('高通滤波所得高频增强图像');

结果如图3.6所示

（a）coins原图 （b）傅里叶变换所得频谱

（c）傅里叶变换去对数频谱 （d）理想高通滤波后图像

（e） 高频增强后图像

图3.6理想高通滤波器增强效果

（2） Butterworth 高通

根据滤波函数，得到如下算法[19]： I1=imread('eight.tif');

f=double(I1);

中国矿业大学2012届本科生毕业论文 第29页

g=fft2(f);

g=fftshift(g); [N1,N2]=size(g); n=2; d0=50;

n1=fix(N1/2); n2=fix(N2/2); for i=1:N1; for j=1:N2;

d=sqrt((i-n1)^2+(j-n2)^2); if d==0;

h=0; else

h=1/(1+(d0/d)^(2*n)); end

result(i,j)=h*g(i,j); end end

result=ifftshift(result); X2=ifft2(result);

X3=uint8(real(X2));

subplot(1,2,1),imshow(I1);

subplot(1,2,2),imshow(X3);

运行结果如图3.7所示：

（a）eight 原图 （b）Butterworth 高通滤波后

图3.7巴特沃斯高通滤波效果

（3） 小结

由以上两个例子我们看出，高通滤波后图像很昏暗，很多细节都看不清。这是由于图像的大部分能量集中在低频区域，而高通滤波使图中各区域的边界得到较明显的增强的同

中国矿业大学2012届本科生毕业论文 第30页

时滤掉了低频分量，使图中原来比较平滑的区域内部的灰度动态范围被压缩，因而整幅图变得比较昏暗。

中国矿业大学2012届本科生毕业论文 第31页

4 结论与对未来展望

增强图像中的有用信息，它可以是一个失真的过程，其目的是要改善图像的视觉效果，针对给定图像的应用场合，有目的地强调图像的整体或局部特性，将原来不清晰的图像变得清晰或强调某些感兴趣的特征，扩大图像中不同物体特征之间的差别，抑制不感兴趣的特征，使之改善图像质量、丰富信息量，加强图像判读和识别效果，满足某些特殊分析的需要。

图像增强可分成两大类：频率域法和空间域法。前者把图像看成一种二维信号，对其进行基于二维傅里叶变换的信号增强。采用低通滤波（即只让低频信号通过）法，可去掉图中的噪声；采用高通滤波法，则可增强边缘等高频信号，使模糊的图片变得清晰。具有代表性的空间域算法有局部求平均值法和中值滤波（取局部邻域中的中间像素值）法等，它们可用于去除或减弱噪声[20]。

图像增强的方法是通过一定手段对原图像附加一些信息或变换数据，有选择地突出图像中感兴趣的特征或者抑制(掩盖)图像中某些不需要的特征，使图像与视觉响应特性相匹配。在图像增强过程中，不分析图像降质的原因，处理后的图像不一定逼近原始图像。图像增强技术根据增强处理过程所在的空间不同，可分为基于空域的算法和基于频域的算法两大类。基于空域的算法处理时直接对图像灰度级做运算基于频域的算法是在图像的某种变换域内对图像的变换系数值进行某种修正，是一种间接增强的算法。

基于空域的算法分为点运算算法和邻域去噪算法。点运算算法即灰度级校正、灰度变换和直方图修正等，目的或使图像成像均匀，或扩大图像动态范围，扩展对比度。邻域增强算法分为图像平滑和锐化两种。平滑一般用于消除图像噪声，但是也容易引起边缘的模糊。常用算法有均值滤波、中值滤波。锐化的目的在于突出物体的边缘轮廓，便于目标识别。常用算法有梯度法、算子、高通滤波、掩模匹配法、统计差值法等。

图像处理内容涉及光学、微电子学、信息学、统计学、数学、计算机科学等领域，是一门综合性很强的交叉学科，其中任何一门学科的发展都将推动图像处理的进一步发展。近年来随着计算机技术和人工智能、视觉心理研究的迅速发展以及处理器硬件上的不断升级，图像处理向更高、更深层次发展，因此带动图像增强技术的长足进步，其应用的需求也随之越来越广泛。

由于对图像质量的要求越来越高，单一的图像增强算法往往难以满足实际需求，因此 几种算法相结合、取长补短、优势互补是图像增强算法发展必然趋势。更深入的研究可从以下几方面着手：

（1）在进一步提高精度的同时着重解决处理速度问题。如， 在航天遥感、气象云图处理方面，巨大的数据量和处理速度仍然是主要矛盾之一，图像处理的发展将围绕HDTV(高清晰度电视)的研制， 开展实时图像处理的理论及技术研究，向着高速、高分辨率、 立体化、多媒体化、 智能化和标准化方向发展。

（2）加强软件研究、开发新的处理方法，特别要注意移植和借鉴其他学科的技术和研究成果，创造新的处理方法。

（3）加强边缘学科的研究工作，促进图像处理技术的发展。如，人的视觉特性、心理学特性等的研究，如果有所突破，将对图像处理技术的发展起到极大的促进作用。

中国矿业大学2012届本科生毕业论文 第32页

（4）加强理论研究， 逐步形成图像处理科学自身的理论体系，新理论与新算法研究。在图像处理领域，近几年来， 引入了一些新的理论并提出了一些新的算法，如小波分析(Wavelet)、分形几何(Fractal)、 形态学(Morphology)、 遗传算法( Genetic Algorithms)、人工神经网络等(Artificial neural networks)。这些理论及建立在其上的算法， 将会成为今后图像处理理论与技术的研究热点[21]。

（5）图像处理领域的标准化，图像的信息量大、数据量大， 因而图像信息的建库、检索和交流是一个重要的问题。就现有的情况看， 软件、硬件种类繁多， 交流和使用极为不便，成为资源共享的严重障碍。应建立图像信息库， 统一存放格式，建立标准子程序，统一检索方法。

数字图像处理经过初创期、发展期、普及期及广泛应用几个阶段，如今已是各个学科竞相研究并在各个领域广泛应用的一门科学。随着科学技术的进步以及人类需求的不断增长，图像处理科学无论是在理论上还是实践上， 均会取得更大的发展。

中国矿业大学2012届本科生毕业论文 第33页

参考文献

[1] 朱秀昌.2011.数字图像处理教程.北京：清华大学出版社.1~20

[2] J. R. Parker,2010, Algorithms for Image Processing and Computer Vision.出版社: Wiley

65~66．

[3] 王宏民，朱莹莹，付源．空域法进行图像增强处理在MATLAB中的仿真与分析．煤矿

机械，2005(5)：56~58． [4] 胡晓军，徐飞.2011.MATLAB应用数字图像处理.西安：西安电子科技大学出版社78~90 [5] Otsu N A．Threshold selection method from gay-level histogram．IEEE

Trsns.SMC.1998,SMC 8：62—66．

[6] 冈萨雷斯.2004.数字图像处理（MATLAB版）.北京：电子工业出版社76~89

[7] ERAMIAN M, MOULD D. Histogram equalization using neighbor hood metrics[A].

Proceedings of the 2nd Canadian Conference on Computer and Robot Vision[C]. Victoria, BS, Canada, 2005. 397-404.

[8] 张倩.2011. 详解MATLAB图像函数及其应用.北京：电子工业出版社116~128

[9] KIM Y T. Contrast enhancement using brightness preserving bi-histogram

equalization[J]. Consumer Electronics, IEEE Transactions on, 1997, 43(1):1-8 [10] 王桥.2009.数字图像处理.北京：科学出版社，200~203

[11] 阮秋琦，阮宇智.2011. 数字图像处理（第三版）.北京：电子工业出版社340~345 [12] SIM K S, TSO C P, TAN Y Y. Recursive sub-image histogram equalization applied to

gray scale images[J]. Pattern Recognition Letters, 2007, 28(10):1209-1221 [13] 章毓晋. 2006. 图像处理（上册）：图像处理. 北京：清华大学出版社.125~201 [14] 雍杨，王敬儒，张启衡．2005.弱小目标低对比度图象增强算法研究． [15] 章毓晋.2006. 图像工程(上):图像处理.北京：清华大学出版社 400~410

[16] CHEN J, PARIS S, DURAND F. Real-time edge aware image processing with the bilateral

grid[J]. ACM Transactions on Graphics (TOG), 2007, 26(3):103 [17] 刘新春.基于局部直方图相关的造影图像边缘检测方法，中国图像图形学报，2000,5（9）

750~754

[18] Richard O. Duda.2009.李宏东等译.模式分类.北京：机械工业出版社.222-225 [19] 朱秀昌，刘峰，胡栋.2008.数字图像处理与图像通信.北京：科学出版社 [20] 张毓晋.2000.图像理解与计算机视觉.北京：清华大学出版社 [21] 马颂德，张正友.1998.计算机视觉.北京：科学出版社 209~211

中国矿业大学2012届本科生毕业论文 第34页

英文原文

Fast Localization of the Optic Disc Using

Projection of Image Features

Ahmed E. Mahfouz and Ahmed S. Fahmy

Abstract-Optic Disc (OD) localization is an important pre-processing step that signi?cantly simpli?es subsequent segmentation of the OD and other retinal structures. Current OD localization techniques suffer from impractically-high computation times (few minutes per image). In this work, we present a fast technique that requires less than a second to localize the OD. The technique is based upon obtaining two projections of certain image features that encode the x- and y- coordinates of the OD. The resulting 1-D projections are then searched to determine the location of the OD. This avoids searching the 2-D image space and, thus, enhances the speed of the OD localization process. Image features such as retinal vessels orientation and the OD brightness are used in the current method. Four publicly-avail-able databases, including STARE and DRIVE, are used to evaluate the proposed technique. The OD was successfully located in 330 images out of 340 images (97%) with an average computation time of 0.65 s. Index Terms—Image features, localization, optic disc, projection.

I. INTRODUCTION

With the new advances in digital modalities for retinal imaging, there is a progressive need of image processing tools that provide fast and reliable segmentation of retinal anatomical structures. The optic disc(OD) is a major retinal structure that usually appears in retinal images as a circular bright object [1]. It is the region where the optic nerve and the retinal and choroidal vessels emerge into the eye [2].A large number of algorithms have been proposed in literature to segment the OD; this includes the use of Hough Transform [3]–[5], active contour models [6], and gradient vector ?ow (GVF) [7] and [8].Nevertheless, the success and ef?ciency of these algorithms depend mainly upon determining a seed point inside the OD, i.e., localization of the OD [6]–[8]. Although manual localization of the OD is suf?cient, the process can be prohibitively cumbersome when dealing with large number of images. This has stimulated several research groups to develop algorithms for automatic localization of the OD [1], [2] and [9]–[12]. OD localization can also be useful for a number of applications. For example, the OD location can serve as a landmark for localizing and segmenting other anatomical structures such as the fovea.

(where the distance between the OD center and the center of the fovea is roughly constant) [2]. The location can also be used to classify left and right eyes in fovea-centered retinal images [13]. In addition, the detection of OD location is sometimes necessary for computing some important diagnostic indices for hypertensive retinopathy based upon vasculature such as central retinal artery equivalent (CRAE) and central retinal vein equivalent (CRVE) [10]. Also, since the OD can be

easily confounded with large exudates and lesions, the detection of its location is important to remove it from a set of candidate lesions [9].

In normal eyes, automatic localization of the OD is simple because it has well-de?ned features. Nevertheless, developing fast and robust methods for automatic localization of the OD could be very challenging due to the presence of retinal pathologies that alter the appearance of the OD signi?cantly and/or have similar properties to the OD [3]. OD localization methods can be classi?ed into two main categories, appearance-basedmethods and model-basedmethods. Appearance-basedmethods identify the location of the OD as the location of the brightest round object within the retinal image. These methods include techniques such as intensity thresholding [4] and [5], highest average variation [14], matched spatial ?lter [12], and principle component analysis[10]. Although these methods are simple and have high success rates in normal images, they fail to correctly localize the OD in diseased retinal images where the pathologies have

中国矿业大学2012届本科生毕业论文 第35页

similar appearance properties to the OD.

Model-based methods depend mainly upon extracting and analyzing the structure of the retinal vessels and de?ning the location of the OD as the point where all the retinal vessels originate [1], [2] and [9].Techniques such as geometrical models [9], convergence of vasculature [1], and Vessel’s direction matched ?lter [2] have a relatively high success rate in diseased images, but they are computationally very expensive because they require segmentation of the retinal vessels as an initial step of the localization process. For example, the geometrical model-based method proposed by Foracchia et al. [9] achieves a success rate of 97.5%, but it takes an average computation time of 2 min to localize the OD in a given image. The Vessel’s direction matched ?lter described by Youssif et al. [2] achieves an accuracy of 98.8%, but it takes an average computation time of 3.5 min per image to correctly locate the OD. The Hausdorff-based template matching technique proposed by Lalonde et al. [11] has a good computation time (1.6 ?0.3s), but the technique relies on three major assumptions: 1) the imaging protocol is known, 2) the OD represents a bright region, and 3) the OD is a circular object. The ?rst assumption requires that the graders mark each image as left, right, or OD centered. The other two assumptions are usually violated in retinal images with pathologies (e.g., [17, images 1, 4–6, 13, 17, 19–21, 25–27, 34, 36, 46, and 219]) leading to low success rate of the OD localization process, i.e., 71.6%in STARE data-base, as reported by Youssef et al. [2]. In this work, a novel fast technique for OD localization is proposed. The new method can be classi?ed as a model-based method in which the OD is considered the region where the main retinal vessels originate in a vertical direction. The computation time of the localization process is signi?cantly enhanced by reducing the problem from one 2-D localization problem to two 1-D problems that does not require segmentation of the retinal vessels. The remaining sections of this manuscript are organized as follows; Section II-A describes the ―easy-to-compute‖ image features that can be used to decompose the image into two 1-D signals. Section II-B contains the methodologies of determining the horizontal and the vertical locations of the OD from the

Fig. 1. Feature_map_1 showing a horizontally sliding window at two different locations, the sliding direction and the projection direction.

resulting two 1-D signals. Section II-C proposes a geometry-based method that can be used to enhance the robustness of the localization process. Section II-D contains the detailed algorithm that can be used to implement the proposed technique and reproduce the results which are displayed in Section III. Section IV contains a discussion of the results and the conclusion, respectively.

中国矿业大学2012届本科生毕业论文 第36页

II. THEORY AND METHODS

A. Projection of Image Features

Searching for the OD location in a 2-D space (image space) renders any localization algorithm highly expensive in terms of computation time. The main objective of this work is to propose a localization algorithm with signi?cantly enhanced speed by converting the typical 2-D localization problem into two 1-D localization problems. This reduction of dimensionality is achieved by projecting certain features from the retinal image onto two orthogonal axes (horizontal and vertical).The resulting two 1-D signals are then searched to determine the horizontal and vertical coordinates of the OD location. The key factor needed for the success of the dimensionality reduction process is to obtain two meaningful 1-D signals that can be used to determine the coordinates of the OD location. In order to produce such 1-D signals, the set of retinal image features to be projected on either axis should be carefully determined.

Two feature maps are used to produce the two 1-Dprojection signals. The ?rst map feature_map_1 is constructed by calculating the difference between the vertical and horizontal edge maps of the retinal image and dividing the result by the intensity map, see Fig. 1. The second map feature_map_2 is constructed by calculating the summation of the vertical and horizontal edge maps and multiplying the result by the intensity map, see Fig. 2. These maps can be used to accurately localize the OD based upon the simple observation that the central retinal artery and vein emerge from the OD mainly in the vertical direction and then progressively branch into the main horizontal vessels, see Fig. 3(a). This vasculature structure of the retina suggests that a vertical window (with height equal to the image height and a proper width) would always be dominated by vertical edges and dark pixels (vertical vessels) when centered at the OD, see location 1 in Fig. 1. Although the window may contain vertical edges at other locations, i.e., corresponding to small vascular branches or lesions, it will always be populated by strong horizontal edges as well, i.e., the edges of the two main horizontal branches of the retinal vessels, see location 2 in Fig. 1. The vasculature structure suggests also that a horizontal window (with height and width equal to the OD diameter) centered at the determined horizontal location would always enclose a large number of strong edges and a large number of bright pixels when centered at the OD, see Fig. 2. This follows the fact that the possibility of having lesions in the regions above or below the OD is very small, because no retinal vessels are present in these regions. It is worth noting that simple gradient operators (the kernel [1 0 1] and its transpose) are used to produce the vertical and horizontal edge maps of the image.

Fig. 2. Feature_map_2 showing a vertically sliding window centered at the OD, the sliding

中国矿业大学2012届本科生毕业论文 第37页

direction and the projection direction.

Fig. 3. (a) Example of a retinal image. (b) Plot of the 1-D signal resulting from projecting feature_map_2 onto the vertical axis. (c) Plot of the 1-D signal resulting from projecting feature_map_1 onto the horizontal axis.

B. OD Localization

In order to localize the OD, the process is split into two steps. In the ?rst step, the image features are projected onto the horizontal axis to determine the horizontal location of the OD. In the second step, the horizontal location, determined from step 1, is searched for the correct vertical location of the OD. The following two sections show these two steps in detail.

It is worth noting that the areas outside the camera aperture (circular region) are excluded using a binary mask generated by thresholding the red component of the image based upon the method described in [3]. Fig. 3(c) shows the 1-D signal resulting from projecting the features encoded in feature_map_1 on the horizontal axis. The value of the signal at each horizontal location is the summation of the values of feature_map_1 inside the vertical window, shown in Fig. 2, when centered at this horizontal location. Notice that the horizontal location of the optic disc is easily identi?ed as the location of the maximum peak of the resulting 1-D signal.

中国矿业大学2012届本科生毕业论文 第38页

Fig. 4. (a) Vertical localization signals corresponding to Peak 1. (b) Retinal image

showing the two candidate OD locations. (c) Vertical localization signals corresponding to Peak 2. (d) Horizontal localization signal.

It is worth Noting that the vessels thickness and the OD diameter are calculated automatically from the image resolution; assuming that the average OD diameter in adults is 1.5mm and the main vessels thickness is 15% of the OD diameter [15].

Assuming that the horizontal location of the OD is successfully identi?ed, the objective now is to determine the vertical location of the OD. This is done by projecting the features encoded in feature_map_2 on the vertical axis, as shown in Fig. 3(b). The value of the signal at each vertical location is the summation of the values of feature_map_2 inside the horizontal window, shown in Fig. 2, when centered at this vertical location. Notice that the vertical location of the optic disc is easily identi?ed as the location of the maximum peak of the resulting 1-D signal.

C. Improving the Robustness

Consider the horizontal signal of the image shown in Fig. 4(b). It can be shown that the true peak corresponding to the OD horizontal location, peak 2 in Fig. 4(d), is not the maximum peak. This is due to the image artifact that appears as a bright spot to the left of the image. If we follow the algorithm described previously, the estimated OD location will be at a point that, by intuition, cannot belong to an optic disc. On the other hand, if the second peak of the horizontal signal is considered a candidate horizontal location for the OD, the estimated OD location will correspond to the correct location of the OD. This observation can be used to enhance the success rate of the technique. That is, instead of considering the maximum peak of the horizontal signal only, a candidate list of possible horizontal OD locations is used. The candidate list contains the locations of the maximum peaks and the vertical localization step is repeated for each candidate horizontal location. This results in possible (2-D) candidate locations of the OD. In order to determine the ?nal location, a set of image features is used to score each candidate. Then, the ?nal location of the OD is taken as the candidate with the maximum scoring index. In this work, the candidate list contains two locations. The scoring index is calculated as the peak strength of the horizontal signal at the candidate location multiplied by a weighing factor.

中国矿业大学2012届本科生毕业论文 第39页

The weighting factor incorporates some a priori knowledge of the typical geometric and appearance properties of the OD.

To calculate the weighting factor, a square window (with edge equal to twice the OD diameter) is centered at the candidate OD location. Then, 10% of the brightest pixels within this window are segmented. If an object (large cluster of bright pixels) exists at the candidate location, the eccentricity, de?ned as the ratio of the object’s minor axis length to the object’s major axis length [16], of this object is calculated. If no object exists, the eccentricity of the candidate location is set to a very small value (e.g., 0.1). Then, the weighting factor of this location is set equal to the calculated eccentricity.

D. Algorithm

Step 1) Get image features

Construct feature_map_1 and feature_map_2:

—feature_map_1(x,y)??Ev?x,y??EH?x,y??/I?x,y?

—feature_map_2?x,y???Ev?x,y??EH?x,y???I?x,y?

Where, Ev and EH are the absolute vertical & horizontal edge images, respectively; and

I?x,y? is the image intensity.

Step 2) Project feature_map_1 on the horizontal axis

1) Define Whrz as a rectangular window of size (image height, 2×main vessel width) centered at the horizontal location x

2) Slide Whrz over feature_map_1 and calculate:

Hproj?x??sum of feature_map_1 values inside Whrz 3) The horizontal location of the OD, ODH, is the location of the maximum peak ofHproj Step 3) Project feature_map_2 on the vertical axis

1) Define Wver as a rectangular window of size (OD diameter, OD diameter) centered at the horizontal location ODH and the vertical location y 2) Slide Wver over feature_map_2 and calculate:

Vproj?y??sum of feature_map_2 inside Wver 3) The vertical location of the OD, ODV, is the location of the maximum peak of Vproj Step 4) Improving the robustness

1) For each candidate (ODH,ODV), select a square region of interest ( ROI ) with edge size of 2×OD diameter

2) Segment the brightest 10% pixels within each ROI 3) Group neighboring pixels into objects

4) Calculate the eccentricity of the largest object:

eccen (ODH,ODV) = ( minor axis length / major axis length )

5) Calculate the Scoring Index ( SI ) of each candidate: SI (ODH,ODV)=

Hproj?ODH??eccen?ODH,ODV?

6) Select the final OD location as the location with the largest SI

III. RESULTS

Four publicly available databases are used to evaluate the accuracy and the computation time of the proposed technique. The four databases are: 1) STARE database (81 images, 605 X700 pixels) [17], 2) DRIVE database (40 images, 565 X 584 pixels) [18],

3) Standard Diabetic Retinopathy Database ―Calibration Level 0‖ (DIARETDB0) (130 images,

中国矿业大学2012届本科生毕业论文 第40页

1500, X1152 pixels) [19].

4) Standard Diabetic Retinopathy Database ―Calibration Level 1‖ (DIARETDB1) (89 images, 1500 X1152 pixels) [19].

The diseased images in the four databases contain signs of DR, such as hard exudates, soft exudates, hemorrhages, and neo-vascularization (NV). The accuracy and computation time results of evaluating the proposed method using these databases are summarized in Table I. The table includes the results of applying the method without constructing the candidate list (the maximum peak of the horizontal localization signal is selected as the correct location) and also the results with the list containing two candidates.

Fig. 5. Selected OD localization results. (a)–(l) Show successful OD localization examples. (m)–(o) Show examples of failure in OD localization. The white ―X‖ indicates the location of the OD as detected by the proposed method.

TABLE I

RESULTS FOR EVALUATING THE PERFORMANCE OF THE PROPOSED METHODIN TERMS OF SUCCESS RATE AND COMPUTATION TIME USING FOUR DATABASES Database STARE DRIVE DIARETDB0 DIARETDB1 Total Normal Images Diseased Images Number of Images Success Rate ( % ) Success Rate ( % ) ( with improvement) Computation Time ( seconds ) 31 50 81 89.0 92.6 0.46 33 7 40 100 100 0.32 20 110 130 94.6 98.5 0.98 5 84 89 96.6 97.8 0.98 89 251 340 94.4 97.0 The proposed method was implemented using Matlab (The MathWorks, Inc.) and the results shown in Table I are acquired by running the developed code on a PC (2.66 Intel Core 2 Due and 4 GB RAM). The detected location of the OD is considered correct if it falls within 60 pixels of a

中国矿业大学2012届本科生毕业论文 第41页

manually identi?ed OD center, as proposed by Hoover et al. in [1], Foracchia et al. in [9], and Youssif et al. in [2]. The center of the OD is manually identi?ed as the point from which the main retinal vessels emerge.

The proposed method achieved a total success rate of 97%, i.e., the OD was correctly located in 330 images out of the 340 images tested. The new method achieved a success rate of 92.6% in STARE with an average computation time of 0.46 s per image and an error of 14 ? 15 pixels (mean ? standard deviation). In addition, the method achieved 100% success rate in DRIVE with an average computation time of 0.32s per image and an error of 11 ? 11 pixels. Fig. 5(a)–(m) show examples of retinal images with the OD correctly localized. Fig. 5(m)–(o) show examples of retinal images in which the proposed method failed to locate the OD.

Table II compares the performance of the proposed method to other OD localization methods, on STARE database, in terms of accuracy and computation time.

Table III illustrates the parameters’ values used in the proposed algorithm to produce the results shown in Table I. It is worth noting that these parameters are calculated automatically from the image resolution and small changes in the windows’ sizes will not alter the results signi?cantly.

IV. DISCUSSION &CONCLUSION

A new method for OD localization is presented. The method achieves fast localization by reducing the dimensionality of the search space from a 2-D space (image space) to two 1-D spaces (two 1-D signals). The process of dimensionality reduction is achieved through the projection of certain image features onto two orthogonal axes (horizontal and vertical). In this work, two features are selected for projection. The ?rst one is the directionality of the retinal vessels (represented by the horizontal and vertical edge maps of the image). The second feature is the intensity pro?le of the OD (a bright circular region with a thin dark slab in the middle). The robustness of the novel technique is guaranteed by evaluating the method using four databases (340 images); most commonly used to evaluate OD localization techniques [1]–[3] and [9]. The method achieved a relatively high success of 97% with all the parameters of algorithm maintained at constant values (except for linear scaling of the sizes of the used projection windows according to the image resolution). That is, there is no need to tailor the algorithm parameters for different databases. TABLE II

COMPARING THE PERFORMANCE OF DIFFERENT OD LOCALIZATIONTECHNIQUES ON STARE DATABASE

Technique Sinthanayothin et al. [14] - Highest Average Variation Lalonde et al. [11] – Template Matching Hoover et al. [1] – Fuzzy Convergence Foracchia et al. [9] – Geometrical Model Youssif et al. [2] – Vessel’s Direction Matched Filter Proposed Method – Projection of Image Features Success Rate* 42.0% 71.6% 89.0% 97.5% 98.8% 92.6% Computation Time** Not reported 1.6 seconds 15 seconds 2.0 minutes 3.5 minutes 0.46 seconds

*The success rates shown are all taken from Youssef et al. [2]. **The

computation times shown are taken from the corresponding references. Note that the computation time reported by Lalonde et al. [11] is for their local database (640 X480 pixels), while STARE is (605 X700 pixels).

Reducing the search space dimensionality resulted in a signi?cantly reduced computation time, compared to other techniques, see Table II. The key factor that greatly enhanced the speed

中国矿业大学2012届本科生毕业论文 第42页

TABLE III

ALGORITHM PARAMETERS

Database OD Diameter Main vessel thickness STARE 100 15 DRIVE 80 12 DIARETDB0,1 215 32 of the algorithm is that the directionality of the retinal vessels, weather vertical or horizontal, is described by the directionality of their corresponding edges in the vertical and horizontal edge maps of the retinal image. The process of convoluting the image with a 3 X 1 gradient mask to produce the edge maps is the most computationally demanding operation in the algorithm. This operation is negligible if compared to the initial step of retinal vessels extraction that is required in all model-based techniques. The latter is usually achieved by applying a 2-D matched ?lter (typically a 10 X 15 mask for the resolution of STARE images) with several orientations (typically at 12 different angles) [20]. The accuracy of the proposed technique is highly dependent upon the accuracy of the horizontal localization process. Hence, accuracy of the proposed technique could be increased if two candidate locations of the OD are estimated. Scoring of these candidates is done by incorporating their appearance and geometric properties. The process of investigating two candidates increases the total success rate of the technique from 94.4% to 97% with no signi?cant computation overhead, given that the process of investigating two candidate locations was applied to all the images.

As shown in Fig. 5, even in the presence of retinal pathologies and/or image artifacts, the selected features were unique to the OD and, thus, allowed proper localization with relatively high accuracy.

ACKNOWLEDGMENT

The authors would like to thank A. Ahmed for the implementing and for providing the computation time of the technique in [1].

REFERENCES

[1] G. A. Hoover and M. Goldbaum, ―Locating the optic nerve in a retinal image using the fuzzy convergence of the blood vessels,‖ IEEE Trans. Med. Imag., vol. 22, no. 8, pp. 951–958, Aug. 2003.

[2] A. Youssif, A. Ghalwash, and A. Ghoneim, ―Optic disc detection from normalized digital fundus images by means of a vessels’ direction matched ?lter,‖ IEEE Trans. Med. Imag., vol. 27, no. 1, pp. 11–18,Jan. 2008.

[3] F. ter Haar, ―Automatic localization of the optic disc in digital color images of the human retina,‖M.S. thesis, Computer Science Dept.,Utrecht Univ., Utrecht, The Netherlands, 2005. [4] S. Tamura, Y. Okamoto, and K. Yanashima, ―Zero crossing interval correction in tracking eye-fundus blood vessels,‖ in Proc. Int. Conf. Pattern Recognit., pp. 227–233.

[5] Z. Liu, O. Chutatape, and S. M. Krishnan, ―Automatic image analysis of fundus photograph,‖

in Proc. Conf. IEEE Eng. Med. Biol. Soc., pp. 524–525.

[6] A. Osareh, M. Mirmehdi, B. Thomas, and R. Markham, ―Colour morphology and snakes for optic disc localisation,‖ in Proc. 6thMed. Image Understand. Anal. Conf., pp. 21–24.

[7] V. Thongnuch and B. Uyyanonvara, ―Automatic optic disk detection from low contrast retinal images of ROP infant using GVF snake,‖ Suranaree J. Sci. Technol., vol. 14, pp. 223–226, Jan./Dec. 2007.

[8] A. Osareh, M. Mitmehdi, and R. Markham, ―Comparison of colour spaces for optic disc

本文来源：https://www.bwwdw.com/article/v8hr.html