图像处理为什么不用dst

2018-10-30 22:19:36 Eastmount 阅读数 8400
  • 再论图像分割

    掌握图像去噪滤波方法及各自特点; 掌握常用的图像边缘检测算子,及Canny算子基本原理 掌握灰度阈值化及大津算法; 掌握常见的图像特征描述方式 了解局部阈值分割、区域生长、分水岭算法及基于轮廓的阈值分割; 通过...

    2537人学习 CSDN就业班
    免费试看

该系列文章是讲解Python OpenCV图像处理知识,前期主要讲解图像入门、OpenCV基础用法,中期讲解图像处理的各种算法,包括图像锐化算子、图像增强技术、图像分割等,后期结合深度学习研究图像识别、图像分类应用。希望文章对您有所帮助,如果有不足之处,还请海涵~

该系列在github所有源代码:https://github.com/eastmountyxz/ImageProcessing-Python
PS:请求帮忙点个Star,哈哈,第一次使用Github,以后会分享更多代码,一起加油。

同时推荐作者的C++图像系列知识:
[数字图像处理] 一.MFC详解显示BMP格式图片
[数字图像处理] 二.MFC单文档分割窗口显示图片
[数字图像处理] 三.MFC实现图像灰度、采样和量化功能详解
[数字图像处理] 四.MFC对话框绘制灰度直方图
[数字图像处理] 五.MFC图像点运算之灰度线性变化、灰度非线性变化、阈值化和均衡化处理详解
[数字图像处理] 六.MFC空间几何变换之图像平移、镜像、旋转、缩放详解
[数字图像处理] 七.MFC图像增强之图像普通平滑、高斯平滑、Laplacian、Sobel、Prewitt锐化详解

前文参考:
[Python图像处理] 一.图像处理基础知识及OpenCV入门函数
[Python图像处理] 二.OpenCV+Numpy库读取与修改像素
[Python图像处理] 三.获取图像属性、兴趣ROI区域及通道处理
[Python图像处理] 四.图像平滑之均值滤波、方框滤波、高斯滤波及中值滤波
[Python图像处理] 五.图像融合、加法运算及图像类型转换
[Python图像处理] 六.图像缩放、图像旋转、图像翻转与图像平移

本篇文章主要讲解Python调用OpenCV实现图像阈值化处理操作,包括二进制阈值化、反二进制阈值化、截断阈值化、反阈值化为0、阈值化为0。全文均是基础知识,希望对您有所帮助。
1.阈值化
2.二进制阈值化
3.反二进制阈值化
4.截断阈值化
5.反阈值化为0
6.阈值化为0

PS:文章参考自己以前系列图像处理文章及OpenCV库函数,同时部分参考网易云视频,推荐大家去学习。同时,本篇文章涉及到《计算机图形学》基础知识,请大家下来补充。

PSS:2019年1~2月作者参加了CSDN2018年博客评选,希望您能投出宝贵的一票。我是59号,Eastmount,杨秀璋。投票地址:https://bss.csdn.net/m/topic/blog_star2018/index

五年来写了314篇博客,12个专栏,是真的热爱分享,热爱CSDN这个平台,也想帮助更多的人,专栏包括Python、数据挖掘、网络爬虫、图像处理、C#、Android等。现在也当了两年老师,更是觉得有义务教好每一个学生,让贵州学子好好写点代码,学点技术,"师者,传到授业解惑也",提前祝大家新年快乐。2019我们携手共进,为爱而生。

一. 阈值化

(注:该部分参考作者的论文《基于苗族服饰的图像锐化和边缘提取技术研究》)

图像的二值化或阈值化(Binarization)旨在提取图像中的目标物体,将背景以及噪声区分开来。通常会设定一个阈值T,通过T将图像的像素划分为两类:大于T的像素群和小于T的像素群。
灰度转换处理后的图像中,每个像素都只有一个灰度值,其大小表示明暗程度。二值化处理可以将图像中的像素划分为两类颜色,常用的二值化算法如公式1所示:

{Y=0gray<TY=255gray>=T\begin{cases} Y=0,gray<T\\ Y=255,gray>=T\\ \end{cases}
当灰度Gray小于阈值T时,其像素设置为0,表示黑色;当灰度Gray大于或等于阈值T时,其Y值为255,表示白色。
Python OpenCV中提供了阈值函数threshold()实现二值化处理,其公式及参数如下图所示:
retval, dst = cv2.threshold(src, thresh, maxval, type)

常用的方法如下表所示,其中函数中的参数Gray表示灰度图,参数127表示对像素值进行分类的阈值,参数255表示像素值高于阈值时应该被赋予的新像素值,最后一个参数对应不同的阈值处理方法。
对应OpenCV提供的五张图如下所示,第一张为原图,后面依次为:二进制阈值化、反二进制阈值化、截断阈值化、反阈值化为0、阈值化为0。
二值化处理广泛应用于各行各业,比如生物学中的细胞图分割、交通领域的车牌设别等。在文化应用领域中,通过二值化处理将所需民族文物图像转换为黑白两色图,从而为后面的图像识别提供更好的支撑作用。下图表示图像经过各种二值化处理算法后的结果,其中“BINARY”是最常见的黑白两色处理。


二. 二进制阈值化

该方法先要选定一个特定的阈值量,比如127。新的阈值产生规则如下:
dst(x,y)={maxValifsrc(x,y)>thresh0otherwise dst(x,y) = \begin{cases} maxVal, if src(x,y)>thresh\\ 0,otherwise\\ \end{cases}
(1) 大于等于127的像素点的灰度值设定为最大值(如8位灰度值最大为255)
(2) 灰度值小于127的像素点的灰度值设定为0
例如,163->255,86->0,102->0,201->255。

关键字为 cv2.THRESH_BINARY,完整代码如下:

#encoding:utf-8
import cv2  
import numpy as np  

#读取图片
src = cv2.imread('test.jpg')

#灰度图像处理
GrayImage = cv2.cvtColor(src,cv2.COLOR_BGR2GRAY)

#二进制阈值化处理
r, b = cv2.threshold(GrayImage, 127, 255, cv2.THRESH_BINARY)
print r

#显示图像
cv2.imshow("src", src)
cv2.imshow("result", b)

#等待显示
cv2.waitKey(0)
cv2.destroyAllWindows()

输出为两个返回值,r为127,b为处理结果(大于127设置为255,小于设置为0)。如下图所示:



三. 反二进制阈值化

该方法与二进制阈值化方法相似,先要选定一个特定的灰度值作为阈值,比如127。新的阈值产生规则如下:
dst(x,y)={0ifsrc(x,y)>threshmaxValotherwise dst(x,y) = \begin{cases} 0, if src(x,y)>thresh\\ maxVal,otherwise\\ \end{cases}
(1) 大于127的像素点的灰度值设定为0(以8位灰度图为例)
(2) 小于该阈值的灰度值设定为255
例如,163->0,86->255,102->255,201->0。

关键字为 cv2.THRESH_BINARY_INV,完整代码如下:

#encoding:utf-8
import cv2  
import numpy as np  

#读取图片
src = cv2.imread('test.jpg')

#灰度图像处理
GrayImage = cv2.cvtColor(src,cv2.COLOR_BGR2GRAY)

#反二进制阈值化处理
r, b = cv2.threshold(GrayImage, 127, 255, cv2.THRESH_BINARY_INV)
print r

#显示图像
cv2.imshow("src", src)
cv2.imshow("result", b)

#等待显示
cv2.waitKey(0)
cv2.destroyAllWindows()

输出结果如下图所示:

该方法得到的结果正好与二进制阈值化方法相反,亮色元素反而处理为黑色,暗色处理为白色。

四. 截断阈值化

该方法需要选定一个阈值,图像中大于该阈值的像素点被设定为该阈值,小于该阈值的保持不变,比如127。新的阈值产生规则如下:
dst(x,y)={thresholdifsrc(x,y)>threshsrc(x,y)otherwise dst(x,y) = \begin{cases} threshold, if src(x,y)>thresh\\ src(x,y),otherwise\\ \end{cases}
(1) 大于等于127的像素点的灰度值设定为该阈值127
(2) 小于该阈值的灰度值不改变
例如,163->127,86->86,102->102,201->127。

关键字为 cv2.THRESH_TRUNC,完整代码如下:

#encoding:utf-8
import cv2  
import numpy as np  

#读取图片
src = cv2.imread('test.jpg')

#灰度图像处理
GrayImage = cv2.cvtColor(src,cv2.COLOR_BGR2GRAY)

#截断阈值化处理
r, b = cv2.threshold(GrayImage, 127, 255, cv2.THRESH_TRUNC)
print r

#显示图像
cv2.imshow("src", src)
cv2.imshow("result", b)

#等待显示
cv2.waitKey(0)
cv2.destroyAllWindows()

输出结果如下图所示::

该处理方法相当于把图像中比较亮(大于127,偏向于白色)的像素值处理为阈值。

五. 反阈值化为0

该方法先选定一个阈值,比如127,接着对图像的灰度值进行如下处理:
dst(x,y)={0ifsrc(x,y)>threshsrc(x,y)otherwise dst(x,y) = \begin{cases} 0, if src(x,y)>thresh\\ src(x,y),otherwise\\ \end{cases}
(1) 大于等于阈值127的像素点变为0
(2) 小于该阈值的像素点值保持不变
例如,163->0,86->86,102->102,201->0。

关键字为 cv2.THRESH_TOZERO_INV,完整代码如下:

#encoding:utf-8
import cv2  
import numpy as np  

#读取图片
src = cv2.imread('test.jpg')

#灰度图像处理
GrayImage = cv2.cvtColor(src,cv2.COLOR_BGR2GRAY)

#反阈值化为0处理
r, b = cv2.threshold(GrayImage, 127, 255, cv2.THRESH_TOZERO_INV)
print r

#显示图像
cv2.imshow("src", src)
cv2.imshow("result", b)

#等待显示
cv2.waitKey(0)
cv2.destroyAllWindows()

输出结果如下图所示:



六. 阈值化为0

该方法先选定一个阈值,比如127,接着对图像的灰度值进行如下处理:
dst(x,y)={src(x,y)ifsrc(x,y)>thresh0otherwise dst(x,y) = \begin{cases} src(x,y), if src(x,y)>thresh\\ 0,otherwise\\ \end{cases}
(1) 大于等于阈值127的像素点,值保持不变
(2) 小于该阈值的像素点值设置为0
例如,163->163,86->0,102->0,201->201。

关键字为 cv2.THRESH_TOZERO,完整代码如下:

#encoding:utf-8
import cv2  
import numpy as np  

#读取图片
src = cv2.imread('test.jpg')

#灰度图像处理
GrayImage = cv2.cvtColor(src,cv2.COLOR_BGR2GRAY)

#阈值化为0处理
r, b = cv2.threshold(GrayImage, 127, 255, cv2.THRESH_TOZERO)
print r

#显示图像
cv2.imshow("src", src)
cv2.imshow("result", b)

#等待显示
cv2.waitKey(0)
cv2.destroyAllWindows()

输出结果如下图所示:

该算法把比较亮的部分不变,比较暗的部分处理为0。

完整五个算法的对比代码如下所示:

#encoding:utf-8
import cv2  
import numpy as np  
import matplotlib.pyplot as plt

#读取图像
img=cv2.imread('test.jpg')
lenna_img = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)
GrayImage=cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)  

#阈值化处理
ret,thresh1=cv2.threshold(GrayImage,127,255,cv2.THRESH_BINARY)  
ret,thresh2=cv2.threshold(GrayImage,127,255,cv2.THRESH_BINARY_INV)  
ret,thresh3=cv2.threshold(GrayImage,127,255,cv2.THRESH_TRUNC)  
ret,thresh4=cv2.threshold(GrayImage,127,255,cv2.THRESH_TOZERO)  
ret,thresh5=cv2.threshold(GrayImage,127,255,cv2.THRESH_TOZERO_INV)

#显示结果
titles = ['Gray Image','BINARY','BINARY_INV','TRUNC','TOZERO','TOZERO_INV']  
images = [GrayImage, thresh1, thresh2, thresh3, thresh4, thresh5]  
for i in xrange(6):  
   plt.subplot(2,3,i+1),plt.imshow(images[i],'gray')  
   plt.title(titles[i])  
   plt.xticks([]),plt.yticks([])  
plt.show()

输出结果如下图所示:

希望文章对大家有所帮助,如果有错误或不足之处,还请海涵。最近经历的事情太多,有喜有悲,关闭了朋友圈,希望通过不断学习和写文章来忘记烦劳,将忧郁转换为动力。哎,总感觉自己在感动这个世界,帮助所有人,而自己却…谁有关心秀璋?晚安。
(By:Eastmount 2018-10-30 晚上10点 https://blog.csdn.net/Eastmount/)

2017-04-15 19:34:01 m0_37264397 阅读数 61094
  • 再论图像分割

    掌握图像去噪滤波方法及各自特点; 掌握常用的图像边缘检测算子,及Canny算子基本原理 掌握灰度阈值化及大津算法; 掌握常见的图像特征描述方式 了解局部阈值分割、区域生长、分水岭算法及基于轮廓的阈值分割; 通过...

    2537人学习 CSDN就业班
    免费试看

    傅立叶变换在图像处理中有非常重要的作用。因为不仅傅立叶分析涉及图像处理很多方面,傅立 叶改进算法,比如离散余弦变换,gabor与小波在图像处理中也有重要的分量。傅立叶变换在图像处理的重要作用:

   1.图像增强与图像去噪

      绝 大部分噪音都是图像的高频分量,通过低通滤波器来滤除高频——噪声; 边缘也是图像的高频分量,可以通过添加高频分量来增强原始图像的边缘;

  2.图像分割之边缘检测

     提 取图像高频分量

  3.图像特征提取:

     形状特征:傅里叶描述子

     纹 理特征:直接通过傅里叶系数来计算纹理特征

     其他特征:将提取的特征值进行傅里叶变 换来使特征具有平移、伸缩、旋转不变性

  4.图像压缩

     可以直接通过傅里叶系数来压缩数据;常用的离散余弦变换是傅立叶变换的实变换;傅立叶变换。

    傅里叶变换是将 时域信号分解为不同频率的正弦信号或余弦函数叠加之和。连续情况下要求原始信号在一个周期内满足绝对可积条件。离散情况下,傅里叶变换一定存在。冈萨雷斯 版<图像处理>里面的解释非常形象:一个恰当的比喻是将傅里叶变换比作一个玻璃棱镜。棱镜是可以将光分解为不同颜色的物理仪器,每个成分的颜 色由波长(或频率)来决定。傅里叶变换可以看作是数学上的棱镜,将函数基于频率分解为不同的成分。当我们考虑光时,讨论它的光谱或频率谱。同样,傅立叶变 换使我们能通过频率成分来分析一个函数。

    傅立叶变换有很多优良的性质。

    如线性, 对称性(可以用在计算信号的傅里叶变换里面);

    时移性:函数在时域中的时移,对 应于其在频率域中附加产生的相移,而幅度频谱则保持不变;

    频移性:函数在时域中乘 以e^jwt,可以使整个频谱搬移w。这个也叫调制定理,通讯里面信号的频分复用需要用到这个特性(将不同的信号调制到不同的频段上同时传输);

   卷积定理:时域卷积等于频域乘积;时域乘积等于频域卷积(附加一个系数)。(图像处理里面 这个是个重点)。

   信号在频率域的表现。

    在频域中,频率越大说明原始信号变化速度越快;频率越小说明原始信号越平缓。当频率为0时,表示直 流信号,没有变化。因此,频率的大小反应了信号的变化快慢。高频分量解释信号的突变部分,而低频分量决定信号的整体形象。

    在图像处理中,频域反应了图像在空域灰度变化剧烈程度,也就是图像灰度的变化速度,也就是图像的梯 度大小。对图像而言,图像的边缘部分是突变部分,变化较快,因此反应在频域上是高频分量;图像的噪声大部分情况下是高频部分;图像平缓变化部分则为低频分 量。也就是说,傅立叶变换提供另外一个角度来观察图像,可以将图像从灰度分布转化到频率分布上来观察图像的特征。书面一点说就是,傅里叶变换提供了一条从 空域到频率自由转换的途径。对图像处理而言,以下概念非常的重要:

    图像高频分量: 图像突变部分;在某些情况下指图像边缘信息,某些情况下指噪声,更多是两者的混合;

    低 频分量:图像变化平缓的部分,也就是图像轮廓信息

    高通滤波器:让图像使低频分量抑 制,高频分量通过

    低通滤波器:与高通相反,让图像使高频分量抑制,低频分量通过

    带通滤波器:使图像在某一部分的频率信息通过,其他过低或过高都抑制

    带阻滤波器,是带通的反。

模板运算与卷积定理

    在时域内做模板运算,实际上就是对图像进行卷积。 模板运算是图像处理一个很重要的处理过程,很多图像处理过程,比如增强/去噪(这两个分不清楚),边缘检测中普遍用到。根据卷积定理,时域卷积等价与频域 乘积。因此,在时域内对图像做模板运算就等效于在频域内对图像做滤波处理。

比如说 一个均值模板,其频域响应为一个低通滤波器;在时域内对图像作均值滤波就等效于在频域内对图像用均值模板的频域响应对图像的频域响应作一个低通滤波。

    图像去噪

    图像去噪 就是压制图像的噪音部分。因此,如果噪音是高频额,从频域的角度来看,就是需要用一个低通滤波器对图像进行处理。通过低通滤波器可以抑制图像的高频分量。 但是这种情况下常常会造成边缘信息的抑制。常见的去噪模板有均值模板,高斯模板等。这两种滤波器都是在局部区域抑制图像的高频分量,模糊图像边缘的同时也 抑制了噪声。还有一种非线性滤波-中值滤波器。中值滤波器对脉冲型噪声有很好的去掉。因为脉冲点都是突变的点,排序以后输出中值,那么那些最大点和最小点 就可以去掉了。中值滤波对高斯噪音效果较差。

    椒盐噪声:对于椒盐采用中值滤波可以 很好的去除。用均值也可以取得一定的效果,但是会引起边缘的模糊。

高斯白噪声:白 噪音在整个频域的都有分布,好像比较困难。

冈萨雷斯版图像处理P185:算术均值 滤波器和几何均值滤波器(尤其是后者)更适合于处理高斯或者均匀的随机噪声。谐波均值滤波器更适合于处理脉冲噪声。

    图像增强

有时候感觉图像增 强与图像去噪是一对矛盾的过程,图像增强经常是需要增强图像的边缘,以获得更好的显示效果,这就需要增加图像的高频分量。而图像去噪是为了消除图像的噪 音,也就是需要抑制高频分量。有时候这两个又是指类似的事情。比如说,消除噪音的同时图像的显示效果显著的提升了,那么,这时候就是同样的意思了。

    常见的图像增强方法有对比度拉伸,直方图均衡化,图像锐化等。前面两个是在空域进行基于像 素点的变换,后面一个是在频域处理。我理解的锐化就是直接在图像上加上图像高通滤波后的分量,也就是图像的边缘效果。对比度拉伸和直方图均衡化都是为了提 高图像的对比度,也就是使图像看起来差异更明显一些,我想,经过这样的处理以后,图像也应该增强了图像的高频分量,使得图像的细节上差异更大。同时也引入 了一些噪音。

    对图像二维傅立叶变换的意义

众所周至,傅立叶变换可以将连续或离散的函数序列从空域映射到频域上,因此,傅立叶变换是信息与信号学中不可获缺的强大工具。但是,由于傅立 叶变换在学习时是以一大堆公式的形式给出的,因此很多人(包括我在内)往往在做了一大堆习题掌握了变换的数学表示却对其变换后的物理意义一无所知,尤其是 自学的时候更是晕头转向。

     这里假设大家对傅立叶变换的数学表示已经很熟悉了,撇开傅立叶变换本身和其在其他领域的应用不谈,只谈图像傅立叶变换前后的对应关系。我们知道傅立叶变换 以前,图像(未压缩的位图)是由对在连续空间(现实空间)上的采样得到一系列点的集合,我们习惯用一个二维矩阵表示空间上各点,则图像可由 z=f(x,y)来表示。由于空间是三维的,图像是二维的,因此空间中物体在另一个维度上的关系就由梯度来表示,这样我们可以通过观察图像得知物体在三维 空间中的对应关系。为什么要提梯度?因为实际上对图像进行二维傅立叶变换得到频谱图,就是图像梯度的分布图,当然频谱图上的各点与图像上各点并不存在一一 对应的关系,即使在不移频的情况下也是没有。傅立叶频谱图上我们看到的明暗不一的亮点,实际上图像上某一点与邻域点差异的强弱,即梯度的大小,也即该点的 频率的大小(可以这么理解,图像中的低频部分指低梯度的点,高频部分相反)。一般来讲,梯度大则该点的亮度强,否则该点亮度弱。这样通过观察傅立叶变换后 的频谱图,也叫功率图(看看频谱图的各点的计算公式就知道为什么叫功率图了:)),我们首先就可以看出,图像的能量分布,如果频谱图中暗的点数更多,那么 实际图像是比较柔和的(因为各点与邻域差异都不大,梯度相对较小),反之,如果频谱图中亮的点数多,那么实际图像一定是尖锐的,边界分明且边界两边像素差 异较大的。对频谱移频到原点以后,可以看出图像的频率分布是以原点为圆心,对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外,还有一个好 处,它可以分离出有周期性规律的干扰信号,比如正玄(sin的正玄,找不到这个字,郁闷)干扰,一副带有正玄干扰,移频到原点的频谱图上可以看出除了中心 以外还存在以某一点为中心,对称分布的亮点集合,这个集合就是干扰噪音产生的,这时可以很直观的通过在该位置放置带阻滤波器消除干扰。



数学公式:

1维的离散序列的DFT变换公式为:

 

2维的离散矩阵的DFT变换公式为:

1.使用模板处理图像相关概念:     

      模板:矩阵方块,其数学含义是一种卷积运算。
      卷积运算:可看作是加权求和的过程,使用到的图像区域中的每个像素分别于卷积核(权矩阵)的每个元素对应相
                乘,所有乘积之和作为区域中心像素的新值。
      卷积核:卷积时使用到的权用一个矩阵表示,该矩阵与使用的图像区域大小相同,其行、列都是奇数,
              是一个权矩阵。
      卷积示例:
              3 * 3 的像素区域R与卷积核G的卷积运算:
              R5(中心像素)=R1G1 + R2G2 + R3G3 + R4G4 + R5G5 + R6G6 + R7G7 + R8G8 + R9G9
            

2.使用模板处理图像的问题:
       边界问题:当处理图像边界像素时,卷积核与图像使用区域不能匹配,卷积核的中心与边界像素点对应,
                 卷积运算将出现问题。
       处理办法:
              A. 忽略边界像素,即处理后的图像将丢掉这些像素。
              B. 保留原边界像素,即copy边界像素到处理后的图像。

3.常用模板:



例子1.:

//【1】以灰度模式读取原始图像并显示
	Mat srcImage = imread("1.jpg", 0);
	if(!srcImage.data ) { printf("读取图片错误,请确定目录下是否有imread函数指定图片存在~! \n"); return false; } 
	imshow("原始图像" , srcImage);   

	//【2】将输入图像延扩到最佳的尺寸,边界用0补充
	int m = getOptimalDFTSize( srcImage.rows );
	int n = getOptimalDFTSize( srcImage.cols ); 
	//将添加的像素初始化为0.
	Mat padded;  
	copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));

	//【3】为傅立叶变换的结果(实部和虚部)分配存储空间。
	//将planes数组组合合并成一个多通道的数组complexI
	Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
	Mat complexI;
	merge(planes, 2, complexI);         

	//【4】进行就地离散傅里叶变换
	dft(complexI, complexI);           

	//【5】将复数转换为幅值,即=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
	split(complexI, planes); // 将多通道数组complexI分离成几个单通道数组,planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
	magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude  
	Mat magnitudeImage = planes[0];

	//【6】进行对数尺度(logarithmic scale)缩放
	magnitudeImage += Scalar::all(1);
	log(magnitudeImage, magnitudeImage);//求自然对数

	//【7】剪切和重分布幅度图象限
	//若有奇数行或奇数列,进行频谱裁剪      
	magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
	//重新排列傅立叶图像中的象限,使得原点位于图像中心  
	int cx = magnitudeImage.cols/2;
	int cy = magnitudeImage.rows/2;
	Mat q0(magnitudeImage, Rect(0, 0, cx, cy));   // ROI区域的左上
	Mat q1(magnitudeImage, Rect(cx, 0, cx, cy));  // ROI区域的右上
	Mat q2(magnitudeImage, Rect(0, cy, cx, cy));  // ROI区域的左下
	Mat q3(magnitudeImage, Rect(cx, cy, cx, cy)); // ROI区域的右下
	//交换象限(左上与右下进行交换)
	Mat tmp;                           
	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);
	//交换象限(右上与左下进行交换)
	q1.copyTo(tmp);                 
	q2.copyTo(q1);
	tmp.copyTo(q2);

	//【8】归一化,用0到1之间的浮点值将矩阵变换为可视的图像格式
	//此句代码的OpenCV2版为:
	//normalize(magnitudeImage, magnitudeImage, 0, 1, CV_MINMAX); 
	//此句代码的OpenCV3版为:
	normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX); 

	//【9】显示效果图
	imshow("频谱幅值", magnitudeImage); 

函数解读:

C++: intgetOptimalDFTSize(int vecsize)

源码解读;

copyMakeBorder

C++: void copyMakeBorder(InputArraysrc, OutputArray dst, int top, int bottom,                                                        int left,int right, intborderType, const Scalar& value=Scalar())

src: 源图像

dst: 目标图像,和源图像有相同的类型,dst.cols=src.cols+left+right; dst.rows=src.rows+dst.top+dst.bottom

top:

bottom:

left:

right: 以上四个参数指定了在src图像周围附加的像素个数。

borderType: 边框类型

value: 当borderType==BORDER_CONSTANT时需要指定该值。


例子2.

  1. int cv::getOptimalDFTSizeint size0 )  
  2. {  
  3.    int a = 0, b = sizeof(optimalDFTSizeTab)/sizeof(optimalDFTSizeTab[0]) -1;  
  4.    if( (unsigned)size0 >= (unsigned)optimalDFTSizeTab[b] )  
  5.        return -1;  
  6.    
  7.    while( a < b )//二分查找合适的size  
  8.     {  
  9.        int c = (a + b) >> 1;  
  10.        if( size0 <= optimalDFTSizeTab[c] )  
  11.            b = c;  
  12.        else  
  13.            a = c+1;  
  14.     }  
  15.    
  16.     returnoptimalDFTSizeTab[b];  
  17. }</span>  

optimalDFTSizeTab定义在namespace cv中,里边的数值为2^x*3^y*5^z

static const int optimalDFTSizeTab[] = {1,2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16,…,                                                                                       2123366400, 2125764000};

 

 


示例代码:

  1. <span style="font-size:18px;">#include <opencv2/core/core.hpp>  
  2. #include<opencv2/highgui/highgui.hpp>  
  3. #include<opencv2/imgproc/imgproc.hpp>  
  4. #include <iostream>  
  5.    
  6. using namespace cv;  
  7. using namespace std;  
  8.    
  9. int main(){  
  10.        Matsrc = imread("fruits.jpg");  
  11.        if(src.empty())  
  12.        {  
  13.               return-1;  
  14.        }  
  15.    
  16.        Matsrc_gray;  
  17.        cvtColor(src,src_gray,CV_RGB2GRAY);//灰度图像做傅里叶变换  
  18.    
  19.        intm = getOptimalDFTSize(src_gray.rows);//2,3,5的倍数有更高效率的傅里叶转换  
  20.        intn = getOptimalDFTSize(src_gray.cols);  
  21.    
  22.        Matdst;  
  23.        ///把灰度图像放在左上角,在右边和下边扩展图像,扩展部分填充为0;  
  24.        copyMakeBorder(src_gray,dst,0,m-src_gray.rows,0,n-src_gray.cols,BORDER_CONSTANT,Scalar::all(0));  
  25.        cout<<dst.size()<<endl;  
  26.    
  27.        //新建一个两页的array,其中第一页用扩展后的图像初始化,第二页初始化为0  
  28.        Matplanes[] = {Mat_<float>(dst), Mat::zeros(dst.size(), CV_32F)};  
  29.        Mat  completeI;  
  30.        merge(planes,2,completeI);//把两页合成一个2通道的mat  
  31.    
  32.        //对上边合成的mat进行傅里叶变换,支持原地操作,傅里叶变换结果为复数。通道1存的是实部,通道2存的是虚部。  
  33.        dft(completeI,completeI);  
  34.    
  35.        split(completeI,planes);//把变换后的结果分割到各个数组的两页中,方便后续操作  
  36.        magnitude(planes[0],planes[1],planes[0]);//求傅里叶变换各频率的幅值,幅值放在第一页中。  
  37.    
  38.        MatmagI = planes[0];  
  39.        //傅立叶变换的幅度值范围大到不适合在屏幕上显示。高值在屏幕上显示为白点,  
  40.        //而低值为黑点,高低值的变化无法有效分辨。为了在屏幕上凸显出高低变化的连续性,我们可以用对数尺度来替换线性尺度:  
  41.        magI+= 1;  
  42.        log(magI,magI);//取对数  
  43.        magI= magI(Rect(0,0,src_gray.cols,src_gray.rows));//前边对原始图像进行了扩展,这里把对原始图像傅里叶变换取出,剔除扩展部分。  
  44.    
  45.    
  46.        //这一步的目的仍然是为了显示。 现在我们有了重分布后的幅度图,  
  47.        //但是幅度值仍然超过可显示范围[0,1] 。我们使用 normalize() 函数将幅度归一化到可显示范围。  
  48.        normalize(magI,magI,0,1,CV_MINMAX);//傅里叶图像进行归一化。  
  49.    
  50.    
  51.        //重新分配象限,使(0,0)移动到图像中心,  
  52.        //在《数字图像处理》中,傅里叶变换之前要对源图像乘以(-1)^(x+y)进行中心化。  
  53.        //这是是对傅里叶变换结果进行中心化  
  54.        intcx = magI.cols/2;  
  55.        intcy = magI.rows/2;  
  56.    
  57.        Mattmp;  
  58.        Matq0(magI,Rect(0,0,cx,cy));  
  59.        Matq1(magI,Rect(cx,0,cx,cy));  
  60.        Matq2(magI,Rect(0,cy,cx,cy));  
  61.        Matq3(magI,Rect(cx,cy,cx,cy));  
  62.    
  63.         
  64.        q0.copyTo(tmp);  
  65.        q3.copyTo(q0);  
  66.        tmp.copyTo(q3);  
  67.    
  68.        q1.copyTo(tmp);  
  69.        q2.copyTo(q1);  
  70.        tmp.copyTo(q2);  
  71.    
  72.         
  73.    
  74.        namedWindow("InputImage");  
  75.        imshow("InputImage",src);  
  76.    
  77.        namedWindow("SpectrumImage");  
  78.        imshow("SpectrumImage",magI);  
  79.    
  80.        waitKey();  
  81.        return0;  
  82. }</span>  


图像卷积原理:

图像处理-线性滤波-1 基础(相关算子、卷积算子、边缘效应)

这里讨论利用输入图像中像素的小邻域来产生输出图像的方法,在信号处理中这种方法称为滤波(filtering)。其中,最常用的是线性滤波:输出像素是输入邻域像素的加权和。

 

1.相关算子(Correlation Operator)

       定义:image image ,其中h称为相关核(Kernel).

        

  步骤:

        1)滑动核,使其中心位于输入图像g的(i,j)像素上

        2)利用上式求和,得到输出图像的(i,j)像素值

        3)充分上面操纵,直到求出输出图像的所有像素值

 

  例:

A = [17  24      15            h = [8     6
     23      14  16                 3     7
         13  20  22                 4     2]
     10  12  19  21             
     11  18  25     9]

计算输出图像的(2,4)元素=image

image

Matlab 函数:imfilter(A,h)

 

2.卷积算子(Convolution)

定义:image image ,其中

   步骤:

        1)将围绕中心旋转180度

        2)滑动核,使其中心位于输入图像g的(i,j)像素上

        3)利用上式求和,得到输出图像的(i,j)像素值

        4)充分上面操纵,直到求出输出图像的所有像素值

       例:计算输出图像的(2,4)元素=image

       image

Matlab 函数:Matlab 函数:imfilter(A,h,'conv')% imfilter默认是相关算子,因此当进行卷积计算时需要传入参数'conv'

3.边缘效应

当对图像边缘的进行滤波时,核的一部分会位于图像边缘外面。

image

常用的策略包括:

1)使用常数填充:imfilter默认用0填充,这会造成处理后的图像边缘是黑色的。

2)复制边缘像素:I3 = imfilter(I,h,'replicate');

image

   

4.常用滤波

fspecial函数可以生成几种定义好的滤波器的相关算子的核。

例:unsharp masking 滤波

1
2
3
4
5
I = imread('moon.tif');
h = fspecial('unsharp');
I2 = imfilter(I,h);
imshow(I), title('Original Image')
figure, imshow(I2), title('Filtered Image')
 
 

图像处理-线性滤波-2 图像微分(1、2阶导数和拉普拉斯算子)

更复杂些的滤波算子一般是先利用高斯滤波来平滑,然后计算其1阶和2阶微分。由于它们滤除高频和低频,因此称为带通滤波器(band-pass filters)。

在介绍具体的带通滤波器前,先介绍必备的图像微分知识。

1 一阶导数

连续函数,其微分可表达为image ,或image                         (1.1)

对于离散情况(图像),其导数必须用差分方差来近似,有

                                   image,前向差分 forward differencing                  (1.2)

                                   image ,中心差分 central differencing                     (1.3)

1)前向差分的Matlab实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
function dimg = mipforwarddiff(img,direction)
% MIPFORWARDDIFF     Finite difference calculations 
%
  DIMG = MIPFORWARDDIFF(IMG,DIRECTION)
%
 Calculates the forward-difference for a given direction
 IMG       : input image
 DIRECTION : 'dx' or 'dy'
 DIMG      : resultant image
%
  See also MIPCENTRALDIFF MIPBACKWARDDIFF MIPSECONDDERIV
  MIPSECONDPARTIALDERIV
  
  Omer Demirkaya, Musa Asyali, Prasana Shaoo, ... 9/1/06
  Medical Image Processing Toolbox
  
imgPad = padarray(img,[1 1],'symmetric','both');%将原图像的边界扩展
[row,col] = size(imgPad);
dimg = zeros(row,col);
switch (direction)   
case 'dx',
   dimg(:,1:col-1) = imgPad(:,2:col)-imgPad(:,1:col-1);%x方向差分计算,
case 'dy',
   dimg(1:row-1,:) = imgPad(2:row,:)-imgPad(1:row-1,:); 
otherwise, disp('Direction is unknown');
end;
dimg = dimg(2:end-1,2:end-1);

2)中心差分的Matlab实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
function dimg = mipcentraldiff(img,direction)
% MIPCENTRALDIFF     Finite difference calculations 
%
  DIMG = MIPCENTRALDIFF(IMG,DIRECTION)
%
 Calculates the central-difference for a given direction
 IMG       : input image
 DIRECTION : 'dx' or 'dy'
 DIMG      : resultant image
%
  See also MIPFORWARDDIFF MIPBACKWARDDIFF MIPSECONDDERIV
  MIPSECONDPARTIALDERIV
  
  Omer Demirkaya, Musa Asyali, Prasana Shaoo, ... 9/1/06
  Medical Image Processing Toolbox
  
img = padarray(img,[1 1],'symmetric','both');
[row,col] = size(img);
dimg = zeros(row,col);
switch (direction)
    case 'dx',
        dimg(:,2:col-1) = (img(:,3:col)-img(:,1:col-2))/2;
    case 'dy',
        dimg(2:row-1,:) = (img(3:row,:)-img(1:row-2,:))/2;
    otherwise,
        disp('Direction is unknown');
end
dimg = dimg(2:end-1,2:end-1);
1   

实例:技术图像x方向导数

1
2
I = imread('coins.png'); figure; imshow(I);
Id = mipforwarddiff(I,'dx'); figure, imshow(Id);

      image image

    原图像                           x方向1阶导数

 

2 图像梯度(Image Gradient)

图像I的梯度定义为image  ,其幅值为image 。出于计算性能考虑,幅值也可用image 来近似。

Matlab函数

1)gradient:梯度计算

2)quiver:以箭头形状绘制梯度。注意放大下面最右侧图可看到箭头,由于这里计算横竖两个方向的梯度,因此箭头方向都是水平或垂直的。

实例:仍采用上面的原始图像

1
2
3
4
5
I = double(imread('coins.png'));
[dx,dy]=gradient(I);
magnitudeI=sqrt(dx.^2+dy.^2);
figure;imagesc(magnitudeI);colormap(gray);%梯度幅值
hold on;quiver(dx,dy);%叠加梯度方向

        image image

                         梯度幅值          梯度幅值+梯度方向

 

3 二阶导数

对于一维函数,其二阶导数image ,即image 。它的差分函数为

                                 image                  (3.1)

 

3.1 普拉斯算子(laplacian operator)

3.1.2 概念

拉普拉斯算子是n维欧式空间的一个二阶微分算子。它定义为两个梯度向量算子的内积

                          image       (3.2)

其在二维空间上的公式为:    image                (3.3)

 

对于1维离散情况,其二阶导数变为二阶差分

1)首先,其一阶差分为image

2)因此,二阶差分为

          image

3)因此,1维拉普拉斯运算可以通过1维卷积核image 实现

 

对于2维离散情况(图像),拉普拉斯算子是2个维上二阶差分的和(见式3.3),其公式为:

image   (3.4)

上式对应的卷积核为

                       image

常用的拉普拉斯核有:

                      image

3.1.2 应用

拉普拉斯算子会突出像素值快速变化的区域,因此常用于边缘检测。

 

 

Matlab里有两个函数

1)del2

计算公式:image image  

2)fspecial:图像处理中一般利用Matlab函数fspecial

h = fspecial('laplacian', alpha) returns a 3-by-3 filter approximating the shape of the two-dimensional Laplacian operator.
The parameter alpha controls the shape of the Laplacian and must be in the range 0.0 to 1.0. The default value for alpha is 0.2.

 

3.1.3 资源

http://fourier.eng.hmc.edu/e161/lectures/gradient/node8.html (非常清晰的Laplacian Operator介绍,本文的主要参考)

http://homepages.inf.ed.ac.uk/rbf/HIPR2/log.htm

 

 
 
 
 

sift算法

 

尺度不变特征转换(Scale-invariant feature transform 或 SIFT)是一种电脑视觉的算法用来侦测与描述影像中的局部性特征,它在空间尺度中寻找极值点,并提取出其位置、尺度、旋转不变量,此算法由 David Lowe 在1999年所发表,2004年完善总结。

Sift算法就是用不同尺度(标准差)的高斯函数对图像进行平滑,然后比较平滑后图像的差别,
差别大的像素就是特征明显的点。

sift可以同时处理亮度,平移,旋转,尺度的变化,利用特征点来提取特征描述符,最后在特征描述符之间寻找匹配


五个步骤

1构建尺度空间,检测极值点,获得尺度不变性

2特征点过滤并进行经确定位,剔除不稳定的特征点

3 在特征点处提取特征描述符,为特征点分配方向直

4声称特征描述子,利用特征描述符寻找匹配点

5计算变换参数

当2幅图像的sift特征向量生成以后,下一步就可以采用关键点特征向量的欧式距离来作为2幅图像中关键点的相似性判定量度


尺度空间:

尺度就是受delta这个参数控制的表示

而不同的L(x,y,delta)就构成了尺度空间,实际上具体计算的时候即使连续的高斯函数,都要被离散为矩阵来和数字图像进行卷积操作

L(x,y,delta)=G(x,y,e)*i(x,y)

尺度空间=原始图像(卷积)一个可变尺度的2维高斯函数G(x,y,e)


G(x,y,e) = [1/2*pi*e^2] * exp[ -(x^2 + y^2)/2e^2] 


为了更有效的在尺度空间检测到稳定的关键点,提出了高斯差分尺度空间,利用不同尺度的高斯差分核与原始图像i(x,y)卷积生成

D(x,y,e)=(G(x,y,ke)-G(x,y,e))*i(x,y)

=L(x,y,ke)-L(x,y,e)

(为避免遍历每个像素点)


高斯卷积:

在组建一组尺度空间后,再组建下一组尺度空间,对上一组尺度空间的最后一幅图像进行二分之一采样,得到下一组尺度空间的第一幅图像,然后进行像建立第一组尺度空间那样的操作,得到第二组尺度空间,公式定义为
         L(x,y,e) = G(x,y,e)*I(x,y)

    图像金字塔的构建:图像金字塔共O组,每组有S层,下一组的图像由上一组图像降采样得到、

高斯差分

    在尺度空间建立完毕后,为了能够找到稳定的关键点,采用高斯差分的方法来检测那些在局部位置的极值点,即采用俩个相邻的尺度中的图像相减,即公式定义为:
        D(x,y,e) = ((G(x,y,ke) - G(x,y,e)) * I(x,y) 
                 = L(x,y,ke) - L(x,y,e)
 咱们再来具体阐述下构造D(x,y,e)的详细步骤:
    1、首先采用不同尺度因子的高斯核对图像进行卷积以得到图像的不同尺度空间,将这一组图像作为金子塔图像的第一层。
    2、接着对第一层图像中的2倍尺度图像(相对于该层第一幅图像的2倍尺度)以2倍像素距离进行下采样来得到金子塔图像的第二层中的第一幅图像,对该图像采用不同尺度因子的高斯核进行卷积,以获得金字塔图像中第二层的一组图像。
    3、再以金字塔图像中第二层中的2倍尺度图像(相对于该层第一幅图像的2倍尺度)以2倍像素距离进行下采样来得到金字塔图像的第三层中的第一幅图像,对该图像采用不同尺度因子的高斯核进行卷积,以获得金字塔图像中第三层的一组图像。这样依次类推,从而获得了金字塔图像的每一层中的一组图像,
 4、对上图得到的每一层相邻的高斯图像相减,就得到了高斯差分图像,如下述第一幅图所示。下述第二幅图中的右列显示了将每组中相邻图像相减所生成的高斯差分图像的结果,限于篇幅,图中只给出了第一层和第二层高斯差分图像的计算
sift算法
 

 

图像处理之卷积概念

 

我们来看一下一维卷积的概念.
连续空间的卷积定义是 f(x)与g(x)的卷积是 f(t-x)g(x) 在t从负无穷到正无穷的积分值.t-x要在f(x)定义域内,所以看上去很大的积分实际上还是在一定范围的.
实际的过程就是f(x) 先做一个Y轴的反转,然后再沿X轴平移t就是f(t-x),然后再把g(x)拿来,两者乘积的值再积分.想象一下如果g(x)或者f(x)是个单位的阶越函数. 那么就是f(t-x)与g(x)相交部分的面积.这就是卷积了.
把积分符号换成求和就是离散空间的卷积定义了.

 

么在图像中卷积卷积地是什么意思呢,就是图像f(x),模板g(x),然后将模版g(x)在模版中移动,每到一个位置,就把f(x)与g(x)的定义域相交的元素进行乘积并且求和,得出新的图像一点,就是被卷积后的图像. 模版又称为卷积核.卷积核做一个矩阵的形状.


卷积定义上是线性系统分析经常用到的.线性系统就是一个系统的输入和输出的关系是线性关系.就是说整个系统可以分解成N多的无关独立变化,整个系统就是这些变化的累加.
如 x1->y1, x2->y2; 那么A*x1 + B*x2 -> A*y1 + B*y2 这就是线性系统. 表示一个线性系统可以用积分的形式 如 Y = Sf(t,x)g(x)dt S表示积分符号,就是f(t,x)表示的是A B之类的线性系数.
看上去很像卷积呀,,对如果f(t,x) = F(t-x) 不就是了吗.从f(t,x)变成F(t-x)实际上是说明f(t,x)是个线性移不变,就是说 变量的差不变化的时候,那么函数的值不变化. 实际上说明一个事情就是说线性移不变系统的输出可以通过输入和表示系统线性特征的函数卷积得到.


http://dept.wyu.edu.cn/dip/DIPPPT2005/����������ϵͳ.ppt


谈起卷积分当然要先说说冲击函数—-这个倒立的小蝌蚪,卷积其实就是为它诞生的。”冲击函数”是狄拉克为了解决一些瞬间作用的物理现象而提出的符号。
古人曰:”说一堆大道理不如举一个好例子”,冲量这一物理现象很能说明”冲击函数”。在t时间内对一物体作用F的力,我们可以让作用时间t很小,作用力F很大,但让Ft的乘积不变,即冲量不变。于是在用t做横坐标、F做纵坐标的坐标系中,就如同一个面积不变的长方形,底边被挤的窄窄的,高度被挤的高高的,在数学中它可以被挤到无限高,但即使它无限瘦、无限高、但它仍然保持面积不变(它没有被挤没!),为了证实它的存在,可以对它进行积分,积分就是求面积嘛!于是”卷积” 这个数学怪物就这样诞生了。说它是数学怪物是因为追求完美的数学家始终在头脑中转不过来弯,一个能瘦到无限小的家伙,竟能在积分中占有一席之地,必须将这个细高挑清除数学界。但物理学家、工程师们确非常喜欢它,因为它解决了很多当时数学家解决不了的实际问题。最终追求完美的数学家终于想通了,数学是来源于实际的,并最终服务于实际才是真。于是,他们为它量身定做了一套运作规律。于是,妈呀!你我都感觉眩晕的卷积分产生了。
例子:
有一个七品县令,喜欢用打板子来惩戒那些市井无赖,而且有个惯例:如果没犯大罪,只打一板,释放回家,以示爱民如子。
有一个无赖,想出人头地却没啥指望,心想:既然扬不了善名,出恶名也成啊。怎么出恶名?炒作呗!怎么炒作?找名人呀!他自然想到了他的行政长官——县令。
无赖于是光天化日之下,站在县衙门前撒了一泡尿,后果是可想而知地,自然被请进大堂挨了一板子,然后昂首挺胸回家,躺了一天,嘿!身上啥事也没有!第二天如法炮制,全然不顾行政长管的仁慈和衙门的体面,第三天、第四天……每天去县衙门领一个板子回来,还喜气洋洋地,坚持一个月之久!这无赖的名气已经和衙门口的臭气一样,传遍八方了!
县令大人噤着鼻子,呆呆地盯着案子上的惊堂木,拧着眉头思考一个问题:这三十个大板子怎么不好使捏?……想当初,本老爷金榜题名时,数学可是得了满分,今天好歹要解决这个问题:
——人(系统!)挨板子(脉冲!)以后,会有什么表现(输出!)?
——费话,疼呗!
——我问的是:会有什么表现?
——看疼到啥程度。像这无赖的体格,每天挨一个板子啥事都不会有,连哼一下都不可能,你也看到他那得意洋洋的嘴脸了(输出0);如果一次连揍他十个板子,他可能会皱皱眉头,咬咬牙,硬挺着不哼
(输出1);揍到二十个板子,他会疼得脸部扭曲,象猪似地哼哼(输出3);揍到三十个板子,他可能会象驴似地嚎叫,一把鼻涕一把泪地求你饶他一命(输出5);揍到四十个板子,他会大小便失禁,勉
强哼出声来(输出1);揍到五十个板子,他连哼一下都不可能(输出0)——死啦!
县令铺开坐标纸,以打板子的个数作为X轴,以哼哼的程度(输出)为Y轴,绘制了一条曲线:
——呜呼呀!这曲线象一座高山,弄不懂弄不懂。为啥那个无赖连挨了三十天大板却不喊绕命呀?
—— 呵呵,你打一次的时间间隔(Δτ=24小时)太长了,所以那个无赖承受的痛苦程度一天一利索,没有叠加,始终是一个常数;如果缩短打板子的时间间隔(建议 Δτ=0.5秒),那他的痛苦程度可就迅速叠加了;等到这无赖挨三十个大板(t=30)时,痛苦程度达到了他能喊叫的极限,会收到最好的惩戒效果,再多打就显示不出您的仁慈了。
——还是不太明白,时间间隔小,为什么痛苦程度会叠加呢?
——这与人(线性时不变系统)对板子(脉冲、输入、激励)的响应有关。什么是响应?人挨一个板子后,疼痛的感觉会在一天(假设的,因人而异)内慢慢消失(衰减),而不可能突然消失。这样一来,只要打板子的时间间隔很小,每一个板子引起的疼痛都来不及完全衰减,都会对最终的痛苦程度有不同的贡献:
t个大板子造成的痛苦程度=Σ(第τ个大板子引起的痛苦*衰减系数)
[衰减系数是(t-τ)的函数,仔细品味]
数学表达为:y(t)=∫T(τ)H(t-τ)
——拿人的痛苦来说卷积的事,太残忍了。除了人以外,其他事物也符合这条规律吗?
——呵呵,县令大人毕竟仁慈。其实除人之外,很多事情也遵循此道。好好想一想,铁丝为什么弯曲一次不折,快速弯曲多次却会轻易折掉呢?
——恩,一时还弄不清,容本官慢慢想来——但有一点是明确地——来人啊,将撒尿的那个无赖抓来,狠打40大板!
卷积及拉普拉斯变换的通俗解释–对于我这类没学过信号系统的人来说太需要了
卷积(convolution, 另一个通用名称是德文的Faltung)的名称由来,是在于当初定义它时,定义成 integ(f1(v)*f2(t-v))dv,积分区间在0到t之间。举个简单的例子,大家可以看到,为什么叫”卷积”了。比方说在(0,100)间积分,用简单的辛普生积分公式,积分区间分成100等分,那么看到的是f1(0)和f2(100)相乘,f1(1)和f2(99)相乘,f1(2)和f2 (98)相乘,……… 等等等等,就象是在坐标轴上回卷一样。所以人们就叫它”回卷积分”,或者”卷积”了。
为了理解”卷积”的物理意义,不妨将那个问题”相当于它的时域的信号与系统的单位脉冲响应的卷积”略作变化。这个变化纯粹是为了方便表达和理解,不影响任何其它方面。将这个问题表述成这样一个问题:一个信号通过一个系统,系统的响应是频率响应或波谱响应,且看如何理解卷积的物理意义。
假设信号函数为f, 响应函数为g。f不仅是时间的函数(信号时有时无),还是频率的函数(就算在某一固定时刻,还有的地方大有的地方小);g也是时间的函数(有时候有反应,有时候没反应),同时也是频率的函数(不同的波长其响应程度不一样)。那我们要看某一时刻 t 的响应信号,该怎么办呢?
这就需要卷积了。
要看某一时刻 t 的响应信号,自然是看下面两点:
1。你信号来的时候正赶上人家”系统”的响应时间段吗?
2。就算赶上系统响应时间段,响应有多少?
响 应不响应主要是看 f 和 g 两个函数有没有交叠;响应强度的大小不仅取决于所给的信号的强弱,还取决于在某频率处对单位强度响应率。响应强度是信号强弱和对单位强度信号响应率的乘积。”交叠”体现在f(t1)和g(t-t1)上,g之所以是”(t-t1)”就是看两个函数错开多少。
由于 f 和 g 两个函数都有一定的带宽分布(假若不用开头提到的”表述变化”就是都有一定的时间带宽分布),这个信号响应是在一定”范围”内广泛响应的。算总的响应信号,当然要把所有可能的响应加起来,实际上就是对所有可能t1积分了。积分范围虽然一般在负无穷到正无穷之间;但在没有信号或者没有响应的地方,积也是白积,结果是0,所以往往积分范围可以缩减。
这就是卷积及其物理意义啊。并成一句话来说,就是看一个时有时无(当然作为特例也可以永恒存在)的信号,跟一个响应函数在某一时刻有多大交叠。
*********拉普拉斯*********
拉普拉斯(1729-1827) 是法国数学家,天文学家,物理学家。他提出拉普拉斯变换(Laplace Transform) 的目的是想要解决他当时研究的牛顿引力场和太阳系的问题中涉及的积分微分方程。
拉普拉斯变换其实是一个数学上的简便算法;想要了解其”物理”意义 — 如果有的话 — 请看我举这样一个例子:
问题:请计算十万乘以一千万。
对于没学过指数的人,就只会直接相乘;对于学过指数的人,知道不过是把乘数和被乘数表达成指数形式后,两个指数相加就行了;如果要问究竟是多少,把指数转回来就是。
“拉 普拉斯变换” 就相当于上述例子中把数转换成”指数” 的过程;进行了拉普拉斯变换之后,复杂的微分方程(对应于上例中”复杂”的乘法) 就变成了简单的代数方程,就象上例中”复杂”的乘法变成了简单的加减法。再把简单的代数方程的解反变换回去(就象把指数重新转换会一般的数一样),就解决了原来那个复杂的微分方程。
所以要说拉普拉斯变换真有” 物理意义”的话,其物理意义就相当于人们把一般的有理数用指数形式表达一样。
另外说两句题外话:
1 。拉普拉斯变换之所以现在在电路中广泛应有,根本原因是电路中也广泛涉及了微分方程。
2。拉普拉斯变换与Z变换当然有紧密联系;其本质区别在于拉氏变换处理的是时间上连续的问题,Z变换处理的是时间上分立的问题。
Signals, Linear Systems, and Convolution
Download from here
 
我们都知道卷积公式,但是它有什么物理意义呢?平时我们用卷积做过很多事情,信号处理时,输出函数是输入函数和系统函数的卷积;在图像处理时,两组幅分辨率不同的图卷积之后得到的互相平滑的图像可以方便处理。卷积甚至可以用在考试作弊中,为了让照片同时像两个人,只要把两人的图像卷积处理即可,这就是一种平滑的过程,可是我们怎么才能真正把公式和实际建立起一种联系呢?生活中就有实例:
     比如说你的老板命令你干活,你却到楼下打台球去了,后来被老板发现,他非常气愤,扇了你一巴掌(注意,这就是输入信号,脉冲),于是你的脸上会渐渐地(贱贱地)鼓起来一个包,你的脸就是一个系统,而鼓起来的包就是你的脸对巴掌的响应。
      好,这样就和信号系统建立起来意义对应的联系。下面还需要一些假设来保证论证的严谨:假定你的脸是线性时不变系统,也就是说,无论什么时候老板打你一巴掌,打在你脸的同一位置(这似乎要求你的脸足够光滑,如果你说你长了很多青春痘,甚至整个脸皮处处连续处处不可导,那难度太大了,我就无话可说了),你的脸上总是会在相同的时间间隔内鼓起来一个相同高度的包来,并且假定以鼓起来的包的大小作为系统输出。好了,那么,下面可以进入核心内容——卷积了!
      如果你每天都到楼下去打台球,那么老板每天都要扇你一巴掌,不过当老板打你一巴掌后,你5分钟就消肿了,所以时间长了,你甚至就适应这种生活了……如果有一天,老板忍无可忍,以0.5秒的间隔开始不间断的扇你的过程,这样问题就来了:第一次扇你鼓起来的包还没消肿,第二个巴掌就来了,你脸上的包就可能鼓起来两倍高,老板不断扇你,脉冲不断作用在你脸上,效果不断叠加了,这样这些效果就可以求和了,结果就是你脸上的包的高度岁时间变化的一个函数了(注意理解)!
      如果老板再狠一点,频率越来越高,以至于你都辨别不清时间间隔了,那么,求和就变成积分了。可以这样理解,在这个过程中的某一固定的时刻,你的脸上的包的鼓起程度和什么有关呢?和之前每次打你都有关!但是各次的贡献是不一样的,越早打的巴掌,贡献越小,这就是说,某一时刻的输出是之前很多次输入乘以各自的衰减系数之后的叠加而形成某一点的输出,然后再把不同时刻的输出点放在一起,形成一个函数,这就是卷积。卷积之后的函数就是你脸上的包的大小随时间变化的函数。本来你的包几分钟就可以消肿,可是如果连续打,几个小时也消不了肿了,这难道不是一种平滑过程么?反映到公式上,f(a)就是第a个巴掌,g(x-a)就是第a个巴掌在x时刻的作用程度,乘起来再叠加就ok了,这就是卷积!
     最后提醒各位,请勿亲身尝试……

卷积的物理意义?
在信号与系统中,两个函数所要表达的物理含义是什么?例如,一个系统,其单位冲激响应为h(t),当输入信号为f(t)时,该系统的输出为y(t)。为什么y(t)是f(t)和h(t)的卷积?(从数学推导我明白,但其物理意义不明白。)y(t)是f(t)和h(t)的卷积表达了一个什么意思?

卷积(convolution, 另一个通用名称是德文的Faltung)的名称由来,是在于当初定义它时,定义成 integ(f1(v)*f2(t-v))dv,积分区间在0到t之间。举个简单的例子,大家可以看到,为什么叫“卷积”了。比方说在(0,100)间积分,用简单的辛普生积分公式,积分区间分成100等分,那么看到的是f1(0)和f2(100)相乘,f1(1)和f2(99)相乘,f1(2)和f2(98)相乘,......... 等等等等,就象是在坐标轴上回卷一样。所以人们就叫它“回卷积分”,或者“卷积”了。
为了理解“卷积”的物理意义,不妨将那个问题“相当于它的时域的信号与系统的单位脉冲响应的卷积”略作变化。这个变化纯粹是为了方便表达和理解,不影响任何其它方面。将这个问题表述成这样一个问题:一个信号通过一个系统,系统的响应是频率响应或波谱响应,且看如何理解卷积的物理意义。
假设信号函数为f, 响应函数为g。f不仅是时间的函数(信号时有时无),还是频率的函数(就算在某一固定时刻,还有的地方大有的地方小);g也是时间的函数(有时候有反应,有时候没反应),同时也是频率的函数(不同的波长其响应程度不一样)。那我们要看某一时刻 t 的响应信号,该怎么办呢?
这就需要卷积了。
其实卷积积分应用广泛用在信号里面,一个是频域一个是时域
 

卷积是个啥?我忽然很想从本质上理解它。于是我从抽屉里翻出自己珍藏了许多年,每每下决心阅读却永远都读不完的《应用傅立叶变换》。
 
3.1 一维卷积的定义
 
函数f(x)与函数h(x)的卷积,由函参量的无穷积分

  定义。这里参量x和积分变量α皆为实数;函数f和h可实可复。
 
定义虽然找到了,但我还是一头雾水。卷积是个无穷积分吗?那它是干啥用的?再往后翻:几何说明、运算举例、基本性质,一堆的公式,就是没有说它是干啥用的。我于是坐在那呆想,忽然第二个困扰我的问题冒了出来:傅立叶变换是个啥?接着就是第三个、第四个、……、第N个问题。
 
傅立叶变换是个啥?听说能将时域上的东东变到频域上分析?哎?是变到频域上还是空间域上来着?到底啥是时域,频域,空间域?
 
上网查傅立叶变换的物理意义,没发现明确答案,倒发现了许多和我一样晕着问问题的人。结果又多出了许多名词,能量?功率谱?图像灰度域?……没办法又去翻那本教材。
 
1.1 一维傅立叶变换的定义与傅立叶积分定理
 
设f(x)是实变量x的函数,该函数可实可复,称积分

为函数f(x)的傅立叶变换。
 
吐血,啥是无穷积分来着?积分是啥来着?还能记起三角函数和差化积、积化和差公式吗?我忽然有种想把高中课本寻来重温的冲动。
 
卷积主要是为了将信号运算从时域转换为频域。
信号的时域的卷积等于频域的乘积。
利用这个性质以及特殊的δ函数可以通过抽样构造简单的调制电路
 
 
我比较赞同卷积的相关性的作用  在通信系统中的接收机部分MF匹配滤波器等就是本质上的相关
匹配滤波器最简单的形式就是原信号反转移位相乘积分得到的近似=相关
相关性越好得到的信号越强   这个我们有一次大作业做的  做地做到呕吐  呵呵
还有解调中一些东西本质就是相关
 

卷积公式  解释  卷积公式是用来求随机变量和的密度函数(pdf)的计算公式。  定义式:  z(t)=x(t)*y(t)= ∫x(m)y(t-m)dm.   已知x,y的pdf,x(t),y(t).现在要求z=x+y的pdf. 我们作变量替显,令  z=x+y,m=x. 雅可比行列式=1.那么,z,m联合密度就是f(z,m)=x(m)y(z-m)*1. 这样,就可以很容易求Z的在(z,m)中边缘分布  即fZ(z)=∫x(m)y(z-m)dm..... 由于这个公式和x(t),y(t)存在一一对应的关系。为了方便,所以记 ∫x(m)y(z-m)dm=x(t)*y(t)   长度为m的向量序列u和长度为n的向量序列v,卷积w的向量序列长度为(m+n-1),   u(n)与v(n)的卷积w(n)定义为: w(n)=u(n)@v(n)=sum(v(m)*u(n-m)),m from 负无穷到正无穷;   当m=n时w(1) = u(1)*v(1)   w(2) = u(1)*v(2)+u(2)*v(1)   w(3) = u(1)*v(3)+u(2)*v(2)+u(3)*v(1)   …   w(n) = u(1)*v(n)+u(2)*v(n-1)+ … +u(n)*v(1)   …   w(2*n-1) = u(n)*v(n)   当m≠n时,应以0补齐阶次低的向量的高位后进行计算  这是数学中常用的一个公式,在概率论中,是个重点也是一个难点。

  卷积公式是用来求随机变量和的密度函数(pdf)的计算公式。
  定义式:
  z(t)=x(t)*y(t)= ∫x(m)y(t-m)dm.
  已知x,y的pdf,x(t),y(t).现在要求z=x+y的pdf. 我们作变量替显,令
  z=x+y,m=x. 雅可比行列式=1.那么,t,m联合密度就是f(z,m)=x(m)y(z-m)*1. 这样,就可以很容易求Z的在(z,m)中边缘分布
  即fZ(z)=∫x(m)y(z-m)dm..... 由于这个公式和x(t),y(t)存在一一对应的关系。为了方便,所以记 ∫x(m)y(z-m)dm=x(t)*y(t)
 
卷积是一种线性运算,图像处理中常见的mask运算都是卷积,广泛应用于图像滤波。castlman的书对卷积讲得很详细。
高斯变换就是用高斯函数对图像进行卷积。高斯算子可以直接从离散高斯函数得到:
for(i=0; i<N; i++)
{
for(j=0; j<N; j++)
{
g[i*N+j]=exp(-((i-(N-1)/2)^2+(j-(N-1)/2)^2))/(2*delta^2));
sum += g[i*N+j];
}
}
再除以 sum 得到归一化算子
N是滤波器的大小,delta自选
首先,再提到卷积之前,必须提到卷积出现的背景。卷积是在信号与线性系统的基础上或背景中出现的,脱离这个背景单独谈卷积是没有任何意义的,除了那个所谓褶反公式上的数学意义和积分(或求和,离散情况下)。
信号与线性系统,讨论的就是信号经过一个线性系统以后发生的变化(就是输入输出和所经过的所谓系统,这三者之间的数学关系)。所谓线性系统的含义,就是,这个所谓的系统,带来的输出信号与输入信号的数学关系式之间是线性的运算关系。
因此,实际上,都是要根据我们需要待处理的信号形式,来设计所谓的系统传递函数,那么这个系统的传递函数和输入信号,在数学上的形式就是所谓的卷积关系。
卷积关系最重要的一种情况,就是在信号与线性系统或数字信号处理中的卷积定理。利用该定理,可以将时间域或空间域中的卷积运算等价为频率域的相乘运算,从而利用FFT等快速算法,实现有效的计算,节省运算代价。

参考:

http://www.cnblogs.com/a-toad/archive/2008/10/24/1318921.html

http://blog.sina.com.cn/s/blog_6d0e97bb01013op2.html

2006-12-21 21:02:00 housisong 阅读数 30504
  • 再论图像分割

    掌握图像去噪滤波方法及各自特点; 掌握常用的图像边缘检测算子,及Canny算子基本原理 掌握灰度阈值化及大津算法; 掌握常见的图像特征描述方式 了解局部阈值分割、区域生长、分水岭算法及基于轮廓的阈值分割; 通过...

    2537人学习 CSDN就业班
    免费试看

        图形图像处理-之-高质量的快速的图像缩放 中篇 二次线性插值和三次卷积插值
                     HouSisong@GMail.com   2006255.12.13

(2015.08.15  PicZoom_ftBilinear_Common更精确的边界公式推导)

(2009.03.07  可以到这里下载缩放算法的完整的可以编译的项目源代码:  http://blog.csdn.net/housisong/archive/2009/03/07/3967270.aspx  )

(2007.11.12 替换了二次线性插值的实现(以前偷懒使用了一个近似公式),改进后在图片边缘的插值效果更好(包括三次卷积插值的边界也更精确);
(2007.09.14 修正三次卷积的MMX版本中表的精度太低(7bit),造成卷积结果误差较大的问题,该版本提高了插值质量,并且速度加快12-25%)
(2007.09.07 PicZoom_ThreeOrder2和PicZoom_ThreeOrder_MMX在缩放的图片宽或高
小于3个像素的时候有一个Bug(边界计算错误);将unsigned long xrIntFloat_16,
yrIntFloat_16的定义改成long xrIntFloat_16,yrIntFloat_16就可以了)
(2007.07.02 ThreeOrder2_Fast一点小的改进,加快14%)
(2007.06.18 优化PicZoom_BilInear_MMX的实现(由138.5fps提高到147.9fps),
            并添加更快的两路展开的实现版本BilInear_MMX_expand2函数;
            补充新的SSE2的实现PicZoom_BilInear_SSE2函数)
(2007.06.06 更新测试数据,编译器由vc6改为vc2005,CPU由赛扬2G改为AMD64x2 4200+(2.1G) )
(2007.03.06 更新)


tag:图像缩放,速度优化,定点数优化,近邻取样插值,二次线性插值,三次线性插值,
   MipMap链,三次卷积插值,MMX,SSE,SSE2,CPU缓存优化

摘要:首先给出一个基本的图像缩放算法,然后一步一步的优化其速度和缩放质量;

高质量的快速的图像缩放 全文 分为:
     上篇 近邻取样插值和其速度优化
     中篇 二次线性插值和三次卷积插值
     下篇 三次线性插值和MipMap链

     补充 使用SSE2优化

正文:
  为了便于讨论,这里只处理32bit的ARGB颜色;
  代码使用C++;涉及到汇编优化的时候假定为x86平台;使用的编译器为vc2005;
  为了代码的可读性,没有加入异常处理代码;
   测试使用的CPU为AMD64x2 4200+(2.37G)  和 Intel Core2 4400(2.00G);


速度测试说明:
  只测试内存数据到内存数据的缩放
  测试图片都是800*600缩放到1024*768; fps表示每秒钟的帧数,值越大表示函数越快


A:近邻取样插值、二次线性插值、三次卷积插值 缩放效果对比

                                   
       原图         近邻取样缩放到0.6倍     近邻取样缩放到1.6倍

                                         
                二次线性插值缩放到0.6倍   二次线性插值缩放到1.6倍

                                       
               三次卷积插值缩放到0.6倍   三次卷积插值缩放到1.6倍

      
 原图 近邻取样缩放到8倍 二次线性插值缩放到8倍 三次卷积插值缩放到8倍 二次线性插值(近似公式)


     近邻取样插值缩放简单、速度快,但很多时候缩放出的图片质量比较差(特别是对于人物、景色等),
图片的缩放有比较明显的锯齿;使用二次或更高次插值有利于改善缩放效果;


B: 首先定义图像数据结构:

#define asm __asm

typedef unsigned 
char TUInt8; // [0..255]
struct TARGB32      //32 bit color
{
    TUInt8  b,g,r,a;          
//a is alpha
};

struct TPicRegion  //一块颜色数据区的描述,便于参数传递
{
    TARGB32
*    pdata;         //颜色数据首地址
    long        byte_width;    //一行数据的物理宽度(字节宽度);
                
//abs(byte_width)有可能大于等于width*sizeof(TARGB32);
    long        width;         //像素宽度
    long        height;        //像素高度
};

//那么访问一个点的函数可以写为:
inline TARGB32& Pixels(const TPicRegion& pic,const long x,const long y)
{
    
return ( (TARGB32*)((TUInt8*)pic.pdata+pic.byte_width*y) )[x];
}



二次线性插值缩放:

C: 二次线性插值缩放原理和公式图示:

        

             缩放后图片                 原图片
            (宽DW,高DH)              (宽SW,高SH)

  缩放映射原理:
  (Sx-0)/(SW-0)=(Dx-0)/(DW-0)   (Sy-0)/(SH-0)=(Dy-0)/(DH-0)
 =>   Sx=Dx*SW/DW                    Sy=Dy*SH/DH

  聚焦看看(Sx,Sy)坐标点(Sx,Sy为浮点数)附近的情况;


        


  对于近邻取样插值的缩放算法,直接取Color0颜色作为缩放后点的颜色;
二次线性插值需要考虑(Sx,Sy)坐标点周围的4个颜色值Color0/Color1/Color2/Color3,
把(Sx,Sy)到A/B/C/D坐标点的距离作为系数来把4个颜色混合出缩放后点的颜色;
( u=Sx-floor(Sx); v=Sy-floor(Sy); 说明:floor函数的返回值为小于等于参数的最大整数 ) 
  二次线性插值公式为:
 tmpColor0=Color0*(1-u) + Color2*u;
 tmpColor1=Color1*(1-u) + Color3*u;
        DstColor =tmpColor0*(1-v) + tmpColor2*v;

  展开公式为:
        pm0=(1-u)*(1-v);
        pm1=v*(1-u);
        pm2=u*(1-v);
        pm3=u*v;
  则颜色混合公式为:
        DstColor = Color0*pm0 + Color1*pm1 + Color2*pm2 + Color3*pm3;

参数函数图示:

     

                                              二次线性插值函数图示

对于上面的公式,它将图片向右下各移动了半个像素,需要对此做一个修正;
  =>   Sx=(Dx+0.5)*SW/DW-0.5; Sy=(Dy+0.5)*SH/DH-0.5;
而实际的程序,还需要考虑到边界(访问源图片可能超界)对于算法的影响,边界的处理可能有各种
方案(不处理边界或边界回绕或边界饱和或边界映射或用背景颜色混合等;文章中默认使用边界饱和来处理超界);
比如:边界饱和函数: 

//访问一个点的函数,(x,y)坐标可能超出图片边界; //边界处理模式:边界饱和
inline TARGB32 Pixels_Bound(const TPicRegion& pic,long x,long y)
{
    
//assert((pic.width>0)&&(pic.height>0));
    bool IsInPic=true;
    
if (x<0) {x=0; IsInPic=false; } else if (x>=pic.width ) {x=pic.width -1; IsInPic=false; }
    
if (y<0) {y=0; IsInPic=false; } else if (y>=pic.height) {y=pic.height-1; IsInPic=false; }
    TARGB32 result
=Pixels(pic,x,y);
    
if (!IsInPic) result.a=0;
    
return result;
}


D: 二次线性插值缩放算法的一个参考实现:PicZoom_BilInear0
  该函数并没有做什么优化,只是一个简单的浮点实现版本;


    inline void Bilinear0(const TPicRegion& pic,float fx,float fy,TARGB32* result)
    {
        
long x=(long)fx; if (x>fx) --x; //x=floor(fx);    
        long y=(long)fy; if (y>fy) --y; //y=floor(fy);
        
        TARGB32 Color0
=Pixels_Bound(pic,x,y);
        TARGB32 Color2
=Pixels_Bound(pic,x+1,y);
        TARGB32 Color1
=Pixels_Bound(pic,x,y+1);
        TARGB32 Color3
=Pixels_Bound(pic,x+1,y+1);

        
float u=fx-x;
        
float v=fy-y;
        
float pm3=u*v;
        
float pm2=u*(1-v);
        
float pm1=v*(1-u);
        
float pm0=(1-u)*(1-v);

        result
->a=(pm0*Color0.a+pm1*Color1.a+pm2*Color2.a+pm3*Color3.a);
        result
->r=(pm0*Color0.r+pm1*Color1.r+pm2*Color2.r+pm3*Color3.r);
        result
->g=(pm0*Color0.g+pm1*Color1.g+pm2*Color2.g+pm3*Color3.g);
        result
->b=(pm0*Color0.b+pm1*Color1.b+pm2*Color2.b+pm3*Color3.b);
    }

void PicZoom_Bilinear0(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;

    unsigned 
long dst_width=Dst.width;
    TARGB32
* pDstLine=Dst.pdata;
    
for (unsigned long y=0;y<Dst.height;++y)
    {
        
float srcy=(y+0.4999999)*Src.height/Dst.height-0.5;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            
float srcx=(x+0.4999999)*Src.width/Dst.width-0.5;
            Bilinear0(Src,srcx,srcy,
&pDstLine[x]);
        }
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_BilInear0      8.3 fps
////////////////////////////////////////////////////////////////////////////////

E: 把PicZoom_BilInear0的浮点计算改写为定点数实现:PicZoom_BilInear1

    inline void Bilinear1(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
    {
        
long x=x_16>>16;
        
long y=y_16>>16;
        TARGB32 Color0
=Pixels_Bound(pic,x,y);
        TARGB32 Color2
=Pixels_Bound(pic,x+1,y);
        TARGB32 Color1
=Pixels_Bound(pic,x,y+1);
        TARGB32 Color3
=Pixels_Bound(pic,x+1,y+1);

        unsigned 
long u_8=(x_16 & 0xFFFF)>>8;
        unsigned 
long v_8=(y_16 & 0xFFFF)>>8;
        unsigned 
long pm3_16=(u_8*v_8);
        unsigned 
long pm2_16=(u_8*(unsigned long)(256-v_8));
        unsigned 
long pm1_16=(v_8*(unsigned long)(256-u_8));
        unsigned 
long pm0_16=((256-u_8)*(256-v_8));

        result
->a=((pm0_16*Color0.a+pm1_16*Color1.a+pm2_16*Color2.a+pm3_16*Color3.a)>>16);
        result
->r=((pm0_16*Color0.r+pm1_16*Color1.r+pm2_16*Color2.r+pm3_16*Color3.r)>>16);
        result
->g=((pm0_16*Color0.g+pm1_16*Color1.g+pm2_16*Color2.g+pm3_16*Color3.g)>>16);
        result
->b=((pm0_16*Color0.b+pm1_16*Color1.b+pm2_16*Color2.b+pm3_16*Color3.b)>>16);
    }

void PicZoom_Bilinear1(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;

    
long xrIntFloat_16=((Src.width)<<16)/Dst.width
    
long yrIntFloat_16=((Src.height)<<16)/Dst.height;
    
const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
    
const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

    unsigned 
long dst_width=Dst.width;

    TARGB32
* pDstLine=Dst.pdata;
    
long srcy_16=csDErrorY;
    
long y;
    
for (y=0;y<Dst.height;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear1(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }


////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_BilInear1     17.7 fps
////////////////////////////////////////////////////////////////////////////////

F: 二次线性插值需要考略边界访问超界的问题,我们可以将边界区域和内部区域分开处理,这样就可以优化内部的插值实现函数了:比如不需要判断访问超界、减少颜色数据复制、减少一些不必要的重复坐标计算等等



    inline void Bilinear2_Fast(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
    {
        unsigned 
long pm3_16=u_8*v_8;
        unsigned 
long pm2_16=(u_8<<8)-pm3_16;
        unsigned 
long pm1_16=(v_8<<8)-pm3_16;
        unsigned 
long pm0_16=(1<<16)-pm1_16-pm2_16-pm3_16;
   
        result
->a=((pm0_16*PColor0[0].a+pm2_16*PColor0[1].a+pm1_16*PColor1[0].a+pm3_16*PColor1[1].a)>>16);
        result
->r=((pm0_16*PColor0[0].r+pm2_16*PColor0[1].r+pm1_16*PColor1[0].r+pm3_16*PColor1[1].r)>>16);
        result
->g=((pm0_16*PColor0[0].g+pm2_16*PColor0[1].g+pm1_16*PColor1[0].g+pm3_16*PColor1[1].g)>>16);
        result
->b=((pm0_16*PColor0[0].b+pm2_16*PColor0[1].b+pm1_16*PColor1[0].b+pm3_16*PColor1[1].b)>>16);
    }

    inline 
void Bilinear2_Border(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
    {
        
long x=(x_16>>16);
        
long y=(y_16>>16);
        unsigned 
long u_16=((unsigned short)(x_16));
        unsigned 
long v_16=((unsigned short)(y_16));

        TARGB32 pixel[
4];
        pixel[
0]=Pixels_Bound(pic,x,y);
        pixel[
1]=Pixels_Bound(pic,x+1,y);
        pixel[
2]=Pixels_Bound(pic,x,y+1);
        pixel[
3]=Pixels_Bound(pic,x+1,y+1);
        
        Bilinear2_Fast(
&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
    }

void PicZoom_Bilinear2(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;

    
long xrIntFloat_16=((Src.width)<<16)/Dst.width
    
long yrIntFloat_16=((Src.height)<<16)/Dst.height;
    
const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
    
const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

    unsigned 
long dst_width=Dst.width;

    
//计算出需要特殊处理的边界
    long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
    if (border_y0>=Dst.height) border_y0=Dst.height;
    
long border_x0=-csDErrorX/xrIntFloat_16+1;     
    
if (border_x0>=Dst.width ) border_x0=Dst.width; 
    
long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1//y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
    if (border_y1<border_y0) border_y1=border_y0;
    
long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1
    
if (border_x1<border_x0) border_x1=border_x0;

    TARGB32
* pDstLine=Dst.pdata;
    
long Src_byte_width=Src.byte_width;
    
long srcy_16=csDErrorY;
    
long y;
    
for (y=0;y<border_y0;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear2_Border(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y0;y<border_y1;++y)
    {
        
long srcx_16=csDErrorX;
        
long x;
        
for (x=0;x<border_x0;++x)
        {
            Bilinear2_Border(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }

        {
            unsigned 
long v_8=(srcy_16 & 0xFFFF)>>8;
            TARGB32
* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
            
for (unsigned long x=border_x0;x<border_x1;++x)
            {
                TARGB32
* PColor0=&PSrcLineColor[srcx_16>>16];
                TARGB32
* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                Bilinear2_Fast(PColor0,PColor1,(srcx_16 
& 0xFFFF)>>8,v_8,&pDstLine[x]);
                srcx_16
+=xrIntFloat_16;
            }
        }

        
for (x=border_x1;x<dst_width;++x)
        {
            Bilinear2_Border(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y1;y<Dst.height;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear2_Border(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_BilInear2     43.4 fps
////////////////////////////////////////////////////////////////////////////////

(F'补充:
  如果不想处理边界访问超界问题,可以考虑扩大源图片的尺寸,加一个边框 (“哨兵”优化);
这样插值算法就不用考虑边界问题了,程序写起来也简单很多! 
  如果对缩放结果的边界像素级精度要求不是太高,我还有一个方案,一个稍微改变的缩放公式:
  Sx=Dx*(SW-1)/DW; Sy=Dy*(SH-1)/DH;  (源图片宽和高:SW>=2;SH>=2)
  证明这个公式不会造成内存访问超界: 


    //fit的基本想法,不超界情况下尽量增大inc的值

    //fit算法推导:
    //寻址不超界要求 (inc*(D-1)>>16)+1 <= S-1
    //即:  inc*(D-1)/(1<<16) <= (S-2)
    //     inc*(D-1) <= (S-2)*(1<<16) +((1<<16)-1)
    //     inc*(D-1) <=  ((S-1)<<16) -1
    //     inc <= (((S-1)<<16) -1)/(D-1)
  比如,按这个公式的一个简单实现: (缩放效果见前面的"二次线性插值(近似公式)"图示)
void PicZoom_ftBilinear_Common(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(2>Src.width)||(2>Src.height)) return;

    
long xrIntFloat_16=((Src.width-1)<<16)/Dst.width; 
    
long yrIntFloat_16=((Src.height-1)<<16)/Dst.height;

    unsigned 
long dst_width=Dst.width;
    
long Src_byte_width=Src.byte_width;
    TARGB32
* pDstLine=Dst.pdata;
    
long srcy_16=0;
    
for (unsigned long y=0;y<Dst.height;++y)
    {
        unsigned 
long v_8=(srcy_16 & 0xFFFF)>>8;
        TARGB32
* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
        
long srcx_16=0;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            TARGB32
* PColor0=&PSrcLineColor[srcx_16>>16];
            Bilinear_Fast_Common(PColor0,(TARGB32
*)((TUInt8*)(PColor0)+Src_byte_width),(srcx_16 & 0xFFFF)>>8,v_8,&pDstLine[x]);
            srcx_16
+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }


)

G:利用单指令多数据处理的MMX指令一般都可以加快颜色的运算;在使用MMX改写之前,利用
32bit寄存器(或变量)来模拟单指令多数据处理;
数据储存原理:一个颜色数据分量只有一个字节,用2个字节来储存单个颜色分量的计算结果,
对于很多颜色计算来说精度就够了;那么一个32bit寄存器(或变量)就可以储存2个计算出的
临时颜色分量;从而达到了单个指令两路数据处理的目的;
单个指令两路数据处理的计算:
  乘法: ((0x00AA*a)<<16) | (0x00BB*a) = 0x00AA00BB * a
    可见只要保证0x00AA*a和0x00BB*a都小于(1<<16)那么乘法可以直接使用无符号数乘法了
  加法: ((0x00AA+0x00CC)<<16) | (0x00BB+0x00DD) = 0x00AA00BB + 0x00CC00DD
    可见只要0x00AA+0x00CC和0x00BB+0x00DD小于(1<<16)那么加法可以直接使用无符号数加法了
  (移位、减法等稍微复杂一点,因为这里没有用到就不推倒运算公式了)


    inline void Bilinear_Fast_Common(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
    {
        unsigned 
long pm3_8=(u_8*v_8)>>8;
        unsigned 
long pm2_8=u_8-pm3_8;
        unsigned 
long pm1_8=v_8-pm3_8;
        unsigned 
long pm0_8=256-pm1_8-pm2_8-pm3_8;

        unsigned 
long Color=*(unsigned long*)(PColor0);
        unsigned 
long BR=(Color & 0x00FF00FF)*pm0_8;
        unsigned 
long GA=((Color & 0xFF00FF00)>>8)*pm0_8;
                      Color
=((unsigned long*)(PColor0))[1];
                      GA
+=((Color & 0xFF00FF00)>>8)*pm2_8;
                      BR
+=(Color & 0x00FF00FF)*pm2_8;
                      Color
=*(unsigned long*)(PColor1);
                      GA
+=((Color & 0xFF00FF00)>>8)*pm1_8;
                      BR
+=(Color & 0x00FF00FF)*pm1_8;
                      Color
=((unsigned long*)(PColor1))[1];
                      GA
+=((Color & 0xFF00FF00)>>8)*pm3_8;
                      BR
+=(Color & 0x00FF00FF)*pm3_8;

        
*(unsigned long*)(result)=(GA & 0xFF00FF00)|((BR & 0xFF00FF00)>>8);
    }

    inline 
void Bilinear_Border_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
    {
        
long x=(x_16>>16);
        
long y=(y_16>>16);
        unsigned 
long u_16=((unsigned short)(x_16));
        unsigned 
long v_16=((unsigned short)(y_16));

        TARGB32 pixel[
4];
        pixel[
0]=Pixels_Bound(pic,x,y);
        pixel[
1]=Pixels_Bound(pic,x+1,y);
        pixel[
2]=Pixels_Bound(pic,x,y+1);
        pixel[
3]=Pixels_Bound(pic,x+1,y+1);
        
        Bilinear_Fast_Common(
&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
    }

void PicZoom_Bilinear_Common(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;

    
long xrIntFloat_16=((Src.width)<<16)/Dst.width
    
long yrIntFloat_16=((Src.height)<<16)/Dst.height;
    
const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
    
const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

    unsigned 
long dst_width=Dst.width;

    
//计算出需要特殊处理的边界
    long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
    if (border_y0>=Dst.height) border_y0=Dst.height;
    
long border_x0=-csDErrorX/xrIntFloat_16+1;     
    
if (border_x0>=Dst.width ) border_x0=Dst.width; 
    
long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1//y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
    if (border_y1<border_y0) border_y1=border_y0;
    
long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1
    
if (border_x1<border_x0) border_x1=border_x0;

    TARGB32
* pDstLine=Dst.pdata;
    
long Src_byte_width=Src.byte_width;
    
long srcy_16=csDErrorY;
    
long y;
    
for (y=0;y<border_y0;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear_Border_Common(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y0;y<border_y1;++y)
    {
        
long srcx_16=csDErrorX;
        
long x;
        
for (x=0;x<border_x0;++x)
        {
            Bilinear_Border_Common(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }

        {
            unsigned 
long v_8=(srcy_16 & 0xFFFF)>>8;
            TARGB32
* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
            
for (unsigned long x=border_x0;x<border_x1;++x)
            {
                TARGB32
* PColor0=&PSrcLineColor[srcx_16>>16];
                TARGB32
* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                Bilinear_Fast_Common(PColor0,PColor1,(srcx_16 
& 0xFFFF)>>8,v_8,&pDstLine[x]);
                srcx_16
+=xrIntFloat_16;
            }
        }

        
for (x=border_x1;x<dst_width;++x)
        {
            Bilinear_Border_Common(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y1;y<Dst.height;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear_Border_Common(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_BilInear_Common   65.3 fps
////////////////////////////////////////////////////////////////////////////////

H:使用MMX指令改写:PicZoom_Bilinear_MMX

    inline void  Bilinear_Fast_MMX(TARGB32* PColor0,TARGB32* PColor1,unsigned long u_8,unsigned long v_8,TARGB32* result)
    {
        asm
        {    
              MOVD      MM6,v_8
              MOVD      MM5,u_8
              mov       edx,PColor0
              mov       eax,PColor1
              PXOR      mm7,mm7

              MOVD         MM2,dword ptr [eax]  
              MOVD         MM0,dword ptr [eax
+4]
              PUNPCKLWD    MM5,MM5
              PUNPCKLWD    MM6,MM6
              MOVD         MM3,dword ptr [edx]  
              MOVD         MM1,dword ptr [edx
+4]
              PUNPCKLDQ    MM5,MM5 
              PUNPCKLBW    MM0,MM7
              PUNPCKLBW    MM1,MM7
              PUNPCKLBW    MM2,MM7
              PUNPCKLBW    MM3,MM7
              PSUBw        MM0,MM2
              PSUBw        MM1,MM3
              PSLLw        MM2,
8
              PSLLw        MM3,
8
              PMULlw       MM0,MM5
              PMULlw       MM1,MM5
              PUNPCKLDQ    MM6,MM6 
              PADDw        MM0,MM2
              PADDw        MM1,MM3

              PSRLw        MM0,
8
              PSRLw        MM1,
8
              PSUBw        MM0,MM1
              PSLLw        MM1,
8
              PMULlw       MM0,MM6
              mov       eax,result
              PADDw        MM0,MM1

              PSRLw        MM0,
8
              PACKUSwb     MM0,MM7
              movd      [eax],MM0 
              
//emms
        }
    }

    
void Bilinear_Border_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
    {
        
long x=(x_16>>16);
        
long y=(y_16>>16);
        unsigned 
long u_16=((unsigned short)(x_16));
        unsigned 
long v_16=((unsigned short)(y_16));

        TARGB32 pixel[
4];
        pixel[
0]=Pixels_Bound(pic,x,y);
        pixel[
1]=Pixels_Bound(pic,x+1,y);
        pixel[
2]=Pixels_Bound(pic,x,y+1);
        pixel[
3]=Pixels_Bound(pic,x+1,y+1);
        
        Bilinear_Fast_MMX(
&pixel[0],&pixel[2],u_16>>8,v_16>>8,result);
    }

void PicZoom_Bilinear_MMX(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;

    
long xrIntFloat_16=((Src.width)<<16)/Dst.width
    
long yrIntFloat_16=((Src.height)<<16)/Dst.height;
    
const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
    
const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

    unsigned 
long dst_width=Dst.width;

    
//计算出需要特殊处理的边界
    long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
    if (border_y0>=Dst.height) border_y0=Dst.height;
    
long border_x0=-csDErrorX/xrIntFloat_16+1;     
    
if (border_x0>=Dst.width ) border_x0=Dst.width; 
    
long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1//y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
    if (border_y1<border_y0) border_y1=border_y0;
    
long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1
    
if (border_x1<border_x0) border_x1=border_x0;

    TARGB32
* pDstLine=Dst.pdata;
    
long Src_byte_width=Src.byte_width;
    
long srcy_16=csDErrorY;
    
long y;
    
for (y=0;y<border_y0;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y0;y<border_y1;++y)
    {
        
long srcx_16=csDErrorX;
        
long x;
        
for (x=0;x<border_x0;++x)
        {
            Bilinear_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }

        {
            unsigned 
long v_8=(srcy_16 & 0xFFFF)>>8;
            TARGB32
* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
            
for (unsigned long x=border_x0;x<border_x1;++x)
            {
                TARGB32
* PColor0=&PSrcLineColor[srcx_16>>16];
                TARGB32
* PColor1=(TARGB32*)((TUInt8*)(PColor0)+Src_byte_width);
                Bilinear_Fast_MMX(PColor0,PColor1,(srcx_16 
& 0xFFFF)>>8,v_8,&pDstLine[x]);
                srcx_16
+=xrIntFloat_16;
            }
        }

        
for (x=border_x1;x<dst_width;++x)
        {
            Bilinear_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y1;y<Dst.height;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    asm emms
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_BilInear_MMX 132.9 fps
////////////////////////////////////////////////////////////////////////////////

 

H' 对BilInear_MMX简单改进:PicZoom_Bilinear_MMX_Ex


void PicZoom_Bilinear_MMX_Ex(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;

    
long xrIntFloat_16=((Src.width)<<16)/Dst.width
    
long yrIntFloat_16=((Src.height)<<16)/Dst.height;
    
const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
    
const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

    unsigned 
long dst_width=Dst.width;

    
//计算出需要特殊处理的边界
    long border_y0=-csDErrorY/yrIntFloat_16+1;              //y0+y*yr>=0; y0=csDErrorY => y>=-csDErrorY/yr
    if (border_y0>=Dst.height) border_y0=Dst.height;
    
long border_x0=-csDErrorX/xrIntFloat_16+1;     
    
if (border_x0>=Dst.width ) border_x0=Dst.width; 
    
long border_y1=(((Src.height-2)<<16)-csDErrorY)/yrIntFloat_16+1//y0+y*yr<=(height-2) => y<=(height-2-csDErrorY)/yr
    if (border_y1<border_y0) border_y1=border_y0;
    
long border_x1=(((Src.width-2)<<16)-csDErrorX)/xrIntFloat_16+1
    
if (border_x1<border_x0) border_x1=border_x0;

    TARGB32
* pDstLine=Dst.pdata;
    
long Src_byte_width=Src.byte_width;
    
long srcy_16=csDErrorY;
    
long y;
    
for (y=0;y<border_y0;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }

    
for (y=border_y0;y<border_y1;++y)
    {
        
long srcx_16=csDErrorX;
        
long x;
        
for (x=0;x<border_x0;++x)
        {
            Bilinear_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }

        {
            
long dst_width_fast=border_x1-border_x0;
            
if (dst_width_fast>0)
            {
                unsigned 
long v_8=(srcy_16 & 0xFFFF)>>8;
                TARGB32
* PSrcLineColor= (TARGB32*)((TUInt8*)(Src.pdata)+Src_byte_width*(srcy_16>>16)) ;
                TARGB32
* PSrcLineColorNext= (TARGB32*)((TUInt8*)(PSrcLineColor)+Src_byte_width) ;
                TARGB32
* pDstLine_Fast=&pDstLine[border_x0];
                asm
                {
                      movd         mm6,v_8
                      pxor         mm7,mm7 
//mm7=0
                      PUNPCKLWD    MM6,MM6
                      PUNPCKLDQ    MM6,MM6
//mm6=v_8
                    
                      mov       esi,PSrcLineColor
                      mov       ecx,PSrcLineColorNext
                      mov       edx,srcx_16
                      mov       ebx,dst_width_fast
                      mov       edi,pDstLine_Fast
                      lea       edi,[edi
+ebx*4]
                      push      ebp
                      mov       ebp,xrIntFloat_16
                      neg       ebx

                loop_start:

                          mov       eax,edx
                          shl       eax,
16
                          shr       eax,
24
                          
//== movzx       eax,dh  //eax=u_8
                          MOVD      MM5,eax
                          mov       eax,edx
                          shr       eax,
16     //srcx_16>>16

                          MOVD         MM2,dword ptr [ecx
+eax*4]  
                          MOVD         MM0,dword ptr [ecx
+eax*4+4]
                          PUNPCKLWD    MM5,MM5
                          MOVD         MM3,dword ptr [esi
+eax*4]  
                          MOVD         MM1,dword ptr [esi
+eax*4+4]
                          PUNPCKLDQ    MM5,MM5 
//mm5=u_8
                          PUNPCKLBW    MM0,MM7
                          PUNPCKLBW    MM1,MM7
                          PUNPCKLBW    MM2,MM7
                          PUNPCKLBW    MM3,MM7
                          PSUBw        MM0,MM2
                          PSUBw        MM1,MM3
                          PSLLw        MM2,
8
                          PSLLw        MM3,
8
                          PMULlw       MM0,MM5
                          PMULlw       MM1,MM5
                          PADDw        MM0,MM2
                          PADDw        MM1,MM3

                          PSRLw        MM0,
8
                          PSRLw        MM1,
8
                          PSUBw        MM0,MM1
                          PSLLw        MM1,
8
                          PMULlw       MM0,MM6
                          PADDw        MM0,MM1

                          PSRLw     MM0,
8
                          PACKUSwb  MM0,MM7
                          MOVd   dword ptr    [edi
+ebx*4],MM0 //write DstColor
                                      
                          add       edx,ebp 
//srcx_16+=xrIntFloat_16
                          inc       ebx
                          jnz       loop_start

                      pop       ebp
                      mov       srcx_16,edx
                }
            }
        }

        
for (x=border_x1;x<dst_width;++x)
        {
            Bilinear_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y1;y<Dst.height;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            Bilinear_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    asm emms
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_Bilinear_MMX_Ex 157.0 fps
////////////////////////////////////////////////////////////////////////////////

I: 把测试成绩放在一起:

////////////////////////////////////////////////////////////////////////////////
//CPU: AMD64x2 4200+(2.37G)  zoom 800*600 to 1024*768
//==============================================================================
// StretchBlt                   232.7 fps  
// PicZoom3_SSE                 711.7 fps
//
// PicZoom_BilInear0              8.3 fps
// PicZoom_BilInear1             17.7 fps
// PicZoom_BilInear2             43.4 fps
// PicZoom_BilInear_Common       65.3 fps
// PicZoom_BilInear_MMX         132.9 fps
// PicZoom_BilInear_MMX_Ex      157.0 fps
////////////////////////////////////////////////////////////////////////////////

补充Intel Core2 4400上的测试成绩:

////////////////////////////////////////////////////////////////////////////////
//CPU: Intel Core2 4400(2.00G)  zoom 800*600 to 1024*768
//==============================================================================
// PicZoom3_SSE                1099.7 fps  
//
// PicZoom_BilInear0             10.7 fps
// PicZoom_BilInear1             24.2 fps
// PicZoom_BilInear2             54.3 fps
// PicZoom_BilInear_Common       59.8 fps
// PicZoom_BilInear_MMX         118.4 fps
// PicZoom_BilInear_MMX_Ex      142.9 fps
////////////////////////////////////////////////////////////////////////////////

 

 

三次卷积插值:

 


J: 二次线性插值缩放出的图片很多时候让人感觉变得模糊(术语叫低通滤波),特别是在放大
的时候;使用三次卷积插值来改善插值结果;三次卷积插值考虑映射点周围16个点(4x4)的颜色来
计算最终的混合颜色,如图;

         
         P(0,0)所在像素为映射的点,加上它周围的15个点,按一定系数混合得到最终输出结果;

         混合公式参见PicZoom_ThreeOrder0的实现;

    插值曲线公式sin(x*PI)/(x*PI),如图:


                             三次卷积插值曲线sin(x*PI)/(x*PI) (其中PI=3.1415926...)

 

K:三次卷积插值缩放算法的一个参考实现:PicZoom_ThreeOrder0
  该函数并没有做过多的优化,只是一个简单的浮点实现版本;



        inline 
double SinXDivX(double x) 
        {
            
//该函数计算插值曲线sin(x*PI)/(x*PI)的值 //PI=3.1415926535897932385; 
            
//下面是它的近似拟合表达式
            const float a = -1//a还可以取 a=-2,-1,-0.75,-0.5等等,起到调节锐化或模糊程度的作用

            
if (x<0) x=-x; //x=abs(x);
            double x2=x*x;
            
double x3=x2*x;
            
if (x<=1)
              
return (a+2)*x3 - (a+3)*x2 + 1;
            
else if (x<=2
              
return a*x3 - (5*a)*x2 + (8*a)*- (4*a);
            
else
              
return 0;
        } 

        inline TUInt8 border_color(
long Color)
        {
            
if (Color<=0)
                
return 0;
            
else if (Color>=255)
                
return 255;
            
else
                
return Color;
        }
        
    
void ThreeOrder0(const TPicRegion& pic,const float fx,const float fy,TARGB32* result)
    {
        
long x0=(long)fx; if (x0>fx) --x0; //x0=floor(fx);    
        long y0=(long)fy; if (y0>fy) --y0; //y0=floor(fy);
        float fu=fx-x0;
        
float fv=fy-y0;

        TARGB32 pixel[
16];
        
long i,j;

        
for (i=0;i<4;++i)
        {
            
for (j=0;j<4;++j)
            {
                
long x=x0-1+j;
                
long y=y0-1+i;
                pixel[i
*4+j]=Pixels_Bound(pic,x,y);
            }
        }

        
float afu[4],afv[4];
        
//
        afu[0]=SinXDivX(1+fu);
        afu[
1]=SinXDivX(fu);
        afu[
2]=SinXDivX(1-fu);
        afu[
3]=SinXDivX(2-fu);
        afv[
0]=SinXDivX(1+fv);
        afv[
1]=SinXDivX(fv);
        afv[
2]=SinXDivX(1-fv);
        afv[
3]=SinXDivX(2-fv);

        
float sR=0,sG=0,sB=0,sA=0;
        
for (i=0;i<4;++i)
        {
            
float aR=0,aG=0,aB=0,aA=0;
            
for (long j=0;j<4;++j)
            {
                aA
+=afu[j]*pixel[i*4+j].a;
                aR
+=afu[j]*pixel[i*4+j].r;
                aG
+=afu[j]*pixel[i*4+j].g;
                aB
+=afu[j]*pixel[i*4+j].b;
            }
            sA
+=aA*afv[i];
            sR
+=aR*afv[i];
            sG
+=aG*afv[i];
            sB
+=aB*afv[i];
        }

        result
->a=border_color((long)(sA+0.5));
        result
->r=border_color((long)(sR+0.5));
        result
->g=border_color((long)(sG+0.5));
        result
->b=border_color((long)(sB+0.5));
    }

void PicZoom_ThreeOrder0(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;


    unsigned 
long dst_width=Dst.width;
    TARGB32
* pDstLine=Dst.pdata;
    
for (unsigned long y=0;y<Dst.height;++y)
    {
        
float srcy=(y+0.4999999)*Src.height/Dst.height-0.5;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            
float srcx=(x+0.4999999)*Src.width/Dst.width-0.5;
            ThreeOrder0(Src,srcx,srcy,
&pDstLine[x]);
        }
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_ThreeOrder0    3.6 fps
////////////////////////////////////////////////////////////////////////////////

L: 使用定点数来优化缩放函数;边界和内部分开处理;对SinXDivX做一个查找表;对border_color做一个查找表;


    static long SinXDivX_Table_8[(2<<8)+1];
    
class _CAutoInti_SinXDivX_Table {
    
private
        
void _Inti_SinXDivX_Table()
        {
            
for (long i=0;i<=(2<<8);++i)
                SinXDivX_Table_8[i]
=long(0.5+256*SinXDivX(i*(1.0/(256))))*1;
        };
    
public:
        _CAutoInti_SinXDivX_Table() { _Inti_SinXDivX_Table(); }
    };
    
static _CAutoInti_SinXDivX_Table __tmp_CAutoInti_SinXDivX_Table;


    
//颜色查表
    static TUInt8 _color_table[256*3];
    
static const TUInt8* color_table=&_color_table[256];
    
class _CAuto_inti_color_table
    {
    
public:
        _CAuto_inti_color_table() {
            
for (int i=0;i<256*3;++i)
                _color_table[i]
=border_color(i-256);
        }
    };
    
static _CAuto_inti_color_table _Auto_inti_color_table;

    
void ThreeOrder_Fast_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
    {
        unsigned long u_8=(unsigned char)((x_16)>>8);
        unsigned 
long v_8=(unsigned char)((y_16)>>8);
        
const TARGB32* pixel=&Pixels(pic,(x_16>>16)-1,(y_16>>16)-1);
        
long pic_byte_width=pic.byte_width;

        
long au_8[4],av_8[4];
        
//
        au_8[0]=SinXDivX_Table_8[(1<<8)+u_8];
        au_8[
1]=SinXDivX_Table_8[u_8];
        au_8[
2]=SinXDivX_Table_8[(1<<8)-u_8];
        au_8[
3]=SinXDivX_Table_8[(2<<8)-u_8];
        av_8[
0]=SinXDivX_Table_8[(1<<8)+v_8];
        av_8[
1]=SinXDivX_Table_8[v_8];
        av_8[
2]=SinXDivX_Table_8[(1<<8)-v_8];
        av_8[
3]=SinXDivX_Table_8[(2<<8)-v_8];

        
long sR=0,sG=0,sB=0,sA=0;
        
for (long i=0;i<4;++i)
        {
            
long aA=au_8[0]*pixel[0].a + au_8[1]*pixel[1].a + au_8[2]*pixel[2].a + au_8[3]*pixel[3].a;
            
long aR=au_8[0]*pixel[0].r + au_8[1]*pixel[1].r + au_8[2]*pixel[2].r + au_8[3]*pixel[3].r;
            
long aG=au_8[0]*pixel[0].g + au_8[1]*pixel[1].g + au_8[2]*pixel[2].g + au_8[3]*pixel[3].g;
            
long aB=au_8[0]*pixel[0].b + au_8[1]*pixel[1].b + au_8[2]*pixel[2].b + au_8[3]*pixel[3].b;
            sA
+=aA*av_8[i];
            sR
+=aR*av_8[i];
            sG
+=aG*av_8[i];
            sB
+=aB*av_8[i];
            ((TUInt8
*&)pixel)+=pic_byte_width;
        }

        result
->a=color_table[sA>>16];
        result
->r=color_table[sR>>16];
        result
->g=color_table[sG>>16];
        result
->b=color_table[sB>>16];
    }

    
void ThreeOrder_Border_Common(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
    {
        
long x0_sub1=(x_16>>16)-1;
        
long y0_sub1=(y_16>>16)-1;
        unsigned 
long u_16_add1=((unsigned short)(x_16))+(1<<16);
        unsigned 
long v_16_add1=((unsigned short)(y_16))+(1<<16);

        TARGB32 pixel[
16];
        
long i;

        
for (i=0;i<4;++i)
        {
            
long y=y0_sub1+i;
            pixel[i
*4+0]=Pixels_Bound(pic,x0_sub1+0,y);
            pixel[i
*4+1]=Pixels_Bound(pic,x0_sub1+1,y);
            pixel[i
*4+2]=Pixels_Bound(pic,x0_sub1+2,y);
            pixel[i
*4+3]=Pixels_Bound(pic,x0_sub1+3,y);
        }
        
        TPicRegion npic;
        npic.pdata     
=&pixel[0];
        npic.byte_width
=4*sizeof(TARGB32);
        
//npic.width     =4;
        
//npic.height    =4;
        ThreeOrder_Fast_Common(npic,u_16_add1,v_16_add1,result);
    }

void PicZoom_ThreeOrder_Common(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;

    
long xrIntFloat_16=((Src.width)<<16)/Dst.width
    
long yrIntFloat_16=((Src.height)<<16)/Dst.height;
    
const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
    
const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

    unsigned 
long dst_width=Dst.width;

    
//计算出需要特殊处理的边界
    long border_y0=((1<<16)-csDErrorY)/yrIntFloat_16+1;//y0+y*yr>=1; y0=csDErrorY => y>=(1-csDErrorY)/yr
    if (border_y0>=Dst.height) border_y0=Dst.height;
    
long border_x0=((1<<16)-csDErrorX)/xrIntFloat_16+1;
    
if (border_x0>=Dst.width ) border_x0=Dst.width;
    
long border_y1=(((Src.height-3)<<16)-csDErrorY)/yrIntFloat_16+1//y0+y*yr<=(height-3) => y<=(height-3-csDErrorY)/yr
    if (border_y1<border_y0) border_y1=border_y0;
    
long border_x1=(((Src.width-3)<<16)-csDErrorX)/xrIntFloat_16+1;; 
    
if (border_x1<border_x0) border_x1=border_x0;

    TARGB32
* pDstLine=Dst.pdata;
    
long srcy_16=csDErrorY;
    
long y;
    
for (y=0;y<border_y0;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            ThreeOrder_Border_Common(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y0;y<border_y1;++y)
    {
        
long srcx_16=csDErrorX;
        
long x;
        
for (x=0;x<border_x0;++x)
        {
            ThreeOrder_Border_Common(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }
        
for (x=border_x0;x<border_x1;++x)
        {
            ThreeOrder_Fast_Common(Src,srcx_16,srcy_16,
&pDstLine[x]);//fast  !
            srcx_16+=xrIntFloat_16;
        }
        
for (x=border_x1;x<dst_width;++x)
        {
            ThreeOrder_Border_Common(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y1;y<Dst.height;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            ThreeOrder_Border_Common(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_ThreeOrder_Common    16.9 fps
////////////////////////////////////////////////////////////////////////////////

 M: 用MMX来优化ThreeOrder_Common函数:ThreeOrder_MMX


    typedef   unsigned long TMMXData32;
    
static TMMXData32 SinXDivX_Table_MMX[(2<<8)+1];
    
class _CAutoInti_SinXDivX_Table_MMX {
    
private
        
void _Inti_SinXDivX_Table_MMX()
        {
            
for (long i=0;i<=(2<<8);++i)
            {
                unsigned 
short t=long(0.5+(1<<14)*SinXDivX(i*(1.0/(256))));
                unsigned 
long tl=| (((unsigned long)t)<<16);
                SinXDivX_Table_MMX[i]
=tl;
            }
        };
    
public:
        _CAutoInti_SinXDivX_Table_MMX() { _Inti_SinXDivX_Table_MMX(); }
    };
    
static _CAutoInti_SinXDivX_Table_MMX __tmp_CAutoInti_SinXDivX_Table_MMX;


    
void __declspec(naked) _private_ThreeOrder_Fast_MMX()
    {
        asm
        {
            movd        mm1,dword ptr [edx]
            movd        mm2,dword ptr [edx
+4]
            movd        mm3,dword ptr [edx
+8]
            movd        mm4,dword ptr [edx
+12]
            movd        mm5,dword ptr [(offset SinXDivX_Table_MMX)
+256*4+eax*4]
            movd        mm6,dword ptr [(offset SinXDivX_Table_MMX)
+eax*4]
            punpcklbw   mm1,mm7
            punpcklbw   mm2,mm7
            punpcklwd   mm5,mm5
            punpcklwd   mm6,mm6
            psllw       mm1,
7
            psllw       mm2,
7
            pmulhw      mm1,mm5
            pmulhw      mm2,mm6
            punpcklbw   mm3,mm7
            punpcklbw   mm4,mm7
            movd        mm5,dword ptr [(offset SinXDivX_Table_MMX)
+256*4+ecx*4]
            movd        mm6,dword ptr [(offset SinXDivX_Table_MMX)
+512*4+ecx*4]
            punpcklwd   mm5,mm5
            punpcklwd   mm6,mm6
            psllw       mm3,
7
            psllw       mm4,
7
            pmulhw      mm3,mm5
            pmulhw      mm4,mm6
            paddsw      mm1,mm2
            paddsw      mm3,mm4
            movd        mm6,dword ptr [ebx] 
//v
            paddsw      mm1,mm3
            punpcklwd   mm6,mm6

            pmulhw      mm1,mm6
            add     edx,esi  
//+pic.byte_width
            paddsw      mm0,mm1

            ret
        }
    }

    inline 
void ThreeOrder_Fast_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
    {
        asm
        {
            mov     ecx,pic
            mov     eax,y_16
            mov     ebx,x_16
            movzx   edi,ah 
//v_8
            mov     edx,[ecx+TPicRegion::pdata]
            shr     eax,
16
            mov     esi,[ecx
+TPicRegion::byte_width]
            dec     eax
            movzx   ecx,bh 
//u_8
            shr     ebx,16
            imul    eax,esi
            lea     edx,[edx
+ebx*4-4]
            add     edx,eax 
//pixel

            mov     eax,ecx
            neg     ecx

            pxor    mm7,mm7  
//0
            
//mov     edx,pixel
            pxor    mm0,mm0  //result=0
            
//lea     eax,auv_7

            lea    ebx,[(offset SinXDivX_Table_MMX)
+256*4+edi*4]
            call  _private_ThreeOrder_Fast_MMX
            lea    ebx,[(offset SinXDivX_Table_MMX)
+edi*4]
            call  _private_ThreeOrder_Fast_MMX
            neg    edi
            lea    ebx,[(offset SinXDivX_Table_MMX)
+256*4+edi*4]
            call  _private_ThreeOrder_Fast_MMX
            lea    ebx,[(offset SinXDivX_Table_MMX)
+512*4+edi*4]
            call  _private_ThreeOrder_Fast_MMX

            psraw     mm0,
3
            mov       eax,result
            packuswb  mm0,mm7
            movd      [eax],mm0
            
//emms
        }
    }

    
void ThreeOrder_Border_MMX(const TPicRegion& pic,const long x_16,const long y_16,TARGB32* result)
    {
        unsigned 
long x0_sub1=(x_16>>16)-1;
        unsigned 
long y0_sub1=(y_16>>16)-1;
        
long u_16_add1=((unsigned short)(x_16))+(1<<16);
        
long v_16_add1=((unsigned short)(y_16))+(1<<16);

        TARGB32 pixel[
16];

        
for (long i=0;i<4;++i)
        {
            
long y=y0_sub1+i;
            pixel[i
*4+0]=Pixels_Bound(pic,x0_sub1  ,y);
            pixel[i
*4+1]=Pixels_Bound(pic,x0_sub1+1,y);
            pixel[i
*4+2]=Pixels_Bound(pic,x0_sub1+2,y);
            pixel[i
*4+3]=Pixels_Bound(pic,x0_sub1+3,y);
        }
        
        TPicRegion npic;
        npic.pdata     
=&pixel[0];
        npic.byte_width
=4*sizeof(TARGB32);
        
//npic.width     =4;
        
//npic.height    =4;
        ThreeOrder_Fast_MMX(npic,u_16_add1,v_16_add1,result);
    }

void PicZoom_ThreeOrder_MMX(const TPicRegion& Dst,const TPicRegion& Src)
{
    
if (  (0==Dst.width)||(0==Dst.height)
        
||(0==Src.width)||(0==Src.height)) return;

    
long xrIntFloat_16=((Src.width)<<16)/Dst.width
    
long yrIntFloat_16=((Src.height)<<16)/Dst.height;
    
const long csDErrorX=-(1<<15)+(xrIntFloat_16>>1);
    
const long csDErrorY=-(1<<15)+(yrIntFloat_16>>1);

    unsigned 
long dst_width=Dst.width;

    
//计算出需要特殊处理的边界
    long border_y0=((1<<16)-csDErrorY)/yrIntFloat_16+1;//y0+y*yr>=1; y0=csDErrorY => y>=(1-csDErrorY)/yr
    if (border_y0>=Dst.height) border_y0=Dst.height;
    
long border_x0=((1<<16)-csDErrorX)/xrIntFloat_16+1;
    
if (border_x0>=Dst.width ) border_x0=Dst.width;
    
long border_y1=(((Src.height-3)<<16)-csDErrorY)/yrIntFloat_16+1//y0+y*yr<=(height-3) => y<=(height-3-csDErrorY)/yr
    if (border_y1<border_y0) border_y1=border_y0;
    
long border_x1=(((Src.width-3)<<16)-csDErrorX)/xrIntFloat_16+1;; 
    
if (border_x1<border_x0) border_x1=border_x0;

    TARGB32
* pDstLine=Dst.pdata;
    
long srcy_16=csDErrorY;
    
long y;
    
for (y=0;y<border_y0;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y0;y<border_y1;++y)
    {
        
long srcx_16=csDErrorX;
        
long x;
        
for (x=0;x<border_x0;++x)
        {
            ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }
        
for (x=border_x0;x<border_x1;++x)
        {
            ThreeOrder_Fast_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]);//fast MMX !
            srcx_16+=xrIntFloat_16;
        }
        
for (x=border_x1;x<dst_width;++x)
        {
            ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]);//border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }
    
for (y=border_y1;y<Dst.height;++y)
    {
        
long srcx_16=csDErrorX;
        
for (unsigned long x=0;x<dst_width;++x)
        {
            ThreeOrder_Border_MMX(Src,srcx_16,srcy_16,
&pDstLine[x]); //border
            srcx_16+=xrIntFloat_16;
        }
        srcy_16
+=yrIntFloat_16;
        ((TUInt8
*&)pDstLine)+=Dst.byte_width;
    }

    asm emms
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicZoom_ThreeOrder_MMX   34.3 fps
////////////////////////////////////////////////////////////////////////////////

 

 N:将测试结果放到一起:

////////////////////////////////////////////////////////////////////////////////
//CPU: AMD64x2 4200+(2.37G)  zoom 800*600 to 1024*768
//==============================================================================
// StretchBlt                   232.7 fps  
// PicZoom3_SSE                 711.7 fps  
// PicZoom_BilInear_MMX_Ex      157.0 fps
//
// PicZoom_ThreeOrder0            3.6 fps
// PicZoom_ThreeOrder_Common     16.9 fps
// PicZoom_ThreeOrder_MMX        34.3 fps
////////////////////////////////////////////////////////////////////////////////

补充Intel Core2 4400上的测试成绩:

////////////////////////////////////////////////////////////////////////////////
//CPU: Intel Core2 4400(2.00G)  zoom 800*600 to 1024*768
//==============================================================================
// PicZoom3_SSE                1099.7 fps  
// PicZoom_BilInear_MMX_Ex      142.9 fps
//
// PicZoom_ThreeOrder0            4.2 fps
// PicZoom_ThreeOrder_Common     17.6 fps
// PicZoom_ThreeOrder_MMX        34.4 fps
////////////////////////////////////////////////////////////////////////////////

 

2019-04-23 16:24:29 Eastmount 阅读数 19610
  • 再论图像分割

    掌握图像去噪滤波方法及各自特点; 掌握常用的图像边缘检测算子,及Canny算子基本原理 掌握灰度阈值化及大津算法; 掌握常见的图像特征描述方式 了解局部阈值分割、区域生长、分水岭算法及基于轮廓的阈值分割; 通过...

    2537人学习 CSDN就业班
    免费试看

该系列文章是讲解Python OpenCV图像处理知识,前期主要讲解图像入门、OpenCV基础用法,中期讲解图像处理的各种算法,包括图像锐化算子、图像增强技术、图像分割等,后期结合深度学习研究图像识别、图像分类应用。希望文章对您有所帮助,如果有不足之处,还请海涵~

该系列在github所有源代码:https://github.com/eastmountyxz/ImageProcessing-Python
PS:请求帮忙点个Star,哈哈,第一次使用Github,以后会分享更多代码,一起加油。

同时推荐作者的C++图像系列知识:
[数字图像处理] 一.MFC详解显示BMP格式图片
[数字图像处理] 二.MFC单文档分割窗口显示图片
[数字图像处理] 三.MFC实现图像灰度、采样和量化功能详解
[数字图像处理] 四.MFC对话框绘制灰度直方图
[数字图像处理] 五.MFC图像点运算之灰度线性变化、灰度非线性变化、阈值化和均衡化处理详解
[数字图像处理] 六.MFC空间几何变换之图像平移、镜像、旋转、缩放详解
[数字图像处理] 七.MFC图像增强之图像普通平滑、高斯平滑、Laplacian、Sobel、Prewitt锐化详解

前文参考:
[Python图像处理] 一.图像处理基础知识及OpenCV入门函数
[Python图像处理] 二.OpenCV+Numpy库读取与修改像素
[Python图像处理] 三.获取图像属性、兴趣ROI区域及通道处理
[Python图像处理] 四.图像平滑之均值滤波、方框滤波、高斯滤波及中值滤波
[Python图像处理] 五.图像融合、加法运算及图像类型转换
[Python图像处理] 六.图像缩放、图像旋转、图像翻转与图像平移
[Python图像处理] 七.图像阈值化处理及算法对比
[Python图像处理] 八.图像腐蚀与图像膨胀
[Python图像处理] 九.形态学之图像开运算、闭运算、梯度运算
[Python图像处理] 十.形态学之图像顶帽运算和黑帽运算
[Python图像处理] 十一.灰度直方图概念及OpenCV绘制直方图
[Python图像处理] 十二.图像几何变换之图像仿射变换、图像透视变换和图像校正
[Python图像处理] 十三.基于灰度三维图的图像顶帽运算和黑帽运算
[Python图像处理] 十四.基于OpenCV和像素处理的图像灰度化处理
[Python图像处理] 十五.图像的灰度线性变换
[Python图像处理] 十六.图像的灰度非线性变换之对数变换、伽马变换
[Python图像处理] 十七.图像锐化与边缘检测之Roberts算子、Prewitt算子、Sobel算子和Laplacian算子
[Python图像处理] 十八.图像锐化与边缘检测之Scharr算子、Canny算子和LOG算子
[Python图像处理] 十九.图像分割之基于K-Means聚类的区域分割
[Python图像处理] 二十.图像量化处理和采样处理及局部马赛克特效
[Python图像处理] 二十一.图像金字塔之图像向下取样和向上取样

前面一篇文章我讲解了Python图像量化、采样处理及图像金字塔。本文主要讲解图像傅里叶变换的相关内容,在数字图像处理中,有两个经典的变换被广泛应用——傅里叶变换和霍夫变换。其中,傅里叶变换主要是将时间域上的信号转变为频率域上的信号,用来进行图像除噪、图像增强等处理。基础性文章,希望对你有所帮助。同时,该部分知识均为杨秀璋查阅资料撰写,转载请署名CSDN+杨秀璋及原地址出处,谢谢!!

1.图像傅里叶变换
2.Numpy实现傅里叶变换
3.Numpy实现傅里叶逆变换
4.OpenCV实现傅里叶变换
5.OpenCV实现傅里叶逆变换


PS:文章参考自己以前系列图像处理文章及OpenCV库函数,同时参考如下文献:
《数字图像处理》(第3版),冈萨雷斯著,阮秋琦译,电子工业出版社,2013年.
《数字图像处理学》(第3版),阮秋琦,电子工业出版社,2008年,北京.
《OpenCV3编程入门》,毛星云,冷雪飞,电子工业出版社,2015,北京.
百度百科-傅里叶变换
网易云课堂-高登教育 Python+OpenCV图像处理
安安zoe-图像的傅里叶变换
daduzimama-图像的傅里叶变换的迷思----频谱居中
tenderwx-数字图像处理-傅里叶变换在图像处理中的应用
小小猫钓小小鱼-深入浅出的讲解傅里叶变换(真正的通俗易懂)


一.图像傅里叶变换原理

傅里叶变换(Fourier Transform,简称FT)常用于数字信号处理,它的目的是将时间域上的信号转变为频率域上的信号。随着域的不同,对同一个事物的了解角度也随之改变,因此在时域中某些不好处理的地方,在频域就可以较为简单的处理。同时,可以从频域里发现一些原先不易察觉的特征。傅里叶定理指出“任何连续周期信号都可以表示成(或者无限逼近)一系列正弦信号的叠加。”

下面引用李老师 “Python+OpenCV图像处理” 中的一个案例,非常推荐同学们去购买学习。如下图所示,他将某饮料的制作过程的时域角度转换为频域角度。

绘制对应的时间图和频率图如下所示:

傅里叶公式如下,其中w表示频率,t表示时间,为复变函数。它将时间域的函数表示为频率域的函数f(t)的积分。

傅里叶变换认为一个周期函数(信号)包含多个频率分量,任意函数(信号)f(t)可通过多个周期函数(或基函数)相加合成。从物理角度理解,傅里叶变换是以一组特殊的函数(三角函数)为正交基,对原函数进行线性变换,物理意义便是原函数在各组基函数的投影。如下图所示,它是由三条正弦曲线组合成。

傅里叶变换可以应用于图像处理中,经过对图像进行变换得到其频谱图。从谱频图里频率高低来表征图像中灰度变化剧烈程度。图像中的边缘信号和噪声信号往往是高频信号,而图像变化频繁的图像轮廓及背景等信号往往是低频信号。这时可以有针对性的对图像进行相关操作,例如图像除噪、图像增强和锐化等。

二维图像的傅里叶变换可以用以下数学公式(15-3)表达,其中f是空间域(Spatial Domain))值,F是频域(Frequency Domain)值

对上面的傅里叶变换有了大致的了解之后,下面通过Numpy和OpenCV分别讲解图像傅里叶变换的算法及操作代码。


二.Numpy实现傅里叶变换

Numpy中的 FFT包提供了函数 np.fft.fft2()可以对信号进行快速傅里叶变换,其函数原型如下所示,该输出结果是一个复数数组(Complex Ndarry)。

fft2(a, s=None, axes=(-2, -1), norm=None)

  • a表示输入图像,阵列状的复杂数组
  • s表示整数序列,可以决定输出数组的大小。输出可选形状(每个转换轴的长度),其中s[0]表示轴0,s[1]表示轴1。对应fit(x,n)函数中的n,沿着每个轴,如果给定的形状小于输入形状,则将剪切输入。如果大于则输入将用零填充。如果未给定’s’,则使用沿’axles’指定的轴的输入形状
  • axes表示整数序列,用于计算FFT的可选轴。如果未给出,则使用最后两个轴。“axes”中的重复索引表示对该轴执行多次转换,一个元素序列意味着执行一维FFT
  • norm包括None和ortho两个选项,规范化模式(请参见numpy.fft)。默认值为无

Numpy中的fft模块有很多函数,相关函数如下:

#计算一维傅里叶变换
numpy.fft.fft(a, n=None, axis=-1, norm=None)
#计算二维的傅里叶变换
numpy.fft.fft2(a, n=None, axis=-1, norm=None)
#计算n维的傅里叶变换
numpy.fft.fftn()
#计算n维实数的傅里叶变换
numpy.fft.rfftn()
#返回傅里叶变换的采样频率
numpy.fft.fftfreq()
#将FFT输出中的直流分量移动到频谱中央
numpy.fft.shift()

下面的代码是通过Numpy库实现傅里叶变换,调用np.fft.fft2()快速傅里叶变换得到频率分布,接着调用np.fft.fftshift()函数将中心位置转移至中间,最终通过Matplotlib显示效果图。

# -*- coding: utf-8 -*-
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

#读取图像
img = cv.imread('test.png', 0)

#快速傅里叶变换算法得到频率分布
f = np.fft.fft2(img)

#默认结果中心点位置是在左上角,
#调用fftshift()函数转移到中间位置
fshift = np.fft.fftshift(f)       

#fft结果是复数, 其绝对值结果是振幅
fimg = np.log(np.abs(fshift))

#展示结果
plt.subplot(121), plt.imshow(img, 'gray'), plt.title('Original Fourier')
plt.axis('off')
plt.subplot(122), plt.imshow(fimg, 'gray'), plt.title('Fourier Fourier')
plt.axis('off')
plt.show()

输出结果如图15-2所示,左边为原始图像,右边为频率分布图谱,其中越靠近中心位置频率越低,越亮(灰度值越高)的位置代表该频率的信号振幅越大。


三.Numpy实现傅里叶逆变换

下面介绍Numpy实现傅里叶逆变换,它是傅里叶变换的逆操作,将频谱图像转换为原始图像的过程。通过傅里叶变换将转换为频谱图,并对高频(边界)和低频(细节)部分进行处理,接着需要通过傅里叶逆变换恢复为原始效果图。频域上对图像的处理会反映在逆变换图像上,从而更好地进行图像处理。

图像傅里叶变化主要使用的函数如下所示:

#实现图像逆傅里叶变换,返回一个复数数组
numpy.fft.ifft2(a, n=None, axis=-1, norm=None)
#fftshit()函数的逆函数,它将频谱图像的中心低频部分移动至左上角
numpy.fft.fftshift()
#将复数转换为0至255范围
iimg = numpy.abs(逆傅里叶变换结果)

下面的代码分别实现了傅里叶变换和傅里叶逆变换。

# -*- coding: utf-8 -*-
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

#读取图像
img = cv.imread('Lena.png', 0)

#傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
res = np.log(np.abs(fshift))

#傅里叶逆变换
ishift = np.fft.ifftshift(fshift)
iimg = np.fft.ifft2(ishift)
iimg = np.abs(iimg)

#展示结果
plt.subplot(131), plt.imshow(img, 'gray'), plt.title('Original Image')
plt.axis('off')
plt.subplot(132), plt.imshow(res, 'gray'), plt.title('Fourier Image')
plt.axis('off')
plt.subplot(133), plt.imshow(iimg, 'gray'), plt.title('Inverse Fourier Image')
plt.axis('off')
plt.show()

输出结果如图15-4所示,从左至右分别为原始图像、频谱图像、逆傅里叶变换转换图像。


四.OpenCV实现傅里叶变换

OpenCV 中相应的函数是cv2.dft()和用Numpy输出的结果一样,但是是双通道的。第一个通道是结果的实数部分,第二个通道是结果的虚数部分,并且输入图像要首先转换成 np.float32 格式。其函数原型如下所示:

dst = cv2.dft(src, dst=None, flags=None, nonzeroRows=None)

  • src表示输入图像,需要通过np.float32转换格式
  • dst表示输出图像,包括输出大小和尺寸
  • flags表示转换标记,其中DFT _INVERSE执行反向一维或二维转换,而不是默认的正向转换;DFT _SCALE表示缩放结果,由阵列元素的数量除以它;DFT _ROWS执行正向或反向变换输入矩阵的每个单独的行,该标志可以同时转换多个矢量,并可用于减少开销以执行3D和更高维度的转换等;DFT _COMPLEX_OUTPUT执行1D或2D实数组的正向转换,这是最快的选择,默认功能;DFT _REAL_OUTPUT执行一维或二维复数阵列的逆变换,结果通常是相同大小的复数数组,但如果输入数组具有共轭复数对称性,则输出为真实数组
  • nonzeroRows表示当参数不为零时,函数假定只有nonzeroRows输入数组的第一行(未设置)或者只有输出数组的第一个(设置)包含非零,因此函数可以处理其余的行更有效率,并节省一些时间;这种技术对计算阵列互相关或使用DFT卷积非常有用

注意,由于输出的频谱结果是一个复数,需要调用cv2.magnitude()函数将傅里叶变换的双通道结果转换为0到255的范围。其函数原型如下:

cv2.magnitude(x, y)

  • x表示浮点型X坐标值,即实部
  • y表示浮点型Y坐标值,即虚部
    最终输出结果为幅值,即:

完整代码如下所示:

# -*- coding: utf-8 -*-
import numpy as np
import cv2
from matplotlib import pyplot as plt

#读取图像
img = cv2.imread('Lena.png', 0)

#傅里叶变换
dft = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)

#将频谱低频从左上角移动至中心位置
dft_shift = np.fft.fftshift(dft)

#频谱图像双通道复数转换为0-255区间
result = 20*np.log(cv2.magnitude(dft_shift[:,:,0], dft_shift[:,:,1]))

#显示图像
plt.subplot(121), plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(result, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

输出结果如图15-5所示,左边为原始“Lena”图,右边为转换后的频谱图像,并且保证低频位于中心位置。


五.OpenCV实现傅里叶逆变换

在OpenCV 中,通过函数cv2.idft()实现傅里叶逆变换,其返回结果取决于原始图像的类型和大小,原始图像可以为实数或复数。其函数原型如下所示:

dst = cv2.idft(src[, dst[, flags[, nonzeroRows]]])

  • src表示输入图像,包括实数或复数
  • dst表示输出图像
  • flags表示转换标记
  • nonzeroRows表示要处理的dst行数,其余行的内容未定义(请参阅dft描述中的卷积示例)

完整代码如下所示:

# -*- coding: utf-8 -*-
import numpy as np
import cv2
from matplotlib import pyplot as plt

#读取图像
img = cv2.imread('Lena.png', 0)

#傅里叶变换
dft = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)
dftshift = np.fft.fftshift(dft)
res1= 20*np.log(cv2.magnitude(dftshift[:,:,0], dftshift[:,:,1]))

#傅里叶逆变换
ishift = np.fft.ifftshift(dftshift)
iimg = cv2.idft(ishift)
res2 = cv2.magnitude(iimg[:,:,0], iimg[:,:,1])

#显示图像
plt.subplot(131), plt.imshow(img, 'gray'), plt.title('Original Image')
plt.axis('off')
plt.subplot(132), plt.imshow(res1, 'gray'), plt.title('Fourier Image')
plt.axis('off')
plt.subplot(133), plt.imshow(res2, 'gray'), plt.title('Inverse Fourier Image')
plt.axis('off')
plt.show()

输出结果如图15-6所示,第一幅图为原始“Lena”图,第二幅图为傅里叶变换后的频谱图像,第三幅图为傅里叶逆变换,频谱图像转换为原始图像的过程。


六.总结

傅里叶变换的目的并不是为了观察图像的频率分布(至少不是最终目的),更多情况下是为了对频率进行过滤,通过修改频率以达到图像增强、图像去噪、边缘检测、特征提取、压缩加密等目的。下一篇文章,作者将结合傅里叶变换和傅里叶逆变换讲解它的应用。

时也,命也。
英语低分数线一分,些许遗憾,但不气馁,更加努力。雄关漫道真如铁,而今迈过从头越,从头越。苍山如海,残阳如血。感谢一路陪伴的人和自己。

无论成败,那段拼搏的日子都很美。结果只会让我更加努力,学好英语。下半年沉下心来好好做科研写文章,西藏之行,课程分享。同时,明天的博士考试加油,虽然裸泳,但也加油!还有春季招考开始准备。

最后补充马刺小石匠精神,当一切都看起来无济于事的时候,我去看一个石匠敲石头.他一连敲了100次,石头仍然纹丝不动。但他敲第101次的时候,石头裂为两半。可我知道,让石头裂开的不是那最后一击,而是前面的一百次敲击的结果。人生路漫漫,不可能一路一帆风顺,暂时的不顺只是磨练自己的必经之路,夜最深的时候也是距黎明最近的时刻,经历过漫漫长夜的打磨,你自身会更加强大。

最后希望这篇基础性文章对您有所帮助,如果有错误或不足之处,请海涵!

(By:Eastmount 2019-04-23 周二下午6点写于花溪 https://blog.csdn.net/Eastmount )

2012-01-06 15:53:38 jia20003 阅读数 27145
  • 再论图像分割

    掌握图像去噪滤波方法及各自特点; 掌握常用的图像边缘检测算子,及Canny算子基本原理 掌握灰度阈值化及大津算法; 掌握常见的图像特征描述方式 了解局部阈值分割、区域生长、分水岭算法及基于轮廓的阈值分割; 通过...

    2537人学习 CSDN就业班
    免费试看

 

图像噪声源于现实世界中数字信号总会受到各种各样的干扰,最终接受的图像和源于的数字信号之间总

是存在一定的差异,对于图像噪声,使用均值滤波和中值滤波来消除图像噪声的做法已经是很常见的图

像消噪手段。

 

一:图像加噪原理

1.     椒盐噪声(Salt And Pepper Noise)

椒盐噪声是一种因为信号脉冲强度引起的噪声,信噪比(Signal NoiseRate)是衡量图像噪声的一个数字指标。

给一副数字图像加上椒盐噪声的处理顺序应该如下:

  1. 指定信噪比 SNR 其取值范围在[0, 1]之间
  2. 计算总像素数目 SP, 得到要加噪的像素数目 NP = SP * (1-SNR)
  3. 随机获取要加噪的每个像素位置P(i, j)
  4. 指定像素值为255或者0。
  5. 重复c, d两个步骤完成所有像素的NP个像素
  6. 输出加噪以后的图像

 

2.     高斯噪声(Gaussian Noise)

高斯噪声的密度取决于公式G(x, sigma) 其中X是代表平均值,sigma代表的标准方差,每个输入像素 Pin, 

一个正常的高斯采样分布公式G(d), 得到输出像素Pout.

       Pout = Pin + XMeans + sigma *G(d)

其中d为一个线性的随机数,G(d)是随机数的高斯分布随机值。

给一副数字图像加上高斯噪声的处理顺序如下:

a.      输入参数sigam 和 X mean

b.      以系统时间为种子产生一个伪随机数

c.      将伪随机数带入G(d)得到高斯随机数

d.      根据输入像素计算出输出像素

e.      重新将像素值防缩在[0 ~ 255]之间

f.       循环所有像素

g.      输出图像

 

二:关键程序解析

1.     椒盐噪声

根据信噪比,获取要加入椒盐噪声的像素数目

int size= (int)(inPixels.length * (1-SNR));

 

随机得到像素,完成椒盐噪声的加入

for(int i=0; i<size; i++) {

int row = (int)(Math.random()* (double)height);

int col = (int)(Math.random()* (double)width);

index= row * width + col;

inPixels[index]= (255 << 24) | (255 << 16) | (255 << 8) | 255;

}

 

2.     高斯噪声

根据标准方差,和伪随机数的范围,首先计算出一个伪随机数d ,根据d得到高斯分布的随机数值,整个代码如下:

    float d = (float)Math.random()*RANDOM_SCOPE - RANDOM_SCOPE/2;

    float sigma2 = sigma*sigma*2;

    float PI2 = (float)Math.PI * 2;

    float sigmaPI2 = (float)Math.sqrt(PI2*sigma);

    float result = (float)Math.exp(-d/sigma2)/sigmaPI2;

伪随机数的范围为[-127~ 127]之间。

 

获取高斯噪声的像素代码如下:

tr = (int)((float)tr + getGaussianValue() + this.means);

tg = (int)((float)tg + getGaussianValue() + this.means);

tb = (int)((float)tb + getGaussianValue() + this.means);

mean是的值为0.

 

三:程序效果如下


加入白色椒盐噪声的图片

 


加入高斯噪声的图片



椒盐噪声的代码如下:

	private BufferedImage addSaltAndPepperNoise(BufferedImage src, BufferedImage dst) {
		int width = src.getWidth();
        int height = src.getHeight();

        if ( dst == null )
            dst = createCompatibleDestImage( src, null );

        int[] inPixels = new int[width*height];
        getRGB( src, 0, 0, width, height, inPixels );
        
        int index = 0;
        int size = (int)(inPixels.length * (1-SNR));

        for(int i=0; i<size; i++) {
        	int row = (int)(Math.random() * (double)height);
        	int col = (int)(Math.random() * (double)width);
        	index = row * width + col;
        	inPixels[index] = (255 << 24) | (255 << 16) | (255 << 8) | 255;
        }

        setRGB( dst, 0, 0, width, height, inPixels );
        return dst;
	}

高斯噪声的代码如下:

private BufferedImage gaussianNoise(BufferedImage src, BufferedImage dst) {
		int width = src.getWidth();
        int height = src.getHeight();

        if ( dst == null )
            dst = createCompatibleDestImage( src, null );

        int[] inPixels = new int[width*height];
        int[][][] tempPixels = new int[height][width][4]; 
        int[] outPixels = new int[width*height];
        getRGB( src, 0, 0, width, height, inPixels );
        int index = 0;
        float inMax = 0;
        float outMax = 0;
        for(int row=0; row<height; row++) {
        	int ta = 0, tr = 0, tg = 0, tb = 0;
        	for(int col=0; col<width; col++) {
        		index = row * width + col;
        		ta = (inPixels[index] >> 24) & 0xff;
                tr = (inPixels[index] >> 16) & 0xff;
                tg = (inPixels[index] >> 8) & 0xff;
                tb = inPixels[index] & 0xff;
                if(inMax < tr) {
                	inMax = tr;
                }
                if(inMax < tg) {
                	inMax = tg;
                }
                if(inMax < tb) {
                	inMax = tb;
                }
                tr = (int)((float)tr + getGaussianValue() + this.means);
                tg = (int)((float)tg + getGaussianValue() + this.means);
                tb = (int)((float)tb + getGaussianValue() + this.means);
                if(outMax < tr) {
                	outMax = tr;
                }
                if(outMax < tg) {
                	outMax = tg;
                }
                if(outMax < tb) {
                	outMax = tb;
                }
                tempPixels[row][col][0] = ta;
                tempPixels[row][col][1] = tr;
                tempPixels[row][col][2] = tg;
                tempPixels[row][col][3] = tb;
        	}
        }

        // Normalization
        index = 0;
        float rate = inMax/outMax;
        for(int row=0; row<height; row++) {
        	int ta = 0, tr = 0, tg = 0, tb = 0;
        	for(int col=0; col<width; col++) {
        		index = row * width + col;
        		ta = tempPixels[row][col][0];
        		tr = tempPixels[row][col][1];
        		tg = tempPixels[row][col][2];
        		tb = tempPixels[row][col][3];

        		tr = (int)((float)tr * rate);
        		tg = (int)((float)tg * rate);
        		tb = (int)((float)tb * rate);
        		outPixels[index] = (ta << 24) | (tr << 16) | (tg << 8) | tb;
        	}
        }
        setRGB( dst, 0, 0, width, height, outPixels );
        return dst;
	}