图像处理的角度

2013-06-03 16:41:12 Augusdi 阅读数 46967

通俗的说:

膨胀算法使图像扩大一圈。
腐蚀算法使二值图像减小一圈。


腐蚀:删除对象边界的某些像素
膨胀:给图像中的对象边界添加像素


算法:


从图像处理角度看,二值图像的腐蚀和膨胀就是将一个小型二值图(结构元素,一般为3*3大小)在一个大的二值图上逐点移动并进行比较,根据比较的结果作出相应处理而已。

膨胀算法:用3X3的结构元素,扫描二值图像的每一个像素,用结构元素与其覆盖的二值图像做“与”运算,如果都为0,结构图像的该像素为0,否则为1.结果:使二值图像扩大一圈。


腐蚀算法:用3X3的结构元素,扫描二值图像的每一个像素,用结构元素与其覆盖的二值图像做“与”运算,如果都为1,结构图像的该像素为1,否则为0.结果:使二值图像减小一圈。


注意:

OpenCV里面的腐蚀膨胀都是针对 白色 目标区域的。

说,膨胀使图像 变大一圈, 那是指 图像中的 白色目标区域 扩大了一圈~~~大笑~~~腐蚀,亦然。

====================================

代码:

#include <stdio.h>
#include <cv.h>
#include <highgui.h>
#include <cvaux.h> //必须引此头文件

 

#pragma comment(lib, "cv.lib") #pragma comment(lib, "cxcore.lib") #pragma comment(lib, "highgui.lib")

int main( int argc, char** argv ) { IplImage *img = cvLoadImage("rice.png",0); if(img == NULL) { printf("img load failed!\n"); return 0; } IplImage *img_erode = cvCreateImage(cvGetSize(img), 8, 1); IplImage *img_dilate = cvCreateImage(cvGetSize(img), 8, 1); cvErode( img,img_erode, NULL,1); //腐蚀 cvDilate( img,img_dilate, NULL,1); //膨胀 cvNamedWindow("img_erode"); cvNamedWindow("img_dilate"); cvShowImage("img_erode",img_erode); cvShowImage("img_dilate",img_dilate); cvWaitKey(-1); cvReleaseImage(&img_dilate); cvReleaseImage(&img_erode); cvDestroyAllWindows(); return 0; }
图示:



附:

Dilate 使用任意结构元素膨胀图像

void cvDilate( const CvArr* src, CvArr* dst, IplConvKernel* element=NULL, int iterations=1 );

src
输入图像.
dst
输出图像.
element
用于膨胀的结构元素。若为 NULL, 则使用 3×3 长方形的结构元素
iterations
膨胀的次数


对彩色图像,每个彩色通道单独处理。


Erode 使用任意结构元素腐蚀图像

void cvErode( const CvArr* src, CvArr* dst, IplConvKernel* element=NULL, int iterations=1 );

src
输入图像.
dst
输出图像.
element
用于腐蚀的结构元素。若为 NULL, 则使用 3×3 长方形的结构元素
iterations
腐蚀的次数


函数 cvErode 对输入图像使用指定的结构元素进行腐蚀,该结构元素决定每个具有最小值象素点的邻域形状:

2017-04-15 19:34:01 m0_37264397 阅读数 61094

    傅立叶变换在图像处理中有非常重要的作用。因为不仅傅立叶分析涉及图像处理很多方面,傅立 叶改进算法,比如离散余弦变换,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

2015-12-29 11:21:24 u013088062 阅读数 32635

  最近有人问我图像处理怎么研究,怎么入门,怎么应用,我竟一时语塞。仔细想想,自己也搞了两年图像方面的研究,做个两个创新项目,发过两篇论文,也算是有点心得,于是总结总结和大家分享,希望能对大家有所帮助。在写这篇教程之前我本想多弄点插图,让文章看起来花哨一点,后来我觉得没必要这样做,大家花时间沉下心来读读文字没什么不好,况且学术和技术本身也不是多么花哨的东西。

  一、图像处理的应用

  这个其实没什么好说的,一种技术的应用价值不是靠嘴上说,而是要看有多少人去搞,很简单的道理。其实我觉得判断一项技术有没有价值、有多大价值的最简单最有效的方法就是看有多少人在研究它。如果大家都在研究,那必然说明它很火,至少目前很火,以后的几年里依然还会火。因此,如果你不确定图像处理是不是有价值,那就查查全国图像处理工程师的数量就行了。

当然这里还是简单提一下。如果你真的只想问“图像处理有什么用?”,相信百度会给出比我更专业的答案。不过作为图像处理的行内人,我还是从几个基本的点来具体说一下。

  1、身份认证

  21世纪是刷脸的时代,这一点无可厚非。首先是银行,据说重庆的银行已经使用了人脸识别的验证系统来进行辅助的认证。其次是门禁系统,以前是指纹,虹膜,现在就是人脸。指纹、虹膜的识别虽然准确,但是有侵犯性,采集的过程中有侵犯性,验证的过程中也有侵犯性,反正让谁天天录指纹(采集指纹信息),用眼睛瞪摄像头(采集虹膜信息),谁都会觉得不舒服的,况且手还会脱皮。相比之下,人脸的识别就方便多了,拍张照片(采集人脸信息)谁都不会介意。最后是监控,一个摄像头所拍的监控能从不同的角度记录成百上千的人(比如车站等密集场所的监控),让警察去辨认将是一个浩大的工程,要是系统能够自动判别人员信息,无疑会给办案带来极大方便。

  2、监控安防

  安防监控可以说是图像处理领域最有潜力的应用领域。如今各个城市都在疯狂的安装监控摄像头,全国时刻都有无数的摄像头在录监控,但是安防监控的后端处理却没跟上。什么是后端处理,简单的说就是监控录像的视频处理。注意这里的视频处理可就不止包含人脸识别了,还有行人检测、异常检测、显著性检测、协同跟踪等。人脸识别之前说了,这里简单说说行人异常检测。图像处理中的行人异常检测在外行人眼中是一个非常有魔力的东西。毕竟让摄像头通过监控视频就能判断出当前画面中谁是好人谁是坏人(当然直接分为好人和坏人还是太武断了),在一般思维中貌似是不可能的。但是请不要忽略这样一个事实,就是目前大部分监控视频的分析处理都是由人工来完成的,警察破案时经常动不动就调出最近几天的监控视频,然后从头看到尾,其工程量之大可想而知。也正是这样才催生了人们对智能监控方面的研究,因为有这方面的实际需求。当然我们的视频分析程序不会直接给出诸如“好人or坏人”这样的武断而片面的判断结果。就当前的技术水平而言,能够统计当下监控画面中的人数(行人检测)、定位他们的人脸(人脸检测)、识别他们的身份(人脸识别)、判别他们的表情(表情识别)、检测他们的动作(异常检测),这就已经够了。这样以后人们就不会再面对长达几十甚至上百个小时的监控视频发呆,而是直接分析计算机给出的数据,当前画面中有多少人、都是谁、谁的动作可疑等。总之,接下来智能监控发展会很迅速,因为需求很迫切。

  3、深度学习(Deep Learning)

  通俗的说一句,“图像处理是深度学习应用的天堂”。深度学习这个概念可能有人不太熟悉,大家可以自行百度,我这里给出一个相对通俗的解释:“如果你不知道什么叫深度学习,就想想《终结者》中的T-800”。当然这句话不是我说的,是出自业界的一位大牛之口。当然这可能有点小片面,不过深度学习确实是公认的新一代人工智能的基础。

  这里举两个例子。一是谷歌的人工大脑项目。谷歌公司可以说是深度学习方面的牵头企业了,其在2012年公布的谷歌大脑项目,动用了1.6万个计算节点,训练数周,得到的人工智能模型能够自主识别猫脸图像,为新一代人工智能开辟了道路,之后微软深度学习研究院、百度深度学习研究院等机构都开始大量投入,各个高校也搞得风声水起,原因很简单,大家都知道它会火。

  第二就是图像识别方面的竞赛。最有权威的就是ImageNet竞赛。大家在一个拥有上千万张,上千类别的图像数据库上训练和测试自己的算法,比拼识别率。近几年来,摘得桂冠的一直都是深度学习模型,确切的说是卷积神经网络。更多有关ImageNet历年的竞赛信息大家可以自行百度。

  说道深度学习在图像处理的应用,不得不提中国的汤晓鸥教授,说他是国内深度学习的领头羊也不为过。他提出的DeepID人脸识别算法(一共分为三代),在一些大规模人脸数据库上的正确率(若LFW库)已经达到了99.75%,单纯从数字上讲的话可以说已经超越了人类的识别率,为此汤教授还开办了公司,开发FaceSDK(虽然还没有公布)。不过拿计算机和人脑相比本身就是不合理的,各有所长嘛。不过可见DeepLearning在图像识别领域的强大威力。至于深度学习与图像处理的关系这里就不用多说了,谷歌大脑识别的是图像,深度学习竞赛用的是图像,DeepID识别的还是图像人脸,虽然深度学习在其他方面诸如语音识别等也有应用,在图像处理依然是其主要的应用领域。

  二、图像处理研究工具

  图像处理的研究分为算法研究和应用两个部分。用到的主要编程语言有Matlab、C/C++、Python等,原因很简单,它们都有很多相应的第三方库,不用我们从零开始编程。

  1、Matlab

  MathWork公司的Matlab软件可以说是算法研究的利器,它的强大之处在于其方便快捷的矩阵运算能力和图形仿真能力,单从简洁性和封装性来说,确实完爆其他语言。但高度封装必然会相应的损失一部分灵活性,况且Matlab严格的讲更像是一个工具,而非一门编程语言。顺便提一句,它在2015年编程语言排行榜中位于第20名,仅次于IOS开发的Objective-C。

  对于算法研究人员(尤其是高校的硕士博士),首选工具自然是matlab,因为它简便快捷,封装性好,更重要的是全世界几乎所有的算法大牛、精英教授都会首先公布对应的Matlab源码,然后在逐步改写成其他语言进行实际应用。所以,如果你想做图像处理方面的研究,Matlab是必须掌握的,而且是熟练掌握。当你有一些想法需要验证时,最好明智的先用matlab编写出来测试。如果你上来就用看似高大上的C++来实验,不仅错误BUG一大堆,到头来可能效果还不佳,就算效果好,时间也会耽搁不少,毕竟算法开发还是要快的,这样才能赶在别人之前发论文。总之,只要是接触图像算法,终究逃不过Matlab,就算你是软件开发的,不研发算法,但总得能看懂别人的Matlab算法吧。

  对于之前没怎么接触过Matlab与图像处理的人,在这里推荐一本相关的书籍《MATLAB图像处理实例详解(附光盘)》。这本书对于Matlab图像处理入门还是很有帮助的。记得我当时刚上研究生时就靠两本书入门的,一是冈萨雷斯的《数字图像处理》,二是这本《MATLAB图像处理实例详解》。不过这里友情提示,在看这类教程(不仅仅是Matlab)时千万不要试图去记忆所有的工具函数,这种做法是十分愚蠢的。正确的做法是根据自己的情况快速翻阅这类工具书,可以找出里面的有实际意义的源码来敲一敲练练手感,至于具体的工具函数,只需要知道Matlab提供了这方面的功能就行了,以后用到了再回来查,或者谷歌百度。我觉得在入门阶段,最重要的不是看了多少书,听了多少课,而是尽快自己敲出一段代码,运行出结果,来建立自信和成就感,这才是支持我们走下去的最实在的动力。记得我当时看了没多久就自己敲了一个蹩脚的车牌检测的Matlab程序,现在看来真是漏洞百出,不过当时我真的很兴奋,很有成就感,觉得自己能干这行,对于初学者来说,这种感受弥足珍贵。

  2、OpenCv

  Opencv是Intel公司开发的C++图像处理工具包,形象的理解为就是C++版的Matlab。当初Intel公司开发这个工具包的初衷也是方便大家共享,希望大家能够在一个共同架构的基础上共同建造摩天大楼,而不是各自在自己的地基上盖平房。与Matlab不同,Opencv是面向开发的,稳定性好,异常处理机制周全,但有一点需要注意,由于Opencv是开源的,那么如果你在项目中直接调用了它的API,那就意味着你的项目也必须开源。因此在真正的产品开发过程中,往往需要从Opencv库里面挖代码,而不是直接调用,幸好Intel公司允许我们看源码,自己编译一把就可以了。

  说道C++和Opencv,有一个问题不得不提,那就是深度学习领域大名鼎鼎的Caffe框架。这是一个典型的基于C++和OpenCv的深度学习框架,由谷歌深度学习团队、“谷歌大脑”负责人之一贾扬清学者编写,并公布了源码。如今各个深度学习机构都在大量使用这个框架进行研究。

  这里同样对推荐两本关于Opencv方面的教程。一本是CSDN博客大牛毛星云写的《OpenCV3编程入门》,这是它根据自己多年的博客整理成的书,很详细,很典型的一本教程,介绍了OpenCv中相对前沿的知识。我翻看过这本教程,中规中矩,里面的代码通俗易懂,尤其适合初学者。当然大家同样要注意不要犯了死读书的毛病,只看它的功能,敲代码练手感即可,不要试图记忆API函数。重要的工具用多了自然会记住,不重要的工具记住了也没用。

  这里推荐的第二本书是《图像识别与项目实践――VC++、MATLAB技术实现》,这本书是一本偏向于工程应用的书,我之所以推荐它是因为它给出了很多有新意、能运行的代码。其中里面有一个项目让我印象很深,是一个车牌检测的实例。简单描述一下:由于车牌中的字符数是固定的,因此它通过判断横向区域笔画的跳变数以及笔画宽度来定位车牌区域。这个想法让人耳目一新,并且它还给出了详细代码,我也亲身试验过,效果还不错。

  这里同样再强调一下,就是一定要尽早入手写程序,建立自信和成就感。我当时学OpenCv正好用它开发了一个人脸性别识别的系统,是一个本科大学生创新计划的需求,效果还可以。

  3、Python

  Python在今年12月份的编程语言排行榜中名列第5,增长迅速。可以说Python已经逐渐成为当下脚本语言的新标准。Python在图像处理算法方面除了其自身简洁的编程优势外,还得益于两个重要的Python类库——Numpy和Theano。

  Numpy是Python的线性代数库,对于矩阵运算能提供很好的支持,并且能够在此基础上进行很多机器学习相关算法的开发仿真,这里推荐一本受到大家广泛认可的书《机器学习实战》,我最近也正在看这本书,里面对好多机器学习领域的经典算法,小到KNN,大到SVM,都给出了详细的介绍以及代码实现(Python版)。Theano是Python的机器学习库,能够方便的实现深度学习(例如卷积神经网络CNN)算法,网上很多对于DeepID算法的复现都是用的这个库。

  人觉得单从图像处理的角度评价的话,Python并没有前面两个工具(Matlab和OpenCv)应用广泛,不过作为通用的脚本语言,我觉得每个程序员都应该去了解了解它,毕竟俗话说没有烂的编程语言,只有烂程序员。我在学Python时第一个自己写的程序就是微信打飞机的小程序,在我的博客中有详细的教程,虽然是参照小甲鱼的《零基础入门学习Python》视频教程写的,但还是蛮有成就感的。

  三、图像处理研究方法

  我觉得,图像处理研究主要可以分为三个部分:基础概念、基本思想、算法研究。

  1、基础概念

  所谓基础概念,就是图像处理里最基本的知识,比如什么是图像?什么是像素?什么是彩色图像等等。没有一个明确的界限来划定什么是基础概念什么是高级知识,因人而异。了解图像处理的基础知识,有一本书是必读的,就是冈萨雷斯编写的、阮秋琦翻译的《数字图像处理》。这本书已经作为图像处理领域的经典教材使用了三十多年,我自己也把这本书看了好几遍,每一遍都会有新的体会。我觉得每一个搞图像的都应该熟读这本书。书中除了有几章内容在讲小波变换、模式识别等相对抽象的内容外,其他内容相对都是很基础的,本科生水平就能看懂。而且我建议要尽早看这本书,如果是研究生尽量在进入课题之前就看一遍,因为这样一本经典的书在进入课题之后可能就没时间看了,以后也顶多是查阅而已。我当初就是在大四的寒假看完了这本书,以后在图像入门的过程中就显得轻松很多。看完这本书,哪怕是只看前几章,明白了什么是图像(二维或者三维矩阵)、什么是像素、彩色图和灰度图、颜色空间、图像滤波、图像噪声、图像频域变换等概念,将来在进行更深一步的研究的话就会方便很多了。

  2、基本思想

  刚开始想把这部分内容命名为“基本算法”,意在介绍图像处理中的一些基本算法,后来仔细想想决定不这么写,因为图像处理是一个非常大的概念,图像处理不等于人脸识别,也不等于模式识别,直接介绍诸如图像处理基本算法之类的内容很容易写成空话,没有什么实际意义。读者有兴趣的话可以直接谷歌百度“图像处理十大经典算法”,上面有我想说的内容。

  万变不离其宗,算法是死的,重在思想。举个例子,我个人是主攻模式识别方向,在这个方向判断一个学生是否入门有一个非常简单的方法,就是“如果你能把图像很自然的想象成高维空间中的一个点”,那就说明在模式识别方面入门了,可以对图像进行分类了。当然标准不是唯一,在其他领域如目标检测也会有其他的判断标准,总之我们要对图像进行处理,那么图像就不再只是图像,它可能会演变成各种不同形式的概念,可能是点,可能是面,还可能是一个坐标空间。在目标跟踪的经典算法粒子滤波中,将一个个的小图像块看做一个个粒子;在子空间理论中,将一系列图像放在一起构建一个成分主空间(例如主成分分析PCA算法等等。,我不会详细介绍这些算法,说多了就显得抽象老套,但我要说的是我们一定要把图像本身理解好,它是一个图像,是一个矩阵,是一个信息的容器,是一种数据的表现形式,图像不一定都必须在视觉上有意义(比如频域的图像)。

  总之图像处理的基本思想还是要立足于图像本身,要深度到图像内部结构中,思维要灵活。我当时做本科毕设时,怎么也不知道图像和高维空间中的点之间有什么对应关系,后来总算有一天,突然就明白了,这也就是所谓的量变产生质变。总之一定要多想,多总结,主动去钻研,才能够真正领悟一些东西。最基本的东西往往蕴藏着深奥的道理,无论你现在多牛多厉害,都不能放掉最本源的东西。多想想图像是什么,有什么本质属性,你可能无法得到准确的答案,但肯定能得到一些有用的感悟(有点像哲学问题了)。

  3、算法研究

  算法研究应该是图像处理的核心工作,尤其是各大高校的博士硕士。这里我并不想谈那些高大上的算法,我更想说的是一些算法研究的一些基础的东西,比如说一些基础课程,比如说矩阵运算。

  研究图像处理的算法,离不开数学。在这里我建议图像处理方面的硕士一定要上两门课:《泛函分析》以及《最优化算法》,有的学校已经将这两门课列为了研究生阶段的必修课程。这两门可可以说是图像处理(至少是模式识别)的基础。我当初没上过最优化算法,但后来也自己补上了,不然真的是寸步难行。至于泛函我当时听课的时候也不是很懂,但是在之后的研究过程中发现很多图像处理的基本知识基本理论都和泛函分析中枯燥的定理如出一辙,没办法,有的东西本身就是枯燥的干货,学着费力,缺它不行。

  其次我想说的是矩阵运算。图像就是矩阵,图像处理就是矩阵运算。大家为什么都喜欢用Matlab,就是因为它的矩阵运算能力实在是太强大,在Matlab的世界中任何变量都是矩阵。同样OpenCv之所以能流行,不仅仅是因为它良好的封装性,也是因为它的矩阵格式,它定义了Mat基础类,允许你对矩阵进行各种操作。Python也不例外,它的Numpy就是一个专门的线性代数库。

  真正在图像编程过程中,那些看着高大上的API函数归根到底都是工具,查查手册就能找到,真正核心还是在算法,算法是由公式编写的,公式的单元是变量,而图像届的变量就是矩阵。所以,熟练去操作矩阵,求秩、求逆、最小二乘,求协方差,都是家常便饭。所以,如果你有幸能上《矩阵分析》这门课,一定要把它看懂,那里面都是干货。

  四、小结

  总之,图像处理就是一个典型的门槛低、厅堂深的领域。不需要太多基础,学过线性代数,会一点编程就够了;但是那些算法却深不可测,是个消耗功夫的活儿。在写这篇教程时我说的很直白,就像和大家对话一样,想到什么说什么。在最后我想说两句题外话,就是不仅仅针对图像处理,对于其他新技术的入门学习也是一样,尽快迈出第一步,尽快去建立自信和成就感,让自己有勇气走下去,然后缺什么补什么就行了。我觉得真正让人望而却步的往往不是技术本身,而是我们对自身的不自信。唯有果断开工,才能战胜心魔。


如果觉得这篇文章对您有所启发,欢迎关注我的公众号,我会尽可能积极和大家交流,谢谢。


2016-02-26 17:48:13 baimafujinji 阅读数 69723

什么是数字图像处理?历史、以及它所研究的内容。

 

说起图像处理,你会想到什么?你是否真的了解这个领域所研究的内容。纵向来说,数字图像处理研究的历史相当悠久;横向来说,数字图像处理研究的话题相当广泛。

数字图像处理的历史可以追溯到近百年以前,大约在1920年的时候,图像首次通过海底电缆从英国伦敦传送到美国纽约。图像处理的首次应用是为了改善伦敦和纽约之间海底电缆发送的图片质量,那时就应用了图像编码,被编码后的图像通过海底电缆传送至目的地,再通过特殊设备进行输出。这是一次历史性的进步,传送一幅图片的时间从原来的一个多星期减少到了3小时。

1950年,美国的麻省理工学院制造出了第一台配有图形显示器的电子计算机——旋风I号(Whirlwind I)。旋风I号的显示器使用一个类似于示波器的阴极射线管(Cathode Ray Tube,CRT)来显示一些简单的图形。1958年美国Calcomp公司研制出了滚筒式绘图仪,GerBer公司把数控机床发展成为平板式绘图仪。在这一时期,电子计算机都主要应用于科学计算,而为这些计算机配置的图形设备也仅仅是作为一种简单的输出设备。

随着计算机技术的进步,数字图像处理技术也得到了很大的发展。1962年,当时还在麻省理工学院攻读博士学位的伊凡·苏泽兰(Ivan Sutherland)成功开发了具有划时代意义的“画板”(Sketchpad)程式。而这正是有史以来第一个交互式绘图系统,同时这也是交互式电脑绘图的开端。从此计算机和图形图像被更加紧密地联系到了一起。鉴于伊凡·苏泽兰为计算机图形学创立所做出的杰出贡献,他于1988年被授予计算机领域最高奖——图灵奖。

1964年,美国加利福尼亚的喷气推进实验室用计算机对“旅行者七号”太空船发回的大批月球照片进行处理,以校正航天器上摄影机中各种类型的图像畸变,收到了明显的效果。在后来的宇航空间技术中,数字图像处理技术都发挥了巨大的作用。

到了20世纪60年代末期,数字图像处理已经形成了比较完善的学科体系,这套理论在20世纪70年代发展得十分迅速,并开始应用于医学影像和天文学等领域。1972年,美国物理学家阿伦·马克利奥德·柯麦科(Allan MacLeodCormack)和英国电机工程师戈弗雷·纽博尔德·豪恩斯弗尔德(Godfrey Newbold Housfield)发明了轴向断层术,并将其用于头颅诊断。世界第一台X射线计算机轴向断层摄影装置由EMI公司研制成功,这也就是人们通常所说的CT(Computer Tomograph)。CT可通过一些算法用感知到的数据去重建通过物体的“切片”图像。这些图像组成了物体内部的再现图像,也就是根据人的头部截面的投影,经计算机处理来进行图像重建。鉴于CT对于医学诊断技术的发展所起到的巨大推动作用,柯麦科和豪恩斯弗尔德于1979年获得了诺贝尔生理或医学奖。

随后在2003年,诺贝尔生理或医学奖的殊荣再次授予了两位在医疗影像设备研究方面做出杰出贡献的科学家——美国化学家保罗·劳特伯尔(Paul Lauterbur)和英国物理学家彼得·曼斯菲尔(Peter Mansfield)。两位获奖者在利用磁共振成像(Magnetic Resonance Imaging,MRI)显示不同结构方面分别取得了开创性成就。瑞典卡罗林斯卡医学院称,这两位科学家在MRI领域的开创性工作,代表了医学诊疗和研究的重大突破。而事实上,核磁共振的成功同样也离不开数字图像处理方面的发展。即使在今天,诸如MRI图像降噪等问题依然是数字图像处理领域的热门研究方向。

说到数字图像的发展历程,还有一项至关重要的成果不得不提,那就是电荷耦合元件(Charge-coupled Device,CCD)。CCD最初是由美国贝尔实验室的科学家维拉德·波义耳(Willard Sterling Boyle)和乔治·史密斯(George Elwood Smith)于1969年发明的。CCD的作用就像胶片一样,它能够把光学影像转化为数字信号。今天人们所广泛使用的数码照相机、数码摄影机和扫描仪都是以CCD为基础发展而来的。换句话说,我们现在所研究的数字图像主要也都是通过CCD设备获取的。由于波义耳和史密斯在CCD研发上所做出的巨大贡献,他们两人共同荣获了2009年度的诺贝尔物理学奖。

数字图像处理在今天是非常热门的技术之一,生活中无处不存在着它的影子,可以说它是一种每时每刻都在改变着人类生活的技术。但长久以来,很多人对数字图像处理存在着较大的曲解,人们总是不自觉地将图像处理和Photoshop联系在一起。大名鼎鼎的Photoshop无疑是当前使用最为广泛的图像处理工具。类似的软件还有Corel公司生产的CorelDRAW等软件。

尽管Photoshop是一款非常优秀的图像处理软件,但它的存在并不代表数字图像处理的全部理论与方法。它所具有的功能仅仅是数字图像处理中的一部分。总的来说,数字图像处理研究的内容主要包括如下几个方面:

  • 1)图像获取和输出
  • 2)图像编码和压缩
  • 3)图像增强与复原
  • 4)图像的频域变换
  • 5)图像的信息安全
  • 6)图像的区域分割
  • 7)图像目标的识别
  • 8)图像的几何变换

但图像处理的研究内容,又不仅限于上述内容!所以说图像处理的研究话题是相当宽泛的。那现在图像处理都应用在哪些领域呢?或许我们可能熟知的例子有(当然,你应该还能举出更多例子):

  • 1)一些专业图像处理软件:Photoshop、CorelDRAW……
  • 2)一些手机APP应用:美图秀秀、玩图……
  • 3)一些医学图像处理应用:MRI、彩超图像处理……
  • 4)一些制造业上的应用:元器件检测、瑕疵检测……
  • 5)一些摄像头、相机上的应用:夜间照片的质量改善……
  • 6)一些电影工业上是应用:换背景、电影特技……

 

什么样的人会去学(或者需要学)图像处理?

 

1)如果你是我上述那些应用领域的从业者,你当然需要掌握图像方面的理论和技术;2)相关专业的研究人员、大专院校的博士生、研究生。

所谓相关专业又是指什么呢?这个答案也可能相当宽泛,例如(但不仅限于此):Computer Science, Software Engineering, Electronic Engineering, Biomedical Engineering, Automation, Control, Applied Mathematics……

 

如何学好图像处理——我的一些箴言

 

1)对于初级入门者

 

一个扎实的基础和对于图像处理理论的完整的、系统的整体认识对于后续的深入研究和实践应用具有非常非常重要的意义。

我经常喜欢拿武侠小说《天龙八部》中的一段情节来向读者说明此中的道理,相信读者对这部曾经被多次搬上银幕的金庸作品已经耳熟能详了。书中讲到有个名叫鸠摩智的番僧一心想练就绝世武学,而且他也算是个相当勤奋的人了。但是,他错就错在太过于急功近利,甚至使用道家的小无相功来催动少林绝技。看上去威力无比,而且可以在短时间内“速成”,但实则后患无穷。最终鸠摩智走火入魔,前功尽废,方才大彻大悟。这个故事其实就告诉我们打牢基础是非常重要的,特别是要取得更长足的发展,就更是要对基本原理刨根问底,力求甚解,从而做到庖丁解牛,游刃有余。

一些看似高深的算法往往是许多基础算法的组合提升。例如,令很多人望而却步的SIFT特征构建过程中,就用到了图像金字塔、直方图、高斯滤波这些非常非常基础的内容。但是,它所涉及的基础技术显然有好几个,如果缺乏对图像处理理论的系统认识,你可能会感觉事倍功半。因为所有的地方好像都是沟沟坎坎。

关于课程——

在这个阶段其实对于数学的要求并不高,你甚至可以从一些感性的角度去形象化的理解图像处理中很多内容(但不包括频域处理方面的内容)。具体到学习的建议,如果有条件(例如你还在高校里读书)你最好能选一门图像处理方面的课程,系统地完整的地去学习一下。这显然是入门的最好办法。如此一来,在建立一个完整的、系统的认知上相当有帮助。如果你没办法在学校里上一门这样的课,网上的一些公开课也可以试试。但现在中文MOOC上还没有这方面的优质课程推荐。英文的课程则有很多,例如美国加州伦斯勒理工学院Rich教授的数字图像处理公开课——https://www.youtube.com/channel/UCaiJlKxXamoODQtlx486qJA?spfreload=10。

关于教材——

显然,只听课其实还不太够,如果能一并读一本书就最好了。其实不用参考很多书,只要一本,你能从头读到尾就很好了。如果你没有条件去上一门课,那读一本来完整的自学一下就更有必要了。这个阶段,去网上到处找博客、看帖子是不行的。因为你特别需要在这个阶段对这门学问建立一个系统的完整的知识体系。东一块、西一块的胡拼乱凑无疑是坑你自己,你的知识体系就像一个气泡,可能看起来很大,但是又脆弱的不堪一击。

现在很多学校采用冈萨雷斯的《数字图像处理》一书作为教材。这是一本非常非常经典的著作。但是我必须要提醒读者:

1)这是一本专门为Electronic Engineering专业学生所写的书。它需要有信号与系统、数字信号处理这两门课作为基础。如果你没有这两门课的基础,你读这本书要么是看热闹,要么就是看不懂。

下面是冈书中的一张插图。对于EE的学生来说,这当然不是问题。但是如果没有我说的那两门课的基础,其实你很难把握其中的精髓。H和h,一个大小一个小写,冈书中有的地方用H,有的地方用h,这都是有很深刻用意的。原作者并没有特别说明它们二者的区别,因为他已经默认你应该知道二者是不同的。事实上,它们一个表示频域信号,一个表示时域信号,这也导致有时候运算是卷积,有时候运算是乘法(当然这跟卷积定理有关)。所以我并不太建议那些没有这方面基础的学生在自学的时候读这本书。

 

2)冈萨雷斯教授的《数字图像处理》第一版是在1977年出版的,到现在已经快40年了;现在国内广泛使用的第二版是2002年出版的(第三版是2007年但是其实二者差异并不大),到现在也有20年左右的时间了。事实上,冈萨雷斯教授退休也有快30年了。所以这本书的内容已经偏于陈旧。数字图像处理这个领域的发展绝对是日新月异,突飞猛进的。特别在最近二三十年里,很多新思路,新方法不断涌现。如果你看了我前面推荐的Rich教授的公开课(这也是当前美国大学正在教学的内容),你一下子就会发现,原来我们的教育还停留在改革开放之前外国的水平上。这其实特别可怕。所以我觉得冈萨雷斯教授的《数字图像处理》作为学习过程中的一个补充还是不错的,但是如果把它作为主参考,那真的就是:国外都洋枪洋炮了,我们还在大刀长矛。

 

那么现在问题来了,对于图像处理学习者而言到底看什么书好呢?我的意见是你可以选择下面两本书中的任何一本《数字图像处理原理与实践(Matlab版)》,以及《数字图像处理:技术详解与Visual C++实践》,当然选择的标准之一就是到底你更擅长使用MATLAB还是C++。

   

 

 

 

2)对于中级水平者

 

纸上得来终觉浅,绝知此事要躬行。对于一个具有一定基础的,想更进一步的中级水平的人来说,这个阶段最重要的就是增强动手实践的能力。

还是说《天龙八部》里面的一个角色——口述武功、叹为观止的王语嫣。王语嫣的脑袋里都是武功秘籍,但问题是她从来都没练过一招一式。结果是,然并卵。所以光说不练肯定不灵啊。特别是,如果你将来想从事这个行业,结果一点代码都不会写,那几乎是不可想象的。学习阶段,最常被用来进行算法开发的工具是Matlab和OpenCV。你可以把这两个东西都理解为一个相当完善的库。当然,在工业中C++用得更多,所以Matlab的应用还是很有限的。前面我们讲到,图像处理研究内容其实包括:图像的获取和编解码,但使用Matlab和OpenCV就会掩盖这部分内容的细节。你当然永远不会知道,JPEG文件到底是如何被解码的。

如果你的应用永远都不会涉及这些话题,那么你一直用Matlab和OpenCV当然无所谓。例如你的研究领域是SIFT、SURF这种特征匹配,可以不必理会编解码方面的内容。但是如果你的研究话题是降噪或者压缩,可能你就绕不开这些内容。最开始学的时候,如果能把这部分内容也自己写写,可能会加深你的理解。以后做高级应用开发时,再调用那些库。所以具体用什么,要不要自己写,是要视你所处的阶段和自己的实际情况而定的。以我个人的经验,在我自学的时候,我就动手写了Magic House,我觉得这个过程为我奠定了一个非常夯实的基础,对于我后续的深入研究很有帮助。

 

下面这个文中,我给出了一些这方面的资源,代码多多,很值得参考学习:图像处理与机器视觉网络资源收罗

http://blog.csdn.net/baimafujinji/article/details/32332079

 

3)对于高级进阶者

 

到了这个程度的读者,编程实现之类的基本功应该不在话下。但是要往深,往高去学习、研究和开发图像处理应用,你最需要的内容就变成了数学。这个是拦在很多处于这个阶段的人面前的一大难题。如果你的专业是应用数学,当然你不会感觉有问题。但如果是其他专业背景的人就会越发感觉痛苦。

如果你的图像处理是不涉及机器学习内容的,例如用Poisson方程来做图像融合,那你就要有PDE数值解方面的知识;如果你要研究KAZE特征,你就必须要知道AOS方面的内容。如果你研究TV降噪,你又要知道泛函分析中的BV空间内容……这些词你可能很多都没听过。总的来说,这块需要的内容包括:复变函数、泛函分析、偏微分方程、变分法、数学物理方法……

如果你要涉足机器视觉方法的内容,一些机器学习和数据挖掘方法的内容就不可或缺。而这部分内容同样需要很强大的数学基础,例如最大似然方法、梯度下降法、欧拉-拉格朗日方程、最小二乘估计、凸函数与詹森不等式……

当然,走到这一步,你也已经脱胎换骨,从小白到大神啦!路漫漫其修远兮,吾将上下而求索。

 

(全文完)

 

 

2014-09-02 13:50:53 gotosola 阅读数 3366
图形图像处理-之-任意角度的高质量的快速的图像旋转
 转自:http://blog.chinaunix.net/uid-20306372-id-1707267.html
 
上篇 纯软件的任意角度的快速旋转
                                HouSisong@GMail.com   2007.04.26

tag:图像旋转,任意角度,图像缩放,速度优化,定点数优化,近邻取样插值,二次线性插值,
   三次卷积插值,MipMap链,三次线性插值,MMX\SSE优化,CPU缓存优化,AlphaBlend,颜色混合,并行

摘要:首先给出一个基本的图像旋转算法,然后一步一步的优化其速度和旋转质量,打破不能软件旋转的神话!

任意角度的高质量的快速的图像旋转 全文 分为:
     上篇 纯软件的任意角度的快速旋转
     中篇 高质量的旋转
     下篇 补充话题(完整AlphaBlend旋转、旋转函数并行化、针对大图片的预读缓冲区优化)

(2007.06.22 优化PicRotary3加快13.6%,优化PicRotarySSE加快6.1%,
            尝试了一下使用SSE2的写缓冲优化MOVNTI)
(2007.04.29 修正一个TRotaryClipData.find_begin的bug)
(2007.05.16 更换测试用电脑和编译器,为了保持测试数据一致和可对比性,更新了测试数据)

正文:
  为了便于讨论,这里只处理32bit的ARGB颜色;
  代码使用C++;涉及到汇编优化的时候假定为x86平台;使用的编译器为vc2005;
  为了代码的可读性,没有加入异常处理代码;
   测试使用的CPU为赛扬2G(新的测试平台的CPU为AMD64x2 4200+(2.37G),测试时使用的单线程执行);
  (一些基础代码和插值原理的详细说明参见作者的《图形图像处理-之-高质量的快速的图像缩放》系列文章)


速度测试说明:
  只测试内存数据到内存数据的缩放
  测试图片都是800*600旋转到1004*1004,测试成绩取各个旋转角度的平均速度值; fps表示每秒钟的帧数,值越大表示函数越快


A:旋转原理和旋转公式:
  推导旋转公式:
             
                       旋转示意图
   有:  tg(b)=y/x                             ----(1) 
         tg(a+b)=y'/x'                         ----(2) 
         x*x + y*y = x'*x' + y'*y'             ----(3)
   有公式:tg(a+b) = ( tg(a)+tg(b) ) / ( 1-tg(a)*tg(b) )  ----(4)
     把(1)代入(4)从而消除参数b;
     tg(a)+y/x = y'/x'*( 1-tg(a)*y/x )                ----(5)
     由(5)可以得x'=y'*(x-y*tg(a))/( x*tg(a)+y )       ----(6)
   把(6)代入(3)从而消除参数x',化简后求得:
     y'=x*sin(a)+y*cos(a);                     ----(7)
   把(7)代入(6),有:
     x'=x*cos(a)-y*sin(a);                     ----(8)

  OK,旋转公式有了,那么来看看在图片旋转中的应用;
  假设对图片上任意点(x,y),绕一个坐标点(rx0,ry0)逆时针旋转RotaryAngle角度后的新的坐标设为(x', y'),有公式: 
  (x平移rx0,y平移ry0,角度a对应-RotaryAngle , 带入方程(7)、(8)后有: )  
  x'= (x - rx0)*cos(RotaryAngle) + (y - ry0)*sin(RotaryAngle) + rx0 ;
  y'=-(x - rx0)*sin(RotaryAngle) + (y - ry0)*cos(RotaryAngle) + ry0 ;

那么,根据新的坐标点求源坐标点的公式为:
  x=(x'- rx0)*cos(RotaryAngle) - (y'- ry0)*sin(RotaryAngle) + rx0 ;
  y=(x'- rx0)*sin(RotaryAngle) + (y'- ry0)*cos(RotaryAngle) + ry0 ;

旋转的时候还可以顺便加入x轴和y轴的缩放和平移,而不影响速度,那么完整的公式为:            
  x=(x'- move_x-rx0)/ZoomX*cos(RotaryAngle) - (y'- move_y-ry0)/ZoomY*sin(RotaryAngle) + rx0 ;
  y=(x'- move_x-rx0)/ZoomX*sin(RotaryAngle) + (y'- move_y-ry0)/ZoomY*cos(RotaryAngle) + ry0 ;
  其中: RotaryAngle为逆时针旋转的角度;
         ZoomX,ZoomY为x轴y轴的缩放系数(支持负的系数,相当于图像翻转);
         move_x,move_y为x轴y轴的平移量;

 一些颜色和图片的数据定义:

#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];
}
//判断一个点是否在图片中

inline bool PixelsIsInPic(const TPicRegion& pic,const long x,const long y)
{
    
return ( (x>=0)&&(x<pic.width) && (y>=0)&&(y<
pic.height) );
}

 

 B:一个简单的浮点实现版本

//////////////////////////////////////////////////////////////////////////////////////////////////////
//函数假设以原图片的中心点坐标为旋转和缩放的中心

void PicRotary0(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double move_y)
{
    
if ( (fabs(ZoomX*Src.width)<1.0e-4|| (fabs(ZoomY*Src.height)<1.0e-4) ) return//太小的缩放比例认为已经不可见

    double rx0=Src.width*0.5;  //(rx0,ry0)为旋转中心 
    double ry0=Src.height*0.5
    
for (long y=0;y<Dst.height;++
y)
    {
        
for (long x=0;x<Dst.width;++
x)
        {
            
long srcx=(long)((x- move_x-rx0)/ZoomX*cos(RotaryAngle) - (y- move_y-ry0)/ZoomY*sin(RotaryAngle) +
 rx0) ;
            
long srcy=(long)((x- move_x-rx0)/ZoomX*sin(RotaryAngle) + (y- move_y-ry0)/ZoomY*cos(RotaryAngle) +
 ry0) ;
            
if
 (PixelsIsInPic(Src,srcx,srcy))
                Pixels(Dst,x,y)
=
Pixels(Src,srcx,srcy);
        }
    }
}

(调用方法比如:
   PicRotary0(ppicDst,ppicSrc,PI/6,0.9,0.9,(dst_wh-ppicSrc.width)*0.5,(dst_wh-ppicSrc.height)*0.5);
   //作用:将图片ppicSrc按0.9的缩放比例旋转PI/6幅度后绘制到图片ppicDst的中心
)

//注:测试图片都是800*600的图片旋转到1004*1004的图片中心 测试成绩取各个旋转角度的平均速度值
////////////////////////////////////////////////////////////////////////////////
//速度测试:                  
//==============================================================================
// PicRotary0              34.9 fps
////////////////////////////////////////////////////////////////////////////////  

旋转结果图示(小图):

   
                 30度                         60度                      90度

   
                 120度                       150度                      180度

  
                 210度                            240度                          270度

    
                 300度                        330度                       360度


C:优化循环内部,化简系数
  1.sin和cos函数是很慢的计算函数,可以在循环前预先计算好sin(RotaryAngle)和cos(RotaryAngle)的值:
    double sinA=sin(RotaryAngle);
    double cosA=cos(RotaryAngle);
  2.可以将除以ZoomX、ZoomY改成乘法,预先计算出倒数:
    double rZoomX=1.0/ZoomX;
    double rZoomY=1.0/ZoomY;
  3.优化内部的旋转公式,将能够预先计算的部分提到循环外(即:拆解公式):
     原:  long srcx=(long)((x- move_x-rx0)/ZoomX*cos(RotaryAngle) - (y- move_y-ry0)/ZoomY*sin(RotaryAngle) + rx0) ;
           long srcy=(long)((x- move_x-rx0)/ZoomX*sin(RotaryAngle) + (y- move_y-ry0)/ZoomY*cos(RotaryAngle) + ry0) ;
     变形为:
           long srcx=(long)( Ax*x + Bx*y +Cx ) ;
           long srcy=(long)( Ay*x + By*y +Cy ) ;
      其中: Ax=(rZoomX*cosA); Bx=(-rZoomY*sinA); Cx=(-(rx0+move_x)*rZoomX*cosA+(ry0+move_y)*rZoomY*sinA+rx0);
            Ay=(rZoomX*sinA); By=(rZoomY*cosA);  Cy=(-(rx0+move_x)*rZoomX*sinA-(ry0+move_y)*rZoomY*cosA+ry0); 
      (提示: Ax,Bx,Cx,Ay,By,Cy都可以在旋转之前预先计算出来)
  改进后的函数为:

void PicRotary1(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double move_y)
{
    
if ( (fabs(ZoomX*Src.width)<1.0e-4|| (fabs(ZoomY*Src.height)<1.0e-4) ) return//太小的缩放比例认为已经不可见

    double rZoomX=1.0/ZoomX;
    
double rZoomY=1.0/
ZoomY;
    
double sinA=
sin(RotaryAngle);
    
double cosA=
cos(RotaryAngle);
    
double Ax=(rZoomX*
cosA); 
    
double Ay=(rZoomX*
sinA); 
    
double Bx=(-rZoomY*
sinA); 
    
double By=(rZoomY*
cosA); 
    
double rx0=Src.width*0.5;  //(rx0,ry0)为旋转中心 

    double ry0=Src.height*0.5
    
double Cx=(-(rx0+move_x)*rZoomX*cosA+(ry0+move_y)*rZoomY*sinA+
rx0);
    
double Cy=(-(rx0+move_x)*rZoomX*sinA-(ry0+move_y)*rZoomY*cosA+
ry0); 

    TARGB32
* pDstLine=
Dst.pdata;
    
double srcx0_f=
(Cx);
    
double srcy0_f=
(Cy);
    
for (long y=0;y<Dst.height;++
y)
    {
        
double srcx_f=
srcx0_f;
        
double srcy_f=
srcy0_f;
        
for (long x=0;x<Dst.width;++
x)
        {
            
long srcx=(long
)(srcx_f);
            
long srcy=(long
)(srcy_f);
            
if
 (PixelsIsInPic(Src,srcx,srcy))
                pDstLine[x]
=
Pixels(Src,srcx,srcy);
            srcx_f
+=
Ax;
            srcy_f
+=
Ay;
        }
        srcx0_f
+=
Bx;
        srcy0_f
+=
By;
        ((TUInt8
*&)pDstLine)+=
Dst.byte_width;
    }
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:                  
//==============================================================================
// PicRotary1               62.0 fps
//////////////////////////////////////////////////////////////////////////////// 

( 在AMD64x2 4200+和VC2005编译下PicRotary1(51.8fps)比PicRotary0(27.1fps)快90%;
  在AMD64x2 4200+和VC6编译下PicRotary1(20.3fps)比PicRotary0(16.1fps)快26%;
  以前在赛扬2G和VC6编译下PicRotary1(8.4fps)反而比PicRotary0(12.7fps)慢50%! :( ) 

D:更深入的优化、定点数优化
  (浮点数到整数的转化也是应该优化的一个地方,这里不再处理,可以参见
    <图形图像处理-之-高质量的快速的图像缩放 上篇 近邻取样插值和其速度优化>中的PicZoom3_float函数)
  1.优化除法: 
    原: double rZoomX=1.0/ZoomX;
         double rZoomY=1.0/ZoomY;
    改写为(优化掉了一次除法):
         double tmprZoomXY=1.0/(ZoomX*ZoomY);  
         double rZoomX=tmprZoomXY*ZoomY;
         double rZoomY=tmprZoomXY*ZoomX;
  2.x86的浮点计算单元(FPU)有一条指令"fsincos"可以同时计算出sin和cos值
    //定义SinCos函数: 同时计算sin(Angle)和cos(Angle)的内嵌x86汇编函数
    void __declspec(naked) __stdcall SinCos(const double Angle,double& sina,double& cosa) 
    { 
        asm
        {
            fld  qword ptr [esp+4]//Angle   
            mov  eax,[esp+12]//&sina
            mov  edx,[esp+16]//&cosa
            fsincos   
            fstp qword ptr [edx]   
            fstp qword ptr [eax]  
            ret 16
        }
    } 
  3.用定点数代替浮点计算

void PicRotary2(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double move_y)
{
    
if ( (fabs(ZoomX*Src.width)<1.0e-4|| (fabs(ZoomY*Src.height)<1.0e-4) ) return//太小的缩放比例认为已经不可见

    double tmprZoomXY=1.0/(ZoomX*ZoomY);  
    
double rZoomX=tmprZoomXY*
ZoomY;
    
double rZoomY=tmprZoomXY*
ZoomX;
    
double
 sinA,cosA;
    SinCos(RotaryAngle,sinA,cosA);
    
long Ax_16=(long)(rZoomX*cosA*(1<<16
)); 
    
long Ay_16=(long)(rZoomX*sinA*(1<<16
)); 
    
long Bx_16=(long)(-rZoomY*sinA*(1<<16
)); 
    
long By_16=(long)(rZoomY*cosA*(1<<16
)); 
    
double rx0=Src.width*0.5;  //(rx0,ry0)为旋转中心 

    double ry0=Src.height*0.5
    
long Cx_16=(long)((-(rx0+move_x)*rZoomX*cosA+(ry0+move_y)*rZoomY*sinA+rx0)*(1<<16
));
    
long Cy_16=(long)((-(rx0+move_x)*rZoomX*sinA-(ry0+move_y)*rZoomY*cosA+ry0)*(1<<16
)); 

    TARGB32
* pDstLine=
Dst.pdata;
    
long srcx0_16=
(Cx_16);
    
long srcy0_16=
(Cy_16);
    
for (long y=0;y<Dst.height;++
y)
    {
        
long srcx_16=
srcx0_16;
        
long srcy_16=
srcy0_16;
        
for (long x=0;x<Dst.width;++
x)
        {
            
long srcx=(srcx_16>>16
);
            
long srcy=(srcy_16>>16
);
            
if
 (PixelsIsInPic(Src,srcx,srcy))
                pDstLine[x]
=
Pixels(Src,srcx,srcy);
            srcx_16
+=
Ax_16;
            srcy_16
+=
Ay_16;
        }
        srcx0_16
+=
Bx_16;
        srcy0_16
+=
By_16;
        ((TUInt8
*&)pDstLine)+=
Dst.byte_width;
    }
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicRotary2             134.2 fps
//////////////////////////////////////////////////////////////////////////////// 


E:优化内部循环的判断函数PixelsIsInPic,将判断展开到内部循环之外,跳过不需要处理的目标像素;
   TRotaryClipData类用于寻找旋转需要处理的边界范围;算法思路是首先寻找原图片中心点对应的;
那条扫描线,然后依次向上和向下寻找边界;  如果想要更快速的边界寻找算法,可以考虑利用像素的直线
  绘制原理来自动生成边界(有机会的时候再来实现它);

                   

                       边界寻找算法图示

struct TRotaryClipData{
public
:
    
long
 src_width;
    
long
 src_height;
    
long
 dst_width;
    
long
 dst_height;
    
long
 Ax_16; 
    
long
 Ay_16; 
    
long
 Bx_16; 
    
long
 By_16; 
    
long
 Cx_16;
    
long
 Cy_16; 
public
:
    
long
 out_dst_up_y;
    
long
 out_dst_down_y;

    
long
 out_src_x0_16;
    
long
 out_src_y0_16;
private
:
    
long
 out_dst_up_x0;
    
long
 out_dst_up_x1;
    
long
 out_dst_down_x0;
    
long
 out_dst_down_x1;
public
:
    inline 
long get_up_x0(){ if (out_dst_up_x0<0return 0;  else return
 out_dst_up_x0; }
    inline 
long get_up_x1(){ if (out_dst_up_x1>=dst_width) return dst_width;  else return
 out_dst_up_x1; }
    inline 
long get_down_x0(){ if (out_dst_down_x0<0return 0;  else return
 out_dst_down_x0; }
    inline 
long get_down_x1(){ if (out_dst_down_x1>=dst_width) return dst_width;  else return
 out_dst_down_x1; }

    inline 
bool is_in_src(long src_x_16,long
 src_y_16)
    {
         
return ( ( (src_x_16>=0)&&((src_x_16>>16)<
src_width) )
               
&& ( (src_y_16>=0)&&((src_y_16>>16)<
src_height) ) );
    }
    
void find_begin_in(long dst_y,long& out_dst_x,long& src_x_16,long&
 src_y_16)
    {
        src_x_16
-=
Ax_16;
        src_y_16
-=
Ay_16;
        
while
 (is_in_src(src_x_16,src_y_16))
        {
            
--
out_dst_x;
            src_x_16
-=
Ax_16;
            src_y_16
-=
Ay_16;
        }
        src_x_16
+=
Ax_16;
        src_y_16
+=
Ay_16;
    }
    
bool find_begin(long dst_y,long& out_dst_x0,long
 dst_x1)
    {
        
long test_dst_x0=out_dst_x0-1
;
        
long src_x_16=Ax_16*test_dst_x0 + Bx_16*dst_y +
 Cx_16;
        
long src_y_16=Ay_16*test_dst_x0 + By_16*dst_y +
 Cy_16;
        
for (long i=test_dst_x0;i<=dst_x1;++
i)
        {
            
if
 (is_in_src(src_x_16,src_y_16))
            {
                out_dst_x0
=
i;
                
if (i==
test_dst_x0)
                    find_begin_in(dst_y,out_dst_x0,src_x_16,src_y_16);
                
if (out_dst_x0<0
)
                {
                    src_x_16
-=(Ax_16*
out_dst_x0);
                    src_y_16
-=(Ay_16*
out_dst_x0);
                }
                out_src_x0_16
=
src_x_16;
                out_src_y0_16
=
src_y_16;
                
return true
;
            }
            
else

            {
                src_x_16
+=Ax_16;
                src_y_16
+=
Ay_16;
            }
        }
        
return false
;
    }
    
void find_end(long dst_y,long dst_x0,long&
 out_dst_x1)
    {
        
long test_dst_x1=
out_dst_x1;
        
if (test_dst_x1<dst_x0) test_dst_x1=
dst_x0;

        
long src_x_16=Ax_16*test_dst_x1 + Bx_16*dst_y +
 Cx_16;
        
long src_y_16=Ay_16*test_dst_x1 + By_16*dst_y +
 Cy_16;
        
if
 (is_in_src(src_x_16,src_y_16))
        {
            
++
test_dst_x1;
            src_x_16
+=
Ax_16;
            src_y_16
+=
Ay_16;
            
while
 (is_in_src(src_x_16,src_y_16))
            {
                
++
test_dst_x1;
                src_x_16
+=
Ax_16;
                src_y_16
+=
Ay_16;
            }
            out_dst_x1
=
test_dst_x1;
        }
        
else

        {
            src_x_16
-=Ax_16;
            src_y_16
-=
Ay_16;
            
while (!
is_in_src(src_x_16,src_y_16))
            {
                
--
test_dst_x1;
                src_x_16
-=
Ax_16;
                src_y_16
-=
Ay_16;
            }
            out_dst_x1
=
test_dst_x1;
        }
    }

    
bool inti_clip(double move_x,double
 move_y)
    {
        
//计算src中心点映射到dst后的坐标

        out_dst_down_y=(long)(src_height*0.5+move_y);
        out_dst_down_x0
=(long)(src_width*0.5+
move_x);
        out_dst_down_x1
=
out_dst_down_x0;
        
//得到初始out_???

        if (find_begin(out_dst_down_y,out_dst_down_x0,out_dst_down_x1))
            find_end(out_dst_down_y,out_dst_down_x0,out_dst_down_x1);
        out_dst_up_y
=
out_dst_down_y;
        out_dst_up_x0
=
out_dst_down_x0;
        out_dst_up_x1
=
out_dst_down_x1;
        
return (out_dst_down_x0<
out_dst_down_x1);
    }
    
bool
 next_clip_line_down()
    {
        
++
out_dst_down_y;
        
if (!find_begin(out_dst_down_y,out_dst_down_x0,out_dst_down_x1)) return false
;
        find_end(out_dst_down_y,out_dst_down_x0,out_dst_down_x1);
        
return (out_dst_down_x0<
out_dst_down_x1);
    }
    
bool
 next_clip_line_up()
    {
        
--
out_dst_up_y;
        
if (!find_begin(out_dst_up_y,out_dst_up_x0,out_dst_up_x1)) return false
;
        find_end(out_dst_up_y,out_dst_up_x0,out_dst_up_x1);
        
return (out_dst_up_x0<
out_dst_up_x1);
    }
};

void PicRotary3_CopyLine(TARGB32* pDstLine,long dstCount,long Ax_16,long Ay_16,
                        
long srcx0_16,long srcy0_16,TARGB32* pSrcLine,long
 src_byte_width)
{
    
for (long x=0;x<dstCount;++
x)
    {
        pDstLine[x]
=Pixels(pSrcLine,src_byte_width,srcx0_16>>16,srcy0_16>>16
);
        srcx0_16
+=
Ax_16;
        srcy0_16
+=
Ay_16;
    }
}

void PicRotary3(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double
 move_y)
{
    
if ( (fabs(ZoomX*Src.width)<1.0e-4|| (fabs(ZoomY*Src.height)<1.0e-4) ) return//太小的缩放比例认为已经不可见

    double tmprZoomXY=1.0/(ZoomX*ZoomY);  
    
double rZoomX=tmprZoomXY*
ZoomY;
    
double rZoomY=tmprZoomXY*
ZoomX;
    
double
 sinA,cosA;
    SinCos(RotaryAngle,sinA,cosA);
    
long Ax_16=(long)(rZoomX*cosA*(1<<16
)); 
    
long Ay_16=(long)(rZoomX*sinA*(1<<16
)); 
    
long Bx_16=(long)(-rZoomY*sinA*(1<<16
)); 
    
long By_16=(long)(rZoomY*cosA*(1<<16
)); 
    
double rx0=Src.width*0.5;  //(rx0,ry0)为旋转中心 

    double ry0=Src.height*0.5
    
long Cx_16=(long)((-(rx0+move_x)*rZoomX*cosA+(ry0+move_y)*rZoomY*sinA+rx0)*(1<<16
));
    
long Cy_16=(long)((-(rx0+move_x)*rZoomX*sinA-(ry0+move_y)*rZoomY*cosA+ry0)*(1<<16
)); 

    TRotaryClipData rcData;
    rcData.Ax_16
=
Ax_16;
    rcData.Bx_16
=
Bx_16;
    rcData.Cx_16
=
Cx_16;
    rcData.Ay_16
=
Ay_16;
    rcData.By_16
=
By_16;
    rcData.Cy_16
=
Cy_16;
    rcData.dst_width
=
Dst.width;
    rcData.dst_height
=
Dst.height;
    rcData.src_width
=
Src.width;
    rcData.src_height
=
Src.height;
    
if (!rcData.inti_clip(move_x,move_y)) return
;

    TARGB32
* pDstLine=
Dst.pdata;
    ((TUInt8
*&)pDstLine)+=(Dst.byte_width*
rcData.out_dst_down_y);
    
while (true//to down

    {
        
long y=
rcData.out_dst_down_y;
        
if (y>=Dst.height) break
;
        
if (y>=0
)
        {
            
long x0=
rcData.get_down_x0();
            PicRotary3_CopyLine(
&pDstLine[x0],rcData.get_down_x1()-
x0,Ax_16,Ay_16,
                rcData.out_src_x0_16,rcData.out_src_y0_16,Src.pdata,Src.byte_width);
        }
        
if (!rcData.next_clip_line_down()) break
;
        ((TUInt8
*&)pDstLine)+=
Dst.byte_width;
    }
   
    pDstLine
=
Dst.pdata;
    ((TUInt8
*&)pDstLine)+=(Dst.byte_width*
rcData.out_dst_up_y);
    
while (rcData.next_clip_line_up()) //to up

    {
        
long y=
rcData.out_dst_up_y;
        
if (y<0break
;
        ((TUInt8
*&)pDstLine)-=
Dst.byte_width;
        
if (y<
Dst.height)
        {
            
long x0=
rcData.get_up_x0();
            PicRotary3_CopyLine(
&pDstLine[x0],rcData.get_up_x1()-
x0,Ax_16,Ay_16,
                rcData.out_src_x0_16,rcData.out_src_y0_16,Src.pdata,Src.byte_width);
        }
    }
}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicRotary3             280.9 fps
//////////////////////////////////////////////////////////////////////////////// 


F:使用SSE的MOVNTQ指令优化CPU写缓冲

 (仅改写了PicRotary3_CopyLine的实现)


void __declspec(naked) __stdcall PicRotarySSE_CopyLine(TARGB32* pDstLine,long dstCount,long Ax_16,long Ay_16,
                        
//                                      [esp+ 4]      [esp+ 8]      [esp+12]   [esp+16]

                        long srcx0_16,long srcy0_16,TARGB32* pSrcLine,long src_byte_width)
                        
//   [esp+20]      [esp+24]          [esp+28]      [esp+32]

{
    
//利用SSE的MOVNTQ指令优化写缓冲的汇编实现

    asm
    {
        push        ebx  
        push        esi  
        push        edi   
        push        ebp  
        
//esp offset 16

        mov         ebx,dword ptr [esp+ 8+16
        mov         esi,dword ptr [esp
+32+16

        mov         edi,dword ptr [esp
+28+16

        mov         eax,dword ptr [esp
+24+16

        mov         ecx,dword ptr [esp
+20+16

        dec         ebx 
        xor         edx,edx 
        test        ebx,ebx 
        mov         dword ptr [esp
+ 8+16
],ebx 
        jle         loop_bound 

        
//
jmp   loop_begin
        
//align 16

    loop_begin:
            mov         ebx,eax 
            add         eax,dword ptr [esp
+16+16

            sar         ebx,
16
 
            imul        ebx,esi 
            add         ebx,edi 
            mov         ebp,ecx 
            add         ecx,dword ptr [esp
+12+16

            sar         ebp,
16
 
            MOVD        MM0,dword ptr [ebx
+ebp*4

            mov         ebx,eax 
            add         eax,dword ptr [esp
+16+16

            sar         ebx,
16
 
            imul        ebx,esi 
            mov         ebp,ecx 
            add         ecx,dword ptr [esp
+12+16

            sar         ebp,
16
 
            add         ebx,edi 
            MOVD        MM1,dword ptr [ebx
+ebp*4

            mov         ebp,dword ptr [esp
+ 4+16

            PUNPCKlDQ   MM0,MM1
            mov         ebx,dword ptr [esp
+ 8+16

            MOVNTQ      qword ptr [ebp
+edx*4
],MM0 
            add         edx,
2
 
            cmp         edx,ebx  
            jl          loop_begin 

            EMMS

    loop_bound:
        cmp         edx,ebx 
        jne         loop_bound_end 
        sar         eax,
16
 
        imul        eax,esi 
        sar         ecx,
16
 
        add         eax,edi 
        mov         eax,dword ptr [eax
+ecx*4

        mov         ecx,dword ptr [esp
+ 4+16

        mov         dword ptr [ecx
+edx*4
],eax 
    loop_bound_end: 
        pop         ebp  
        pop         edi  
        pop         esi  
        pop         ebx  
        ret         
32
  
    }
}

void PicRotarySSE(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double
 move_y)
{
    
if ( (fabs(ZoomX*Src.width)<1.0e-4|| (fabs(ZoomY*Src.height)<1.0e-4) ) return//太小的缩放比例认为已经不可见

    double tmprZoomXY=1.0/(ZoomX*ZoomY);  
    
double rZoomX=tmprZoomXY*
ZoomY;
    
double rZoomY=tmprZoomXY*
ZoomX;
    
double
 sinA,cosA;
    SinCos(RotaryAngle,sinA,cosA);
    
long Ax_16=(long)(rZoomX*cosA*(1<<16
)); 
    
long Ay_16=(long)(rZoomX*sinA*(1<<16
)); 
    
long Bx_16=(long)(-rZoomY*sinA*(1<<16
)); 
    
long By_16=(long)(rZoomY*cosA*(1<<16
)); 
    
double rx0=Src.width*0.5;  //(rx0,ry0)为旋转中心 

    double ry0=Src.height*0.5
    
long Cx_16=(long)((-(rx0+move_x)*rZoomX*cosA+(ry0+move_y)*rZoomY*sinA+rx0)*(1<<16
));
    
long Cy_16=(long)((-(rx0+move_x)*rZoomX*sinA-(ry0+move_y)*rZoomY*cosA+ry0)*(1<<16
)); 

    TRotaryClipData rcData;
    rcData.Ax_16
=
Ax_16;
    rcData.Bx_16
=
Bx_16;
    rcData.Cx_16
=
Cx_16;
    rcData.Ay_16
=
Ay_16;
    rcData.By_16
=
By_16;
    rcData.Cy_16
=
Cy_16;
    rcData.dst_width
=
Dst.width;
    rcData.dst_height
=
Dst.height;
    rcData.src_width
=
Src.width;
    rcData.src_height
=
Src.height;
    
if (!rcData.inti_clip(move_x,move_y)) return
;

    TARGB32
* pDstLine=
Dst.pdata;
    ((TUInt8
*&)pDstLine)+=(Dst.byte_width*
rcData.out_dst_down_y);
    
while (true//to down

    {
        
long y=
rcData.out_dst_down_y;
        
if (y>=Dst.height) break
;
        
if (y>=0
)
        {
            
long x0=
rcData.get_down_x0();
            PicRotarySSE_CopyLine(
&pDstLine[x0],rcData.get_down_x1()-
x0,Ax_16,Ay_16,
                rcData.out_src_x0_16,rcData.out_src_y0_16,Src.pdata,Src.byte_width);
        }
        
if (!rcData.next_clip_line_down()) break
;
        ((TUInt8
*&)pDstLine)+=
Dst.byte_width;
    }
   
    pDstLine
=
Dst.pdata;
    ((TUInt8
*&)pDstLine)+=(Dst.byte_width*
rcData.out_dst_up_y);
    
while (rcData.next_clip_line_up()) //to up 

    {
        
long y=
rcData.out_dst_up_y;
        
if (y<0break
;
        ((TUInt8
*&)pDstLine)-=
Dst.byte_width;
        
if (y<
Dst.height)
        {
            
long x0=
rcData.get_up_x0();
            PicRotarySSE_CopyLine(
&pDstLine[x0],rcData.get_up_x1()-
x0,Ax_16,Ay_16,
                rcData.out_src_x0_16,rcData.out_src_y0_16,Src.pdata,Src.byte_width);
        }
    }

    asm  sfence 
//刷新写入

}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicRotarySEE           306.3 fps
////////////////////////////////////////////////////////////////////////////////

F':尝试利用SSE2新增的MOVNTI指令优化CPU写缓冲


void __declspec(naked) __stdcall PicRotarySSE2_CopyLine(TARGB32* pDstLine,long dstCount,long Ax_16,long Ay_16,
                        
//                                       [esp+ 4]      [esp+ 8]      [esp+12]   [esp+16]

                        long srcx0_16,long srcy0_16,TARGB32* pSrcLine,long src_byte_width)
                        
//   [esp+20]      [esp+24]          [esp+28]      [esp+32]

{
    
//利用SSE2的MOVNTI指令优化写缓冲的汇编实现

    asm
    {
        push        ebx  
        push        esi  
        push        edi   
        push        ebp  
        
//esp offset 16

        mov         ebx,dword ptr [esp+ 8+16
        mov         esi,dword ptr [esp
+32+16

        mov         edi,dword ptr [esp
+28+16

        mov         eax,dword ptr [esp
+24+16

        mov         ecx,dword ptr [esp
+20+16

        dec         ebx 
        xor         edx,edx 
        test        ebx,ebx 
        mov         dword ptr [esp
+ 8+16
],ebx 
        jle         loop_bound 

        jmp   loop_begin
        align 
16

    loop_begin:
        mov         ebx,eax 
        add         eax,dword ptr [esp
+16+16
        sar         ebx,
16
 
        imul        ebx,esi 
        add         ebx,edi 
        mov         ebp,ecx 
        add         ecx,dword ptr [esp
+12+16

        sar         ebp,
16
 
        mov         ebx,dword ptr [ebx
+ebp*4

        mov         ebp,dword ptr [esp
+ 4+16

        MOVNTI         dword ptr [ebp
+edx*4
],ebx 
        mov         ebx,eax 
        add         eax,dword ptr [esp
+16+16

        sar         ebx,
16
 
        imul        ebx,esi 
        mov         ebp,ecx 
        add         ecx,dword ptr [esp
+12+16

        sar         ebp,
16
 
        add         ebx,edi 
        mov         ebx,dword ptr [ebx
+ebp*4

        mov         ebp,dword ptr [esp
+4+16

        MOVNTI         dword ptr [ebp
+edx*4+4
],ebx 
        mov         ebx,dword ptr [esp
+ 8+16

        add         edx,
2
 
        cmp         edx,ebx 
        jl          loop_begin 
    loop_bound:
        cmp         edx,ebx 
        jne         loop_bound_end 
        sar         eax,
16
 
        imul        eax,esi 
        sar         ecx,
16
 
        add         eax,edi 
        mov         eax,dword ptr [eax
+ecx*4

        mov         ecx,dword ptr [esp
+ 4+16

        mov         dword ptr [ecx
+edx*4
],eax 
    loop_bound_end: 
        pop         ebp  
        pop         edi  
        pop         esi  
        pop         ebx  
        ret         
32
  
    }
}

void PicRotarySSE2(const TPicRegion& Dst,const TPicRegion& Src,double RotaryAngle,double ZoomX,double ZoomY,double move_x,double
 move_y)
{
    
if ( (fabs(ZoomX*Src.width)<1.0e-4|| (fabs(ZoomY*Src.height)<1.0e-4) ) return//太小的缩放比例认为已经不可见

    double tmprZoomXY=1.0/(ZoomX*ZoomY);  
    
double rZoomX=tmprZoomXY*
ZoomY;
    
double rZoomY=tmprZoomXY*
ZoomX;
    
double
 sinA,cosA;
    SinCos(RotaryAngle,sinA,cosA);
    
long Ax_16=(long)(rZoomX*cosA*(1<<16
)); 
    
long Ay_16=(long)(rZoomX*sinA*(1<<16
)); 
    
long Bx_16=(long)(-rZoomY*sinA*(1<<16
)); 
    
long By_16=(long)(rZoomY*cosA*(1<<16
)); 
    
double rx0=Src.width*0.5;  //(rx0,ry0)为旋转中心 

    double ry0=Src.height*0.5
    
long Cx_16=(long)((-(rx0+move_x)*rZoomX*cosA+(ry0+move_y)*rZoomY*sinA+rx0)*(1<<16
));
    
long Cy_16=(long)((-(rx0+move_x)*rZoomX*sinA-(ry0+move_y)*rZoomY*cosA+ry0)*(1<<16
)); 

    TRotaryClipData rcData;
    rcData.Ax_16
=
Ax_16;
    rcData.Bx_16
=
Bx_16;
    rcData.Cx_16
=
Cx_16;
    rcData.Ay_16
=
Ay_16;
    rcData.By_16
=
By_16;
    rcData.Cy_16
=
Cy_16;
    rcData.dst_width
=
Dst.width;
    rcData.dst_height
=
Dst.height;
    rcData.src_width
=
Src.width;
    rcData.src_height
=
Src.height;
    
if (!rcData.inti_clip(move_x,move_y)) return
;

    TARGB32
* pDstLine=
Dst.pdata;
    ((TUInt8
*&)pDstLine)+=(Dst.byte_width*
rcData.out_dst_down_y);
    
while (true//to down

    {
        
long y=
rcData.out_dst_down_y;
        
if (y>=Dst.height) break
;
        
if (y>=0
)
        {
            
long x0=
rcData.get_down_x0();
            PicRotarySSE2_CopyLine(
&pDstLine[x0],rcData.get_down_x1()-
x0,Ax_16,Ay_16,
                rcData.out_src_x0_16,rcData.out_src_y0_16,Src.pdata,Src.byte_width);
        }
        
if (!rcData.next_clip_line_down()) break
;
        ((TUInt8
*&)pDstLine)+=
Dst.byte_width;
    }
   
    pDstLine
=
Dst.pdata;
    ((TUInt8
*&)pDstLine)+=(Dst.byte_width*
rcData.out_dst_up_y);
    
while (rcData.next_clip_line_up()) //to up 

    {
        
long y=
rcData.out_dst_up_y;
        
if (y<0break
;
        ((TUInt8
*&)pDstLine)-=
Dst.byte_width;
        
if (y<
Dst.height)
        {
            
long x0=
rcData.get_up_x0();
            PicRotarySSE2_CopyLine(
&pDstLine[x0],rcData.get_up_x1()-
x0,Ax_16,Ay_16,
                rcData.out_src_x0_16,rcData.out_src_y0_16,Src.pdata,Src.byte_width);
        }
    }

    asm  sfence 
//刷新写入

}

////////////////////////////////////////////////////////////////////////////////
//速度测试:
//==============================================================================
// PicRotarySEE2          304.2 fps
////////////////////////////////////////////////////////////////////////////////

一张效果图:
 //程序使用的调用参数:
    const long testcount=2000;
    long dst_wh=1004;
    for (int i=0;i<testcount;++i)
    {
        double zoom=rand()*(1.0/RAND_MAX)+0.5;
        PicRotarySSE(ppicDst,ppicSrc,rand()*(PI*2/RAND_MAX),zoom,zoom,((dst_wh+ppicSrc.width)*rand()*(1.0/RAND_MAX)-ppicSrc.width),(dst_wh+ppicSrc.height)*rand()*(1.0/RAND_MAX)-ppicSrc.height);
    }

 //ps:如果很多时候源图片绘制时可能落在目标区域的外面,那么需要写一个剪切算法快速排除不必要的绘制

 

一张测试函数速度的时候生成的图像:


G:旋转测试的结果放到一起:

//注:测试图片都是800*600的图片旋转到1004*1004的图片中心,测试成绩取各个旋转角度的平均速度值  
////////////////////////////////////////////////////////////////////////////////
//速度测试:  (测试CPU为AMD64x2 4200+(2.37G),单线程)
//==============================================================================
// PicRotary0              34.9 fps
// PicRotary1              62.0 fps
// PicRotary2             134.2 fps
// PicRotary3             280.9 fps
// PicRotarySEE           306.3 fps
// PicRotarySEE2          304.2 fps
//(PicRotarySSE2_Block    316.6 fps (参见《下篇 补充话题》))
////////////////////////////////////////////////////////////////////////////////

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

////////////////////////////////////////////////////////////////////////////////
//速度测试:  (测试CPU为Intel Core2 4400(2.00G)单线程)
//==============================================================================
// PicRotary0              58.6 fps
// PicRotary1              82.1 fps
// PicRotary2             167.9 fps
// PicRotary3             334.9 fps
// PicRotarySEE           463.1 fps
// PicRotarySEE2          449.3 fps
//(PicRotarySSE2_Block    351.3 fps (参见《下篇 补充话题》))
////////////////////////////////////////////////////////////////////////////////