图像处理中常用的矩阵知识

2016-07-04 10:38:28 qq_26499769 阅读数 23102


2.3.3 赋范空间


每个实数或复数,都有相对应的绝对值或者模,每一个n维矢量,也都可以定义其长度。如果把“长度”的概念推广到一般抽象空间中的元素上,就可以得到范数这个概念。



本节完。



2.3.6 希尔伯特空间


定义:在由内积所定义的范数意义下完备的内积空间称为希尔伯特(Hilbert)空间
希尔伯特空间是一类性质非常好的线性赋范空间,在工程上有着非常广泛的应用,而且在希尔伯特空间中最佳逼近问题可以得到比较完满的解决。




傅立叶变换以高等数学(微积分)中的傅立叶级数为基础发展而来,它是信号处理(特别是图像处理)中非常重要的一种时频变换手段,具有重要应用。在图像编码、压缩、降噪、数字水印方面都有重要意义。此外,快速傅立叶变换算法还位列20世纪十大算法之列,它是“动态规划”策略在算法设计中的杰出代表。本文将详细介绍图像中的傅立叶变换及其快速算法。对于傅立叶变换的数学原理还不是很理解的同学,建议参考本系列前面已经发布的傅立叶级数相关内容,争取彻底搞懂相关数学原理。一知半解、不求甚解,都是自欺欺人的表现。


6.1.2   数字图像的傅立叶变换

为了在科学计算和数字信号处理等领域使用计算机进行傅立叶变换,必须将函数f(t)定义在离散点而非连续域内,且须满足有限性或周期性条件。这种情况下,使用离散傅立叶变换。将连续函数f(t)等间隔采样就得到一个离散序列f(x),假设采样N次,则这个离散序列可以表示为{f(0),f(1),f(2),...,f(N-1)}。如果令x为离散实变量,u为离散频率变量,则一维离散傅立叶变换的正变换定义为



图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在平面空间上的梯度。从傅立叶频谱图上看到的明暗不一的亮点,实际上图像上某一点与邻域点差异的强弱,即梯度的大小,也即该点的频率的大小(可以这么理解,图像中的低频部分指低梯度的点,高频部分相反)。通常,梯度大则该点的亮度强,否则该点亮度弱。这样通过观察傅立叶变换后的频谱图,也叫功率图。在功率图中我们可以看出图像的能量分布,如果频谱图中暗的点数更多,那么实际图像是比较柔和的(因为各点与邻域差异都不大,梯度相对较小),反之,若频谱图中亮的点数多,那么实际图像一定是尖锐的,边界分明且边界两边像素差异较大的。对频谱移频到原点以后,可以看出图像的频率分布是以原点为圆心,对称分布的。变换最慢的频率成分(u = v = 0)对应一幅图像的平均灰度级。当从变换的原点移开时,低频对应着图像的慢变换分量,较高的频率开始对应图像中变化越来越快的灰度级。这些是物体的边缘和由灰度级的突发改变(如噪声)标志的图像成分。通常在进行傅立叶变换之前用(-1)^(x+y)乘以输入的图像函数,这样便可将傅立叶变换的原点(0,0)移到(M/2,N/2)上。


6.1.3  快速傅立叶变换的算法


离散傅立叶变换(DFT)已经成为数字信号处理和图像处理的一种重要手段,但是DFT的计算量太大,速度太慢,这令其实用性大打折扣。1965年,Cooley和Tukey提出了一种快速傅立叶变换算法(Fast Fourier Transform,FFT),极大地提供了傅立叶变换的速度。正是FFT的出现,才使得傅立叶变换得以广泛地应用。


FFT并不是一种新的变换,它只是傅立叶变换算法实现过程的一种改进。FFT中比较常用的是蝶形算法。蝶形算法主要是利用傅立叶变换的可分性、对称性和周期性来简化DFT的运算量。下面就来介绍一下蝶形算法的基本思想。


由于二维离散傅立叶变换具有可分离性, 即它可由两次一维离散傅立叶变换计算得到,因此仅研究一维离散傅立叶变换的快速算法即可。一维离散傅立叶变换的公式为





上述FFT是将f(x)序列按x的奇偶进行分组计算的,称之为时间抽选FFT。如果将频域序列的F(u)按u的奇偶进行分组计算,也可实现快速傅立叶计算,这称为频率抽选FFT。

通过对图6-6的观察可以发现,蝶形算法的频率域是按照正常顺序排列的,而空间域是按照一种叫做“码位倒序”的方式排列的。这个倒序的过程可以采用下面的方法来实现:将十进制的数转化成二进制,然后将二进制的序列倒序重排,最后再把颠倒顺序后的二进制数转换成十进制数。倒序重排的程序是一段经典程序,它以巧妙的构思、简单的语句用完成了倒序重排的功能。表6-1给出了倒序重排的示例。






6.4.3 主成分变换的实现


本小节通过一个算例验证一下之前的推导。在前面给出的例子中,各点在原始的

由于方程是齐次的,所以不独立。因为系数矩阵有零行列式,所以方程有非无效解。从两个方程的任何一个可见

现在考虑该结论该如何解释。特征向量g1和g2是在原坐标系中用来定义主成分轴的向量,如图6-20所示,其中,e1和e2分别是水平和垂直的方向向量。显而易见,这些数据在新坐标系中是非相关的。该新坐标系是原坐标系的旋转,出于这种原因,可以将主成分变换理解为旋转变换(即使在高维空间上亦是如此)。

6.4.4 基于K-L变换的图像压缩

从图像压缩的角度出发,我们必然希望变换系数协方差矩阵Σ中除对角线外的所有协方差均为零,成为对角线矩阵,即原来像素间的相关性经变换后全部解除,或者至少大部分协方差要等于或接近于零。为此,需要选择适当的变换矩阵,它作用于Σx 后使其变成对角线型。通过前面的分析和推导,可知这样的变换矩阵是存在的。如果用协方差矩阵Σx 的特征向量作变换的基向量,即由Σx 的特征向量作为正交变换的变换矩阵,就可以得到对角线型的变换域协方差矩阵Σy 。K-L变换就是采用这种矩阵进行变换的正交变换,它可以在变换域完全解除相关性,因此是理论上的最佳变换。同时,换一个角度也可以证明,K-L变换是均方误差最小准则下的最佳变换,即当压缩比确定的情况下,采用K-L变换后,重建图像的均方误差比采用任何其他正交变换的都小。

但是回顾之前进行的K-L变换,哪个步骤可以称为图像压缩的切入点呢?一幅大小为M×N的图像,它的协方差矩阵Σx大小为MN×MN。由上述K-L变换理论可知,对X进行K-L变换的变换矩阵就是Σx的特征向量矩阵,该矩阵大小亦为MN×MN,其大小远远大于原始图像数据矩阵。而且要在解码时恢复原图像,不但需要变换后的系数矩阵Y,还需要知道逆变换矩阵(也就是变换矩阵的转置)。如果不经过任何处理就这样直接将K-L变换用于数字图像的压缩编码,不但达不到任何数据压缩的效果,还极大的增加了数据量。即使仅保留一个最大的特征值,变换矩阵中和该特征值对应的特征向量为M×N维,系数矩阵 Y 保留的元素为一个。要重建图像数据,需要保留的元素个数为仍大于原矩阵,所以达不到压缩的目的。另外,求一个矩阵的协方差矩阵和特征向量矩阵,都是非常复杂的运算过程,需要大量的计算。当X比较大时,运算时间的耗用可能是非常大的。有时甚至会出现因为过于复杂而导致Σx和变换矩阵无法求解的情况。


要解决上述问题,可以考虑将图像分成若干个小块,然后对每个小块分别进行K-L变换(这与本章前面的处理方式基本保持一致)。这样使得Σx和变换矩阵都比较小,计算机处理起来比较容易而且速度快。这里仍然将图像划分为多个不重叠的8×8小块(当图像垂直和水平方向的像素数不是8的倍数时补0,使之均为8的倍数)。然后再分别对每一个小块执行K-L变换,变换矩阵的数目为K个,每个矩阵大小为64×64,仅变换矩阵就要记录K×64×64个数据,还是远远大于原始数据的个数M×N。是否可以让变换矩阵的数量变得少些,最好只保留一个变换矩阵。回忆前面做K-L变换的例子,变换矩阵的大小与输入矩阵的维度有关,而与样本数量无关,据此可以将每个8×8小块变成一个行向量(也就是一个64维的数组),原图中的每一个小方块都是一个64维的样本。所以最后只需要一个64×64的变换矩阵即可,它对于原图像的任意一个数据块都适用。这样的处理方式并不是完全意义上的K-L 变换,因为采用分块的处理方式,各个数据块之间的相关性是没有消除的。但实验亦表明,这样的K-L 变换虽然不能完全消除图像各像素点之间的相关性,也能达到很好的去相关效果,在去相关性性能上优于离散余弦变换。


图像数据经K-L变换后,得到的系数矩阵Y大部分系数都很小,接近于零。只有很少的几个系数的数值比较大,这正是K-L变换所起到的去除像素间的相关性,把能量集中分布在较少的变换系数上的作用的结果。据此,在图像数据压缩时,系数矩阵Y保留M个分量,其余分量则舍去。在实际编程开发中,经K-L变换后的系数矩阵中的数值都是按从大到小的顺序排列的,所以直接舍去后面的64M个分量即可。通过这一步的处理,便可动态地调节压缩编码系统的压缩比和重建图像的质量。解码时,首先做K-L逆变换,然后将上述过程逆转,可以得到重建后的图像数据矩阵。


我们在MATLAB中编写的示例程序演示了运用上述方法对图像实施基于K-L变换的压缩处理的过程。最后可以通过编程实现基于K-L变换的图像压缩算法并测试其压缩效果,所得之测试结果如图6-22所示。该程序验证了三种不同的压缩比,即舍去排序后的系数矩阵中的32/64(对应压缩比50%)、48/64(对应压缩比75%)以及56/64(对应压缩比87.5%)。相关测试程序源码读者可以从本书的在线支持资源中下载得到。


最后需要补充说明的是尽管K-L变换可以将数据之间的相关性完全去除,所以理论上是一种最理想的数据压缩方案,但它在实际应用过程中仍然受到很大局限。例如,它没有快速算法,不同的图像所对应的变换矩阵也不同,从这个角度来说,单纯将K-L变换直接应用于图像数据压缩的理论价值要大于实际价值。它的存在意义,一方面是可以作为理论验证的参考模板,另一方面就是需要对原始算法加以改进后再付诸应用。




6.4.2 主成分变换的推导

前面提到的一国经济增长与城市化水平关系的问题是典型二维问题,而协方差也只能处理二维问题,那维数多了自然就需要计算多个协方差,所以自然会想到使用矩阵来组织这些数据。为了帮助读者理解上面给出的协方差矩阵定义,在此举一个简单的三维的例子,假设数据集有 {x,y,z} 三个维度,则协方差矩阵为


可见,协方差矩阵是一个对称的矩阵,而且对角线是各个维度上的方差。下面通过一个例子来尝试演算协方差矩阵(很多数学软件都为该操作提供了支持)。需要提醒读者注意的是,协方差矩阵计算的是不同维度之间的协方差,而不是不同样本之间的。例如有一个样本容量为 9 的三维数据,如下


根据公式,计算协方差需要计算均值,那是按行计算均值还是按列呢,前面也特别强调了,协方差矩阵是计算不同维度间的协方差,要时刻牢记这一点。样本矩阵的每行是一个样本,每列为一个维度,所以要按列计算均值。经过计算,不难得到上述数据对应的协方差矩阵如下


众所周知,为了描述一个点在直角坐标系中的位置,至少需要两个分量。图6-17所示是两个二维数组,其中左图显示的各个点之间相关性微乎其微,而右图所示的各个点之间则高度相关,显然数据散布在一定角度内较为集中。对于右图而言,只要知道某个点一维分量的大小就可以大致确定其位置,两个分量中任一分量的增加或者减少都能引起另一分量相应的增减。相反,左图中的情况却不是这样。


对之前给出的协方差矩阵定义式稍加改写,以使其获得计算上更为直观的便利。则有在X矢量空间(或坐标系下),协方差矩阵Σx的无偏计算公式为


表6-2给出了对于图6-17中左图所示的6个样本点的集合,以及经计算后求得的样本集协方差矩阵和相关矩阵的结果。应当注意,协方差矩阵和相关矩阵二者都是沿对角线对称的。从相关矩阵来看,各个数据分量间存在不相关关系的明显事实就是协方差矩阵(以及相关矩阵)中非对角线元素都是零。

最终计算可得





1.1.3 函数的极限


本小节介绍两个重要的函数极限,并讨论它们的应用。

重要极限1:


此外,该重要极限的另一种形式也常常被用到,即


综上,结论得证。
由此,也很容易推出如下结论,证明从略,有兴趣的读者可以自行尝试推导




1.3.2 内积与外积


因为cos(π/2)=0。当然,这也是众多教科书上介绍向量内积最开始时常常用到的一种定义方式。但必须明确,这种表示方式仅仅是一种非常狭隘的定义。如果从这个定义出发来介绍向量内积,其实是本末倒置的。因为对于高维向量而言,夹角的意义是不明确的。例如,在三维坐标空间中,再引入一维时间坐标,形成一个四维空间,那么时间向量与空间向量的夹角该如何解释呢?所以读者务必明确,首先应该是给出如本小节最开始时给出的内积定义,然后才能由此给出二维或三维空间下的夹角定义。在此基础上,我们来证明余弦定律。



若根据a·b = |a||b|cosθ这个定义,因为0<=cosθ<=1,显然柯西-施瓦茨不等式是成立的。但是这样的证明方式同样又犯了本末倒置的错误。柯西-施瓦茨不等式并没有限定向量的维度,换言之它对于任意维度的向量都是成立的,这时夹角的定义是不明确的。正确的思路同样应该从本小节最开始的定义出发来证明柯西-施瓦茨不等式,因为存在这样一个不等式关系,然后我们才会想到内积与向量模的乘积之间存在一个介于0和1之间的系数,然后我们才用cosθ来表述这个系数,于是才会得到a·b = |a||b|cosθ这个表达式。下面就来证明柯西-施瓦茨不等式。


证明:

与内积类似,向量a,b的外积也可以狭义地定义为






1.4.5   卷积定理及其证明


卷积定理是傅立叶变换满足的一个重要性质。卷积定理指出,函数卷积的傅立叶变换是函数傅立叶变换的乘积。换言之,一个域中的卷积对应于另一个域中的乘积,例如,时域中的卷积对应于频域中的乘积。


这一定理对拉普拉斯变换、Z变换等各种傅立叶变换的变体同样成立。需要注意的是,以上写法只对特定形式的变换正确,因为变换可能由其它方式正规化,从而使得上面的关系式中出现其它的常数因子。
下面我们来证明时域卷积定理,频域卷积定理的证明与此类似,读者可以自行证明。
证明:将卷积的定义


傅立叶变换的作用在频域对信号进行分析,我们可以把时域的信号看做是若干正弦波的线性叠加,傅立叶变换的作用正是求得这些信号的幅值和相位。既然固定的时域信号是若干固定正弦信号的叠加,在不改变幅值的情况下,在时间轴上移动信号,也就相当于同时移动若干正弦信号,这些正弦信号的相位改变、但幅值不变,反映在频域上就是傅立叶变换结果的模不变、而相位改变。所以,时移性质其实就表明当一个信号沿时间轴平移后,各频率成份的大小不发生改变,但相位发生变化。
既然这里提到了傅立叶变换的性质,这里我们还将补充一些关于帕塞瓦尔定理的有关内容。该定理最早是由法国数学家帕塞瓦尔(Marc-Antoine Parseval)在1799年推导出的一个关于级数的理论,该定理随后被应用于傅立叶级数。帕塞瓦尔定理的表述是这样的:




综上所述,原结论得证。
前面我们也介绍过复数形式的傅立叶级数,下面我们来推导与复数形式傅立叶变换相对应的帕塞瓦尔等式。这里再次给出傅立叶级数的复数形式表达式,具体推导过程请读者参阅前文



帕塞瓦尔定理把一个信号的能量或功率的计算和频谱函数或频谱联系起来了,它表明一个信号所含有的能量(功率)恒等于此信号在完备正交函数集中各分量能量(功率)之和。换言之,能量信号的总能量等于各个频率分量单独贡献出来的能量的连续和;而周期性功率信号的平均功率等于各个频率分量单独贡献出来的功率之和。


1.1.2 级数的敛散


关于上面这个级数敛散性的讨论,在数学史上曾经是一个非常有名的问题。大数学家莱布尼兹曾经在惠更斯的指导下对级数的敛散性进行过研究。后来莱布尼兹的学生伯努利兄弟(雅各·伯努利和约翰·伯努利)从他们老师的某些研究成果出发,最终证明了调和级数的发散性,以及几何级数的收敛性。但是几何级数最终收敛到多少这个问题却一直困扰着他们。最终,雅各布也不得不带着几分绝望的恳求宣告了他的失败:“如果有人能够发现并告知我们迄今为止尚未解出的难题的答案,我们将不胜感谢。”所幸的是,几何级数到底等于多少这个难题最终被约翰·伯努利的学生欧拉所破解。欧拉使用了一种极其巧妙的方法得出





1.1 极限及其应用


极限的概念是微积分理论赖以建立的基础。在研究极限的过程中,我们一方面会证明许多在图像处理中将要用到的公式,另一方面还会得到所谓的自然常数(或称纳皮尔常数)。图像处理技术中的很多地方都会遇到它,例如用来对图像进行模糊降噪的高斯函数,以及泊松噪声中都会有自然常数出现。而且在本文稍往后的内容还会讲到欧拉公式,届时自然常数还将会再次出现。

1.1.1 数列的极限






1.3.7 曲面积分





关于这部分内容的讨论,既阐明了第二类曲面积分的实际意义,其实也明确了两类曲面积分之间的关联。需要说明的是,在后面的介绍中,我们将更多地采用通量这个提法来替代此前所用的流量。通量是更广义的说法,如果考虑的向量场是流速场的话,那么通量就是流量,如果考虑的是电场或者磁场的话,那么通量就是电通量或者磁通量。






在泛函分析中,索伯列夫空间并不像 巴拿赫空间或者希尔伯特空间那么引入注意。但是在图像处理中,索伯列夫空间在介绍BV空间(有界变差函数空间)时,会被提到。而BV函数空间对于理解TV算法(偏微分方程在图像处理中的重要内容)至关重要!所以我特别在“图像处理中的数学原理详解”系列文章中留出一个小节来对索伯列夫空间进行必要的介绍。


2.3.7 索伯列夫空间



由广义导数的定义可以看出,这种导数不是关于函数的个别点处局部性质反映,因为它是通过在整个区间上积分的极限来确定的,而积分是一种关于函数的整体性质的概念。但也应该指出,广义导数其实是对通常意义下导数概念的推广。如果函数本身是通常意义下可微的,则其导函数与广义导数是一致的。






2.3.5 内积空间

      前面我们已经讨论过关于内积的话题,此处以公理化的形式给出内积的定义。





2.3.2 距离空间
        尽管在线性空间上我们已经可以完成简单的线性运算,但这仍然不能满足我们的需求。为了保证数学刻画的精确性,还必须引入距离的概念。本文最初是从极限开始讲起的,它是因此微积分的必备要素之一,而极限的概念显然也是基于距离上无限接近这样一种角度来描述的。



       由此,在距离空间中,可以引入“任意逼近”的概念,即极限概念。一般来说,一个集合如果能够在其中确切地引入任意逼近的概念,就称之为“拓扑空间”。而距离空间是一种最常用的拓扑空间。





2.3  泛函与抽象空间

牛顿说:“把简单的问题看得复杂,可以发现新领域;把复杂的问题看得简单,可以发现新规律。”而从历史的角度来看,一个学科的发展也亦是如此。随着学科的发展,最开始的一个主干方向会不断衍生出各自相对独立的分支,这也就是所谓“把简单的问题看得复杂”的过程。然而,一旦学科发展到一定程度之后,某些分支学科又开始被抽象综合起来,这也就是所谓“把复杂的问题看得简单”的过程。例如,在很长一段时间里,物理学家们都把电和磁看成是两种独立的物理现象在研究,当学科研究积累到一定程度时,麦克斯韦就创立了电磁学从而完成了物理学中的一次大综合。而在数学发展的历史中,几何与代数也曾经在很长的一段时间里是彼此独立的。直到笛卡尔引入了直角坐标系的概念之后,人们才开始建立了一种代数与几何之间的联系,也就是所谓的解析几何。泛函分析也是对以往许多数学问题或者领域进行高度抽象和综合的结果,其主要研究对象之一是抽象空间。其实在学习线性代数的过程中,读者已经建立了一种从矩阵到线性方程组之间的一种联系。而在泛函分析中,实数系、矩阵、多项式以及函数族这些看似关联不大的概念都可以抽成空间。由于泛函分析是一门比较晦涩抽象的学问,读者应该注意联系以往学习中比较熟悉的一些已知的、具体的概念,从而帮助自己理解那些全新的、抽象的概念。此外,需要说明的是本部分内容的重点在于有关定义或者概念的介绍,希望读者能够努力领会这些定义或者概念。

2.3.1  线性空间



2018-07-12 17:49:37 dfdfdsfdfdfdf 阅读数 5807

fishing-panhttps://blog.csdn.net/u013921430转载请注明出处】

前言

       Hessian Matrix(海森矩阵)在图像处理中有广泛的应用,比如边缘检测、特征点检测等。而海森矩阵本身也包含了大量的数学知识,例如泰勒展开、多元函数求导、矩阵、特征值等。写这篇博客的目的,就是想从原理和实现上讲一讲Hessian Matrix,肯定有不足的地方,希望大家批评指正。

泰勒展开及海森矩阵

      将一个一元函数f(x)在x0处进行泰勒展开,可以得到以下公式。


       其中余项为皮亚诺余项

       其中二阶导数的部分映射到二维以及多维空间就是Hessian Matrix。在二维图像中,假设图像像素值关于坐标(x, y)的函数是f(x, y),那么将f(x+dx,y+dy)在f(x0, y0)处展开,得到如下式子;


     如果将这个式子用矩阵表示,并且舍去余项,则式子会是下面这个样子。


      上面等式右边的第三项中的第二个矩阵就是二维空间中的海森矩阵了;从而有了一个结论,海森矩阵就是空间中一点处的二阶导数。进而推广开来,多维空间中的海森矩阵可以表示为。


 

海森矩阵的意义

     众所周知,二阶导数表示的导数的变化规律,如果函数是一条曲线,且曲线存在二阶导数,那么二阶导数表示的是曲线的曲率,曲率越大,曲线越是弯曲。以此类推,多维空间中的一个点的二阶导数就表示该点梯度下降的快慢。以二维图像为例,一阶导数是图像灰度变化即灰度梯度,二阶导数就是灰度梯度变化程度,二阶导数越大灰度变化越不具有线性性(这里有一点绕口了,意思就是这里灰度梯度改变越大,不是线性的梯度)。

      但是在二维图像中,海森矩阵是二维正定矩阵,有两个特征值和对应的两个特征向量。两个特征值表示出了图像在两个特征向量所指方向上图像变化的各向异性。如果利用特征向量与特征值构成一个椭圆,那么这个椭圆就标注出了图像变化的各向异性。那么在二维图像中,什么样的结构最具各项同性,又是什么样的结构的各向异性更强呢?很显然,圆具有最强的各项同性,线性越强的结构越具有各向异性。如下图;


注:图中箭头的方向不一定正确,我只是随意标注

       且特征值应该具有如下特性;

λ1

λ2

图像特征

-High

-High

斑点结构(前景为亮)

+High

+High

斑点结构(前景为暗)

Low

-High

线性结构(前景为亮)

Low

+High

线性结构(前景为暗) 

海森矩阵的应用

       海森矩阵的应用很多,我在此不一一列举。上文中提到了矩阵的特征值与特征向量构成的椭圆表现出了图像的各向异性,这种各项异性图像处理就得到了应用。以二维图像为例,图像中的点性结构具有各项同性,而线性结构具有各向异性。因此我们可以利用海森矩阵对图像中的线性结构进行增强,滤去点状的结构和噪声点。同样,也可以用于找出图像中的点状结构,滤除其他信息。

       当然,我们在使用海森矩阵时,不需要把图像进行泰勒展开,我们只需要直接求取矩阵中的元素即可。一般,对数字图像进行二阶求导使用的是以下方法;

       

      但是这种方法鲁棒性很差,容易受到图像中局部信号的干扰,甚至可以说,这种求导方式是存在争议的。因为这一点的二阶导数也可以采用如下方法表示;

                                              

       也可以采用其他方式表示,而且往往不同的方法求得的值不同,因为这种方法只把包含其自身在内的三个点的信息囊括进去,信息量不足。因此,提倡大家换一种方法。根据线性尺度空间理论(LOG),对一个函数求导,等于函数与高斯函数导数的卷积。式子如下;


       由于高斯模板可以将周围一矩形范围内所有的点的信息都包含进来,这样就不会有误差。所以利用图像求取hessian矩阵中的元素时,将图像与高斯函数的二阶导数做卷积即可。在此,为大家提供高斯函数的二阶偏导。



      在编写程序时,我们只需要事先将图像分别于三个模板进行卷积,生成三种偏导数的“图”,然后每次根据需要索引对应位置的偏导数即可。

代码

        需要说明的是,我在代码中运用了一些OpenCV的函数;

  1. ///---------------------------///
  2. //---海森矩阵二维图像增强
  3. //---潘正宇---2018.03.27
  4. #include <iostream>
  5. #include <vector>
  6. #include<opencv2/opencv.hpp>
  7. #include<opencv2/highgui/highgui.hpp>
  8. #include<opencv2/imgproc/imgproc.hpp>
  9. #include <map>
  10. #define STEP 6
  11. #define ABS(X) ((X)>0? X:(-(X)))
  12. #define PI 3.1415926
  13. using namespace std;
  14. using namespace cv;
  15. int main()
  16. {
  17. Mat srcImage = imread("G:\\博客\\图像处理\\hessian\\hessian_matrix\\Multiscale vessel enhancement filtering1.bmp");
  18. if (srcImage.empty())
  19. {
  20. cout << "图像未被读入";
  21. system("pause");
  22. return 0;
  23. }
  24. if (srcImage.channels() != 1)
  25. {
  26. cvtColor(srcImage, srcImage, CV_RGB2GRAY);
  27. }
  28. int width = srcImage.cols;
  29. int height = srcImage.rows;
  30. Mat outImage(height, width, CV_8UC1,Scalar::all(0));
  31. int W = 5;
  32. float sigma = 01;
  33. Mat xxGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));
  34. Mat xyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));
  35. Mat yyGauKernel(2 * W + 1, 2 * W + 1, CV_32FC1, Scalar::all(0));
  36. //构建高斯二阶偏导数模板
  37. for (int i = -W; i <= W;i++)
  38. {
  39. for (int j = -W; j <= W; j++)
  40. {
  41. xxGauKernel.at<float>(i + W, j + W) = (1 - (i*i) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * PI*pow(sigma, 4)));
  42. yyGauKernel.at<float>(i + W, j + W) = (1 - (j*j) / (sigma*sigma))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(-1 / (2 * PI*pow(sigma, 4)));
  43. xyGauKernel.at<float>(i + W, j + W) = ((i*j))*exp(-1 * (i*i + j*j) / (2 * sigma*sigma))*(1 / (2 * PI*pow(sigma, 6)));
  44. }
  45. }
  46. for (int i = 0; i < (2 * W + 1); i++)
  47. {
  48. for (int j = 0; j < (2 * W + 1); j++)
  49. {
  50. cout << xxGauKernel.at<float>(i, j) << " ";
  51. }
  52. cout << endl;
  53. }
  54. Mat xxDerivae(height, width, CV_32FC1, Scalar::all(0));
  55. Mat yyDerivae(height, width, CV_32FC1, Scalar::all(0));
  56. Mat xyDerivae(height, width, CV_32FC1, Scalar::all(0));
  57. //图像与高斯二阶偏导数模板进行卷积
  58. filter2D(srcImage, xxDerivae, xxDerivae.depth(), xxGauKernel);
  59. filter2D(srcImage, yyDerivae, yyDerivae.depth(), yyGauKernel);
  60. filter2D(srcImage, xyDerivae, xyDerivae.depth(), xyGauKernel);
  61. for (int h = 0; h < height; h++)
  62. {
  63. for (int w = 0; w < width; w++)
  64. {
  65. //map<int, float> best_step;
  66. /* int HLx = h - STEP; if (HLx < 0){ HLx = 0; }
  67. int HUx = h + STEP; if (HUx >= height){ HUx = height - 1; }
  68. int WLy = w - STEP; if (WLy < 0){ WLy = 0; }
  69. int WUy = w + STEP; if (WUy >= width){ WUy = width - 1; }
  70. float fxx = srcImage.at<uchar>(h, WUy) + srcImage.at<uchar>(h, WLy) - 2 * srcImage.at<uchar>(h, w);
  71. float fyy = srcImage.at<uchar>(HLx, w) + srcImage.at<uchar>(HUx, w) - 2 * srcImage.at<uchar>(h, w);
  72. float fxy = 0.25*(srcImage.at<uchar>(HUx, WUy) + srcImage.at<uchar>(HLx, WLy) - srcImage.at<uchar>(HUx, WLy) - srcImage.at<uchar>(HLx, WUy));*/
  73. float fxx = xxDerivae.at<float>(h, w);
  74. float fyy = yyDerivae.at<float>(h, w);
  75. float fxy = xyDerivae.at<float>(h, w);
  76. float myArray[2][2] = { { fxx, fxy }, { fxy, fyy } }; //构建矩阵,求取特征值
  77. Mat Array(2, 2, CV_32FC1, myArray);
  78. Mat eValue;
  79. Mat eVector;
  80. eigen(Array, eValue, eVector); //矩阵是降序排列的
  81. float a1 = eValue.at<float>(0, 0);
  82. float a2 = eValue.at<float>(1, 0);
  83. if ((a1>0) && (ABS(a1)>(1+ ABS(a2)))) //根据特征向量判断线性结构
  84. {
  85. outImage.at<uchar>(h, w) = pow((ABS(a1) - ABS(a2)), 4);
  86. //outImage.at<uchar>(h, w) = pow((ABS(a1) / ABS(a2))*(ABS(a1) - ABS(a2)), 1.5);
  87. }
  88. }
  89. }
  90. //----------做一个闭操作
  91. Mat element = getStructuringElement(MORPH_RECT, Size(3, 2));
  92. morphologyEx(outImage, outImage, MORPH_CLOSE, element);
  93. imwrite("temp.bmp", outImage);
  94. imshow("[原始图]", outImage);
  95. waitKey(0);
  96. system("pause");
  97. return 0;
  98. }

输出结果

      左侧是原图,中间是增强后的结果,右侧是将增强后的结果做了二值化的结果图;

                                        

结果分析

      从结果来看,图像中的大部分的线性结构都被增强了,但是有一些细微的结构并未被增强太多,而且有些粗的结构中出现了空洞,其实这都与求导窗口的大小有关,求导窗口太小,很多粗的结构会出现中空的现象,因为中心区域被认为是点结构了;求导窗口太大,就容易出现细微结构丢失的情况。此外,高斯模板的方差选取也影响了偏导数的大小。其实,这种方式是使用一个方差为s 的高斯核的二阶导数生成一个探测核,用于测量导数方向范围内(-s,s)内外区域之间的对比度。


      但是同一图像中,线性结构的粗细肯定是不同的,同样的窗口大小是无法全部适用的,针对上面的问题,有人提出了多模板的方法。即对一个点用多种尺度的高斯模板进行卷积,然后选择各向异性最强的结果作为该点的输出。值得一试。

      回过头来看,根据海森矩阵,还可以确定一张图像中的角点部分,即前面表格中提到的两个特征值的绝对值都较大的情况。其实这就是Harris 角点检测的主要思想。

最后

      这篇博客准备了一段时间了,主要是公式的原因,我想使用markdown编辑器以及LaTex公式编辑器编辑公式,但是没搞定,还需要琢磨一下,等我琢磨好了,就把这篇文章中的图片格式的公式换掉。


参考文献:

     Frangi A.F., Niessen W.J., Vincken K.L., Viergever M.A. (1998) Multiscale vessel enhancement filtering. In: Wells W.M., Colchester A., Delp S. (eds) Medical Image Computing and Computer-Assisted Intervention — MICCAI’98. MICCAI 1998. Lecture Notes in Computer Science, vol 1496. Springer, Berlin, Heidelberg

2016-04-21 20:44:32 yangdashi888 阅读数 2784

6.5 矩阵的运算及其运算规则

一、矩阵的加法与减法


  1、运算规则
  设矩阵
  则
     
  简言之,两个矩阵相加减,即它们相同位置的元素相加减!
  注意:只有对于两个行数、列数分别相等的矩阵(即同型矩阵),加减法运算才有意义,即加减运算是可行的.

  2、 运算性质 (假设运算都是可行的)
  满足交换律和结合律
  交换律 
  结合律 

二、矩阵与数的乘法


  1、 运算规则
  乘矩阵A,就是将数乘矩阵A中的每一个元素,记为
  特别地,称称为的负矩阵.
  2、 运算性质
  满足结合律和分配律
  结合律: (λμ)A=λ(μA) ; (λ+μ)A =λA+μA
  分配律: λ (A+B)=λA+λB

  典型例题
  例6.5.1 已知两个矩阵
  满足矩阵方程,求未知矩阵
   由已知条件知
    
    

三、矩阵与矩阵的乘法


  1、 运算规则
  设,则A与B的乘积是这样一个矩阵:
  (1) 行数与(左矩阵)A相同,列数与(右矩阵)B相同,即
  (2) C的第行第列的元素由A的第行元素与B的第列元素对应相乘,再取乘积之和.

  典型例题
  例6.5.2 设矩阵
  计算
   的矩阵.设它为
    

    
  想一想:设列矩阵,行矩阵的行数和列数分别是多少呢
  是3×3的矩阵,是1×1的矩阵,即只有一个元素.

  课堂练习
  1、设,求
  2、在第1道练习题中,两个矩阵相乘的顺序是A在左边,B在右边,称为A左乘B或B右乘A.如果交换顺序,让B在左边,A在右边,即A右乘B,运算还能进行吗?请算算试试看.并由此思考:两个矩阵应当满足什么条件,才能够做乘法运算.
  3、设列矩阵,行矩阵,求,比较两个计算结果,能得出什么结论吗?
  4、设三阶方阵,三阶单位阵为,试求,并将计算结果与A比较,看有什么样的结论.

  解:
  第1题
  第2题
  对于
  求是有意义的,而是无意义的.

  结论1 只有在下列情况下,两个矩阵的乘法才有意义,或说乘法运算是可行的:左矩阵的列数=右矩阵的行数.
  第3题
  矩阵,的矩阵.
         
           
    结论2 在矩阵的乘法中,必须注意相乘的顺序.即使在均有意义时,也未必有=成立.可见矩阵乘法不满足交换律.
  第4题
  计算得:
  结论3 方阵A和它同阶的单位阵作乘积,结果仍为A,即
  单位阵在矩阵乘法中的作用相当于数1在我们普通乘法中的作用.

  典型例题
  例6.5.3 设,试计算
   
      
      
    
      
      
    结论4 两个非零矩阵的乘积可以是零矩阵.由此若,不能得出的结论.

  例6.5.4 利用矩阵的乘法,三元线性方程组
  可以写成矩阵的形式
  若记系数、未知量和常数项构成的三个矩阵分别为
  则线性方程组又可以简写为矩阵方程的形式:

  2、 运算性质(假设运算都是可行的)
  (1) 结合律 
  (2) 分配律 (左分配律);
         (右分配律).
  (3) 
   3、 方阵的幂
 
定义:设A是方阵,是一个正整数,规定
显然,记号表示个A的连乘积.

四、矩阵的转置


  1、 定义
 
定义:将矩阵A的行换成同序号的列所得到的新矩阵称为矩阵A的转置矩阵,记作
  例如,矩阵的转置矩阵为
  2、运算性质(假设运算都是可行的)
  (1) 
  (2) 
  (3) 
  (4) 是常数.

  典型例题
  例6.5.5  利用矩阵
  验证运算性质:
   
  而
    
  所以
   

 
定义:如果方阵满足,即,则称A为对称矩阵
  对称矩阵的特点是:它的元素以主对角线为对称轴对应相等.

五、方阵的行列式


  1、定义
 
定义:由方阵A的元素所构成的行列式(各元素的位置不变),称为方阵A的行列式,记作

  2 、运算性质
  (1) (行列式的性质)
  (2) ,特别地:
  (3) 是常数,A的阶数为n)
  思考:设A为阶方阵,那么的行列式与A的行列式之间的关系为什么不是,而是

  不妨自行设计一个二阶方阵,计算一下
  例如,则
  于是,而
  思考:,有几种方法可以求
    方法一:先求矩阵乘法,得到一个二阶方阵,再求其行列式.
    方法二:先分别求行列式,再取它们的乘积.

3、逆矩阵
     逆矩阵: 设A是数域上的一个n阶方阵,若在相同数域上存在另一个n阶矩阵B,使得: AB=BA=E。 则我们称B是A的逆矩阵,而A则被称为可逆矩阵。其用A^(-1)表示。
   逆矩阵的性质为:
    1 矩阵A可逆的充要条件是A的行列式不等于0。
   2 可逆矩阵一定是方阵
3 如果矩阵A是可逆的,A的逆矩阵是唯一的。
4 可逆矩阵也被称为非奇异矩阵、满秩矩阵。
5 两个可逆矩阵的乘积依然可逆。
6 可逆矩阵的转置矩阵也可逆。
7 矩阵可逆当且仅当它是满秩矩阵

 其中的方阵:就是行和列数相等的矩阵。
行列式:其就是求一个矩阵的数值,一般使用上三角或者下三角。其如下:


满秩矩阵:一个判断线性方程是否有解的关键条件,其是其说的就是要线性无关,其可分为列满秩和行满秩。其概念是:
满秩矩阵(non-singular matrix): 设A是n阶矩阵, 若r(A) = n, 则称A为满秩矩阵。但满秩不局限于n阶矩阵。若矩阵秩等于行数,称为行满秩;若矩阵秩等于列数,称为列满秩。既是行满秩又是列满秩则为n阶矩阵即n阶方阵。其中的R(A)指的就是求其秩。其中的行满秩就是行向量线性无关;列满秩就是列向量线性无关;其示意图如下:



   3、 列(行)满秩阵与方程关系:
 
这也是为什么相机标定的时候每个视场只能提前四个有用的点,这时因为第五个点肯定跟前面的四个点其中一个线性相关,意思就是(x,y)(X,Y)这个肯定跟前面的一个点成倍数关系,这样子在求秩的时候,有成倍数的行会被化成零,则这就成了降秩了,就不是满秩了。








2018-03-14 10:34:53 qq_26614295 阅读数 2188

在学习Opencv中的矩阵掩码时候遇到一个关于掩码矩阵表示的紧凑形式的问题,这里本人找了很多关于矩阵卷积的内容,大概说下比较通俗的理解:

1首先上一个叫做卷积核的矩阵

在Opencv官方文档的掩码矩阵中,我们可以看到会出现这样一个矩阵,先说说对掩码矩阵的公式的理解,对于图中一个像素点,视觉效果会收到相邻的像素点的影响,分别是上面的像素点,下面的像素点,左边的像素点和右边的像素点,我们把该像素点乘以5倍,然后分别减去这4个相邻的像素点的值,会发生以下情况,如果改像素点的值比较大,周围的像素点值比较小,那么相减之后的结果必然比该像素点的值大;反之,如果该像素点值比较小,周围像素点值比较大,那么相减之后的结果必然比该像素点值小,从而达到亮的点更加亮,暗的点更加暗,也就是锐化的效果。

那么这个卷积核矩阵就是为了实现这个功能,下面会对这个矩阵的原理做详细分析


2再看看我设定的图像输入矩阵,这里假设图像是单通道


3矩阵掩码操作

操作步骤大概就是先固定住卷积核矩阵不动,也就是上面灰色的矩阵,然后确定图像输入矩阵待求的像素点如(0,0),也就是上面图像输入矩阵的“1”;然后把这个待求点对准卷积核矩阵的中心点,也就是上图卷积核矩阵中的“5”,然后对这两个矩阵求卷积。下面先上图:

看到了吗,“1”和“5”已经对准,下面就是求卷积,求卷积的方法我大概说下我的理解,就是把图像矩阵重叠的部分与卷积核矩阵重叠部分相乘,然后不重叠的部分就是用0乘以对应卷积核矩阵的像素值,然后相加得出的结果,也就是为什么卷积核矩阵中心“5”周围的点要设置为“-1”的原因,就是为了用这些“-1”乘以图像像素点周围4个相邻像素点,不明白的话可以看我下面的计算,当然图中知识以图像矩阵中(0,0)也就是“1”为例,对于图像矩阵,我们必须把其余8个位置的像素值都算出来。

y(0,0) = 0*0 + 0*-1 + 0*0 +0*-1 +1*5 + 2*-1 + 0*0 + 4*-1 + 5*0  

从上面的计算中可以看出,对于图像矩阵中“1”这个点,和他相邻的只有“2”和“4”两个,所以计算中体现出了2*-1和4*-1

剩下图像矩阵中的8个点也是按照这个原理做卷积运算,例如我要求图像矩阵中的“5”点,就把这个“5”点和卷积核矩阵中心的“5”对齐,然后做卷积运算。






2018-05-04 09:49:58 weixin_42026802 阅读数 856

   数学不行就先拿别人的看看吧,比较全:

滤波:FName1为原始数据文件名,FName2为滤波后数据文件名。procedure TForm.Smooth(FName1,FName2:string);var  L1 : longint;  I, Icount, Fp1, Fp2: integer;  pre:array [1..14] of Byte;  Last:array [1..3] of Byte;  X, Y ,Z: buffer;  F:file of byte;begin  assignfile(F,FName1);  reset(F);  datalength:=filesize(F);  closefile(f);  Fp1:=fileopen(Fname1,0);  Icount:=fileread(Fp1,X,14);  L1:=Icount;  assignfile(F,fname2);  rewrite(F);  closefile(F);  Fp2:=fileopen(Fname2,1);  for I:=1 to 14 do    pre[I]:=X[I];  Y[1]:=X[1];  Y[2]:=(X[2]+2*X[1])div 3;  Y[3]:=(x[3]+2*x[2]+3*x[1])div 6;  Y[4]:=(x[4]+2*x[3]+3*x[2]+4*x[1])div 10;  Y[5]:=(x[5]+2*x[4]+3*x[3]+4*x[2]+5*x[1])div 15;  Y[6]:=(x[6]+2*x[5]+3*x[4]+4*x[3]+5*x[2]+6*x[1])div 21;  Y[7]:=(x[7]+2*x[6]+3*x[5]+4*x[4]+5*x[3]+6*x[2]+7*x[1])div 28;  Y[8]:=(x[8]+2*x[7]+3*x[6]+4*x[5]+5*x[4]+6*x[3]+7*x[2]+8*x[1])div 36;  Y[9]:=(x[9]+2*x[8]+3*x[7]+4*x[6]+5*x[5]+6*x[4]+7*x[3]+8*x[2]+7*x[1])div 43;  Y[10]:=(x[10]+2*x[9]+3*x[8]+4*x[7]+5*x[6]+6*x[5]+7*x[4]+8*x[3]+7*x[2]+6*x[1])div 49;  Y[11]:=(x[11]+2*x[10]+3*x[9]+4*x[8]+5*x[7]+6*x[6]+7*x[5]+8*x[4]+7*x[3]+6*x[2]+5*x[1])div 54;  Y[12]:=(x[12]+2*x[11]+3*x[10]+4*x[9]+5*x[8]+6*x[7]+7*x[6]+8*x[5]+7*x[4]+6*x[3]+5*x[2]+4*x[1])div 58;  Y[13]:=(x[13]+2*x[12]+3*x[11]+4*x[10]+5*x[9]+6*x[8]+7*x[7]+8*x[6]+7*x[5]+6*x[4]+5*x[3]+4*x[2]+3*x[1])div 61;  Y[14]:=(x[14]+2*x[13]+3*x[12]+4*x[11]+5*x[10]+6*x[9]+7*x[8]+8*x[7]+7*x[6]+6*x[5]+5*x[4]+4*x[3]+3*x[2]+2*x[1])div 63;  Filewrite(fp2,Y,14);  last[1]:=Y[12];  last[2]:=Y[13];  last[3]:=Y[14];  while L1<dataLength do    begin     Icount:=fileread(Fp1,X,N);     L1:=L1+Icount;     Y[1]:=(x[1]+2*pre[14]+3*pre[13]+4*pre[12]+5*pre[11]+6*pre[10]+7*pre[9]+8*pre[8]+7*pre[7]+6*pre[6]+5*pre[6]+4*pre[4]+3*pre[3]+2*pre[2]+1*pre[1])div 64;     Y[2]:=(x[2]+2*X[1]+3*pre[14]+4*pre[13]+5*pre[12]+6*pre[11]+7*pre[10]+8*pre[9]+7*pre[8]+6*pre[7]+5*pre[6]+4*pre[6]+3*pre[4]+2*pre[3]+1*pre[2])div 64;     Y[3]:=(X[3]+2*x[2]+3*X[1]+4*pre[14]+5*pre[13]+6*pre[12]+7*pre[11]+8*pre[10]+7*pre[9]+6*pre[8]+5*pre[7]+4*pre[6]+3*pre[6]+2*pre[4]+1*pre[3])div 64;     Y[4]:=(X[4]+2*X[3]+3*x[2]+4*X[1]+5*pre[14]+6*pre[13]+7*pre[12]+8*pre[11]+7*pre[10]+6*pre[9]+5*pre[8]+4*pre[7]+3*pre[6]+2*pre[5]+1*pre[4])div 64;     Y[5]:=(X[5]+2*X[4]+3*X[3]+4*x[2]+5*X[1]+6*pre[14]+7*pre[13]+8*pre[12]+7*pre[11]+6*pre[10]+5*pre[9]+4*pre[8]+3*pre[7]+2*pre[6]+1*pre[5])div 64;     Y[6]:=(X[6]+2*X[5]+3*X[4]+4*X[3]+5*x[2]+6*X[1]+7*pre[14]+8*pre[13]+7*pre[12]+6*pre[11]+5*pre[10]+4*pre[9]+3*pre[8]+2*pre[7]+1*pre[6])div 64;     Y[7]:=(X[7]+2*X[6]+3*X[5]+4*X[4]+5*X[3]+6*x[2]+7*X[1]+8*pre[14]+7*pre[13]+6*pre[12]+5*pre[11]+4*pre[10]+3*pre[9]+2*pre[8]+1*pre[7])div 64;     Y[8]:=(X[8]+2*X[7]+3*X[6]+4*X[5]+5*X[4]+6*X[3]+7*x[2]+8*X[1]+7*pre[14]+6*pre[13]+5*pre[12]+4*pre[11]+3*pre[10]+2*pre[9]+1*pre[8])div 64;     Y[9]:=(X[9]+2*X[8]+3*X[7]+4*X[6]+5*X[5]+6*X[4]+7*X[3]+8*x[2]+7*X[1]+6*pre[14]+5*pre[13]+4*pre[12]+3*pre[11]+2*pre[10]+1*pre[9])div 64;     Y[10]:=(X[10]+2*X[9]+3*X[8]+4*X[7]+5*X[6]+6*X[5]+7*X[4]+8*X[3]+7*x[2]+6*X[1]+5*pre[14]+4*pre[13]+3*pre[12]+2*pre[11]+1*pre[10])div 64;     Y[11]:=(X[11]+2*X[10]+3*X[9]+4*X[8]+5*X[7]+6*X[6]+7*X[5]+8*X[4]+7*X[3]+6*x[2]+5*X[1]+4*pre[14]+3*pre[13]+2*pre[12]+1*pre[11])div 64;     Y[12]:=(X[12]+2*X[11]+3*X[10]+4*X[9]+5*X[8]+6*X[7]+7*X[6]+8*X[5]+7*X[4]+6*X[3]+5*x[2]+4*X[1]+3*pre[14]+2*pre[13]+1*pre[12])div 64;     Y[13]:=(X[13]+2*X[12]+3*X[11]+4*X[10]+5*X[9]+6*X[8]+7*X[7]+8*X[6]+7*X[5]+6*X[4]+5*X[3]+4*x[2]+3*X[1]+2*pre[14]+1*pre[13])div 64;     Y[14]:=(X[14]+2*X[13]+3*X[12]+4*X[11]+5*X[10]+6*X[9]+7*X[8]+8*X[7]+7*X[6]+6*X[5]+5*X[4]+4*X[3]+3*x[2]+2*X[1]+1*pre[14])div 64;     Z[1]:=(Y[1]+LAST[3]+LAST[2]+LAST[1]) div 4;     Z[2]:=(Y[2]+Y[1]+LAST[3]+LAST[2]) div 4;     Z[3]:=(Y[3]+Y[2]+Y[1]+LAST[3]) div 4;     for I :=15 to ICount do      Y[I]:=(X[I]+2*X[I-1]+3*X[I-2]+4*X[I-3]+5*X[I-4]+6*X[I-5]+7*X[I-6]+8*X[I-7]+7*X[I-8]+6*X[I-9]+5*X[I-10]+4*X[I-11]+3*x[I-12]+2*X[I-13]+1*X[I-14]) div 64;     for I:= 4 to Icount do      Z[I]:=(Y[I]+Y[I-1]+Y[I-2]+Y[I-3]) div 4;     for I:=1 to 14 do      pre[I]:=X[Icount-I+1];      last[1]:=Y[Icount-2];      Last[2]:=Y[Icount-1];      Last[3]:=Y[Icount];      filewrite(fp2,Z,Icount);    end;      fileclose(fp2);      fileclose(fp1);end;

unit hhx_effectex;interfaceuses hhx_Effects;var mxEmbossColor:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:        (( 0, 0, 0, 0, 0, 0, 0),        ( 0, 0, 0, 0, 0, 0, 0),        ( 0, 0,-1,-1,-1, 0, 0),        ( 0, 0, 0, 1, 0, 0, 0),        ( 0, 0, 1, 1, 1, 0, 0),        ( 0, 0, 0, 0, 0, 0, 0),        ( 0, 0, 0, 0, 0, 0, 0));        Divisor:1;        Bias:0;        FilterName:'彩色浮雕';);var mxEmbossLight:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0,-1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));        Divisor:1;        Bias:192;        FilterName:'高亮浮雕';);var mxEmbossMedium:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));        Divisor:1;        Bias:192;        FilterName:'中值浮雕';);var mxEmbossDark:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));        Divisor:1;        Bias:128;        FilterName:'黑色浮雕';);var mxEdgeEnhance:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;        Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0,-2,16,-2, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));        Divisor:4;        Bias:0;        FilterName:'边缘增强 (线性锐化)';);var mxBlurBartlett:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx7;        Matrix:      (( 1, 2, 3, 4, 3, 2, 1),      ( 2, 4, 6, 8, 6, 4, 2),      ( 3, 6, 9,12, 9, 6, 3),      ( 4, 8,12,16,12, 8, 4),      ( 3, 6, 9,12, 9, 6, 3),      ( 2, 4, 6, 8, 6, 4, 2),      ( 1, 2, 3, 4, 3, 2, 1));        Divisor:256;        Bias:0;        FilterName:'Bartlett模糊 (线性模糊)';);var mxBlurGaussian:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx7;        Matrix:     (( 1, 4, 8, 10, 8, 4, 1),      ( 4,12,25,29,25,12, 4),      ( 8,25,49,58,49,25, 8),      (10,29,58,67,58,29,10),      ( 8,25,49,58,49,25, 8),      ( 4,12,25,29,25,12, 4),      ( 1, 4, 8, 10, 8, 4, 1));      Divisor:999;      Bias:0;      FilterName:'高斯模糊 (线性模糊)';);var mxNegative:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0,-1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:255;      FilterName:'反相 (线性效果)';);var mxAverage:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:9;      Bias:0;      FilterName:'平均值滤波 (线性模糊)';);var mxBlur:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 2, 4, 2, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:16;      Bias:0;      FilterName:'模糊 Blur';);var mxBlurSoftly:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 3, 1, 0, 0),      ( 0, 0, 3,16, 3, 0, 0),      ( 0, 0, 1, 3, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:32;      Bias:0;      FilterName:'轻度模糊 Blur softly';);var mxBlurMore:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx5;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 1, 4, 6, 4, 1, 0),      ( 0, 2, 6, 8, 6, 2, 0),      ( 0, 1, 4, 6, 4, 1, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:64;      Bias:0;      FilterName:'进一步模糊 Blur more';);var mxPrewitt:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 1,-2, 1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'Prewitt边缘 (线性边缘检测)';);var mxTraceContour:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:      (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-6,-2,-6, 0, 0),      ( 0, 0,-1,32,-1, 0, 0),      ( 0, 0,-6,-2,-6, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:240;      FilterName:'轮廓描绘 (线性边缘检测)';);//      FilterName:'Trace contour (Edge detect linear)';);var mxSharpen:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0,-1,16,-1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:8;      Bias:0;      FilterName:'锐化 (线性锐化)';);var mxSharpenMore:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0,-1,12,-1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:4;      Bias:0;      FilterName:'进一步锐化 (线性锐化)';);var mxSharpenLess:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0,-1,24,-1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:16;      Bias:0;      FilterName:'轻度锐化 (线性锐化)';);var mxUnSharpMask:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0,-2,16,-2, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:4;      Bias:0;      FilterName:'非锐化模板 (线性锐化)';);//      FilterName:'非锐化模板Unsharp mask (Sharpen linear)';);

var mxEdgesStrong:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 3, 1, 0, 0),      ( 0, 0, 3,-16,3, 0, 0),      ( 0, 0, 1, 3, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'强化边缘 (线性边缘检测)';);var mxEdgesWeak:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 1, 0, 0, 0),      ( 0, 0, 1,-4, 1, 0, 0),      ( 0, 0, 0, 1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'边缘弱化 (线性边缘检测)';);var mxEtch:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-6, 2,-6, 0, 0),      ( 0, 0,-1,32,-1, 0, 0),      ( 0, 0,-6,-2,-6, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:240;      FilterName:'腐蚀';);//      FilterName:'Etch (Effects linear)';);var mxLaplacianHV:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0,-1, 0, 0, 0),      ( 0, 0,-1, 4,-1, 0, 0),      ( 0, 0, 0,-1, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'Laplacian水平/纵向边缘 (线性边缘检测)';);//      FilterName:'Laplacian horz./vert. (Edge detect linear)';);var mxLaplacianOmni:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0,-1, 8,-1, 0, 0),      ( 0, 0,-1,-1,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'Laplacian所有方向边缘?(线性边缘检测)';);//      FilterName:'Laplacian omnidir? (Edge detect linear)';);var mxSharpenDirectional:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-3,-3,-3, 0, 0),      ( 0, 0, 0,16, 0, 0, 0),      ( 0, 0, 1, 1, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:10;      Bias:0;      FilterName:'有方向的锐化';);//      FilterName:'Sharpen directional (Sharpen linear)';);var mxSobelPass:TGraphicFilter      =(FilterType:ftLinear;MatrixSize:mx3;      Matrix:     (( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 1, 2, 1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0,-1,-2,-1, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0),      ( 0, 0, 0, 0, 0, 0, 0));      Divisor:1;      Bias:0;      FilterName:'Sobel边缘 (线性边缘检测)';);//      FilterName:'Sobel pass (Edge detect linear)';);//******************************************************************************////  Two-Pass filters////******************************************************************************var nmxLaplacianInvert:TMultiPassGraphicFilter=      (FilterType:ftMultiPass;       Filters:(@mxLaplacianOmni,@mxNegative,@mxZero,@mxZero);       Functions:(gaNone,gaNone,gaNone);       FilterName:'Laplacian Negative');implementationend.