精华内容
下载资源
问答
  • 图像处理设计主要有以下几种处理:图像增强(灰度变换、直方图增强、空间域滤波、频率域滤波)、图像分割、图像配准等等。 图像增强图像增强作为基本的图像处理技术,目的在于通过对图像进行加工使其比原始图像...

    图像处理设计主要有以下几种处理:图像增强(灰度变换、直方图增强、空间域滤波、频率域滤波)、图像分割、图像配准等等。

    1. 图像增强:

    图像增强作为基本的图像处理技术,目的在于通过对图像进行加工使其比原始图像更适合于特定应用,即图像灰度增强是根据特定需要有目的进行。医学图像由于成像设备和获取条件等影响,可能会出现图像质量的退化:另外影像医生希望获得对比度高、细节丰富、可读性好的图像以降低阅片强度便于诊断。通过图像增强改善图像的视觉质量,让观察者能够看到更加直接、清晰、适于分析的信息。传统的图像灰度增强方法可分为空域法和频域法两大类。空域法图像灰度增强直接对图像中像素灰度值进行运算处理,如线性灰度变换、非线性灰度变换、直方图均衡化处理等。频域法图像灰度增强首先对图像进行频域变换,然后对各频谱成分进行相应操作,最后经过频域逆变换获得所需结果。任何一种图像灰度增强算法都只是在特定的场合下才可以达到较为满意的增强效果。

      • 灰度变换:
    1. 线性灰度变换:

    是将原图像的灰度动态范围按线性关系扩展到指定范围或整个动态范围。

    a、b为原图像所占用灰度级别的最小值和最大值,c、d为增强后的图像所占用灰度级别的最小值和最大值。

     

    该方法实现简单,算法复杂度低,对于曝光不足或过度、成像设备的非线性、记录设备的动态范围太突等原因引起的图像对比度过低的情况,线性灰度增强将取得较好的灰度增强效果。

     

     

    1. 对数灰度变换:

    对数变换的一般表达式为

     

    其中,c为一个常数,r小于或等于0,。此种变换使一窄带低灰度输入图像值映射为一个宽带输出值。相对应的是输入灰度的高调整值。可以利用这种变换来扩展被压缩的高值图像中的暗像素。相对应的是反对数变换的调整值。一般对数函数的所有曲线都能完成图像灰度的扩散/压缩。事实上,对数函数有其重要特征,即在很大程度上压缩了图像像素值的动态范围,其应用的一个典型例子是傅里叶频谱,它的像素值具有很大的动态范围。

      • 直方图增强:

    在图像处理中,一种最简单且实用的工具是图像的灰度直方图。通过图像灰度直方图的分布情况,可以大致判断一幅图像的质量。如果一幅图像的灰度直方图挤压在一个较小的灰度范围内,图像的灰度动态范围就小,图像的对比度就差,图像的质量也就不好。反之,图像的灰度动态范围大,图像的对比度就好。要改善图像的灰度动态范围小的问题,一个直观的想法就是修改图像的直方图。常用的修改直方图的方法主要有:灰度变换和直方图增强。灰度变换又称为对比度扩展与调整,它是一种逐像素点对图像进行变换的增强方法,一般是通过线性或非线性函数对图像的灰度进行逐点修改来实现图像增强的。直方图增强技术是一种通过改变图像的全部或局部对比度进行图像增强的技术,其典型处理方法之一就是直方图均衡化。直方图均衡化的思想是把原始图像中的像素灰度做某种映射变换,使变换后图像灰度的概率密度为均匀分布,即变换后的图像灰度级均匀,这意味着图像灰度的动态范围的增加,提高图像的对比度。

    直方图均衡化方法有以下两个特点:

    (a) 根据各灰度级别出现频率的大小,对各个灰度级别进行相应程度的增强,即各个级别之间的间距相应增大。因此,这种方法对于对比度较弱的图像进行处理是很有效的。

    (b) 因为直方图是近似的概率密度函数,所以用离散灰度级作变换一般得不到完全平坦的结果。另外,变换后的灰度级可能会減少,这种现象叫做“简并”现象。由于简并现象的存在,处理后的灰度级总是要减少的,这是像素灰度有限的必然结果。由于上述原因,数字图像的直方图均衡只是近似的。

      • 空间域滤波:

    空域滤波是在图像空间中借助模板进行邻域操作完成的,各种空域滤波器根据功能主要分成平滑滤波器和锐化滤波器两类。图像平滑的目的主要是消除图像中的噪声,而图像的锐化则是为了增强被模糊的细节,如图像的边缘等。因此,在图像增强中主要用到两类空间滤波器。

        1. 平滑滤波器:

    平滑滤波器主要用来減弱或消除图像中的噪声成分,从而提高图像的信噪比。平滑滤波器类似于后面讲到的频域中的低通滤波器。因为高频分量对应图像中的区域边缘与噪声等灰度值具有较大、较快变化的部分,滤波器将噪声减弱或消除的同时,也会减弱图像的边缘信息。

    1. 均值滤波法:

    这是一种在空间域对图像进行简单平滑处理,其原理就是用某像素邻域内的各点灰度值的平均值代替该像素的原值。这种处理可以减小图像灰度的“尖锐”变化。由于图像噪声一般为“尖锐”变换的白噪声,所以均值滤波一般用来处理图像中的噪声。一个典型的均值滤波器是它的各元素值相等,且各元素的和为 1。 另外一种均值滤波器,它采取加权平均的方式,即不同的掩膜元素具有不同的权值,从而突出了一些像素的重要性,处于掩膜中心位置的像素较距离掩膜中心较远的其他像素更重要。这样做可以在降低蝶声的同时,减轻平滑处理所带来的边缘信息模糊效应。均值滤波法的优点是容易实现对噪声的抑制,缺点是容易使目标轮廓变得模糊,而且会减弱有用的细节信息。

    1. 中值滤波法:

    中值滤波法是一种非线性滤波。这种方法运算简单,对孤立噪声的平滑效果比均值滤波方法好,而且它能较好地保护图像边界,但是这种算法会使图像失掉细线和小块的目标区域。中值滤波实际上就是用一个含有奇数个像素的滑动窗口,将窗口的正中点的灰度值用窗口内各点的中值代替。例如,若窗口长度为 5,窗口中像素的灰度值分别为 80、90、 200、110、120。按从小到大排列,结果为 80、90、110、120、200,其中间位置上的值 为 110。于是原来窗口正中的灰度值 200 就由 110 代替。如果灰度值 200 是一个噪声的尖峰,则将被滤波;如果它是一个有用的信号,那么此方法处理的结果将会造成信号的损失。中间值取法如下:当邻域内的像素数为奇数时,取排序后的中间像素的灰度值,当邻域内的像素数为偶数时,取排序后的中间两像素的灰度值的平均值。当模板中心像素点处在边缘像素时,模板会超出图像范围,在这种情况下,图像边缘上的像素点保留原值即可。中值滤波的突出优点是在消除噪声 的同时,还能防止边缘模糊。对图像进行中值滤波时,通常选择的滤波窗口是方形的(具有奇数行和列)。在实际使用窗口时,窗口的尺寸一般先取3,再取5,依次增大,直到滤波效果满意为止:对于有缓变的较长轮廓线目标的图像,采用方形或圆形窗口较合适:对于包含尖顶角目标的图像,采用十字形窗口较合适。使用二维中值滤波方法最需要注意的是保持图像中有效的细线状结果。

        1. 锐化滤波器:

    锐化滤波器主要用来通过增强图像的边缘信息,突显图像中感兴趣区域的轮廓。平滑滤波器类似于后面讲到的频域中高通滤波器。由于图像中的边缘信息与噪声都处在高频成分,因此,锐化滤波器将图像边缘锐化的同时,也会降低图像的信噪比。尽管这些滤波器的功能不同,但在空域实现这些功能的方式都是相似的,即都是利用模板与图像做卷积运算。在图像处理中,把消减图像模糊、突出目标边界与图像细节的增强方法称为图像锐化,因此锐化技术常用于加强图像中的目标边界和图像细节。图像模糊是常见的图像降质问题。在图像获取、传输及处理过程中有许多因素会使图像变模糊。造成图像模糊的原因很多,例如,在重建过程中由卷积反投影法、迭代算法等引起的模糊,CT系统扫描过程中投影数据不完全也会造成图像模糊,X线的散射、 环境噪声、部分容积效应等都会使图像变模糊。由于图像模糊,使得图像中的一些重要信息如图像的边缘细节等难以被识别,这样会给图像判读造成很大困难。大量研究表明,各种图像变模糊的物理过程的数字模型含有求和、平均或积分运算。在某些场合应用中,可以不必深究图像变模糊具体的物理过程及其数学模型,而根据各种图像变模糊过程都有相加或积分运算这一共同点,在空间域中运用微分运算突出变化以增强图像,基本的方法是在原图像加上一个其微商后的图像,以抵消积分运算带来的图像模糊。由傅里叶变换的徽分性质可知,通过在频域中用提升信号高频分量的方法增强图像,也可以处理图像的模糊问题。另外,即使图像可能没有什么模糊失真,在人或机器分析图像时也常常需要突出目标边界和灰度细节。因此,图像锐化的目的是加强图像轮廓,使图像看起来比较清时。图像轮廓是灰度陡然变化的部分,包含丰富的空间高频成分,高频分量相对突出,则轮廓清晰。

      • 频率域滤波:

    对图像做傅里叶变换就可以得到它的频谱。在频率域中,零频率分量相当于图像的平均灰度,低频率分量对应于平滑的图像信号,较高频率的分量对应于图像中的细节和边界。通常,我们认为噪声的频率也处于高频分量中,因此,可以通过处理图像的高频部分来平滑图像。反之,去掉低频部分,可实现对图像的锐化和边缘提取。在频域中进行增强的主要步骤有:计算需要增强图像的傅里叶变换;将其与一个根据需要设计的转移函数相乘;再将结果进行傅里叶反变换,得到增强图像。常用的频域增强方法有低通滤波和高通滤波。

        1. 低通滤波器:

    图像中的边缘和噪声都对应着图像傅里叶变换中的高频部分,所以,在频率域中,通过滤波器转移函数衰减图像的高频信息,而使低频信息畅通无阻地保留下来的过程称为低通滤波器。低通滤波器抑制了反应图像边界特征的高频信息以及包括在高频中的孤立点噪声,起到平滑图像去噪声的增强作用。理想低通滤波器的截止频率D0决定于通过滤波器的能量比,D0越小,通过的能量越小。一般来说,D0应使得图像中我们感兴趣的细节能量大部分通过,而截断大部分不感兴趣的部分。然而,这本身就存在一定的矛盾,截断部分不仅包括不感兴趣的部分,还包含感兴趣的部分。一个有效的办法就是通过控制通过的能量比来控制图像的 滤波程度。理想低通滤波器只能在数学上清楚地定义出来,而在实际的电子器件中是不能实现的。理想的低通滤波器的转移函数为:

     

    其中,D(u,v) 是频率平面点(u,v)到频率平面原点(0,0)的距离;D0是一个规定的非负整数,称为截止频率。传递函数表明,以D0为半径的圆域内所有频率分量无损地通过,圆域外的所有频率分量被完全滤除。

     

     

        1. 高通滤波器:

    高通滤波器是为了衰减或抑制低频分量,而保留高频分量的滤波形式。因为边缘及灰度急剧变化部分与高频分量相关联,在频率域中进行高通滤波将使图像得到锐化处理。与理想低通滤波器一样,理想高通滤波器也无法用实际的电子器件硬件来实现。理想的高通滤波器(IHPF)的传递函数为:

     

    其中,D(u,v) 是频率平面点(u,v)到频率平面原点(0,0)的距离;D0为频率平面上从原点算起的截至距离,称为截止频率。

    1. 图像分割:

    医学图像分割是提取影像中特殊组织从而进行定量和定性分析及可视化必不可少的环节。在医学图像处理中图像分割有着广泛的应用,如组织结构分析、运动分析、三维可视化、图像引导手术、肿瘤放射治疗、治疗评估等研究中都是假设己对图像进行了准确分割,或以图像分割为基础的。目前国内外有关医学图像分割算法的研究成果非常丰富,但鉴于医学图像本身的复杂性和多样性以及临床应用对医学图像分割的准确度和处理速度要求较高,目前的分割算法还远末达到要求,所以医学图像分割算法的研究仍是医学图像处理与分析的热点。图像分割是根据图像的基些特征或特征集合的相似性对图像像素进行聚类,把图像平面划分成若干个具有某些一致性的不重登区域,使得同一区域中的像素特征具有一致性,而不同区域间像素的特征具有非一致性。这里的特性可以是像素的灰度、颜色、纹理等。

        1. 边缘检测:

    边缘检测技术是最早研究的图像分割技术之一,其依据是区域边缘上的像素灰度值变化剧烈,通过检测不同均匀区域之间的边缘来解决图像分割问题。边缘属于局部概念,它是一组位于两个区域边界上相连像素的集合。实际上,由于图像采集技术不完善使获得的图像边缘存在模糊,所以用“类斜面”的数字模型模拟边缘更为准确。因此我们可以定义“边缘宽度”的概念,即从初始灰度级跃变到最终灰度级的斜坡长度,此长度取决于斜坡的斜度并受边缘模糊程度的影响。分别利用一阶导数和二阶导数处理上述斜坡数字边缘模型。一阶导数可用来检测图像中某一像素是否是边缘像素;二阶导数可以用来判断某个边缘像素是在边缘暗的一边还是在边缘亮的一边。对图像的每条边缘,二阶导数生成一个正值和一个负值,连接这两个值的虚构直线将在边缘中点附近穿过零点。注意,二阶导数的这个过零点性质对于确定粗边缘的中心非常有用。因此,边缘是图像灰度值不连续的结果,这种不连续性可以利用导数检测到。在数字图像中利用差分来近似微分进行求解,借助空域微分算子卷积实现边缘检测。

    1. 梯度算子:

    梯度对应一阶导数,即梯度算子是一阶导数算子。对一个连续函数f(x,y),它在位置(x,y)的梯度可表示为一个矢量:

     

    从向量分析中可知,梯度向量指向坐标为(x,y)的f最大变化率方向。这个矢量的幅度(常简称为梯度)和方向角分别为:

     

    其中偏导数需对每个像素计算,在实际中常用小区域模板卷积来近似计算,用绝对值对梯度进行近似。对Gx和Gy各用一个模板,也就是将两个模板组合起来构成一个梯度算子。根据模板的大小和其中元素值的不同,有多种不同的算子。

    Roberts算子定义为:

     

     

    Prewitt算子定义为:

     

    Sobel算子定义为:

    算子运行时采取类似卷积的方式,将模板在图像上移动并在每个位置上计算对应中心像素的梯度值,即对一幅灰度图像求梯度所得的结果是一幅梯度图Robert算子基于可通过任意一对互相垂直方向上的差分来计算梯度的原理,采用两对角线方向相邻像素之差近似梯度幅值来检测边缘。它检测斜向边缘的效果好于水平和垂直边缘,具有计算简单、定位精度高、对蛛声敏感等诸多特点。Sobel算子、Prewitt算子比Robert算子的抗虽声能力要强一些,这是因为 Sobel算子和Prewitt算子是 8 邻域算子,而Robert算子是一个 4 邻域算子,并且 Sobel算子对噪声具有一定的平滑作用,能提供较为精确的边缘方向信息,但它同时也会检测出许多的伪边缘,边缘定位精度不够高。

     

     

    1. 拉普拉斯算子:

    拉普拉斯(Laplacian)算子是一种二阶微分算子,对一个连续函数f(x, y),它在位置(x,y)的拉普拉斯值定义如下:

     

    函数的拉普拉斯值可借助模板求得。对模板的基本要求是,对应中心像素的系数为正,中心像素的邻近像素系数为负,且它们的和为零。实际中,常根据二阶微分算子过零点的性质来确定其边缘位置。拉普拉斯算子很少直接用于边隐检测,主要是因为Laplacian算子对图像噪声非常敏感,容易产生双边缘并且 Laplacian 算子不能检测边缘的方向。拉普拉斯算子易受噪点影响,为减少噪点敏感度,可使用高斯型拉普拉斯(Laplacian of Gaussian,LOG)算子。

     LOG 算子首先对图像进行高斯平滑,然后进行拉普拉斯运算。因为高斯函数可减少噪点,因此该算子对喝点不太敏感,且拉普拉斯模板可使检测到假边缘的概率減到最小。用于卷积的 LOG 函数可定义为

     

    1. Canny算子:

    Canny算子是在综合考虑了前述的传统边缘检测算子后提出的,因为这些边缘检测算子都有一些共同要求:好的检测效果,也就是说对边缘的错误检测要尽可能的低,在图像的边缘出现的地方检测结果中必须有,而且不能出现虚假边缘,在所有使用边缘检测做更深入的研究工作的系统中,它的性能都依赖于边缘检测的误差;对边缘的定位要准确,标记出来的边缘位置要和图像的真正边缘的中心位置充分的接近;对同一个边缘要有低的响应次数,即单个边缘产生多个响应的概率要尽可能的低,而且对虚假边缘的响应要尽可能的抑制。在明确了这三个要求后,Canny算子应运而生,其实现边缘检测的步骤如下:用高斯滤波器平滑图像;计算平滑后的图像的梯度幅值和方向;对梯度幅值采用非极大值抑制,其过程为找出图像梯度中的局部极大值点,把其他非极大值点置零而得到细化的边缘;用双阈值算法检测和连接边缘。

        1. 阈值分割:

    阈值分割是一种简单有效的图像分割技术。所谓的阈值法就是选用一个或几个阈值将图像的灰度级分为几个部分,认为属于同一个部分的像素属于同一目标。假设一幅图像由亮目标和暗背景两部分组成,选取阈值可将目标和背景分开,将灰度值大于阈值的像素点归为目标,其余的像素点归为背景。在阈值分割中确定阈值很关键,合适的阈值可以方便地将图像分割开来。

    1. 全局阈值法(双峰法、最大类间方差法):

    双峰法:假设由目标和背景组成的图像其灰度各具有单峰分布特性,在目标或背景内的相邻像素的灰度高度相关,但目标和背景交界处的像素灰度有很大的差别。则图像灰度直方图就可以看作是由分别对应目标和背景的两个单峰直方图混合而成。此时,若这两个分布数量接近,均值差相距足瞬远且均方差也足的小,则这个直方图就是双峰的。对于这类的图像最简单的方法就是利用双峰法来进行分割。双峰法的实现过程就是通过对图像进行逐个像素扫描,并用数组记录每个灰度级像素点的个数,求出两个单峰的峰顶, 接着在两峰之间找出可作为谷底的像素值作为阈值。找到阈值以后,把大于阈值的像素分成一类,把小于阈值的像素分成一类。

    最大类间方差(Otsu)法:最大类间方差法是一种自适应的阈值确定方法,简称Otsu。它根据灰度特性将图像分成目标和背景两部分。目标和背景之间的类间方差越大,说明构成图像的两部分差别越大,部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小。因此,使类间方差最大的分割意味着错分概率最小。对于图像I(x,y),目标和背景的分割阈值记作T,属于目标的像素点数占整幅图像的比例记为,其平均灰度为,背景像素点数占整幅图像的比 例为,其平均灰度为。图像的总平均灰度记为,类间方差记为g假设图像大小为,背景较暗,像素灰度值小于阈值T的像素个数记作,像素灰度大于阈值T的像素个数记作,则有:

    采用遍历的方法得到使类间方差最大的阈值,即为所求的最佳阈值。Otsu 法属于一种单阈值的分割方法,当图像中目标相比背景而言所占比例很小时,该方法分割结果可能不好。

     

     

    1. 局部阈值法(直方图变换法):

    如果图像的直方图波峰很尖、很空且具有对称性,而波谷很深,就容易找到一个很好的阈值。但是在实际应用中,图像常受到噪声等因素的影响使原本分离的峰间波谷被填充,要检测两峰间的波谷就很困难。为解决这类问题可采用将像素自身性质和一些像素间邻域的局部性质相结合的方法,常用的有直方图变换法。直方图变换的基本思想是利用一些像素邻域的局部性质变换原来的直方图以得到一个新的直方图。这个新直方图与原直方图相比,或者峰间的谷更深,或者谷转变成峰从而更易检测。常用的像素邻域局部性质是像素的梯度值。变换的直方图一般可分为两类:利用具有低梯度值像素的直方图;利用具有高梯度值像素的直方图。已知目标和背景内部的像素具有较低的梯度值,而它们边界上的像素具有较高的梯度值。如果作出仅具有低梯度值像素的直方图,那么这个新直方图中对应内部点的峰基本不变,但谷应比原直方图要深。更一般地,可计算一个加权的直方图,其中吐给具有低灰度值的像素权重大一些。在这样的直方图中,边界点贡献小而内部点贡献大,峰基本不变而谷变深,所以峰谷差距加大。利用具有高梯度值像素的直方图,可作出仅具有高梯度值像素的直方图。直方图在对应目标和背景的边界像素灰度级处有一个峰。这个峰主要由边界像素构成,对应这个 峰的灰度值就可先作分割用的阁值。更一般地,可以计算一个加权的直方图,赋给高梯度值的像素权重大一些。这样在统计直方图时梯度值为爱的像素就不必考虑,而具有大梯度值的像素将得到较大的权重。

    1. 动态阈值法:

    用与坐标相关的一系列阈值来对图像进行分割的方法叫做动态阈值法。它的基本思想:首先,将图像分解成一系列子图像,这些子图像可以互相重复也可以相邻。如果分解成的子图像比较小,则由阴影或对比度空间变化等问题对分割结果造成的影响就会比较小; 然后,对每个子图计算一个阈值,此时的阈值可根据图像具体情况采用上述的阈值选取法。通过对这些子图获得的阈值进行插值就可得到对图像中每个像素进行分割所需的阈值:最后,将图像中每个像素与对应的阈值相比较即可实现分割。总地来讲,阅值分割法计算简单,分割速度快,但它忽略了图像的空间信息,这导致了闻值分割法对噪声的灰度不均匀性很敏感。所以对医学图像进行具体分割时,常结合图像的其他信息,比如,图像的梯度、纹理等局部统计信息。

    1. 图像配准:

    图像配准是指对一幅图像进行一定的几何变换映射到另一幅图像中,使得两幅图像中的相关点达到空间上的一致。医学图像配准是医学图像处理过程中的一项基本任务,临床上通常需要对同一个病人进行多种模式或同一种模式的多次成像,即同时从几幅图像获得信息,进行综合分析。利用图像配准技术,继而将上述图像融合起来,在同一幅。图像上同时表达人体的多方面信息,从医学影像上反应人体的内部结构和功能状态,更加直接提供人体解剖和生理病理信息。图像配准可以定义为两幅图像在空间和灰度上的一一映射的过程,也就是说,将两幅图像中对应于空间同一位置的点联系起来,这里的映射一般称为变换。数学语言描述为:若用给定尺寸的二维矩阵表示两幅图像,分别表示相应位置上的灰度值,则图像间的配准关系可用下面的公式来表示:     

     

    其中,g则表示一个一维亮度变换;f代表一个二维空间坐标变换。

    图像配准问题的任务包括找到最优的空间和亮度变换关系,使得图像相对于失配源得到匹配。通常情况下亮度变换g在图像预处理阶段就得到了纠正,所以寻找两幅图像之间的空间变换或几何变换是配准的关键。

    医学图像配准技术是特征空间、搜索空间、搜索算法和相似性测度四个不同方法的组合:特征空间是指从参考图像和浮动图像中提取的可用于图像配准的特征;搜索空间进行变换的方式以及变换范围;搜索算法决定下一步的具体方向以得到最优的变换参数;相似性测度是用来度量图像间相似性的一种标准。

    一般配准的基本步骤如下:首先根据待配准数据的特性确定配准模型, 包括选择适合的特征空间和变换搜索空间,并根据特征空间的具体形式定义图像之间的相似性测度函数;每次对浮动图像施加变换,计算变换后的浮动图像和参考图像所能达到的相似度。不断改变变换参数,使得相似性测度函数达到最优(这个过程需要选择有效的搜索算法实现,即配准实际上可以转化为多参数的最优化问题);由配准模型求解出配准变换参数后,将其作用于浮动图像。变换后的浮动图像被 认为和参考图像达到了空间上的匹配,即两幅图像对应的各点位置已经一一配准。

    医学图像配准方法可以分为两大类: 一类是直接利用图像本身信息进行配准处理,包括基于特征点的配准方法、基于表面的配准方法、基于像素的配准方法;另一类是变换图像,根据变换后的信息对图像进行配准,包括基于小波变换的配准、基于傅里叶变换的配准等。

    1. 基于特征点的图像配准方法:特征点有外部特征点和内部特征点之分。

    外部特征点是在受试者颅骨中嵌入的螺针,玻璃珠、铬合金珠、明胶球等。内部特征点是指解剖结构上容易定位的,或者是几何上的极值点,如耳蜗尖端拐点处、两个线形结构的交点、血管的分叉或相交处、某一表面上的特定拓扑属性、一个沟回的可识别部分等。基于特征的图像配准算法的步骤可分为: 特征提取、特征匹配、模型参数估计、图像变换和灰度插值。

    1. 基于表面的配准方法:首先提取两幅图像中对应的曲线或曲面,然后根

    据这些对应的曲线或曲面决定几何变换。其中最典型的就是头帽算法。具体实现方法是从一幅图像轮廓提取点集,该点集称作帽子(“hat”),从 另一幅图像轮廓提取表面模型,该模型叫做头(“head”)。一般用体积较大的患者图像,或在图像体积大小差不多时用分辨率较高的图像来产生头表面模型。利用头帽法不仅可实现头颅等三维刚形体图像的配准,而且可用于三维弹性图像的配准。但是,这类方法的最大缺陷是配准精度受限于分割精度,配准前要求勾画相互对应的表面轮预,而对于边界模糊的功能成像图,如SPECT,它们的表面轮廓不易提取出来,不易使用此法。

    1. 基于像素的图像配准方法:直接利用图像中的灰度信息。由于这类方法

    不需要提取图像的解剖特征,不需要对图像进行分割或数据缩减,而且极大地利用了图像信息,近些年成为人们最感兴趣和重视的研究方法。这类配准方法可分为两种:一种是利用图像的统计信息,典型方法是基于矩和主轴法。该方法对数据缺失较敏感,配准结果不太精确,但算法自动、快速易实现,主要被用作预配准,以减少后续精确配准时优化算法的搜索区间和计算时间。另一种是利用图像中的所有灰度信息。

    1. 基于傅里叶变换的方法:这是最为典型的变换域配准方法之一。傅里叶

    变换用于图像配准的基本理论是:图像的平移不影响傅里叶变换的幅值:旋转图像在频域相当于对其傅里叶变换作相同角度的旋转。相位相关法就是一种基于傅里叶功率谱的频域相关技术。该技术利用傅里叶变换的平移特性,解决了仅存在平移关系的图像间的配准问题。

    1. 小波变换:这是一种信号的时间一频率分析方法,具有多分辨率分析的

    特点:在低频部分具有较高的频率分新率和较低的时间分辨率,在高频部分具有较高的时间分辨率和较低的频率分辨率。对两幅图像的伸缩、旋转、平移等配准问题,可以转化为对其作小波分解后两幅图像近似分量的伸缩、旋转、平移配准问题。两幅原图像的平移量若为(2x,2y),则它们分别分解后的两幅图像近似分量的平移量为(x,y),而旋转角度与小波分解以前相等。基于上述原理,将待配准的两幅图像进行 3 层小波分解,先对最低分辨率图像进行配准,然后利用配准结果确定前一层的搜索范围,在确定的范围内进行较高分辨率图像的配准,有效地减少了运算量。由于采用多分辨率小波分解的过程本身也是一个图像平滑的过程,所以在用小波做配准时,为了避免陷入局部极值,常常在最低分辨率的图像中采用优化算法。

    展开全文
  • F=fft2(f,p,q) p,q为填充的尺寸(滤波需要) 傅里叶谱(元素幅度)S=abs(F) imshow(S,[]) 原点移到中心fc=fftshift(F) (频域滤波不需要) 计算反正切(相角)phi=atan2(I,R) 空域滤波器获得频域滤波器 %频域滤波 ...

    傅里叶变换

    F=fft2(f,p,q) p,q为填充的尺寸(滤波需要)

    傅里叶谱(元素幅度)S=abs(F)  imshow(S,[])

    原点移到中心fc=fftshift(F)  (频域滤波不需要)

    计算反正切(相角)phi=atan2(I,R)

    空域滤波器获得频域滤波器

    %频域滤波
    f=imread('d.tif');
    imshow(f,[])
    F=fft2(f);
    S=fftshift(log(1+abs(F)));
    figure,imshow(S,[])
    h=fspecial('sobel');
    freqz2(h) %三维图
    PQ=paddedsize(size(f)); %尺寸
    H=freqz2(h,PQ(1),PQ(2));
    H1=ifftshift(H);
    imshow(abs(H),[])
    figure,imshow(abs(H1),[])
    
    %滤波后的空间域
    gs=imfilter(f,h);%空间滤波
    gf=dftfilt(f,H1)%频率滤波
    figure,imshow(abs(gs)>0.2*abs(max(gs(:))))
    figure,imshow(abs(gf)>0.2*abs(max(gf(:))))
    

     

    频域直接生成滤波器

     

    低通滤波(平滑)

    %高斯低通滤波
    f=imread('e.tif');
    [f,revertclass]=tofloat(f);
    PQ=paddedsize(size(f));
    [u,v]=dftuv(PQ(1),PQ(2));
    D=hypot(u,v); %计算到原点0,0距离 原点为中心用fftshift
    D0=0.05*PQ(2);
    F=fft2(f,PQ(1),PQ(2));%频域傅里叶
    H=exp(-(D.^2)/(2*(D0^2)));%高斯低通滤波器
    g=dftfilt(f,H);%滤波
    g=revertclass(g);
    figure,imshow(fftshift(H))
    figure,imshow(log(1+abs(fftshift(F))),[])
    figure,imshow(g)

     

    线框及表面绘制

     

    3D线框:mesh(H) mesh(H(1:k:end,1:k:end))

    颜色colormap

    网格grid off

    轴axis off

    观察点view(方位角,仰角)

     

    高通滤波(锐化)

     理想高通滤波器IHPF
     巴特沃思高通滤波器BHBF
     高斯高通滤波器(GHPF的结果比BHBF和IHPF的结果更平滑)
     频率域的拉普拉斯算子
     钝化模板、高频提升滤波和高频加强 滤波

     

    %高通滤波
    f=imread('e.tif');
    [f,revertclass]=tofloat(f);
    PQ=paddedsize(size(f));
    [u,v]=dftuv(PQ(1),PQ(2));
    D=hypot(u,v); %计算到原点0,0距离 原点为中心用fftshift
    D0=0.05*PQ(2);
    H=hpfilter('gaussian',PQ(1),PQ(2),D0);
    g=dftfilt(f,H);%滤波
    figure,imshow(g)

     

     

     

    高频强调滤波(高通滤波叠加平均灰度)

    %高通滤波
    f=imread('f.tif');
    figure,imshow(f)
    %[f,revertclass]=tofloat(f);
    PQ=paddedsize(size(f));
    [u,v]=dftuv(PQ(1),PQ(2));
    D=hypot(u,v); %计算到原点0,0距离 原点为中心用fftshift
    D0=0.05*PQ(2);
    HBW=hpfilter('btw',PQ(1),PQ(2),D0,2);
    H=0.5+2*HBW;
    gbw=dftfilt(f,HBW);
    gbw=gscale(gbw);%将灰度值自动缩放到0~255之间 
    figure,subplot(221),imshow(gbw,[])
    title('高通滤波图像'); 
    ghw=histeq(gbw,256);
    subplot(222),imshow(ghw,[])
    title('高通均衡图像'); 
    gbf=dftfilt(f,H);
    gbf=gscale(gbf);
    subplot(223),imshow(gbf,[])
    title('强调后图像'); 
    ghe=histeq(gbf,256);%直方图均衡化  
    subplot(224),imshow(ghe,[])
    title('强调均衡图像'); 

     

     

     

    选择性滤波

    带阻带通

    陷波带阻 陷波带通

     

    展开全文
  • 一)空间域滤波与频率域滤波  1)空间域滤波  空间域滤波是指在图像空间中借助模板对图像领域进行操作,处理图像每一个像素值。主要分为线性滤波和非线性滤波两类,根据功能可分为平滑滤波器和锐化滤波器。...

    一)空间域滤波与频率域滤波

     1)空间域滤波

          空间域滤波是指在图像空间中借助模板对图像领域进行操作,处理图像每一个像素值。主要分为线性滤波和非线性滤波两类,根据功能可分为平滑滤波器和锐化滤波器。平滑可通过低通来实现,平滑的目的有两类,一是模糊,目的是在提取较大的目标前去除太小的细节或将目标内的小尖端连接起来;二是去噪。锐化则可用高通滤波来实现,锐化的目的是为了增强被模糊的细节。

          在matlab中实现空间域滤波,有很多类型,如均值、中值、索贝尔、高斯、拉普拉斯、高斯-拉普拉斯等,但各有差异。下面是用matlab实现的代码:


    [html] view plain copy
    1. %空间域滤波  
    2. clc;close all;  
    3. I=imread('1.tif');  
    4. w1=fspecial('average',[3 3]);  
    5. w2=fspecial('sobel');  
    6. w3=fspecial('gaussian',[3 3],0.5);  
    7. w4=fspecial('laplacian',0.2);  
    8. w5=fspecial('log',[5 5],0.5);  
    9. g1=imfilter(I,w1,'replicate');  
    10. g2=imfilter(I,w2,'replicate');  
    11. g3=imfilter(I,w3,'replicate');  
    12. g4=imfilter(I,w4,'replicate');  
    13. g5=imfilter(I,w5,'replicate');  
    14. g6=medfilt2(I);  
    15. subplot(3,3,1);imshow(I);title('原图');  
    16. subplot(3,3,2);imshow(g1);title('均值滤波');  
    17. subplot(3,3,3);imshow(g2);title('索贝尔滤波');  
    18. subplot(3,3,4);imshow(g3);title('高斯滤波');  
    19. subplot(3,3,5);imshow(g4);title('拉普拉斯滤波');  
    20. subplot(3,3,6);imshow(g5);title('高斯-拉普拉斯滤波');  
    21. subplot(3,3,7);imshow(g6);title('中值滤波');  

    其运行结果如下:


        空间域滤波均采用matlab自带的函数,参数也采用默认的值,可见中值滤波的效果最好。

        其中各滤波器的原理和优缺点如下:

       A)均值滤波:由fspecial函数生成的w1是一个大小为3*3的矩形平均滤波器,再用imfilter这个函数使这个掩模的中心逐个滑过图像的每个像素,输出为模板限定的相应领域像素与滤波器系数乘积结果的累加和。由处理结果可见均值滤波器的效果使每个点的像素都平均到它的领域去了,噪声明显减少了很多,效果较好。

       B)索贝尔滤波:w2是一个大小为3*3的sobel滤波器sv,用来近似计算垂直梯度,在图像的任何一点使用此算子,将会产生对应的梯度矢量或是其法矢量。但是Sobel导数并不是真正的导数,这是因为Sobel算子定义于一个离散空间之上,它真正表示的是多项式拟合,用较大的核的话会在更多像素上进行拟合,会更加正确。而较小的核对噪声会更加敏感,此时用sobel算子近似计算导数的缺点精度比较低,这种不精确性在试图估计图像的方向导数 (使用y/x滤波器响应的反正切得到的图像梯度的方向)。比如对于3*3的Sobel滤波器,梯度角度接近水平或者垂直方向的时候, 这样的不准确性会比较明显。由滤波效果可见到图像的边缘凸显了出来,sobel算子主要用于边缘检测。

       C)高斯滤波:高斯滤波器是平滑线性滤波器的一种,线性滤波器很适合于去除高斯噪声。而非线性滤波则很适合用于去除脉冲噪声,中值滤波就是非线性滤波的一种。高斯滤波就是对整幅图像进行加权平均的过程,每一个像素点的值,都由其本身和邻域内的其他像素值经过加权平均后得到。高斯滤波器是带有权重的平均值,即加权平均,中心的权重比邻近像素的权重更大,这样就可以克服边界效应。高斯滤波如果采用3×3掩模的具体公式如下:

    g(x,y)={f(x-1,y-1)+f(x-1,y+1)+f(x+1,y-1)+f(x+1,y+1)+[f(x-1,y)+f(x,y-1)+f(x+1,y)+f(x,y+1)]*2+f(x,y)*4}/16

    其中,f(x,y)为原图像中(x,y)像素点的灰度值,g(x,y)为经过高斯滤波和的值。由处理效果可看出高斯滤波的减噪能力较好。

     

       D)拉普拉斯滤波:拉普拉斯算子是n维欧式空间的一个二阶微分算子。拉普拉斯算子会突出像素值快速变化的区域,因此常用于边缘检测。由效果可见图像的边界得到了增强。

       E)中值滤波:中值滤波法是一种非线性平滑技术,它将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术,中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替,让周围的像素值接近的真实值,从而消除孤立的噪声点。方法是用某种结构的二维滑动模板,将板内像素按照像素值的大小进行排序,生成单调上升(或下降)的为二维数据序列。二维中值滤波输出为g(x,y)=med{f(x-k,y-l),(k,l∈W)} ,其中,f(x,y),g(x,y)分别为原始图像和处理后图像。W为二维模板,通常为3*3,5*5区域,也可以是不同的的形状,如线状,圆形,十字形,圆环形等。中值滤波对于斑点噪声和椒盐噪声来说尤其有用。保存边缘的特性使它在不希望出现边缘模糊的场合也很有用。由上图效果可见中值滤波的效果最好。

         由以上分析可知,各种滤波器各有优劣,适用情况也不尽相同,线性滤波器很适合于去除高斯噪声,而非线性滤波则很适合用于去除脉冲噪声,如中值滤波很适合去除椒盐噪声。使用起来要视具体实际情况而定。


    2)频率域滤波

      

         开始时看书看了些原理想自己用算法实现发现有困难,就去百度了一下,写了下面的一些代码:

        

    [html] view plain copy
    1. %频率域滤波  
    2. clc;close all;  
    3. f=imread('1.tif');  
    4. f=im2double(f);  
    5.   
    6. F=fft2(double(f));%傅里叶变换  
    7.   
    8. F=fftshift(F);%将变换的原点移到频率矩形的中心  
    9. [M,N]=size(f);  
    10.   
    11. %理想低通滤波  
    12.   
    13. D0=input('输入截止频率');  
    14. h1=zeros(M,N);  
    15. for i=1:M  
    16.     for j=i:N  
    17.         if(sqrt(((i-M/2)^2+(j-N/2)^2))<D0)  
    18.             h1(i,j)=1;  
    19.           
    20.         end  
    21.     end  
    22. end  
    23. G1=F.*h1;  
    24. G1=ifftshift(G1);  
    25. g1=real(ifft2(G1));  
    26.   
    27. %巴特沃斯低通滤波  
    28. n=input('巴特沃斯滤波器的阶数 n=');    
    29. n1=fix(M/2);  
    30. n2=fix(N/2);  
    31. h2=zeros(M,N);  
    32. for i=1:M  
    33.     for j=1:N  
    34.         d=sqrt((i-n1)^2+(j-n2)^2);  
    35.         h2=1./(1+(d./D0).^(2*n));  
    36.          
    37.     end  
    38. end  
    39. G2=F.*h2;  
    40. G2=ifftshift(G2);  
    41. g2=real(ifft(G2));  
    42.   
    43. %高斯低通滤波  
    44. h3=zeros(M,N);  
    45. for i=1:M  
    46.     for j=1:N  
    47.         h3=exp(-(d.^2)./(2*(D0^2)));  
    48.     end  
    49. end  
    50. G3=F.*h3;  
    51. G3=ifftshift(G3);  
    52. g3=real(ifft(G3));  
    53.   
    54. subplot(2,3,1);imshow(f);title('原图');  
    55. subplot(2,3,2);imshow(g1);title('理想低通滤波');  
    56. subplot(2,3,3);imshow(g2);title('巴特沃斯低通滤波');  
    57. subplot(2,3,4);imshow(g3);title('高斯低通滤波');  

    其运行效果如下:


     

       可见出现了问题,理想低通滤波的效果并不好,巴特沃斯和高斯滤波器的输出图像成了黑色,后经反复检查,调试,并没找出问题所在,于是打算尝试另一种方法。代码和结果如下:

    [html] view plain copy
    1. %巴特沃斯低通滤波  
    2. clc;close all;  
    3. f=imread('1.tif');  
    4. f=im2double(f);  
    5. M=2*size(I,1);  
    6. N=2*size(I,2);                        %滤波器的行列数    
    7. u=-M/2:(M/2-1);  
    8. v=-N/2:(N/2-1);  
    9. [U,V]=meshgrid(u,v);  
    10. D=sqrt(U.^2+V.^2);  
    11. D0=50;                                %截止频率  
    12. n=6;  
    13. H=1./(1+(D./D0).^(2*n));              %设计巴特沃斯滤波器  
    14. F=fftshift(fft2(I,size(H,1),size(H,2)));%傅里叶变换  
    15. G=F.*H;  
    16. L=ifft2(fftshift(G));                   %傅里叶反变换  
    17. L=L(1:size(I,1),1:size(I,2));  
    18. subplot(121);imshow(f);  
    19. subplot(122);imshow(L);  

    [html] view plain copy
    1. %高斯低通滤波  
    2. clc;close all;  
    3. I=imread('1.tif');    
    4. I=im2double(I);    
    5. M=2*size(I,1);  
    6. N=2*size(I,2);                        %滤波器的行列数    
    7. u=-M/2:(M/2-1);  
    8. v=-N/2:(N/2-1);  
    9. [U,V]=meshgrid(u,v);  
    10. D=sqrt(U.^2+V.^2);  
    11. D0=20;  
    12. H=exp(-(D.^2)./(2*(D0^2)));          %设计高斯滤波器  
    13. J=fftshift(fft2(I,size(H,1),size(H,2)));  
    14. G=J.*H;  
    15. L=ifft2(fftshift(G));  
    16. L=L(1:size(I,1),1:size(I,2));  
    17. figure;  
    18. subplot(121);imshow(I);  
    19. subplot(122);imshow(L);   


        通过这种方法运行出的结果是正常的。低通滤波滤掉了图像频谱中的高频成分,仅让低频部分通过,即变化剧烈的成分减少了,结果是使图像便模糊。比较可见高斯滤波的效果最好。

       频率域高通滤波和低通滤波实现的原理差不多,只是在设计滤波器时的公式有些差异,我就没有再重复实现了。

    展开全文
  • 频率域滤波

    2019-06-12 12:07:00
    1、频率域与空间域之间的关系  在频率域滤波基础——...而后再将频域滤波处理后的图像反变换回空间域,从而达到图像增强的目的,这样做的一个最主要的吸引力再域频率域滤波具有直观性的特点。  根据著名的卷...

     1、频率域与空间域之间的关系

      在频率域滤波基础——傅里叶变换计算及应用基础中,我们已经知道了如何将图像转换到频率域,以及如何将频率域图像通过傅里叶逆变换转换回图像,这样一来,可以利用空域图像与频谱之间的对应关系,尝试将空域卷积滤波变换为频域滤波,而后再将频域滤波处理后的图像反变换回空间域,从而达到图像增强的目的,这样做的一个最主要的吸引力再域频率域滤波具有直观性的特点。

      根据著名的卷积定律,两个二维连续函数再空间域的卷积可由其相应的两个傅里叶变换乘积的反变换而得到,反之,在频率域中的卷积可由在空间域中乘积的傅里叶变换而得到。即

    $$
    \begin{array}{l}{f_{1}(t) \leftrightarrow F_{1}(\omega), f_{2}(t) \leftrightarrow F_{2}(\omega)} \\ {f_{1}(t) * f_{2}(t) \leftrightarrow F_{1}(\omega) \cdot F_{2}(\omega)}\end{array}
    $$

    2、频率域滤波的基本步骤

      本文重点就来讲一讲如何通过频率域滤波来实现图像的平滑和锐化。首先将频率域滤波的步骤来做一个小结如下:

    • 给定一幅大小为MXN的输入图像f(x,y),选择填充参数P=2M.Q=2N.
    • 对f(x,y)添加必要数量的0,形成大小为PXQ的填充图像$f_{p}(x, y)$。
    • 用$(-1)^{x+y}$乘以$f_{p}(x, y)$,移到其变换中心。
    • 计算上图的DFT,得到F(u,v).
    • 生成一个实的,对称的滤波函数$H(u, v)$,其大小为PXQ,中心在$(P / 2, Q / 2)$处。用阵列相乘得到乘积$G(u, v)=H(u, v) F(u, v)$;即$G(i, k)=H(i, k) F(i, k)$
    • 得到处理后的图像:

      $$
      \mathrm{g}_{p}(x, y)=\left\{\text {real}\left[\mathfrak{J}^{-1}[G(u, v)]\right\}(-1)^{x+y}\right.
      $$

    其中,为忽略由于计算不准确导致的寄生复成分,选择实部,下标p指出我们处理的是填充后的图像。

    • 从$g_{p}(x, y)$的左上象限提取MXN区域,得到最终处理结果$g(x, y)$

      由上述叙述可知,滤波的关键取决于滤波函数$H(u, v)$,常常称为滤波器,或滤波器传递函数,因为它在滤波中抑制或除去了频谱中的某些分量,而保留其他的一些频率不受影响,从而达到滤波的目的。而本节将讲解一些常用的滤波器。

    3、使用频率域滤波平滑图像

      在空间域我们已经讲过图像平滑的目的及空间域平滑滤波,现在在频率域滤波中我们首先讲解三种频率平滑滤波器:理想滤波器、巴特沃斯滤波器、高斯滤波器。这三种滤波器涵盖了从非常急剧(理想)到非常平滑的滤波范围。

    1、理想低通滤波器(ILPF)

       首先最容易想到得衰减高频成分得方法就是在一个称为截止频率得位置截断所有的高频成分,将图像频谱中所有高于这一截止频率的频谱成分设为0,低于截止频率的成分保持不变。如图所示(左图是其透视图,右图是其图像表示):

     

      

      能够达到这种用效果的滤波器我们称之为理想滤波器,用公式表示为:

    $$
    H(u, v)=\left\{\begin{array}{ll}{1} & {\text { if } D(u, v) \leq D_{0}} \\ {0} & {\text { if } D(u, v)>D_{0}}\end{array}\right.
    $$

       其中,D0表示通带的半径。D(u,v)的计算方式也就是两点间的距离,很简单就能得到。

    $$
    D(u, v)=\sqrt{\left(u-\frac{P}{2}\right)^{2}+\left(v-\frac{Q}{2}\right)^{2}}
    $$

     

    算法在OpenCV中的实现如下所示:

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    
    
    int main()
    {
    	//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
    	//这是在ubuntu上运行的,p[1]为控制台给出的参数,即图片路径
       //如果不知道怎么传入参数,可以直接改为
    	Mat input=imread("2222.jpg",CV_LOAD_IMAGE_GRAYSCALE);
    	//image.jpg必须放在当前目录下,image.jpg即为输入图片
    	imshow("input", input);
    	int w = getOptimalDFTSize(input.cols);
    	int h = getOptimalDFTSize(input.rows);
    	Mat padded;
    	copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
    	padded.convertTo(padded, CV_32FC1);
    	imshow("padded", padded);
    	for (int i = 0; i < padded.rows; i++)//中心化操作,其余操作和上一篇博客的介绍一样
    	{
    		float *ptr = padded.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)	ptr[j] *= pow(-1, i + j);
    	}
    	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
    	Mat complexImg;
    	merge(plane, 2, complexImg);
    	dft(complexImg, complexImg);
    	//****************************************************
    	Mat mag = complexImg;
    	mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//这里为什么&上-2具体查看opencv文档
    	//其实是为了把行和列变成偶数 -2的二进制是11111111.......10 最后一位是0
    	//获取中心点坐标
    	int cx = mag.cols / 2;
    	int cy = mag.rows / 2;
    
    	float D0 = 20;
    	for (int y = 0; y < mag.rows; y++)
    	{
    		double* data = mag.ptr<double>(y);
    		for (int x = 0; x < mag.cols; x++)
    		{
    			double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
    			if (d <= D0) {
    				//小于截止频率不做操作
    			}
    			else {
    				data[x] = 0;
    			}
    		}
    	}
    
    	//******************************************************************
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("mag", plane[0]);
    
    	//******************************************************************
    	//*************************idft*************************************
    	idft(mag, mag);
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("idft-mag", plane[0]);
    	waitKey();
    	return 0;
    }
    

    从左至右依次是:原图,优化过后的图、傅里叶变换图、理想低通滤波器处理过后的高斯变换图、理想低通滤波器处理过后的效果图:

     

     

     

     

     

     

      

     2、巴特沃思低通滤波器(BLPF)

      截止频率位于距离原点$D_{0}$处的n阶巴特沃斯滤波器的传递函数定义为:

    $$
    H(u, v)=\frac{1}{1+\left[D(u, v) / D_{0}\right]^{2 n}}
    $$

      式中

    $$
    D(u, v)=\sqrt{\left(u-\frac{P}{2}\right)^{2}+\left(v-\frac{Q}{2}\right)^{2}}
    $$

      下图即为巴特沃斯滤波器的透视图及其图像显示:

      

     

       与理想低通滤波不同,巴特沃斯低通滤波器并没有再通过频率和滤除频率之间给出明显截止的急剧不确定性。对于具有平滑传递函数的滤波器,可在这样一点上定义截止频率,即使H(u,v)下降为其最大值的某个百分比的点。

    同理,我们可以在OpenCV下实现相关算法:

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    
    
    int main()
    {
    	//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
    	//这是在ubuntu上运行的,p[1]为控制台给出的参数,即图片路径
       //如果不知道怎么传入参数,可以直接改为
    	Mat input=imread("2222.jpg",CV_LOAD_IMAGE_GRAYSCALE);
    	//image.jpg必须放在当前目录下,image.jpg即为输入图片
    	imshow("input", input);
    	//Butterworth_Low_Paass_Filter(input);
    	int w = getOptimalDFTSize(input.cols);
    	int h = getOptimalDFTSize(input.rows);
    	Mat padded;
    	copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
    	padded.convertTo(padded, CV_32FC1);
    	imshow("padded", padded);
    	for (int i = 0; i < padded.rows; i++)//中心化操作,其余操作和上一篇博客的介绍一样
    	{
    		float *ptr = padded.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)	ptr[j] *= pow(-1, i + j);
    	}
    	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
    	Mat complexImg;
    	merge(plane, 2, complexImg);
    	//****************************************************
    	dft(complexImg, complexImg);
    	split(complexImg, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("complexImg", plane[0]);
    	//****************************************************
    	Mat mag = complexImg;
    	mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//这里为什么&上-2具体查看opencv文档
    	//其实是为了把行和列变成偶数 -2的二进制是11111111.......10 最后一位是0
    	//获取中心点坐标
    	int cx = mag.cols / 2;
    	int cy = mag.rows / 2;
    
    
    	int n = 10;//表示巴特沃斯滤波器的次数
    	//H = 1 / (1+(D/D0)^2n)
    
    	double D0 = 10;
    
    	for (int y = 0; y < mag.rows; y++) {
    		double* data = mag.ptr<double>(y);
    		for (int x = 0; x < mag.cols; x++) {
    			//cout << data[x] << endl;
    			double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
    			//cout << d << endl;
    			double h = 1.0 / (1 + pow(d / D0, 2 * n));
    			if (h <= 0.5) {
    				data[x] = 0;
    			}
    			else {
    				//data[x] = data[x]*0.5;
    				//cout << h << endl;
    			}
    
    			//cout << data[x] << endl;
    		}
    	}
    
    
    	//******************************************************************
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("mag", plane[0]);
    
    	//******************************************************************
    	//*************************idft*************************************
    	idft(mag, mag);
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("idft-mag", plane[0]);
    	waitKey();
    	return 0;
    }
    

      从左至右依次是:原图,优化过后的图、傅里叶变换图、巴特沃斯滤波器平滑后的高斯变换图、巴特沃斯滤波器平滑后的效果图:

      

     

    3、高斯低通滤波器(GLPF)

         高斯低通滤波的频率域二维形式由下式给出:

    $$
    H(u, v)=e^{-D^{2}(u, v) / 2 D_{0}^{2}}
    $$

      

      式中

    $$
    D(u, v)=\sqrt{\left(u-\frac{P}{2}\right)^{2}+\left(v-\frac{Q}{2}\right)^{2}}
    $$

      是离频率矩阵中心的距离。$D_{0}$是截止频率,当$D(u, v)=D_{0}$时,GLPF下降到其最大值的0.607处。

      高斯滤波器的透视图及其图像显示如下图所示:

      

    在OpenCV下实现高斯低通滤波算法如下:

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    
    
    
    
    int main()
    {
    	//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
    	//这是在ubuntu上运行的,p[1]为控制台给出的参数,即图片路径
       //如果不知道怎么传入参数,可以直接改为
    	Mat input=imread("2222.jpg",CV_LOAD_IMAGE_GRAYSCALE);
    	//image.jpg必须放在当前目录下,image.jpg即为输入图片
    	imshow("input", input);
    	int w = getOptimalDFTSize(input.cols);
    	int h = getOptimalDFTSize(input.rows);
    	Mat padded;
    	copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
    	padded.convertTo(padded, CV_32FC1);
    	imshow("padded", padded);
    	for (int i = 0; i < padded.rows; i++)//中心化操作,其余操作和上一篇博客的介绍一样
    	{
    		float *ptr = padded.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)	ptr[j] *= pow(-1, i + j);
    	}
    	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
    	Mat complexImg;
    	merge(plane, 2, complexImg);
    	dft(complexImg, complexImg);
    	//************************gaussian****************************
    	Mat gaussianBlur(padded.size(), CV_32FC2);
    	float D0 = 2 * 10 * 10;
    	for (int i = 0; i < padded.rows; i++)
    	{
    		float*p = gaussianBlur.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)
    		{
    			float d = pow(i - padded.rows / 2, 2) + pow(j - padded.cols / 2, 2);
    			p[2 * j] = expf(-d / D0);
    			p[2 * j + 1] = expf(-d / D0);
    		   
    		}
    	}
    	multiply(complexImg, gaussianBlur, gaussianBlur);//矩阵元素对应相乘法,注意,和矩阵相乘区分
    	//*****************************************************************	
    	split(complexImg, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("dft", plane[0]);
    	//******************************************************************
    	split(gaussianBlur, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("gaussianBlur", plane[0]);
    
    	//******************************************************************
    	//*************************idft*************************************
    	idft(gaussianBlur, gaussianBlur);
    	split(gaussianBlur, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("idft-gaussianBlur", plane[0]);
    	waitKey();
    	return 0;
    }
    

    从左至右依次是:原图,优化过后的图、傅里叶变换图、高斯低通滤波后的高斯变换图、高斯平滑过后的效果图:

        

     

    4、使用频率域滤波锐化图像

      前面已经讲过了通过衰减图像的傅里叶变换的高频成分可以平滑图像。因为边缘和其他灰度的急剧变化与高频成分有关,所以图像的锐化可在频率域通过滤波来实现,高通滤波器会衰减傅里叶变换中的低频成分而不扰乱高频信息。故我们可以考虑用高通滤波器来进行图像的锐化,一个高频滤波器是从给定低通滤波器用下式得到的:

    $$
    H_{\mathrm{HP}}(u, v)=1-H_{\mathrm{LP}}(u, v)
    $$

      与低通滤波器对应,这里介绍理想高通滤波器,巴特沃斯高通滤波器,高斯高通滤波器三种高通滤波器。

    1、理想高通滤波器 (IHPF)

      与低通对应,理想高通滤波器定义为:

    $$
    H(u, v)=\left\{\begin{array}{ll}{0} & {\text { if } D(u, v) \leq D_{0}} \\ {1} & {\text { if } D(u, v)>D_{0}}\end{array}\right.
    $$

       透视图和图像显示如下:

      

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    
    
    int main()
    {
    	//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
    	//这是在ubuntu上运行的,p[1]为控制台给出的参数,即图片路径
       //如果不知道怎么传入参数,可以直接改为
    	Mat input = imread("2222.jpg", CV_LOAD_IMAGE_GRAYSCALE);
    	//image.jpg必须放在当前目录下,image.jpg即为输入图片
    	imshow("input", input);
    	int w = getOptimalDFTSize(input.cols);
    	int h = getOptimalDFTSize(input.rows);
    	Mat padded;
    	copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
    	padded.convertTo(padded, CV_32FC1);
    	imshow("padded", padded);
    	for (int i = 0; i < padded.rows; i++)//中心化操作,其余操作和上一篇博客的介绍一样
    	{
    		float *ptr = padded.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)    ptr[j] *= pow(-1, i + j);
    	}
    	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
    	Mat complexImg;
    	merge(plane, 2, complexImg);
    	dft(complexImg, complexImg);
    	//****************************************************
    	Mat mag = complexImg;
    	mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//这里为什么&上-2具体查看opencv文档
    	//其实是为了把行和列变成偶数 -2的二进制是11111111.......10 最后一位是0
    	//获取中心点坐标
    	int cx = mag.cols / 2;
    	int cy = mag.rows / 2;
    
    	float D0 = 20;
    	for (int y = 0; y < mag.rows; y++)
    	{
    		double* data = mag.ptr<double>(y);
    		for (int x = 0; x < mag.cols; x++)
    		{
    			double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
    			if (d > D0) {
    				//小于截止频率不做操作
    			}
    			else {
    				data[x] = 0;
    			}
    		}
    	}
    
    	//******************************************************************
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("mag", plane[0]);
    
    	//******************************************************************
    	//*************************idft*************************************
    	idft(mag, mag);
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("idft-mag", plane[0]);
    	waitKey();
    	return 0;
    }
    

      结果如下:

      

    2、巴特沃思高通滤波器 (BHPF)

       同理n阶巴特沃斯高通滤波器的定义为

    $$
    H(u, v)=\frac{1}{1+\left[D_{0} / D(u, v)\right]^{2 n}}
    $$

     透视图和图像显示如下:

     

      

     

     

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    
    
    int main()
    {
    	//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
    	//这是在ubuntu上运行的,p[1]为控制台给出的参数,即图片路径
       //如果不知道怎么传入参数,可以直接改为
    	Mat input = imread("2222.jpg", CV_LOAD_IMAGE_GRAYSCALE);
    	//image.jpg必须放在当前目录下,image.jpg即为输入图片
    	imshow("input", input);
    	//Butterworth_Low_Paass_Filter(input);
    	int w = getOptimalDFTSize(input.cols);
    	int h = getOptimalDFTSize(input.rows);
    	Mat padded;
    	copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
    	padded.convertTo(padded, CV_32FC1);
    	imshow("padded", padded);
    	for (int i = 0; i < padded.rows; i++)//中心化操作,其余操作和上一篇博客的介绍一样
    	{
    		float *ptr = padded.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)    ptr[j] *= pow(-1, i + j);
    	}
    	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
    	Mat complexImg;
    	merge(plane, 2, complexImg);
    	//****************************************************
    	dft(complexImg, complexImg);
    	split(complexImg, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("complexImg", plane[0]);
    	//****************************************************
    	Mat mag = complexImg;
    	mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//这里为什么&上-2具体查看opencv文档
    	//其实是为了把行和列变成偶数 -2的二进制是11111111.......10 最后一位是0
    	//获取中心点坐标
    	int cx = mag.cols / 2;
    	int cy = mag.rows / 2;
    
    
    	int n = 1;//表示巴特沃斯滤波器的次数
    	//H = 1 / (1+(D/D0)^2n)
    
    	double D0 = 30;
    
    	for (int y = 0; y < mag.rows; y++) {
    		double* data = mag.ptr<double>(y);
    		for (int x = 0; x < mag.cols; x++) {
    			//cout << data[x] << endl;
    			double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
    			//cout << d << endl;
    			double h = 1.0 / (1 + pow(d / D0, 2 * n));
    			if (h > 0.5) {//低通变高通
    				data[x] = 0;
    			}
    			else {
    				//data[x] = data[x]*0.5;
    				//cout << h << endl;
    			}
    
    			//cout << data[x] << endl;
    		}
    	}
    
    
    	//******************************************************************
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("mag", plane[0]);
    
    	//******************************************************************
    	//*************************idft*************************************
    	idft(mag, mag);
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("idft-mag", plane[0]);
    	waitKey();
    	return 0;
    }
    

      结果如下:

        

     

    3、高斯高通滤波器(GHPF)

       高斯高通滤波器定义为:

    $$
    H(u, v)=1-\mathrm{e}^{-D^{2}(u, v) / 2 D_{0}^{2}}
    $$

     透视图和图像显示如下:

         

     

    OpenCV实现为:

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    
    int main()
    {
    	//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
    	//这是在ubuntu上运行的,p[1]为控制台给出的参数,即图片路径
       //如果不知道怎么传入参数,可以直接改为
    	Mat input=imread("2222.jpg",CV_LOAD_IMAGE_GRAYSCALE);
    	//image.jpg必须放在当前目录下,image.jpg即为输入图片
    	imshow("input", input);
    	int w = getOptimalDFTSize(input.cols);
    	int h = getOptimalDFTSize(input.rows);
    	Mat padded;
    	copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
    	padded.convertTo(padded, CV_32FC1);
    	imshow("padded", padded);
    	for (int i = 0; i < padded.rows; i++)//中心化操作,其余操作和上一篇博客的介绍一样
    	{
    		float *ptr = padded.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)	ptr[j] *= pow(-1, i + j);
    	}
    	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
    	Mat complexImg;
    	merge(plane, 2, complexImg);
    	dft(complexImg, complexImg);
    	//************************gaussian****************************
    	Mat gaussianSharpen(padded.size(), CV_32FC2);
    	float D0 = 2 * 10 * 10;
    	for (int i = 0; i < padded.rows; i++)
    	{
    		float*q = gaussianSharpen.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)
    		{
    			float d = pow(i - padded.rows / 2, 2) + pow(j - padded.cols / 2, 2);
    
    			q[2 * j] = 1 - expf(-d / D0);
    			q[2 * j + 1] = 1 - expf(-d / D0);
    		}
    	}
    	multiply(complexImg, gaussianSharpen, gaussianSharpen);
    	//*****************************************************************	
    	split(complexImg, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("dft", plane[0]);
    	//******************************************************************
    
    	split(gaussianSharpen, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("gaussianSharpen", plane[0]);
    	//******************************************************************
    	//*************************idft*************************************
    	idft(gaussianSharpen, gaussianSharpen);
    
    	split(gaussianSharpen, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("idft-gaussianSharpen", plane[0]);
    	waitKey();
    	return 0;
    }
    

      从左至右依次是:原图,优化过后的图、傅里叶变换图、高斯锐化后的高斯变换图、高斯锐化过后的效果图:

     

    6、选择率滤波器

       在很多应用中,图像处理的目的是处理指定频段或频率矩阵的小区域,此时我们需要使用选择滤波器对其进行处理,选择滤波器分为带通滤波器和带阻滤波器和陷波器。

    6.1、带阻滤波器、带通滤波器

      所谓的带通滤波器即是将低频滤波和高频滤波相结合,最后我们可以利用带阻滤波器滤掉周期性出现的噪声,由前面可知,

    理想带通滤波器公式为

    $$
    H(\mathrm{u}, v)=\left\{\begin{array}{l}{0 i f D_{0}-\frac{W}{2} \leq D \leq D_{0}+\frac{W}{2}} \\ {1}\end{array}\right.
    $$

     巴特沃斯带通滤波器公式为

    $$
    H(u, v)=\frac{1}{1+\left[\frac{D W}{D^{2}-D_{0}^{2}}\right]^{2 n}}
    $$

    高斯带通滤波器公式为:

    $$
    H(u, v)=1-\mathrm{e}^{-\left[\frac{D^{2}-D_{0}^{2}}{D W}\right]^{2}}
    $$

    与由低通滤波器得到高通滤波器相同的方法,由带阻滤波器我们可以得到带通滤波器

    $$
    H_{\mathrm{BP}}(u, v)=1-H_{\mathrm{BR}}(u, v)
    $$

     示例

    以下就是理想带通滤波器的OpenCV实现及效果图

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    
    
    int main()
    {
    	//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
    	//这是在ubuntu上运行的,p[1]为控制台给出的参数,即图片路径
       //如果不知道怎么传入参数,可以直接改为
    	Mat input = imread("2222.jpg", CV_LOAD_IMAGE_GRAYSCALE);
    	//image.jpg必须放在当前目录下,image.jpg即为输入图片
    	imshow("input", input);
    	int w = getOptimalDFTSize(input.cols);
    	int h = getOptimalDFTSize(input.rows);
    	Mat padded;
    	copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
    	padded.convertTo(padded, CV_32FC1);
    	imshow("padded", padded);
    	for (int i = 0; i < padded.rows; i++)//中心化操作,其余操作和上一篇博客的介绍一样
    	{
    		float *ptr = padded.ptr<float>(i);
    		for (int j = 0; j < padded.cols; j++)    ptr[j] *= pow(-1, i + j);
    	}
    	Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
    	Mat complexImg;
    	merge(plane, 2, complexImg);
    	dft(complexImg, complexImg);
    	//****************************************************
    	Mat mag = complexImg;
    	mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//这里为什么&上-2具体查看opencv文档
    	//其实是为了把行和列变成偶数 -2的二进制是11111111.......10 最后一位是0
    	//获取中心点坐标
    	int cx = mag.cols / 2;
    	int cy = mag.rows / 2;
    
    	float D0 = 20;
    	for (int y = 0; y < mag.rows; y++)
    	{
    		double* data = mag.ptr<double>(y);
    		for (int x = 0; x < mag.cols; x++)
    		{
    			double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
    			if ( D0<= d&& d <= D0+30 ) {
    			
    			}
    			else {
    				data[x] = 0;
    			}
    		}
    	}
    
    	//******************************************************************
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	plane[0] += Scalar::all(1);
    	log(plane[0], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("mag", plane[0]);
    
    	//******************************************************************
    	//*************************idft*************************************
    	idft(mag, mag);
    	split(mag, plane);
    	magnitude(plane[0], plane[1], plane[0]);
    	normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
    	imshow("idft-mag", plane[0]);
    	waitKey();
    	return 0;
    }
    

        

     

    6.2、陷波滤波器

       陷波器是一个相对于带通滤波更有用的选择滤波器,其拒绝事先定义好的一个关于频率矩阵中心的一个领域的频域。零相移滤波器必须是关于远点对称的,因此,一个中心位于$\left(u_{0}, v_{0}\right)$的陷波在位置$\left(-u_{0},-v_{0}\right)$必须有一个对应的陷波。陷波带阻滤波器可以用中心已被平移到陷波滤波器中心的高通滤波器的乘积来构造。一般形式为:

    $$
    H_{N R}(u, v)=\prod_{k=1}^{Q} H_{k}(u, v) H_{-k}(u, v)
    $$

     其中,$H_{k}(u, v) H_{-k}(u, v)$是高通滤波器,他们的中心分别位于$\left(u_{k}, v_{k}\right)$和$\left(-u_{k},-v_{k}\right)$处些中心是根据频率矩形的中心(M/2,N/2)确定的。对于每个滤波器,距离计算由下式执行:

    $$
    D_{k}(u, v)=\left[\left(u-M / 2-u_{k}\right)^{2}+\left(v-N / 2-v_{k}\right)^{2}\right]^{1 / 2}
    $$

    $$
    D_{-k}(u, v)=\left[\left(u-M / 2+u_{k}\right)^{2}+\left(v-N / 2+v_{k}\right)^{2}\right]^{1 / 2}
    $$

     

     

     关于摩尔纹,简单来说是差拍原理的一种表现。从数学上讲,两个频率接近的等幅正弦波叠加,合成信号的幅度将按照两个频率之差变化。如果感光元件CCD(CMOS)像素的空间频率与影像中条纹的空间频率接近,就会产生摩尔纹。我们手机对着电脑拍照的时候看到 的波纹即是摩尔纹,现在手机都自带摩尔纹消除了,所以这里无法用手机得到带摩尔纹的图片,这里给出数字图像处理的摩尔纹图片链接,带摩尔纹图片下载的链接为:http://www.imageprocessingplace.com/DIP-3E/dip3e_book_images_downloads.htm第四章的Fig0464(a)(car_75DPI_Moire).tif

      下面通过OpenCV来完成摩尔纹消除实验

    示例:

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    
    typedef cv::Mat Mat;
    
    Mat image_add_border(Mat &src)
    {
    	int w = 2 * src.cols;
    	int h = 2 * src.rows;
    	std::cout << "src: " << src.cols << "*" << src.rows << std::endl;
    
    	cv::Mat padded;
    	copyMakeBorder(src, padded, 0, h - src.rows, 0, w - src.cols,
    		cv::BORDER_CONSTANT, cv::Scalar::all(0));
    	padded.convertTo(padded, CV_32FC1);
    	std::cout << "opt: " << padded.cols << "*" << padded.rows << std::endl;
    	return padded;
    }
    
    //transform to center 中心化
    void center_transform(cv::Mat &src)
    {
    	for (int i = 0; i < src.rows; i++) {
    		float *p = src.ptr<float>(i);
    		for (int j = 0; j < src.cols; j++) {
    			p[j] = p[j] * pow(-1, i + j);
    		}
    	}
    }
    
    //对角线交换内容
    void zero_to_center(cv::Mat &freq_plane)
    {
    	//    freq_plane = freq_plane(Rect(0, 0, freq_plane.cols & -2, freq_plane.rows & -2));
    		//这里为什么&上-2具体查看opencv文档
    		//其实是为了把行和列变成偶数 -2的二进制是11111111.......10 最后一位是0
    	int cx = freq_plane.cols / 2; int cy = freq_plane.rows / 2;//以下的操作是移动图像  (零频移到中心)
    	cv::Mat part1_r(freq_plane, cv::Rect(0, 0, cx, cy));//元素坐标表示为(cx,cy)
    	cv::Mat part2_r(freq_plane, cv::Rect(cx, 0, cx, cy));
    	cv::Mat part3_r(freq_plane, cv::Rect(0, cy, cx, cy));
    	cv::Mat part4_r(freq_plane, cv::Rect(cx, cy, cx, cy));
    
    	cv::Mat tmp;
    	part1_r.copyTo(tmp);//左上与右下交换位置(实部)
    	part4_r.copyTo(part1_r);
    	tmp.copyTo(part4_r);
    	part2_r.copyTo(tmp);//右上与左下交换位置(实部)
    	part3_r.copyTo(part2_r);
    	tmp.copyTo(part3_r);
    }
    
    
    void show_spectrum(cv::Mat &complexI)
    {
    	cv::Mat temp[] = { cv::Mat::zeros(complexI.size(),CV_32FC1),
    					  cv::Mat::zeros(complexI.size(),CV_32FC1) };
    	//显示频谱图
    	cv::split(complexI, temp);
    	cv::Mat aa;
    	cv::magnitude(temp[0], temp[1], aa);
    	//    zero_to_center(aa);
    	cv::divide(aa, aa.cols*aa.rows, aa);
    	double min, max;
    	cv::Point min_loc, max_loc;
    	cv::minMaxLoc(aa, &min, &max,&min_loc, &max_loc);
    	std::cout << "min: " << min << " max: " << max << std::endl;
    	std::cout << "min_loc: " << min_loc << " max_loc: " << max_loc << std::endl;
    	cv::circle(aa, min_loc, 20, cv::Scalar::all(1), 3);
    	cv::circle(aa, max_loc, 20, cv::Scalar::all(1), 3);
    
    	cv::imshow("src_img_spectrum", aa);
    }
    
    //频率域滤波
    cv::Mat frequency_filter(cv::Mat &padded, cv::Mat &blur)
    {
    	cv::Mat plane[] = { padded, cv::Mat::zeros(padded.size(), CV_32FC1) };
    	cv::Mat complexIm;
    
    	cv::merge(plane, 2, complexIm);
    	cv::dft(complexIm, complexIm);//fourior transform
    	show_spectrum(complexIm);
    
    	cv::multiply(complexIm, blur, complexIm);
    	cv::idft(complexIm, complexIm, CV_DXT_INVERSE);//idft
    	cv::Mat dst_plane[2];
    	cv::split(complexIm,dst_plane);
    	center_transform(dst_plane[0]);
    	//    center_transform(dst_plane[1]);
    
    	cv::magnitude(dst_plane[0], dst_plane[1], dst_plane[0]);//求幅值(模)
    //    center_transform(dst_plane[0]);        //center transform
    
    	return dst_plane[0];
    }
    
    //陷波滤波器
    cv::Mat notch_kernel(cv::Mat &scr, std::vector<cv::Point> &notch_pot, float D0)
    {
    	cv::Mat notch_pass(scr.size(), CV_32FC2);
    	int row_num = scr.rows;
    	int col_num = scr.cols;
    	int n = 4;
    	for (int i = 0; i < row_num; i++) {
    		float *p = notch_pass.ptr<float>(i);
    		for (int j = 0; j < col_num; j++) {
    			float h_nr = 1.0;
    			for (unsigned int k = 0; k < notch_pot.size(); k++) {
    				int u_k = notch_pot.at(k).y;
    				int v_k = notch_pot.at(k).x;
    				float dk = sqrt(pow((i - u_k), 2) +
    					pow((j - v_k), 2));
    				float d_k = sqrt(pow((i + u_k), 2) +
    					pow((j + v_k), 2));
    				float dst_dk = 1.0 / (1.0 + pow(D0 / dk, 2 * n));
    				float dst_d_k = 1.0 / (1.0 + pow(D0 / d_k, 2 * n));
    				h_nr = dst_dk * dst_d_k * h_nr;
    				//                std::cout <<  "notch_pot: " << notch_pot.at(k) << std::endl;
    			}
    			p[2 * j] = h_nr;
    			p[2 * j + 1] = h_nr;
    
    		}
    	}
    
    	cv::Mat temp[] = { cv::Mat::zeros(scr.size(), CV_32FC1),
    					   cv::Mat::zeros(scr.size(), CV_32FC1) };
    	cv::split(notch_pass, temp);
    
    	std::string name = "notch滤波器d0=" + std::to_string(D0);
    	cv::Mat show;
    	cv::normalize(temp[0], show, 1, 0, CV_MINMAX);
    	cv::imshow(name, show);
    	return notch_pass;
    }
    
    std::string name_win("Notch_filter");
    cv::Rect g_rectangle;
    bool g_bDrawingBox = false;//是否进行绘制
    cv::RNG g_rng(12345);
    
    std::vector<cv::Point>notch_point;
    
    void on_MouseHandle(int event, int x, int y, int flags, void*param);
    void DrawRectangle(cv::Mat& img, cv::Rect box);
    
    int main()
    {
    	
    
    	Mat srcImage = cv::imread("Fig0464(a)(car_75DPI_Moire).tif", cv::IMREAD_GRAYSCALE);
    	if (srcImage.empty())
    		return -1;
    	imshow("src_img", srcImage);
    	Mat padded = image_add_border(srcImage);
    	center_transform(padded);
    	Mat plane[] = { padded, cv::Mat::zeros(padded.size(), CV_32FC1) };
    	Mat complexIm;
    
    	merge(plane, 2, complexIm);
    	cv::dft(complexIm, complexIm);//fourior transform
    	Mat temp[] = { cv::Mat::zeros(complexIm.size(),CV_32FC1),
    					  cv::Mat::zeros(complexIm.size(),CV_32FC1) };
    	//显示频谱图
    	split(complexIm, temp);
    	Mat aa;
    	magnitude(temp[0], temp[1], aa);
    	divide(aa, aa.cols*aa.rows, aa);
    
    
    	cv::namedWindow(name_win);
    	cv::imshow(name_win, aa);
    
    	cv::setMouseCallback(name_win, on_MouseHandle, (void*)&aa);
    	Mat tempImage = aa.clone();
    	int key_value = -1;
    	while (1) {
    		key_value = cv::waitKey(10);
    		if (key_value == 27) {//按下ESC键查看滤波后的效果
    			break;
    		}
    
    	}
    	if (!notch_point.empty()) {
    		for (unsigned int i = 0; i < notch_point.size(); i++) {
    			cv::circle(tempImage, notch_point.at(i), 3, cv::Scalar(1), 2);
    			std::cout << notch_point.at(i) << std::endl;
    		}
    	}
    	cv::imshow(name_win, tempImage);
    
    	Mat ker = notch_kernel(complexIm, notch_point, 30);
    	cv::multiply(complexIm, ker, complexIm);
    
    	split(complexIm, temp);
    	magnitude(temp[0], temp[1], aa);
    	divide(aa, aa.cols*aa.rows, aa);
    	imshow("aa", aa);
    	cv::idft(complexIm, complexIm, CV_DXT_INVERSE);//idft
    	cv::Mat dst_plane[2];
    	cv::split(complexIm, dst_plane);
    	center_transform(dst_plane[0]);
     //center_transform(dst_plane[1]);
    
      //   cv::magnitude(dst_plane[0],dst_plane[1],dst_plane[0]);  //求幅值(模)
    
    	cv::normalize(dst_plane[0], dst_plane[0], 1, 0, CV_MINMAX);
    	imshow("dest", dst_plane[0]);
    	cv::waitKey(0);
    
    	return 1;
    }
    
    
    void on_MouseHandle(int event, int x, int y, int falgs, void* param)
    {
    	Mat& image = *(cv::Mat*)param;
    
    	switch (event) {//鼠标移动消息
    	case cv::EVENT_MOUSEMOVE: {
    		if (g_bDrawingBox) {//标识符为真,则记录下长和宽到Rect型变量中
    
    			g_rectangle.width = x - g_rectangle.x;
    			g_rectangle.height = y - g_rectangle.y;
    		}
    	}
    							  break;
    
    	case cv::EVENT_LBUTTONDOWN: {
    		g_bDrawingBox = true;
    		g_rectangle = cv::Rect(x, y, 0, 0);//记录起点
    		std::cout << "start point( " << g_rectangle.x << "," << g_rectangle.y << ")" << std::endl;
    	}
    								break;
    
    	case cv::EVENT_LBUTTONUP: {
    		g_bDrawingBox = false;
    		bool w_less_0 = false, h_less_0 = false;
    
    		if (g_rectangle.width < 0) {//对宽高小于0的处理
    			g_rectangle.x += g_rectangle.width;
    			g_rectangle.width *= -1;
    			w_less_0 = true;
    
    		}
    		if (g_rectangle.height < 0) {
    			g_rectangle.y += g_rectangle.height;
    			g_rectangle.height *= -1;
    			h_less_0 = true;
    		}
    		std::cout << "end Rect[ " << g_rectangle.x << "," << g_rectangle.y << "," << g_rectangle.width << ","
    			<< g_rectangle.height << "]" << std::endl;
    
    		if (g_rectangle.height > 0 && g_rectangle.width > 0) {
    			Mat imageROI = image(g_rectangle).clone();
    			double min, max;
    			cv::Point min_loc, max_loc;
    			cv::minMaxLoc(imageROI, &min, &max, &min_loc, &max_loc);
    			cv::circle(imageROI, max_loc, 10, 1);
    			max_loc.x += g_rectangle.x;
    			max_loc.y += g_rectangle.y;
    
    			notch_point.push_back(max_loc);
    
    			cv::circle(image, max_loc, 10, 1);
    			//            cv::imshow( "ROI", imageROI );
    			cv::imshow("src", image);
    		}
    
    	}
    							  break;
    	}
    }

    运行程序后,用鼠标在频谱图上画圈,将四个陷波点标记后,在命令框内按下ESC键即可看到陷波结果:

     

     

    7、同态滤波

       图像的形成是由光源的入射和反射到成像系统中来的到物体的外观特征的一个过程,假如我们用形如f(x,y)的二维函数来表示图像,则f(x,y)可以由两个分量来表示:(1)入射到被观察场景的光源照射总量(2)场景中物体所反射光的总量。这两个分量分别被称为入射分量和反射分量,可以分别表示为i(x,y)和r(x,y)。两个函数的乘积为f(x,y)

    $$
    f(x, y)=i(x, y) r(x, y)
    $$

    其中$0<i(x, y)<\infty ,0<r(x, y)<1$,由于两个函数的乘积的傅里叶变换是不可分离的,我们不能轻易地使用上式分别对照明和反射的频率成分进行操作,即

    $$
    \mathcal{F}(F(x, y)) \neq \mathcal{F}(i(x, y)) \mathcal{F}(r(x, y))
    $$

    但是,假设我们定义

    $$
    \begin{aligned} z(x, y) &=\ln F(x, y) \\ &=\ln i(x, y)+\ln r(x, y) \end{aligned}
    $$

    则有

    $$
    \begin{aligned} \mathcal{F}(z(x, y)) &=\mathcal{F}(\ln F(x, y)) \\ &=\mathcal{F}(\ln i(x, y))+\mathcal{F}(\ln r(x, y)) \end{aligned}
    $$

    $$
    Z(\omega, \nu)=I(\omega, \nu)+R(\omega, \nu)
    $$

    其中Z,I,R分别是z,$\ln i(x, y)$和$\ln r(x, y)$的傅里叶变换。函数Z表示低频照明图像和高频反射图像两个图像的傅立叶变换,其传递函数如下所示:

    如果我们现在应用具有抑制低频分量和增强高频分量的传递函数的滤波器,那么我们可以抑制照明分量并增强反射分量。从而有

    $$
    \begin{aligned} S(\omega, \nu) &=H(\omega, \nu) Z(\omega, \nu) \\ &=H(\omega, \nu) I(\omega, \nu)+H(\omega, \nu) R(\omega, \nu) \end{aligned}
    $$

    其中S是傅里叶变换的结果。在空间领域

    $$
    \begin{aligned} s(x, y) &=\mathcal{F}^{-1}(S(\omega, \nu)) \\ &=\mathcal{F}^{-1}(H(\omega, \nu) I(\omega, \nu))+\mathcal{F}^{-1}(H(\omega, \nu) R(\omega, \nu)) \end{aligned}
    $$

    由定义:

    $$
    i^{\prime}(x, y)=\mathcal{F}^{-1}(H(\omega, \nu) I(\omega, \nu))
    $$

    $$
    r^{\prime}(x, y)=\mathcal{F}^{-1}(H(\omega, \nu) R(\omega, \nu))
    $$

     我们有

    $$
    s(x, y)=i^{\prime}(x, y)+r^{\prime}(x, y)
    $$

    最后,因为z(x,y)是通过去输入图像的对数形成的,我们可以通过取滤波后的而结果的指数这一反处理来形成输出图像:

    $$
    \begin{aligned} \hat{F}(x, y) &=\exp [s(x, y)] \\ &=\exp \left[i^{\prime}(x, y)\right] \exp \left[r^{\prime}(x, y)\right] \\ &=i_{0}(x, y) r_{0}(x, y) \end{aligned}
    $$

    以上便是同态滤波的整个流程,可以总结如下:

    \ begin {figure} \ par \ centerline {\ psfig {figure = figure55.ps,width = 5in}} \ par \ end {figure}

      图像的照射分量通常由慢的空间变化来表征,而反射分量往往引起突变,特别是在不同物体的连接部分。这些特性导致图像取对数后的傅里叶变换的低频成分与照射分量相联系,而高频成分与反射分量相联系。虽然这些联系只是粗略的近似,但它们用在图像滤波中是有益的。

      使用同态滤波器可更好的控制照射成分和反射成分。这种控制需要指定一个滤波器H(u,v),它可用不同的可控方法影响傅里叶变换的低频和高频成分。例如,采用形式稍微变化一下的高斯高通滤波器可得到函数:

    $$
    H(u, v)=\left(\gamma_{H}-\gamma_{L}\right)\left[1-e^{-c\left[D^{2}(u, v) / D_{0}^{2}\right]}\right]+\gamma_{L}
    $$

      其中$D(u, v)=\sqrt{\left(u-\frac{P}{2}\right)^{2}+\left(v-\frac{Q}{2}\right)^{2}}$,常数c控制函数边坡的锐利度。

    利用OpenCV的相关实现如下所示:

    示例:

     

    #include "stdafx.h"
    #include <math.h>
    #include <random>
    #include <iostream>
    #include <opencv2\core\core.hpp>
    #include <opencv2\highgui\highgui.hpp>
    #include <opencv2\imgproc\imgproc.hpp>
    
    using namespace cv;
    using namespace std;
    Mat Butterworth_Homomorphic_Filter(Size sz, int D, int n, float high_h_v_TB, float low_h_v_TB, Mat& realIm)
    {
    	Mat single(sz.height, sz.width, CV_32F);
    	cout << sz.width << " " << sz.height << endl;
    	Point centre = Point(sz.height / 2, sz.width / 2);
    	double radius;
    	float upper = (high_h_v_TB * 0.1);
    	float lower = (low_h_v_TB * 0.1);
    	long dpow = D * D;
    	float W = (upper - lower);
    
    	for (int i = 0; i < sz.height; i++)
    	{
    		for (int j = 0; j < sz.width; j++)
    		{
    			radius = pow((float)(i - centre.x), 2) + pow((float)(j - centre.y), 2);
    			float r = exp(-n * radius / dpow);
    			if (radius < 0)
    				single.at<float>(i, j) = upper;
    			else
    				single.at<float>(i, j) = W * (1 - r) + lower;
    		}
    	}
    
    	single.copyTo(realIm);
    	Mat butterworth_complex;
    	//make two channels to match complex
    	Mat butterworth_channels[] = { Mat_<float>(single), Mat::zeros(sz, CV_32F) };
    	merge(butterworth_channels, 2, butterworth_complex);
    
    	return butterworth_complex;
    }
    //DFT 返回功率谱Power
    Mat Fourier_Transform(Mat frame_bw, Mat& image_complex, Mat &image_phase, Mat &image_mag)
    {
    	Mat frame_log;
    	frame_bw.convertTo(frame_log, CV_32F);
    	frame_log = frame_log / 255;
    	/*Take the natural log of the input (compute log(1 + Mag)*/
    	frame_log += 1;
    	log(frame_log, frame_log); // log(1 + Mag)
    
    	/*2. Expand the image to an optimal size
    	The performance of the DFT depends of the image size. It tends to be the fastest for image sizes that are multiple of 2, 3 or 5.
    	We can use the copyMakeBorder() function to expand the borders of an image.*/
    	Mat padded;
    	int M = getOptimalDFTSize(frame_log.rows);
    	int N = getOptimalDFTSize(frame_log.cols);
    	copyMakeBorder(frame_log, padded, 0, M - frame_log.rows, 0, N - frame_log.cols, BORDER_CONSTANT, Scalar::all(0));
    
    	/*Make place for both the complex and real values
    	The result of the DFT is a complex. Then the result is 2 images (Imaginary + Real), and the frequency domains range is much larger than the spatial one. Therefore we need to store in float !
    	That's why we will convert our input image "padded" to float and expand it to another channel to hold the complex values.
    	Planes is an arrow of 2 matrix (planes[0] = Real part, planes[1] = Imaginary part)*/
    	Mat image_planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
    	/*Creates one multichannel array out of several single-channel ones.*/
    	merge(image_planes, 2, image_complex);
    
    	/*Make the DFT
    	The result of thee DFT is a complex image : "image_complex"*/
    	dft(image_complex, image_complex);
    
    	/***************************/
    	//Create spectrum magnitude//
    	/***************************/
    	/*Transform the real and complex values to magnitude
    	NB: We separe Real part to Imaginary part*/
    	split(image_complex, image_planes);
    	//Starting with this part we have the real part of the image in planes[0] and the imaginary in planes[1]
    	phase(image_planes[0], image_planes[1], image_phase);
    	magnitude(image_planes[0], image_planes[1], image_mag);
    
    	//Power
    	pow(image_planes[0], 2, image_planes[0]);
    	pow(image_planes[1], 2, image_planes[1]);
    
    	Mat Power = image_planes[0] + image_planes[1];
    
    	return Power;
    }
    void Inv_Fourier_Transform(Mat input, Mat& inverseTransform)
    {
    	/*Make the IDFT*/
    	Mat result;
    	idft(input, result, DFT_SCALE);
    	/*Take the exponential*/
    	exp(result, result);
    
    	vector<Mat> planes;
    	split(result, planes);
    	magnitude(planes[0], planes[1], planes[0]);
    	planes[0] = planes[0] - 1.0;
    	normalize(planes[0], planes[0], 0, 255, CV_MINMAX);
    
    	planes[0].convertTo(inverseTransform, CV_8U);
    }
    //SHIFT
    void Shifting_DFT(Mat &fImage)
    {
    	//For visualization purposes we may also rearrange the quadrants of the result, so that the origin (0,0), corresponds to the image center.
    	Mat tmp, q0, q1, q2, q3;
    
    	/*First crop the image, if it has an odd number of rows or columns.
    	Operator & bit to bit by -2 (two's complement : -2 = 111111111....10) to eliminate the first bit 2^0 (In case of odd number on row or col, we take the even number in below)*/
    	fImage = fImage(Rect(0, 0, fImage.cols & -2, fImage.rows & -2));
    	int cx = fImage.cols / 2;
    	int cy = fImage.rows / 2;
    
    	/*Rearrange the quadrants of Fourier image so that the origin is at the image center*/
    	q0 = fImage(Rect(0, 0, cx, cy));
    	q1 = fImage(Rect(cx, 0, cx, cy));
    	q2 = fImage(Rect(0, cy, cx, cy));
    	q3 = fImage(Rect(cx, cy, cx, cy));
    
    	/*We reverse each quadrant of the frame with its other quadrant diagonally opposite*/
    	/*We reverse q0 and q3*/
    	q0.copyTo(tmp);
    	q3.copyTo(q0);
    	tmp.copyTo(q3);
    
    	/*We reverse q1 and q2*/
    	q1.copyTo(tmp);
    	q2.copyTo(q1);
    	tmp.copyTo(q2);
    }
    void homomorphicFiltering(Mat src, Mat& dst)
    {
    	Mat img;
    	Mat imgHls;
    	vector<Mat> vHls;
    
    	if (src.channels() == 3)
    	{
    		cvtColor(src, imgHls, CV_BGR2HSV);
    		split(imgHls, vHls);
    		vHls[2].copyTo(img);
    	}
    	else
    		src.copyTo(img);
    
    	//DFT
    	//cout<<"DFT "<<endl;
    	Mat img_complex, img_mag, img_phase;
    	Mat fpower = Fourier_Transform(img, img_complex, img_phase, img_mag);
    	Shifting_DFT(img_complex);
    	Shifting_DFT(fpower);
    	//int D0 = getRadius(fpower,0.15);
    	int D0 = 10;
    	int n = 2;
    	int w = img_complex.cols;
    	int h = img_complex.rows;
    	//BHPF
    //  Mat filter,filter_complex;
    //  filter = BHPF(h,w,D0,n);
    //  Mat planes[] = {Mat_<float>(filter), Mat::zeros(filter.size(), CV_32F)};
    //  merge(planes,2,filter_complex);
    
    	int rH = 150;
    	int rL = 50;
    	Mat filter, filter_complex;
    	filter_complex = Butterworth_Homomorphic_Filter(Size(w, h), D0, n, rH, rL, filter);
    
    	//do mulSpectrums on the full dft
    	mulSpectrums(img_complex, filter_complex, filter_complex, 0);
    
    	Shifting_DFT(filter_complex);
    	//IDFT
    	Mat result;
    	Inv_Fourier_Transform(filter_complex, result);
    
    	if (src.channels() == 3)
    	{
    		vHls.at(2) = result(Rect(0, 0, src.cols, src.rows));
    		merge(vHls, imgHls);
    		cvtColor(imgHls, dst, CV_HSV2BGR);
    	}
    	else
    		result.copyTo(dst);
    
    }
    int main()
    {
    	Mat img = imread("test.jpg");
    	imshow("rogin img", img);
    	Mat dst;
    	homomorphicFiltering(img, dst);
    	imshow("homoMorphic filter", (img+dst)/2);
    	waitKey();
    	return 0;
    }
    

      运行结果如下所示:

     

     

     

     

     

     

    OpenCV 频率域处理:离散傅里叶变换、频率域滤波

    CS425 Lab: Frequency Domain Processing

    基于opencv的理想低通滤波器和巴特沃斯低通滤波器

    opencv 频率域滤波实例

    OpenCV 频率域实现钝化模板、高提升滤波和高频强调滤波 C++

    【数字图像处理】4.10:灰度图像-频域滤波 概论

    【数字图像处理】4.11:灰度图像-频域滤波 滤波器

    【数字图像处理】4.12:灰度图像-频域滤波 同态滤波

    基于OpenCV的同态滤波

    计算机视觉(六):频率域滤波器

    OpenCV使用陷波滤波器减少摩尔波纹 C++

    转载于:https://www.cnblogs.com/noticeable/p/10941851.html

    展开全文
  • 本节书摘来自异步社区《精通Matlab数字图像处理与识别》一书中的第6章,第6.1节,...精通Matlab数字图像处理与识别在很多情况下,频率域滤波和空间域滤波可以视为对于同一个图像增强问题的殊途同归的两种解决方式。...
  • 常用空间域滤波算法: 均值滤波 ...常用频率域滤波算法: 理想低通滤波 高通滤波 巴特沃斯低通滤波 DCT滤波 双边滤波 均值滤波算法: ##################################################...
  • 目录频率域滤波基础频率域的其他特性频率域滤波基础知识频率域滤波步骤小结空间域和频率域滤波之间的对应关系使用低通频率域滤波器平滑图像理想低通滤波器(ILPF)高斯低通滤波器(GLPF)巴特沃斯低通滤波器低通滤波的...
  •  空间域滤波是指在图像空间中借助模板对图像领域进行操作,处理图像每一个像素值。主要分为线性滤波和非线性滤波两类,根据功能可分为平滑滤波器和锐化滤波器。平滑可通过低通来实现,平滑的目的有两类,一是模糊,...
  • 1、频率域平滑滤波器(1) 理想低通滤波器(ILPF)(2) Butterworth低通滤波器(BLPF)(3) 高斯低通滤波器(GLPF)(常用)2、频率域锐化滤波器(1) 理想高通滤波器(IHPF)(2) Butterworth高通滤波器(BHPF)(3) 高斯高通滤波器...
  • 空间域滤波和频率域滤波 1.空间域滤波 空间域滤波是指在图像空间中借助模板对图像领域进行操作,处理图像每一个像素值。主要分为线性滤波和非线性滤波两类,根据功能可分为平滑滤波器和锐化滤波器。平滑可通过低通来...
  • 很多情况下,频率域滤波和空间域滤波可视为对于同一图像增强问题的殊途同归的两种解决方式。一些在空间域困难的增强任务,在频率域中变得非常普通。 傅里叶级数:任何周期函数可表示为不同频域的正弦和/或余弦之和的...
  • 文章目录前言一、钝化滤波以及高提升滤波1.计算步骤2.MATLAB代码二、高频强调滤波1.高频强调滤波器2.MATLAB代码 ...除此之外,还有同态滤波以及拉普拉斯滤波频率域图像增强方法。 一、钝化滤波以及高提升滤波 1.
  • 多数图像滤波例子是处理图像增强问题。
  • 数字图像处理(Digital Image Processing)是通过计算机对图像进行去除噪声、增强、复原、分割、提取特征等处理的方法和技术。本专栏将以学习笔记形式对数字图像处理的重点基础知识进行总结整理,欢迎大家一起学习交流...
  • 目录使用高通滤波器锐化图像由低通滤波器得到理想、高斯和巴特沃斯高通滤波器指纹增强频域中的拉普拉斯钝化掩蔽、高提升滤波和高频强调滤波同态滤波 使用高通滤波器锐化图像 由低通滤波器得到理想、高斯和巴特沃斯...
  • 频率域图像增强:以修改图像的傅里叶变换为基础的方法; 高、低通滤波,同态滤波。 另外还有彩色图像的增强。 1.空间域图像增强 1.1基于灰度变换的图像增强 属于点处理技术,这部分内容参考之前的笔记点运算部分...
  • 频率域图像增强

    千次阅读 2016-12-27 20:32:27
    频率域图像增强 高通滤波器和低通滤波器 本章的典型案例分析 利用频域滤波消除周期噪声 频域滤波基础频域滤波与空域滤波的关系 傅立叶变换可以将图像从空域变换到频域,而傅立叶反变换则可以将图像的频谱逆变换为...
  • 频率域图像增强

    千次阅读 2018-08-03 17:59:26
    1频率域滤波与空间域滤波殊途同归 空间域图像增强与频率域图像增强是两种截然不同的技术,实际上在相当程度上说它们是在不同的领域做相同的事情,是殊途同归的,,只是有些滤波更适合在空间域完成,而有些则更适合...
  • 频率域滤波与空间域滤波殊途同归,空间域图像增强与频率域图像增强是两种截然不同的技术,实际上在相当程度上说它们是在不同的领域做相同的事情,只是有些滤波更适合在空间域完成,而有些则更适合在频率域中完成。...
  • 为什么要在频率域研究图像增强 1.可以利用频率成分和图像外表之间的对应关系。一些在空间域表述困难的增强任务,在频率域中变得非常普通 2.滤波在频率域更为直观,它可以解释空间域滤波的某些性质 3.可以在频率域...
  • 1 频率域滤波 在频率域中,傅里叶变换可以做到转换过程不丢失任何信息。 2 傅里叶变换 傅里叶级数:满足狄利赫里的正弦函数都可以用正弦函数和余弦函数构成的无穷级数。 傅里叶变换:傅里叶变换在这里不做介绍,数字...
  • 频率域图像增强技术

    2019-09-27 12:38:39
    1、在图像中,像元的灰度值随位置变化的频繁程度可以用频率来表示,这是一种随位置变化的空间频率。 ...是指连续像元的灰度值的最高值与最低值的差。... 滤波是指在图像空间域(x,y)或者频率域...
  • 三、频率域滤波 四、二维离散变换的性质 五、实现 六、相关资源 序言:对图像而言,空间频率是指单位长度内亮度作为周期性变化的次数! 一、傅里叶变换基础 一维傅里叶变换数学推导 首先,我们知道傅里叶级数...
  • 1. 平滑的频率域滤波器 理想低通滤波器ILPF:半径为D0的圆内所有频率无衰减,之外全部衰减。 巴特沃斯滤波器BLPF:阶数越高越接近理想低通滤波器,振铃效应越大。一阶没有振铃效应,二阶微小,是有效的低通滤波...
  • (4)MATLAB 频率域图像增强

    千次阅读 2013-11-26 13:41:29
    频域滤波可以用来消除周期噪声。  傅里叶变换  图像较平滑,低频部分对应的幅值较大,图像灰度变化越剧烈,其频谱高频分量较强I1 = imread('cell.tif'); %读入原图像 fcoef = fft2(double(I1)); %做fft变换 ...
  • 卷积(Convolutions)滤波是通过消除特定的空间频率来使图像增强,根据增强类型(低频、中频和高频)不同可分为低通滤波、带通滤波和高通滤波。此外还可以增强图像某些方向特征的方向滤波等。它们的核心部分是卷积核...

空空如也

空空如也

1 2 3 4 5 6
收藏数 116
精华内容 46
关键字:

频率域滤波图像增强