二次插值算法
- 格式:doc
- 大小:235.50 KB
- 文档页数:8
常见的插值方法及其原理这一节无可避免要接触一些数学知识,为了让本文通俗易懂,我们尽量绕开讨厌的公式等。
为了进一步的简化难度,我们把讨论从二维图像降到一维上。
首先来看看最简单的‘最临近像素插值’。
A,B是原图上已经有的点,现在我们要知道其中间X位置处的像素值。
我们找出X位置和A,B位置之间的距离d1,d2,如图,d2要小于d1,所以我们就认为X处像素值的大小就等于B处像素值的大小。
显然,这种方法是非常苯的,同时会带来明显的失真。
在A,B中点处的像素值会突然出现一个跳跃,这就是为什么会出现马赛克和锯齿等明显走样的原因。
最临近插值法唯一的优点就是速度快。
图10,最临近法插值原理接下来是稍微复杂点的‘线性插值’(Linear)线性插值也很好理解,AB两点的像素值之间,我们认为是直线变化的,要求X点处的值,只需要找到对应位置直线上的一点即可。
换句话说,A,B间任意一点的值只跟A,B有关。
由于插值的结果是连续的,所以视觉上会比最小临近法要好一些。
线性插值速度稍微要慢一点,但是效果要好不少。
如果讲究速度,这是个不错的折衷。
图11,线性插值原理其他插值方法立方插值,样条插值等等,他们的目的是试图让插值的曲线显得更平滑,为了达到这个目的,他们不得不利用到周围若干范围内的点,这里的数学原理就不再详述了。
图12,高级的插值原理如图,要求B,C之间X的值,需要利用B,C周围A,B,C,D四个点的像素值,通过某种计算,得到光滑的曲线,从而算出X的值来。
计算量显然要比前两种大许多。
好了,以上就是基本知识。
所谓两次线性和两次立方实际上就是把刚才的分析拓展到二维空间上,在宽和高方向上作两次插值的意思。
在以上的基础上,有的软件还发展了更复杂的改进的插值方式譬如S-SPline, Turbo Photo等。
他们的目的是使边缘的表现更完美。
插值(Interpolation),有时也称为“重置样本”,是在不生成像素的情况下增加图像像素大小的一种方法,在周围像素色彩的基础上用数学公式计算丢失像素的色彩。
lipschitz常数求法Lipschitz常数是函数分析和数学优化领域中的一个重要概念,它描述了函数的导数变化的上限。
Lipschitz常数的计算对于优化算法、最优控制和机器学习等许多领域都具有重要意义。
在本文中,我们将介绍Lipschitz常数的定义以及两种常用的求解方法。
1. Lipschitz常数的定义Lipschitz常数用来衡量函数的导数变化的上限。
对于一个实数域上的函数 f(x),如果存在一个常数 L,使得对于任意的 x1 和 x2,都有: |f(x1) - f(x2)| ≤ L |x1 - x2|则称函数 f(x) 是Lipschitz连续的,而常数 L 就是该函数的Lipschitz常数。
2. 间隔法(Interval Method)间隔法是一种较为直观的Lipschitz常数求解方法。
它通过计算一个函数在给定区间内的变化率来估计Lipschitz常数。
假设我们要求解函数 f(x) 在区间 [a, b] 上的Lipschitz常数,首先可以在该区间上随机选择两个点 x1 和 x2,并计算其对应函数值之间的差值和 x1、x2 之间的差值的比值。
然后,通过遍历区间内的每一对点,取所有比值的最大值作为Lipschitz常数的估计。
即可得到:L = max(|(f(x1)-f(x2))/(x1-x2)|)3. 二次插值法(Quadratic Interpolation Method)二次插值法是一种更精确的Lipschitz常数求解方法,它利用了函数的二阶导数信息。
假设 f(x) 是定义在 [a,b] 区间上的二次可微函数,并且其二阶导数满足 Lipschitz 条件,即存在 L > 0,对于所有 x∈[a,b],都有|f''(x)| ≤ L。
在该条件下,可以得到 Lipschitz 常数的上限估计 L = L_f =max(|f''(x)| : x∈[a,b])。
halcon曲线插值算法
Halcon是一个强大的机器视觉软件库,提供了多种曲线插值算法。
以下是Halcon 中常用的曲线插值算法:
线性插值:线性插值是最简单的插值方法,通过构造两条直线来连接已知的点。
虽然线性插值的计算速度较快,但对于非线性变化的曲线,其拟合效果可能不够理想。
二次插值:二次插值比线性插值更加复杂,通过构造二次多项式来逼近已知的点。
二次插值能够更好地处理非线性变化的曲线,但计算量较大。
样条插值:样条插值是一种更复杂的插值方法,通过构造样条函数(如三次样条函数)来逼近已知的点。
样条插值能够得到更加平滑的曲线,适合处理具有突变点的数据。
多项式插值:多项式插值是通过构造多项式来逼近已知的点。
多项式插值的拟合效果较好,但计算量较大,且容易受到噪声的影响。
立方插值:立方插值是一种较新的插值方法,通过构造立方多项式来逼近已知的点。
立方插值的拟合效果非常好,能够得到非常平滑的曲线,但计算量较大。
在选择合适的曲线插值算法时,需要考虑数据的特点、计算资源和精度要求等因素。
对于要求较高精度的应用,可以考虑使用立方插值或样条插值;对于计算资源有限的应用,可以考虑使用线性插值或二次插值。
二次插值法亦是用于一元函数在确定的初始区间搜索极小点的一种方法。
它属于曲线拟合方法的畴。
一、基本原理在求解一元函数的极小点时,常常利用一个低次插值多项式来逼近原目标函数,然后求该多项式的极小点(低次多项式的极小点比较容易计算),并以此作为目标函数的近似极小点。
如果其近似的程度尚未达到所要求的精度时,可以反复使用此法,逐次拟合,直到满足给定的精度时为止。
常用的插值多项式为二次或三次多项式,分别称为二次插值法和三次插值法。
这里我们主要介绍二次插值法的计算公式。
假定目标函数在初始搜索区间中有三点、和,其函数值分别为、和(图1},且满足,,即满足函数值为两头大中间小的性质。
利用这三点及相应的函数值作一条二次曲线,其函数为一个二次多项式(1)式中、、为待定系数。
图1根据插值条件,插值函数与原函数在插值结点、、处函数值相等,得(2)为求插值多项式的极小点,可令其一阶导数为零,即(3)解式(3)即求得插值函数的极小点(4)式(4)中要确定的系数可在方程组(2)中利用相邻两个方程消去而得:(5)(6)将式(5)、(6)代入式(4)便得插值函数极小值点的计算公式:(7)把取作区间的另一个计算点,比较与两点函数值的大小,在保持两头大中间小的前提下缩短搜索区间,从而构成新的三点搜索区间,再继续按上述方法进行三点二次插值运算,直到满足规定的精度要求为止,把得到的最后的作为的近似极小值点。
上述求极值点的方法称为三点二次插值法。
为便于计算,可将式(7)改写为(8)式中:(9)(10)二、迭代过程及算法框图(1)确定初始插值结点通常取初始搜索区间的两端点及中点为,,。
计算函数值,,,构成三个初始插值结点、、。
(2)计算二次插值函数极小点按式(8)计算,并将记作点,计算。
若本步骤为对初始搜索区间的第一次插值或点仍为初始给定点时,则进行下一步(3);否则转步骤(4)(3)缩短搜索区间缩短搜索区间的原则是:比较函数值、,取其小者所对应的点作为新的点,并以此点左右两邻点分别取作新的和,构成缩短后的新搜索区间。
经典数值算法及其maple实现经典数值算法是计算机科学中常用的一种算法,用于解决数值计算问题。
这些算法被广泛应用于科学计算、工程计算、金融计算等领域。
下面列举了10个经典数值算法及其Maple实现。
1. 二分法(Bisection Method)二分法是一种求解方程根的迭代算法。
通过将区间不断地二分,确定方程在给定区间内的根的近似值。
具体实现如下:```Maplebisection := proc(f, a, b, tol)local c, fc;while abs(b - a) > tol doc := (a + b) / 2;fc := evalf(f(c));if f(a) * fc < 0 thenb := c;elsea := c;end if;end do;return (a + b) / 2;end proc;```2. 牛顿法(Newton's Method)牛顿法是一种求解方程根的迭代算法。
通过利用函数的切线逼近方程的根,求得根的近似值。
具体实现如下:```Maplenewton := proc(f, x0, tol)local x, fx, dfx;x := x0;repeatfx := evalf(f(x));dfx := evalf(D(f)(x));x := x - fx / dfx;until abs(fx) < tol;return x;end proc;```3. 高斯消元法(Gaussian Elimination)高斯消元法是一种求解线性方程组的算法。
通过将线性方程组转化为阶梯形矩阵,再利用回代法求解方程组的解。
具体实现如下: ```MaplegaussianElimination := proc(A, b)local n, i, j, k, factor;n := RowDimension(A);for k from 1 to n-1 dofor i from k+1 to n dofactor := A[i, k] / A[k, k];for j from k+1 to n doA[i, j] := A[i, j] - factor * A[k, j];end do;b[i] := b[i] - factor * b[k];end do;end do;return A, b;end proc;```4. 欧拉方法(Euler's Method)欧拉方法是一种求解常微分方程初值问题的算法。
二次插值算法范文二次插值算法是一种在离散数据点之间进行数据估计的方法,它通过利用已知数据点及其相邻数据点的信息,来计算出两个数据点之间的值。
二次插值算法可以应用于图像处理、信号处理、数值分析等领域,在实际应用中具有很高的效率和准确性。
二次插值算法的基本原理是通过已知的数据点来构造一个二次函数,然后利用这个函数来估计两个数据点之间的值。
具体而言,假设我们有一组数据点{(x1, y1), (x2, y2), ..., (xn, yn)},我们希望估计在一些位置x处的值y。
首先,我们选择x所在的两个已知数据点,假设它们分别为(xi, yi)和(xi+1, yi+1),然后我们构造一个二次函数f(x) = ax^2 + bx + c,满足条件f(xi) = yi和f(xi+1) = yi+1、通过解这个方程组,我们可以得到二次函数的系数a、b和c,然后通过计算f(x)即可得到在位置x处的估计值y。
下面我们来具体讨论二次插值算法的实现。
首先,我们需要定义一个函数quadratic_interpolation,它接受三个参数:x, data和index。
其中x为需要估计的位置,data为已知的数据点集合,index为x所在的两个数据点的索引。
```pythondef quadratic_interpolation(x, data, index):x0, y0 = data[index]x1, y1 = data[index + 1]a=(y0-2*y1+y2)/((x0-x1)*(x0-x2))b=(y1-y0)/(x1-x0)-a*(x0+x1)c=y0-a*x0**2-b*x0return a*x**2 + b*x + c```在这段代码中,我们首先通过传入的index参数找到x所在的两个数据点,然后根据这两个数据点构造二次函数的系数a、b和c,并最终返回在位置x处的估计值。
接下来,我们可以使用这个二次插值算法来对一组数据进行插值。
interp2d 二阶插值函数什么是二阶插值函数interp2d?interp2d是一种二阶插值函数,用于插值二维平面上的离散数据点以获得连续函数的近似值。
它可以用来处理不规则的二维数据,如图像数据、地理信息数据等。
在插值算法中,二阶插值是相对简单且对于大多数应用来说已经足够精确的一种插值方法。
在程序中,interp2d是用Python实现的。
为什么需要二阶插值函数interp2d?在实际生活和工作中,我们常常会遇到只有有限离散数据点的情况。
例如,一个地面测量的高程图数据,仅包含有限的离散数据点。
但是,当我们需要得到该地区的完整高程图时,需要利用插值函数对这些离散数据点进行插值,获得空间上连续的高程函数。
此时,interp2d即可派上用场。
interp2d可以将已知数据点组成的离散网格,拟合一个连续函数,并且可以对该函数在任意给定位置的值进行估计。
这就在一定程度上解决了离散数据无法直接处理的问题,同时保证了生成的函数的平滑性和准确性。
如何使用二阶插值函数interp2d?使用interp2d进行二维数据插值通常需要以下步骤:1. 准备好离散数据点的坐标和对应的函数值。
通常情况下,这些数据点以网格的形式呈现,即横坐标和纵坐标都是等间距的。
在Python中,可以使用numpy库生成这样的离散数据点。
2. 创建interp2d对象并调用。
在Python中,可以使用如下代码创建interp2d对象:from scipy.interpolate import interp2df = interp2d(x, y, z, kind='cubic')其中,x和y是数据点的横纵坐标,z是数据点在对应坐标上的函数值,kind参数指定了插值所使用的方法。
在这里我们指定为'cubic',即使用三次样条插值方法,一般而言效果最好。
3. 调用生成的函数计算插值结果。
在生成了interp2d对象之后,可以通过调用生成的函数来计算插值结果。
数值分析实验报告线性插值和二次插值计算ln0.54的近似值数值分析实验报告线性插值和二次插值计算ln0.54的近似值篇一:数值分析-用线性插值及二次插值计算数值分析上机报告习题:给出f(x)?lnx的数值表,用线性插值及二次插值计算ln0.54的近似值。
解:(1)用线性插值计算 Matla b程序 x=0.54; a=[0.5,0.6];b=[-0.693147,-0.510826]; l1=b (1)*((x-a(2))/(a(1)-a (2))); l2=b(2)*((x-a(1))/(a(2)-a(1))); y=l1+l2 y = -0.6202(2)用抛物插值计算 Ma tlab程序 x=0.54; a=[0.4,0.5,0.6]; b=[-0.916291,-0.693147,-0.510826]; A=b(1)*(x-a(2))*(x-a(3))/((a (1)-a(2))*(a(1)-a(3))); B=b(2)*(x-a (1))*(x-a(3))/((a(2)-a(1))*(a(2)-a(3))); C=b(3)*(x-a(1))*(x-a(2))/((a(3)-a(1))*(a(3)-a(2)));y=A+B+C y= -0.6153篇二:数值分析上机实验报告二实验报告二题目:如何求解插值函数摘要:在工程测量和科学实验中,所得到的数据通常都是离散的,如果要得到这些离散点意外的其他点的数值,就需要根据这些已知数据进行插值。
这里我们将采用多种插值方法。
前言:(目的和意义)掌握Lagrange,Netn,Hermi te,线性,三次样条插值法的原理及应用,并能求解相应问题。
数学原理:主要的插值法有:多项式插值法、拉格朗日插值法、线性插值法、牛顿插值法,H ermite插值法三次样条插值法等。
高质量的快速的图像缩放——二次线性插值和三次卷积插值限制条件:为了便于讨论,这里只处理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 __asmtypedef 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_BilInear1inline 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)(255-v_8));unsigned long pm1_16=(v_8*(unsigned long)(255-u_8));unsigned long pm0_16=((255-u_8)*(255-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+1;long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;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]); //bordersrcx_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+1;long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;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/yrif (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)/yrif (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]); //bordersrcx_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]);//bordersrcx_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]);//bordersrcx_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]); //bordersrcx_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)证明这个公式不会造成内存访问超界:要求Dx=DW-1时: sx+1=int( (dw-1)/dw*(dw-1) ) +1 <= (sw-1)有: int( (sw-1)*(dw-1)/dw ) <=sw-2(sw-1)*(dw-1)/dw <(sw-1)(dw-1) /dw<1(dw-1) <dw比如,按这个公式的一个简单实现: (缩放效果见前面的"二次线性插值(近似公式)"图示)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+1;long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;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/yrif (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)/yrif (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]); //bordersrcx_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]);//bordersrcx_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]);//bordersrcx_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]); //bordersrcx_16+=xrIntFloat_16;}srcy_16+=yrIntFloat_16;((TUInt8*&)pDstLine)+=Dst.byte_width;}}//////////////////////////////////////////////////////////////////////////////////速度测试://============================================================================== // PicZoom_BilInear_Common 65.3 fps//////////////////////////////////////////////////////////////////////////////// H: 使用MMX指令改写:PicZoom_Bilinear_MMXinline void Bilinear_Fast_MMX(TARGB32* PColor0,TARGB32* PColor1, unsigned long u_8,unsigned long v_8,TARGB32* result){asm{MOVD MM6,v_8MOVD MM5,u_8mov edx,PColor0mov eax,PColor1PXOR mm7,mm7MOVD MM2,dword ptr [eax]MOVD MM0,dword ptr [eax+4]PUNPCKLWD MM5,MM5PUNPCKLWD MM6,MM6MOVD MM3,dword ptr [edx]MOVD MM1,dword ptr [edx+4]PUNPCKLDQ MM5,MM5PUNPCKLBW MM0,MM7PUNPCKLBW MM1,MM7PUNPCKLBW MM2,MM7PUNPCKLBW MM3,MM7PSUBw MM0,MM2PSUBw MM1,MM3PSLLw MM2,8PSLLw MM3,8PMULlw MM0,MM5PMULlw MM1,MM5PUNPCKLDQ MM6,MM6PADDw MM0,MM2PADDw MM1,MM3PSRLw MM0,8PSRLw MM1,8PSUBw MM0,MM1PSLLw MM1,8PMULlw MM0,MM6mov eax,resultPADDw MM0,MM1PSRLw MM0,8PACKUSwb MM0,MM7movd [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+1;long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;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/yrif (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)/yrif (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]);//bordersrcx_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]);//bordersrcx_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]); //bordersrcx_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_Exvoid 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+1;long yrIntFloat_16=((Src.height)<<16)/Dst.height+1;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/yrif (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)/yrif (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]); //bordersrcx_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]);//bordersrcx_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_8pxor mm7,mm7 //mm7=0PUNPCKLWD MM6,MM6PUNPCKLDQ MM6,MM6//mm6=v_8mov esi,PSrcLineColormov ecx,PSrcLineColorNextmov edx,srcx_16mov ebx,dst_width_fastmov edi,pDstLine_Fastlea edi,[edi+ebx*4]push ebpmov ebp,xrIntFloat_16neg ebxloop_start:mov eax,edxshl eax,16shr eax,24//== movzx eax,dh //eax=u_8MOVD MM5,eaxmov eax,edxshr eax,16 //srcx_16>>16MOVD MM2,dword ptr [ecx+eax*4]MOVD MM0,dword ptr [ecx+eax*4+4]PUNPCKLWD MM5,MM5MOVD MM3,dword ptr [esi+eax*4]MOVD MM1,dword ptr [esi+eax*4+4]PUNPCKLDQ MM5,MM5 //mm5=u_8PUNPCKLBW MM0,MM7PUNPCKLBW MM1,MM7PUNPCKLBW MM2,MM7PUNPCKLBW MM3,MM7PSUBw MM0,MM2PSUBw MM1,MM3PSLLw MM2,8PSLLw MM3,8PMULlw MM0,MM5PMULlw MM1,MM5PADDw MM0,MM2PADDw MM1,MM3PSRLw MM0,8PSRLw MM1,8PSUBw MM0,MM1PSLLw MM1,8PMULlw MM0,MM6PADDw MM0,MM1PSRLw MM0,8PACKUSwb MM0,MM7MOVd dword ptr [edi+ebx*4],MM0 //write DstColor add edx,ebp //srcx_16+=xrIntFloat_16inc ebxjnz loop_startpop ebpmov srcx_16,edx}}}for (x=border_x1;x<dst_width;++x){Bilinear_Border_MMX(Src,srcx_16,srcy_16,&pDstLine[x]);//bordersrcx_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]); //bordersrcx_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_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)*x - (4*a);elsereturn 0;}inline TUInt8 border_color(long Color){if (Color<=0)return 0;else if (Color>=255)return 255;elsereturn 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);}}。
二次插值法亦是用于一元函数在确定的初始区间内搜索极小点的一种方法。
它属于曲线拟合方法的范畴。
一、基本原理
在求解一元函数的极小点时,常常利用一个低次插值多项式来逼近原目标函数,
然后求该多项式的极小点(低次多项式的极小点比较容易计算),并以此作为目标函数
的近似极小点。
如果其近似的程度尚未达到所要求的精度时,可以反复使用此法,逐次拟合,直到满足给定的精度时为止。
常用的插值多项式为二次或三次多项式,分别称为二次插值法和三次插值法。
这里我们主要介绍二次插值法的计算公式。
假定目标函数在初始搜索区间中有三点、和
,其函数值分别为、和(图1},且满足,
,即满足函数值为两头大中间小的性质。
利用这三点及相应的函数值作一条二次曲
线,其函数为一个二次多项式
(1)
式中、、为待定系数。
图1
根据插值条件,插值函数与原函数在插值结点、、处函数值相等,得
(2)
为求插值多项式的极小点,可令其一阶导数为零,即
(3)
解式(3)即求得插值函数的极小点(4)
式(4)中要确定的系数可在方程组(2)中利用相邻两个方程消去而得:
(5)
(6)
将式(5)、(6)代入式(4)便得插值函数极小值点的计算公式:
(7)
把取作区间内的另一个计算点,比较与两点函数值的大小,在保持两头大中间小的前提下缩短搜索区间,从而构成新的三点搜索区间,再继续按上述方
法进行三点二次插值运算,直到满足规定的精度要求为止,把得到的最后的作为的近似极小值点。
上述求极值点的方法称为三点二次插值法。
为便于计算,可将式(7)改写为
(8)
式中:
(9)
(10)
二、迭代过程及算法框图
(1)确定初始插值结点
通常取初始搜索区间的两端点及中点为,,。
计算函数值,,,构成三个初始插值结点、、。
(2)计算二次插值函数极小点
按式(8)计算,并将记作点,计算。
若本步骤为对初始搜索区间的第一次插值或点仍为初始给定点时,则进行下一步(3);否则转步骤(4)
(3)缩短搜索区间
缩短搜索区间的原则是:比较函数值、,取其小者所对应的点作为新的点,并以此点左右两邻点分别取作新的和,构成缩短后的新搜索区间。
其具体方法则如图2所示,根据原区间中和的相对位置以及函数值和之比较有a、b、c、d四种情况,图中阴影线部分表示丢去的区间。
在对新区间三个新点的代号作依次、
、的一般化处理后,计算其函数值,并令,,,返回步骤(2)。
图2(a)
图2(b)
图2(c)
图2(d)
(4)判断迭代终止条件
在一般情况下,因是前一次插值函数的极小值点,是本次插值函数的极小值点,若和的距离足够小时,即满足,或和两者原函数值已很接
近,即满足,则停止迭代,这时,若,输出极小值点,
极小值;否则,即时,输出极小值点,极小值。
如不满足上述迭代终止条件,则返回步骤(3),再次缩短搜索区间,直至最后满足终止条件。
按上述步骤设计的二次插值法算法框图见图3。
图3
算法框图中有几点需作些说明。
1.判别框?若成立,按式(9)和式(10)则有
说明三个插值结点、、在一条直线上;
2.判别框?若不成立,说明落在区间之外。
上述两种情况只是在区间已缩得很小,由于三个插值结点已十分接近,计算机的舍入误差才可能使其发生。
此时取和作为最优解应是合理的。
3.在初始搜索区间第一次插值或仍为初始给定点时,和并不代表前后二次插值函数极小点,因而判别式并不能确切地反映该不该终止迭代,这时应进行步骤(3)缩短搜索区间,直至初始点第一次由代替,使用判别式?进行终止判别才具意义。
为此,算法框图中设置开关K=0和K=1分别表示初始点第一次由代替前和后的状态。