第六章图像插值
第六章插值法
(6.1)
则称 y ( x )为函数 f ( x ) 的插值多项式, 称
[ a , b ] 为插值区间,
条件(6.1)称为插值条件, x i 称为插值节点。 几何意义如图所示。由插值函数的定义可知求插值多项 式 y ( x )，即是使曲线 y y ( x ) 与曲线 y f ( x ) 在平面上有 n+1个交点 ( x , f ( x ))， i 0,1, 2, , n
解
1. 线性插值
L1 ( x ) x9 49 2 x4 94 3 2 5 ( x 9) 3 5 ( x 4)
7 L1 ( 7 )
2 5
(7 9)
3 5
(7 4)
13 5
2 .6
2. 抛物插值
7 L 2 (7 ) ( 7 9 )( 7 16 ) ( 4 9 )( 4 16 ) 2 ( 7 4 )( 7 16 ) ( 9 4 )( 9 16 ) 3 ( 7 4 )( 7 9 ) (16 4 )( 16 9 ) 4
4
该线性方程组的系数行列式是Vandermonde行列式
1 1 1 x0 x1 xn x0 x1 xn
2 2 2
x0 x1
n n
0i j n
( x j xi点是相异的节点,故该系数行列式不等于0。由Cramer 法则，方程组有一组唯一的解 a 0 , a1 , , a n 求出方程组的解就得到了插值多项式，称之为待定系数法。但 当节点数较多时，求解方程组的工作量和产生的误差都比较大。 由插值多项式的唯一性，可以利用更简便实用并可行的方法来 构造插值多项式。 5
第六章 插值法
二、Newton 插值多项式(n次)
f ( x) f ( x0 ) f ( x0 , x1 )( x x0 ) f ( x0 , x1 , x2 )( x x0 )( x x1 )
f ( x0 , x1 ,, xn )( x x0 )( x x1 )( x xn1 )
f ( x0 , x1 )
f ( x1 , x2 ) f ( x0 , x1 , x2 ) f ( x2 , x3 ) f ( x1 , x2 , x3 ) f ( x 3 , x 4 ) f ( x 2 , x 3 , x4 )
f ( x k 1 , x k ) f ( x k 2 , x k 1 , x k )
x
1 1 2
y
一阶差商 二阶差商
2 3 1
1 2 2
5 6
N 2 x f ( x0 ) f ( x0 , x1 ) x x0 f ( x0 , x1 , x2 ) x x0 x x1
1 x 1 5 x 1 x 1 2 2 6
差分表6-2
xi x0 x1 x2 f ( xi ) f0 f1 f1 f2 f 2 x3 x4 xn f3 f4 f n1 2 f 2 f f 0 2 f 3 f 4 f n f
f0
f 0
前插的系数
2 ff 0 2 0 2 f1 0 3 ff 0 f 3 f1 3
解：设 0 1, x1 1, x2 2 x
x x1 x x2 x 1 x 2 1 2 l0 x ( x 3 x 2) x0 x1 x0 x2 1 1 1 2 6 1 2 x x0 x x2 x 1 x 2 l1 x ( x x 2) x1 x0 x1 x2 1 11 2 2 x x0 x x1 x 1 x 1 1 2 l2 x ( x 1) x2 x0 x2 x1 2 12 1 3
常见的插值方法及其原理
常见的插值方法及其原理这一节无可避免要接触一些数学知识，为了让本文通俗易懂，我们尽量绕开讨厌的公式等。
为了进一步的简化难度，我们把讨论从二维图像降到一维上。
首先来看看最简单的‘最临近像素插值’。
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)，有时也称为“重置样本”，是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩。
图像的几何运算
图像的⼏何运算@⽬录图像的⼏何运算是指引起图像⼏何形状发⽣改变的变换。
与点运算不同的是，⼏何运算可以看成是像素在图像内的移动过程，该移动过程可以改变图像中物体对象之间的空间关系。
1.图像的插值图像插值是指利⽤已知邻近像素点的灰度值来产⽣位置像素点的灰度值，以便由原始图像再⽣成具有更⾼分辨率的图像。
插值是在不⽣成新的像素的情况下对原图像的像素重新分布，从⽽改变像素数量的⼀种⽅法。
在图像放⼤过程中，像素也相应的增加，增加的过程就是‘插值’发⽣作⽤的过程，‘’插值程序⾃动选择信息较好的像素作为增加、弥补空⽩像素的空间，⽽并⾮只使⽤近邻的像素，所以在放⼤图像时，图像看上去会⽐较平滑、⼲净。
⽆论使⽤何种插值⽅法，⾸先都需要找到与输出图像像素相对应的输⼊图像点，然后再通过计算该点附近某⼀像素集合的权平均值来指定输出像素的灰度值。
像素的权是根据像素到点的距离来⽽定的，不同插值⽅法的区别就在于考虑的像素集合不同。
最常见的插值⽅法如下：（1）向前映射法：通过输⼊图像像素的位置，计算输出图像对应像素的位置，将该位置像素的灰度值按某种⽅式分配到输出图像相邻的四个像素。
（2）向后映射法：通过输出图像像素位置，计算输⼊图像对应像素的位置，根据输⼊图像相邻四个像素的灰度值计算该位置像素的灰度值。
（3）最近邻插值：表⽰输出像素将被指定为像素点所在位置处的像素值。
（4）双线性插值：表⽰输出像素值是像素2×2邻域内的平均值。
（5）双三次插值：表⽰输出像素值是像素4×4邻域内的权平均值。
在MATLAB中，interp2函数⽤于对图像进⾏插值处理，该函数的调⽤⽅法如下：A=interp2(X,Y,Z,IX,IY):Z为要插值的原始图像，IX和IY为图像的新⾏和新列clear allclose allclcI2=imread('eight.tif');subplot(231)imshow(I2)title('原始图像')Z1=interp2(double(I2),2,'nearest');%最近邻插值法Z1=uint8(Z1);subplot(232)imshow(Z1)title('最近邻插值')Z2=interp2(double(I2),2,'linear');%线性插值法Z2=uint8(Z2);subplot(232)imshow(Z2)title('线性插值法')Z3=interp2(double(I2),2,'spline');%三次样条插值法Z3=uint8(Z3);subplot(234)imshow(Z3);title('三次样条插值');Z4=interp2(double(I2),2,'cubic');%⽴⽅插值法Z4=uint8(Z4);subplot(235);imshow(Z4);title('⽴⽅插值')2.旋转与平移变换旋转变换的表达式为⽤齐次矩阵表⽰为在MATLAB中，使⽤imrotate函数来旋转⼀幅图像，调⽤格式如下：B=imrotate(A,ANGLE,METHOD,BBOX)其中，A是需要旋转的图像；ANGLE是旋转的⾓度，正值为逆时针；METHOD是插值⽅法；BBOX表⽰旋转后的显⽰⽅式。
第6章 函数逼近与函数插值
第六章 函数逼近与函数插值本章介绍函数逼近与插值的有关理论和算法. 函数逼近问题与插值问题两者既有联系又有区别，它们都是用较简单的函数来近似未知的、或表达式较复杂的函数. 一般来说，函数逼近是要在整个区间、或一系列离散点上整体逼近被近似函数，而在进行插值时，则须保证在若干自变量点上的函数值与被近似函数相等.6.1 函数逼近的基本概念进行函数逼近一般是在较简单的函数类Φ中找一个函数p(x)来近似给定的函数f(x)，以使得在某种度量意义下误差函数p (x )−f(x)最小. 被逼近函数f(x)可能是较复杂的连续函数，也可能是只在一些离散点上定义的表格函数，而函数类Φ可以是多项式、分段多项式、三角函数、有理函数，等等. 函数逼近问题中度量误差的手段主要是函数空间的范数，下面先介绍函数空间的范数、内积等有关概念，然后讨论函数逼近问题的不同类型.6.1.1 函数空间线性空间的概念大家都很熟悉，其定义中包括一个元素集合和一个数域，以及满足一定运算规则的“加法”和“数乘”运算. 简单说，若这个元素集合对于“加法”和“数乘”运算封闭，则为一线性空间. 线性空间的元素之间存在线性相关和线性无关两种关系，进而又有空间的基和维数的概念.在这里我们先考虑连续函数形成的线性空间. 例如C [a,b ]按函数加法、以及函数与实数乘法，构成一个线性空间. 对于[a,b]区间上所有k 阶导数连续的函数全体C k [a,b ]，也类似地构成一个线性空间. 我们一般讨论实数函数，因此对应的是实数域ℝ，若讨论复数函数，则相应的是复数域ℂ. 另外，与线性代数中讨论的向量空间ℝn 不同，连续函数空间是无限维的.对线性空间可以定义范数的概念（见3.1.2节）. 针对实连续函数空间C [a,b ]，与向量空间类似，可定义如下三种函数的范数(function norm)：1) ∞-范数 设f (x )∈C [a,b ]，则‖f (x )‖∞=max x∈[a,b ]|f (x )| .其几何意义如图6-1所示，即函数值绝对值的最大值.2) 1-范数‖f (x )‖1=∫|f (x )|dx b a .其几何意义如图6-2所示，即函数曲线与横轴之间的面积总和.3) 2-范数‖f (x )‖2=[∫f 2(x )dx b a ]1/2. 2-范数也常称为平方范数，其几何意义与1-范数类似. 线性空间还有一个重要概念是内积，它定义了空间中两个元素的一种运算. 下面给出一般的复数域上线性空间内积的定义.定义6.1：设S为实数域ℝ上的线性空间，∀u,v∈S，定义值域为ℝ的二元运算〈u,v〉，若满足1)〈u,v〉=〈v,u〉, (可交换性)2)〈αu,v〉=α〈u,v〉, ∀α∈ℂ(线性性1)3)〈u+v,w〉=〈u,w〉+〈v,w〉, ∀w∈S(线性性2)4)〈u,u〉≥0，当且仅当u=O时①，〈u,u〉=0, (非负性)则称〈u,v〉为一种实内积运算(inner product). 定义了内积的线性空间称为实内积空间.应说明的是，将定义6.1加以扩展可在更一般的实数域ℂ上定义内积，区别只是将第1条性质改为共轭可交换性：〈u,v〉=〈v,u〉 .例如复向量的内积为: 〈u,v〉=u T v̅，可以验证它满足上述共轭可交换性. 下面只考虑实内积，但得到的结果都可以类似地推广到复内积空间. 另外，定义6.1的条件2还说明零元素与任意元素的内积均等于0.根据内积的线性性可推出：〈α1u1+α2u2,v〉=α1〈u1,v〉+α2〈u2,v〉,∀α1,α2∈ℂ,(6.1) 更一般地有：〈∑αj u j nj=1,v〉=∑αj〈u j,v〉nj=1,∀α1,⋯,αn∈ℂ.(6.2)这里主要考虑函数空间，则(6.2)式表明，线性组合函数(与另一函数作)内积等于(相应各个函数)内积的线性组合.可以规定一种依赖于内积运算的范数：‖u‖≡√〈u,u〉 .易知这种内积导出的范数满足范数定义的三个条件（见3.1.2节），详细证明过程留给读者思考. 应注意，在向量空间中，由内积导出的范数等同于向量的2-范数. 在实函数空间C[a,b]中，一般定义内积为〈u(x),v(x)〉=∫u(x)v(x)dxba,(6.3) 因此，由它导出的范数也等同于函数空间的2-范数.下面介绍与内积有关的两个重要定理.定理6.1：设S为实内积空间，∀u,v∈S，有：|〈u,v〉|2≤〈u,u〉∙〈v,v〉 .(6.4) 这是著名的柯西-施瓦茨不等式(Cauchy-Schwarz inequality).定理6.1的证明留给读者思考，若u,v为三维向量，也请思考该定理有什么几何含义？定理6.2：设S为实内积空间，u1,…,u n∈S，则格莱姆矩阵（Gram matrix）G=[〈u1,u1〉〈u2,u1〉⋯〈u n,u1〉〈u1,u2〉〈u2,u2〉⋯〈u n,u2〉⋮⋮⋱⋮〈u1,u n〉〈u2,u n〉⋯〈u n,u n〉](6.5)非奇异的充要条件是u1,…,u n线性无关.[证明] 首先要用到线性代数中的一个基本结论：矩阵G非奇异⟺det(G)≠0⟺齐次线性方程组Ga=0只有全零解.设向量a=[a1,…,a n]T，则方程Ga=0可写成：①这里用正体的字母O表示线性空间的零元素.∑a j 〈u j ,u k 〉nj=1=0,k =1,2,⋯,n (6.6)下面证明方程组(6.6)只有恒零解的充分必要条件是u 1,…,u n 线性无关. 先证必要性，即已知方程组(6.6)只有恒零解，要证u 1,…,u n 线性无关. 采用反证法，若u 1,…,u n 线性相关，即存在不全为0的一组系数{αj ,j =1,⋯,n}使∑αj u j n j=1=O ，则∑αj 〈u j ,u k 〉n j=1=〈∑αj u j nj=1,u k 〉=〈O,u k 〉=0,(k =1,…,n ),即这组{αj }是方程组(6.6)的解，与已知条件矛盾！再证明充分性，即已知u 1,…,u n 线性无关，要证方程组(6.6)只有全零解. 仍采用反证法，若方程组(6.6)存在不全为零的一组解{αj }，则∑αj 〈u j ,u k 〉n j=1=〈∑αj u j nj=1,u k 〉=0,k =1,…,n将上述方程中第k 个方程乘以αk ，累加所有方程得到，〈∑αj u j n j=1,∑αj u j nj=1〉=0 ,根据内积的定义，必有∑αj u j n j=1=O , 也就是说存在不全为0的一组{αj }j=1n 使∑αj u j n j=1=O ，这与u 1,…,u n 线性无关的已知条件矛盾！综上所述，完成了定理的证明.应注意，格莱姆矩阵是实对称矩阵，并且当u 1,…,u n 线性无关时，它是对称正定矩阵. 针对实函数空间C[a, b]，常常有权函数、加权内积的概念.定义6.2：若函数ρ(x )≥0,∀x ∈[a,b]，且满足1) ∫x k ρ(x )dx ba 存在，(k =0,1,…),2) 对非负连续函数g (x )，若∫g (x )ρ(x )dx =0b a 可推出g (x )≡0,则称ρ(x)为区间[a,b]上的权函数(weight function).关于权函数的定义，说明几点：● 定义中对连续性没有要求，即ρ(x )可能不是连续函数；第1个条件要求的是ρ(x )与多项式乘积为可积函数.● 定义中第2条件的意义不是很直观，较直观的一种等价形式为：不存在子区间(c,d )⊆[a,b]，使ρ(x )=0,∀x ∈(c,d )，即“权函数在[a,b]中任一子区间不恒为零”. ● 一般遇到的C [a,b ]中非负函数（一定有界、可积），若不在某一子区间恒为零，则都可作权函数.定义6.3：若ρ(x )为区间[a,b]上的权函数，则可定义C [a,b ]上的内积为：〈u (x ),v (x )〉=∫ρ(x )u (x )v (x )dx b a ,(6.7)并称其为加权内积(weighted inner product).容易验证加权内积满足一般内积的定义，并且常用的函数内积(6.3)式是加权内积的特例，其对应于权函数ρ(x )≡1的情况. 根据加权内积，也可以导出范数，这种范数可看成是广义的2-范数，其公式为：‖f(x)‖=[∫ρ(x )f 2(x )dx b a ]12⁄ .6.1.2 函数逼近的不同类型在函数逼近问题中，用简单函数p(x)来近似f(x)，并要求误差最小. 这里度量误差大小的标准是范数，采用不同范数时其问题的性质是不同的. 下面分两种情况作些讨论.1) ∞-范数考虑误差函数p (x )−f (x )的∞-范数，假设函数的定义域为[a, b]，则可设ε=‖p (x )−f (x )‖∞=max x∈[a,b ]|p (x )−f (x )| , 因此有−ε≤p (x )−f (x )≤ε,∀x ∈[a,b ]，即p (x )−ε≤f (x )≤p (x )+ε, ∀x ∈[a,b ]图6-3显示了函数p (x ),f (x ), 以及‖p (x )−f (x )‖∞之间的关系，从中可以看出，在∞-范数意义下的逼近要求使ε尽量小，也就是要p (x )在整个区间上“一致地”接近f (x ). 因此，采用∞-范数的函数逼近问题常称为最佳一致逼近.2) 1-范数和2-范数先看看误差函数p (x )−f (x )的1-范数，‖p (x )−f (x )‖1=∫|p (x )−f (x )|dx ba令A =‖p (x )−f (x )‖1，则它表示p (x )和f (x )两个函数曲线之间的面积（如图6-4所示）. 在1-范数意义下的逼近，要求使A 尽量小，也就是要p (x )与f (x )曲线之间的总面积尽量小，反映出这种逼近有整个区间上“平均”误差尽量小的含义（在某个子区间上误差可能很大）.2-范数的意义与1-范数大体上类似，由于它更容易处理，在实际的逼近问题中一般采用图6-3 函数p (x ),f (x ), 以及‖p (x )−f (x )‖∞之间的关系.图6-4 函数p (x ),f (x ), 以及‖p (x )−f (x )‖1之间的关系.2-范数. 这种逼近称为最佳平方逼近或最小二乘逼近(least squares fitting).从直观上看，采用∞-范数的最佳一致逼近效果更好一些，而最佳平方逼近具有平均误差最小的含义.除了度量误差函数可采用不同的范数，被逼近函数也可分为连续函数和表格函数两种情况. 表格函数就是仅在一系列离散自变量点上已知函数值的函数，可通过函数值组成的向量来刻画，有关逼近问题的求解有特殊的处理方法. 而在逼近函数类方面，多项式函数是最常用的一种. 下面给出魏尔斯特拉斯定理（Weierstrass Theorem ），它是用多项式函数进行逼近的一个重要依据.定理6.3：设f (x )∈C[a,b]，则对任何ϵ>0，总存在一个多项式P (x )，使‖P (x )−f (x )‖∞<ϵ在[a, b]上一致成立.该定理的证明已超出了本书的要求，因此不做讨论. 值得一提的是，若f (x )∈C[0,1]，伯恩斯坦多项式(Bernstein polynomial)②B n (f,x )=∑f (k )Q k (x )nk=0 , 其中Q k (x )=(n k)x k (1−x )n−k , 就是满足定理要求的多项式P (x ). 注意B n (f,x )为n 次多项式，并且可以证明，lim n→∞B n (f,x )=f(x)在[0, 1]上一致成立. 因此，C[0,1]中的任意函数都可以用伯恩斯坦多项式（一致）逼近到任意好的程度. 应注意，它一般不是多项式函数类ℙn 中的最佳一致逼近.最后说明一点，求最佳一致逼近多项式的方法比较复杂，感兴趣的读者请参考[4, 9]. 本章后面主要介绍求最佳平方逼近的方法，它有很广泛的应用.6.2 连续函数的最佳平方逼近为了记号的方便，在6.2节和6.3节的介绍中记函数的自变量为t.6.2.1 一般的法方程方法一. 问题描述假设对f (t )∈C [a,b ]进行函数逼近，逼近函数类Φ应是形式简单的函数类，比如多项式函数、三角函数、有理函数，等等，并且它是有限维的线性子空间. 设Φ=span {φ1(t ),…,φn (t )}，则Φ的任一元素可表示为：S (t )=Σj=1n x j φj (t ), (6.8)其中x 1,…,x n ∈ℝ.连续函数的最佳平方逼近问题就是求S (t )∈Φ，使 ‖S (t )−f (t )‖2达到最小值. 利用公式(6.8)以及2-范数的定义，上述问题等价于最小化F =‖S (t )−f (t )‖22=∫[Σj=1n x j φj (t )−f (t )]2dt b a .(6.9)F 是关于实系数x 1,x 2,…,x n 的多元函数，需求出F 的最小值对应的那组系数x 1,x 2,…,x n .二. 法方程方法下面推导如何求(6.9)式的最小值点. 为了记号简便，省略函数记号中的“(t )”，即直接② 由原苏联数学家伯恩斯坦（1880—1968）于1912年提出.f ̃=f (3)=f (2)−2v 2T f (2)v 2T v 2v 2=[ −4.2061330.399807−0.004750130.0009512830.00195269], 此时矩阵A 经变换为： R =A (3)=[ −2.236068−3.35410200.790569000000] . 根据算法6.3，需求解方程R 1x =b ，其中R 1=[−2.236068−3.35410200.790569],b =[−4.2061330.399807]. 解得：x =[1.12250.5057]T ，即拟合公式为y ̃=1.1225+0.5057t ，它与例6.6, 6.7得到的结果是一样的.根据表格函数与其函数值向量的对应关系可证明，算法6.3与通过Gram-Schmidt 正交化过程求最佳逼近函数的方法在数学上是等价的. 不同之处在于：前者不涉及正交函数族，直接得到原基函数对应的拟合系数；前者的主要计算是矩阵的QR 分解，它可通过Householder 变换或Givens 旋转变换等不同方法实现. 由于算法6.3直接利用矩阵的QR 分解的特点，它更易于实现和应用，而且稳定性比算法6.2好. 最后说明一点，若初始的表格函数φ1(t ),…,φn (t )线性相关，矩阵A 不是列满秩的，QR 分解也能进行，但得到的上三角阵R 1奇异. 可以证明，这种情况下有无穷多个最小二乘解，详细的讨论请参考[6].一. 问题背景1945年7月16日，美国科学家在新墨西哥州Los Alamos沙漠试爆了世界上第一颗原子弹，这一事件令全球震惊. 但在当时有关原子弹爆炸的任何资料都是保密的，而很多其他国家的科学家非常想知道这次爆炸的威力有多大.两年之后，美国政府首次公开了这次爆炸的录像带，而其他数据和资料仍然不被外界所知. 英国物理学家G. I. Taylor(1886 ~ 1975)通过研究原子弹爆炸的录像带，建立数学模型对爆炸所释放出的能量进行了估计，得到估计值与若干年后正式公布的爆炸能量21 kt 相当接近(1 kt 为1千吨TNT 炸药的爆炸能量). Taylor 是如何根据爆炸录像估计的呢？主要是通过测量爆炸形成的“蘑菇云”半径来进行估计的（如图(A)）. 因为爆炸产生的冲击波从中心点向外传播，爆炸的能量越大，在相同时间内冲击波传播得越远、蘑菇云的半径就越大. Taylor 通过图(A) 原子弹爆炸的蘑菇云.*t 的单位为ms, r 的单位为m.然后通过量纲分析法建立了蘑菇云半径r 与时间t 和爆炸能量E 的关系式，利用上述数据最后求出了爆炸的能量.二. 数学模型考虑到原子弹爆炸在极短的时间内释放出巨大的能量，蘑菇云半径r 主要与时间t 、爆炸能量E 、以及空气密度ρ等几个参数有关. 通过仔细分析这几个量的单位，采用量纲分析法得到如下的蘑菇云半径的近似表达式：r =(t 2E )15. 其中r , t , E 的单位分别为米(m), 秒(s)和焦耳，而空气密度ρ的值为1.25 (kg m 3⁄). 对这次原子弹爆炸来说，E 为一固定值，因此r 与t 2成正比. 图(B)是根据蘑菇云半径与对应时刻的数据画出的散点图，它大体反映了这个趋势. 接下来的问题是如何求未知的参数E .三. 求解过程首先，改写蘑菇云半径的公式为r =at b 的形式，通过测量数据拟合出参数a 和b ，来验证量纲分析法得到的公式. 要作线性最小二乘拟合，进一步改写公式为：lnr =lna +blnt . 根据测量数据我们得到lnr 和lnt 的数据，将它们的函数关系拟合为1次多项式，得到系数b =0.4094，其值与前面分析的结果2/5非常接近，从而验证了量纲分析得到的公式.为了更为准确地计算爆炸能量E ，将蘑菇云半径公式改写为：5lnr −2lnt =ln (E ) . 此时可根据测量数据得到5lnr −2lnt 对应的一组数据，将它拟合为0次多项式（常数），设得到拟合系数为c ，则E ≈ρ∙e c .根据此方法算出E ≈8.6418×1013，单位为焦耳，查表得知1kt=4.184×1012焦耳，因此爆炸能量约等于20.65 kt.6.4函数插值与拉格朗日插值法函数插值可看作一种“特殊”的函数逼近问题，其逼近采用的“度量”准则是要求在插值节点处误差函数的值为0. 本节先介绍关于插值(interpolation)的一些基本概念，然后讨论最简单的一种多项式插值——拉格朗日插值法.图(B) 蘑菇云半径与对应时刻的数据 rt个节点：x 0<x 1<⋯<x n 进行插值，只需将B −k k (x ),B −k+1k (x ),⋯,B n−1k (x )这n+k 个k 次B-样条函数进行组合. 可以证明，它们在区间[x 0,x n ]上的部分组成n+k 个线性无关的基函数. 因此，对于满足额外边界条件的[x 0,x n ]上的k 次样条函数，可唯一地用这些基函数的线性组合表示. 感兴趣地读者可以推导B i 3(x )的表达式，然后利用插值条件和边界条件列方程求这些基函数对应的系数，进而推导出三次样条插值函数的表达式. 这个计算过程将与上一小节的方法得到相同的结果.利用B-样条基函数，可得到确定和计算各阶样条插值的有效而稳定的方法. 此外，它在计算机图形学、几何建模，以及数值求解微分方程等领域都有广泛的应用.评述关于多项式逼近和插值问题的研究历史悠久，应用面也很广. 本章只讨论了一元函数的最佳平方逼近，更多的相关内容，包括多元函数的逼近、正交多项式等，可参考下述文献：● P . J. Davis, Interpolation and Approximation , Dover, 1975.● W. Cheney, Introduction to Approximation Theory , AMS Chelsea Publishing, 2nd edition,1998.● G. A. Baker, and P . R. Graves-Morris, Pade Approximations , Cambridge University Press,2nd edition, 1996.● W. Gautschi, “Orthogonal polynomials: Applications and computation,” Acta Numerica ,Vol. 5, pp. 45-119, 1996.最佳平方逼近的法方程方法在1795年由高斯提出. 格莱姆－斯密特正交化方法在1883年由格莱姆提出，1907年斯密特给出了现代算法. 在求解最小二乘问题中使用QR 分解方法，特别是使用Householder 变换的方法是在1965年由G. Golub ⑥提出的. 最小二乘方法是统计学的重要工具，也称为回归分析，很多常用的数据处理软件（比如微软公司的Excel 软件）都具有这个功能. 本章讨论的线性最小二乘问题实际上是一种最简化的形式，即假设待逼近函数是基函数的线性组合. 在实际应用中还常遇到非线性最小二乘问题，它属于非线性优化问题，见参考文献[6]及其中给出的更多文献. 另外，若考虑所有参量都带有随机误差的情形，则成为完全最小二乘问题，有关详细讨论见文献：● S. Van Huffel and J. Vandewalle, The Total Least Squares Problem , SIAM Press, 1991. 本章也没有讨论拟合的基函数可能线性相关的情况，这在实际中可能由于拟合模型的不合理或数值误差造成，它使得矩阵A 列不满秩. 此时最佳平方逼近解不唯一，要得到实际有用的一个逼近解，需采用列重排的QR 分解等技术，更多讨论参见文献[6]及其他文献.多项式插值问题历史非常悠久，牛顿、拉格朗日等都在这方法做出了很多贡献. 除了将函数值作为条件的插值问题，插值条件中包括各阶导数值的情况也常见于各种工程应用中. 目前，常用的文档编辑软件都已使用保形分段插值来绘制曲线，例如微软公司的Word 和Power Point 软件. 样条函数是1946年由Schoenberg 首先提出的，本章只讨论了一维数据的样条插值和B-样条函数，实际问题中还有高维的插值问题，尤其在计算机图形学中二维B-样条是一个重要的工具. 关于样条的参考文献主要有：● C. de Boor, A Practical Guide to Splines , Springer-Verlag, 2nd edition, 1984.● E. V. Shikin and A. I. Plis, Handbook on Splines for the User , CRC Press, 1995.最后，列表说明Matlab 中与本章讨论的函数逼近与插值有关的命令和功能.⑥ Gene H. Golub (1932-2007), 美国斯坦福大学计算机系教授，美国科学院、工程院、艺术与科学院三院院士，著名的数值计算专家，1996年出版的著作”Matrix Computations ” [21]被奉为矩阵计算领域的经典.线拟合与样条插值的功能.[本章知识点]: 连续函数的范数；内积及其性质；内积空间的格莱姆矩阵、及其非奇异的充要条件；权函数与加权内积；最佳一致逼近与最佳平方逼近的概念；法方程方法求连续函数的最佳平方逼近；最佳平方逼近的误差；正交函数族与Gram-Schimdit正交化过程；勒让德多项式；用正交函数族作最佳平方逼近；曲线拟合的线性最小二乘问题；线性最小二乘问题的矩阵描述；法方程方法解线性最小二乘问题；表格函数的线性无关性与相关性；利用矩阵的QR分解解线性最小二乘问题；插值的基本概念；范德蒙矩阵与多项式插值的存在唯一性；拉格朗日插值公式；拉格朗日插值余项公式；牛顿插值公式；差商的计算；牛顿插值余项公式；高次多项式插值的问题；分段线性插值；埃尔米特插值；分段三次埃尔米特插值；保形分段插值；三次样条插值及边界条件；三次样条插值的构造方法；三弯矩方程；几种插值的比较；B-样条函数的基本概念与性质.算法背后的历史：拉格朗日与插值法约瑟夫·路易斯·拉格朗日（Joseph-Louis Lagrange，1736年1月25日—1813年4月10日）是法国数学家、物理学家. 他在数学、力学和天文学三个领域中都有巨大的贡献，其中尤以数学方面的成就最为突出. 拉格朗日与同时代的勒让德（Legendre）、拉普拉斯（Laplace）并称为法国的3L.拉格朗日于1736年生于意大利西北部的都灵. 17岁时，开始专攻当时迅速发展的数学分析. 1756年，受欧拉的举荐，拉格朗日被任命为普鲁士科学院通讯院士. 1766年赴柏林任普鲁士科学院数学部主任，居住柏林达20年之久，这是他一生科学研究的鼎盛时期. 在此期间，他完成了著作《分析力学》. 1786年加入了巴黎科学院成立的研究法国度量衡统一问题的委员会，并出任法国米制委员会主任. 1795年建立了法国最高学术机构——法兰西研究院后，拉格朗。
数值分析第六章_数值插值方法
M n1 (n 1)!
n1 ( x)
说明：
n=1时，
R1 ( x)
1 2
f
( )2 (x)
1 2
f
( )(x
x0 )(x
x1)
n=2时，
( [x0 , x1])
R2 (x)
1 6
f
( )(x
x0 )(x
x1)(x
x2 )
( [x0 , x2 ])
,
x1,
Hale Waihona Puke xn)1
x1
x12
x1n
n
( xi
ni j1
xj)
1 xn xn2 xnn
因 xi x j (i j) 故上式不为0。
据Cramer法则，方程组解存在且唯一。 故Pn (x)存在且唯一。虽然直接求解上述方程组 可求得插值多项式，但繁琐复杂，一般不用。
得关于a0,a1,…,an的n+1阶线性方程组
a0 a1x0 a0 a1x1
an x0n an x1n
y0 y1
a0 a1xn an xnn yn
其系数行列式是Vandermonde行列式
1 x0 x02 x0n
V
( x0
jk jk
(j,k=0,1)
称l0 (x)及l1 (x)为线性插值基函数。
2. 抛物插值：n=2情形
假定插值节点为x0, x1, x2 ，求二次插值多项式 L2 (x)，使 L2(xj)=yj (j=0,1,2) y= L2 (x)的几何意义就是过 (x0, y0)，(x1, y1) ， (x2, y2)三点的抛物线。 采用基函数方法，设
图像插值算法总结
图像插值算法总结插值指的是利⽤已知数据去预测未知数据，图像插值则是给定⼀个像素点，根据它周围像素点的信息来对该像素点的值进⾏预测。
当我们调整图⽚尺⼨或者对图⽚变形的时候常会⽤到图⽚插值。
⽐如说我们想把⼀个4x4的图⽚，就会产⽣⼀些新的像素点（如下图红点所⽰），如何给这些值赋值，就是图像插值所要解决的问题, 图⽚来源常见的插值算法可以分为两类：⾃适应和⾮⾃适应。
⾃适应的⽅法可以根据插值的内容来改变（尖锐的边缘或者是平滑的纹理），⾮⾃适应的⽅法对所有的像素点都进⾏同样的处理。
⾮⾃适应算法包括：最近邻，双线性，双三次，样条，sinc，lanczos等。
由于其复杂度, 这些插值的时候使⽤从0 to 256 (or more) 邻近像素。
包含越多的邻近像素，他们越精确，但是花费的时间也越长。
这些算法可以⽤来扭曲和缩放照⽚。
⾃适应算法包括许可软件中的许多专有算法，例如：Qimage，PhotoZoom Pro和正版Fractals。
这篇博客通过opencv中cv.resize()函数介绍⼀些⾮⾃适应性插值算法cv2.resize(src, dsize[, dst[, fx[, fy[, interpolation]]]]) → dst其中interpolation的选项包括，图⽚来源我们主要介绍最近邻，线性插值，双三次插值三种插值⽅式，下图是对双三次插值与⼀些⼀维和⼆维插值的⽐较。
⿊⾊和红⾊/黄⾊/绿⾊/蓝⾊点分别对应于插值点和相邻样本。
点的⾼度与其值相对应。
图⽚来源于最近邻顾名思义最近邻插值就是选取离⽬标点最近的点的值(⿊点，原来就存在的点)作为新的插⼊点的值，⽤opencv进⾏图像处理时，根据srcX = dstX* (srcWidth/dstWidth)srcY = dstY * (srcHeight/dstHeight)得到来计算⽬标像素在源图像中的位置，dstY代表输出图Y的坐标，srcY代表原图Y的坐标，srcX、srcY同理。
多项式插值理论
第六章 多项式插值理论一、区间[a , b ]上的一般插值理论 (从有限维子空间出发的逼近方法)① 对无限维函数空间的一个元素f (x ) 进行逼近，关于f (x ) 的情况仅知道一部分(1、若干点的函数值或导数值已知； 2、满足一些控制方程)② 选择一个由固定基函数张成的有限维函数子空间 基函数性质: ⎪⎩⎪⎨⎧、线性无关、完备的条件、满足基本的函数已知321③ 选择n X 中的元素)()(~x P x f n 或，在一定的约束条件下，使)(~x f 良好的逼近()x f ,即 令)(~x f = n n c c c φφφ+++ 2211关于()x f 在插值区间上有不大的误差（包括一定的光滑性逼近）。
④ 良好逼近的判断ε<-f f ~e .g . Tchebycheff 范数，|| f || = |)(|max x f bx a ≤≤ 称为一致逼近。
⑤ 约束条件: （依据对()x f 的了解来确定）i / 插值约束)()(~i i x f x f = 1n i ≤≤ i x ∈ (a , b ) 且i x 互不相同；ii / 插值与光滑性混合约束(1)、 )()(~i i x f x f = 1k i ≤≤ i x ∈(a , b ) 且i x 互不相同(2)、 )()(~i i x f x f '=' 1k i ≤≤ i x ∈(a,b) 且互不相同(3)、 )(~x f 的二阶导数存在iii / 变分约束 （以下两种约束不再具有严格的插值含义，这里可能仅知道被插函数满足某些控制方程）依据|| f -f ~||在n X 中为最小的条件，即确定常数n c c c 21, 使f ~的解由下列形式的极小化问题得到：|| f -0f ~|| = min{|| f -f ~||：f ~n X ∈}Note ：这里的||·|| 不局限于切比雪夫范数和2－范数，可能是某种内积定义的范数；这也是固体力学求近似解的基本方法（如，有限元就是能量的变分）。
曲线插值
2 曲线的参数表示
曲线的参数方程为
x x(t)
y
y(t)
z z(t)
归一化处理：为了方便起见，可以将参数t的范围区 间规范化成[0，1]。
参数化表示比显式、隐式有更多的优点！ 参数化表示方式易于用矢量和矩阵运算，对曲率、 斜率等的计算也有别于传统方式。
3
3 插值与拟合
设已知某个函数关系 y f (x)在某些离散点上的函 数值：
11
4.1 拉格朗日插值
构造各个插值节点上的基函数 li(x)(i=0,1,…,n) 满足如下条件
xi
x0
x1
x2
L
xn
l0 (x)
1
0
0
L
0
l1 ( x)
0
L
L
1
0
L
L
L
0
LL
ln (x)
0
0
0
L
1
12
4.1 拉格朗日插值
求n次多项式lk(x)(i=0,1,…,n)， k = 0, 1,…, n
29
5.1 五点光滑法
设每两个数据点的三次多项式为：
y ax3 bx2 cx d
插值条件为：
yi f (xi )
yi1 f ( xi1)
y(
xi
)
ti
y( xi1) ti1
由此可以对每一个原始等高线的折线段求出一组a,b,c,d四个 参数，构建一个三次多项式，实现光滑。
30
5.1 五点光滑法
19
4.2 埃尔米特插值
一般来说，给定 m+1 个插值条件，就可以构造出一 个 m 次 Hermite 插值多项式
两个典型的 Hermite 插值
什么叫插值,图像插值
什么叫插值,图像插值 插值（Interpolation/resampling）是一种图像处理方法，它可以为数码图像增加或减少象素的数目。
某些数码运用插值的方法创造出象素比传感器实际能产生象素多的图像，或创造数码变焦产生的图像。
实际上，几乎所有的图像处理软件支持一种或以上插值方法。
图像放大后锯齿现象的强弱直接反映了图像处理器插值运算的成熟程度。
下面的例子是一幅106*40的图像放大成450%的效果：最接近原则插值（Nearest Neighbor Interpolation） 最接近原则插值是最简单的插值方法，它的本质就是放大象素。
新图像的象素颜色是原图像中与创造的象素位置最接近象素的颜色。
如果把原图像放大200%，1个象素就会被放大成（2*2）4个与原象素颜色相同的象素。
多数的图像浏览和编辑软件都会使用这种插值方法放大数码图像，因为这不会改变原图像的颜色信息，并且不会产生防锯齿效果。
同理，在实际放大照片中这种方法并不合适，因为这种插值会增加图像的可见锯齿。
双线性插值（Bilinear Interpolation） 在双线性插值中，新创造的象素值，是由原图像位置在它附近的（2 x -2）4个邻近象素的值通过加权平均计算得出的。
这种平均算法具有放锯齿效果，创造出来的图像拥有平滑的边缘，锯齿难以察觉。
双三次插值（Bicubic interpolation） 双三次插值是一种更加复杂的插值方式，它能创造出比双线性插值更平滑的图像边缘。
请读者留意下图中的眼睫毛部分，在这个地方，软件通过双三次插值创造了一个象素，而这个象素的象素值是由它附近的（4 x 4）个邻近象素值推算出来的，因此精确度较高。
双三次插值方法通常运用在一部分图像处理软件、打印机驱动程序和数码相机中，对原图像或原图像的某些区域进行放大。
Adobe Photoshop CS 更为用户提供了两种不同的双三次插值方法：双三次插值平滑化和双三次插值锐化。
第六章 插值法
返回
前进
y n (x)
y f (x)
y0 x0 x1 y1 x2 y2 yn x xn
因此，所谓插值，即是在x0,x1,…,xn中任意插入一个x， 要求对应的f (x)，具体做法是按上述方法构造n(x)以 n(x)近似f (x)。 -
返回
前进
插值法是求函数值的一种逼近方法，是数值分析中的 基本方法之一，作为基础，后面微分，积分，微分方程在 进行离散化处理时，要用到，作为一种逼近方法，本身也 有广泛的应用价值，如在拱桥建设中，拱轴，拱腹的设计 节点与具体施工设计点常常可能不重合。如图5-2所示。
n
(5 -1)
( 5 - 2)
使其满足在给定点处与f(x)相同，即满足插值条件：
(i 0,1,2,, n)
n(x)称为插值多项式，xi(i=0,1,2,…,n)称为插值节点， [a,b]称为插值区间。
从几何上看（如图5-1所示），n次多项式插值就是 过n+1个点yi = f (xi)(i=0,1,…,n)，作一条多项式曲线 y = (x)近似曲线y = f (x) :
a
积分I是很困难的，构造近似函数使积分容易计算，并且 使之离散化能上机计算求出积分I，都要用到插值逼近。
返回
代数插值
前进
解决上述问题的方法有两类：一类是对于一组离 散点(xi,f (xi)) (i = 0,1,2,…,n)，选定一个便于计算的函 数形式(x)，如多项式，分段线性函数，有理式，三 角函数等，要求(x)通过点(xi)=f (xi) (i = 0,12,…,n), 由此确定函数(x)作为f (x)的近似。这就是插值法。 另一类方法在选定近似函数的形式后，不要求近似 函数过已知样点，只要求在某种意义下它在这些点 上的总偏差最小。这类方法称为曲线（数据）拟合 法，将在下一章介绍。 本章主要讨论构造插值多项式的几种常用的方法及 其误差 用插值法求函数的近似表达式时，首先要选定 函数的形式。可供选择的函数很多，常用的是多项式 函数。因为多项式函数计算简便，只需用加、减、乘 等运算，便于上机计算，而且其导数与积分仍为多项式。
常用三种图像插值算法
常见图像插值算法只有3种么？电脑摄像头最高只有130万像素的，800万是通过软件修改的。
何为数码插值(软件插值)插值(Interpolation)，有时也称为“重置样本”，是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩。
简单地说，插值是根据中心像素点的颜色参数模拟出周边像素值的方法，是数码相机特有的放大数码照片的软件手段。
一、认识插值的算法“插值”最初是电脑术语，后来引用到数码图像上来。
图像放大时，像素也相应地增加，但这些增加的像素从何而来？这时插值就派上用场了。
插值就是在不生成像素的情况下增加图像像素大小的一种方法，在周围像素色彩的基础上用数学公式计算丢失像素的色彩(也有些相机使用插值，人为地增加图像的分辨率)。
所以在放大图像时，图像看上去会比较平滑、干净。
但必须注意的是插值并不能增加图像信息。
以图1为原图(见图1)，以下是经过不同插值算法处理的图片。
1.最近像素插值算法最近像素插值算法(Nearest Neighbour Interpolation)是最简单的一种插值算法，当图片放大时，缺少的像素通过直接使用与之最接近的原有像素的颜色生成，也就是说照搬旁边的像素，这样做的结果是产生了明显可见的锯齿(见图2)。
2.双线性插值算法双线性插值算法(Bilinear Interpolation)输出的图像的每个像素都是原图中四个像素(2×2)运算的结果，这种算法极大程度上消除了锯齿现象(见图3)。
3.双三次插值算法双三次插值算法(Bicubic Interpolation)是上一种算法的改进算法，它输出图像的每个像素都是原图16个像素(4×4)运算的结果(见图4)。
这种算法是一种很常见的算法，普遍用在图像编辑软件、打印机驱动和数码相机上。
4.分形算法分形算法(Fractal Interpolation)是Altamira Group提出的一种算法，这种算法得到的图像跟其他算法相比更清晰、更锐利(见图5)。
图像的基本运算
图像的基本运算图像的基本运算包括以下几类：图像的点运算；图像的代数运算；图像的几何运算；图像的逻辑运算和图像的插值。
下面将依次介绍这几种运算。
一、点运算点运算是指对一幅图像中每个像素点的灰度值进行计算的方法。
点运算通过对图像中每个像素值进行计算，改善图像显示效果的操作，也称对比度增强，对比度拉伸，灰度变换，可以表示为B(x,y)=f(A(x,y))。
这是一种像素的逐点运算，是原始图像与目标图像之间的映射关系，不改变图像像素的空间关系。
可以提高图像的对比度，增加轮廓线等。
可分为：（1）线性点运算：输出灰度级与输入灰度级之间呈线性关系。
（2）非线性点运算：输出灰度级与输入灰度级之间呈非线性关系。
二、代数运算代数运算是指将两幅或多幅图像通过对应像素之间的加、减、乘、除运算得到输出图像的方法。
对于相加和相乘的情形，可能不止有两幅图像参加运算。
如果记A(x,y)和B(x,y)为输入图像，C(x,y)为输出图像。
那么，四种代数运算的数学表达式如下：（1）C(x,y)=A(x,y)+B(x,y)加法运算可以实现以下两个目的：1.1去除叠加性随机噪声；1.2生成图像叠加效果。
（2）C(x,y)=A(x,y)-B(x,y)减法运算可以实现以下两个目的：2.1消除背景影响；2.2检查同一场景两幅图像之间的变化。
（3）C(x,y)=A(x,y)*B(x,y)乘法运算可以实现以下两个目的：3.1图像的局部显示；3.2图像的局部增强。
（4）C(x,y)=A(x,y)/B(x,y)乘法运算可以实现以下三个目的：4.1遥感图像的处理中；4.2消除图像数字化设备随空间变化的影响。
4.3校正成像设备的非线性影响。
还可以通过适当的组合形成涉及几幅图像的复合代数运算。
三、几何运算几何运算就是改变图像中物体对象（像素）之间的空间关系。
从变换性质来分，几何变换可以分为图像的位置变换（平移、镜像、旋转）、形状变换（放大、缩小）以及图像的复合变换等。
三边滤波的图像插值
三边滤波的图像插值Junhua LiuInstitute for Signals and Information ProcessingLanzhou University, Lanzhou, China 730000文摘—本文中，在传统的图像插值方法处理问题会出现轮廓锯齿和模糊文物问题，所以提出了一种新的图像插值方法称为三边滤波插值法。
该方法利用平滑不同自然图像轮廓的行为，并试图通过抑制边缘凸角落像素克服轮廓锯齿问题。
首先给出了一个初步估计带有适当的重量未知像素的四个最近的邻域，这是由相应的边缘信息和空间距离决定。
然后，给出以锐利的边缘，基于相似性进一步改进的权重的介绍，在此基础上，最后估计未知像素获得。
仿真结果表明相比传统方法该方案产生更少的模糊文物和锐利的边缘，仍然保留简单和效率，方案也更先进。
1介绍图像插值是从低分辨率图像产生高分辨率图像的过程。
它有许多应用，例如在医疗成像和取证，图像马赛克，图像的几何变换，视频尺寸转换，等等。
一个插值方案的性能关键因素在于锐利的边缘和自由的文物，计算效率的另一个重要因素[ 2]。
许多插值方案被提出，它们大体可以分为线性和非线性的方法。
线性方法估算值未知像素的线性组合，他们的邻居是已知的。
这些方法包括最近邻插值，双线性插值，三次样条插值，插值内核等，并被广泛地用于计算效率[ 1]，[ 4]，[ 5]，[ 11]。
然而，由于他们是基于空间不变的功能，他们无法捕捉快速演变的统计边缘，因而往往产生模糊和恼人的锯齿状的插值图像[ 6]。
非线性插值计划用来提高主观质量的插值图像，他们通常是根据当地的特点给出的图像。
[ 2]，[ 7]使用局部边缘方向信息估计未知的像素。
[ 9]采用了插值内核，以便于适应当地方向等照度线。
[ 6]和[ 12]模型图像的像素值的分段自回归（杆）的过程，并利用标准参数捕获图像局部结构。
不同的是，[ 6 ]估计一个像素每次在[ 12]估计像素块。
数值分析第六章插值法PPT学习教案
j0 jk
n
( xk
xj)
n
j0 jk
x xj xk x j
j0
jk
称 lk (x) 为关于基点 xi 的n次插值基函数(i=0,1,…,n)
第10页/共81页
以n+1个n次基本插值多项式 lk (x)(k 0,1,, n) 为基础,就能直接写出满足插值条件
P(xi ) f (xi ) (i 0,1,2,, n) 的n次代数插值多项式。
l0 (x) 9l1(x) 23l2 (x) 3l3 (x)
11 x3 45 x 2 1 x 1
第六章 函数插值
第六章 函数插值实践中常有这样的问题：由实验得到某一函数y = f (x )在一系列点x 0, x 1,…, x n 处的值y 0, y i ,…, y n ，其函数的解析表达式是未知的，需要构造一个简单函数P (x )作为y = f (x )的近似表达式；或者y = f (x )虽有解析式，但计算复杂，不便于使用，需要用一个比较简单且易于计算的函数P (x )去近似代替它；本章所介绍的插值法就是建立这种近似公式的基本方法。
§1 代数插值 设已知某个函数关系y = f (x )在某些离散点上的函数值：nn y y y y yx x x x x 21210 （6.1）插值问题就是根据这些已知数据来构造函数y = f (x )的一种简单的近似表达式，以便于计算点i x x 的函数值)(x f ，或计算函数的一阶、二阶导数值。
一种常用的方法就是从多项式中选一个P n (x )，使得n i y x P i i n ,,2,1,0,)(（6.2）作为f (x )的近似。
因为多项式求值方便，且还有直到n 阶的导数。
我们称满足关系（6.2）的函数P n (x )为f (x )的一个插值函数，称x 0, x 1,…, x n 为插值节点，并称关系（6.2）为插值原则。
这种用代数多项式作为工具来研究插值的方法叫做代数插值。
设 x 0 < x 1< …< x n记a = x 0, b = x n ，则 [a, b] 为插值区间。
插值多项式存在的唯一性： 设所要构造的插值多项式为：n n n x a x a x a a x P 2210)(由插值条件n i y x P ii n ,,1,0)(得到如下线性代数方程组：n n n n n n nn nya x a x a y a x a x a y a x a x a101111000100111 此方程组的系数行列式为ni j j in nnnnn x xx x x x x x x x x D 021211020)(111此为范得蒙行列式，在线性代数课中，已经证明当j i x x ，;,2,1n i n j ,2,1 时，D 0，因此，P n (x )由a 0, a 1,…, a n 唯一确定。
图像插值算法
图像插值算法
常用的有以下三者：
1、最邻近元法
这是最简单的一种插值方法，不需要计算，在待求象素的四邻象素中，将距离待求象素最近的邻象素灰度赋给待求象素。
设i+u，j+v （i，j为正整数，u，v为大于零小于1的小数，下同）为待求象素坐标，则待求象素灰度的值f（i+u，j+v）。
2、双线性内插法
双线性内插法是利用待求象素四个邻象素的灰度在两个方向上作线性内插。
3、三次内插法
该方法利用三次多项式S（x）求逼近理论上最佳插值函数sin （x）/x，待求像素（x，y）的灰度值由其周围16个灰度值加权内插得到。
图像插值算法及其实现
图像插值算法及其实现sensor、codec、display device都是基于pixel的，⾼分辨率图像能呈现更多的detail，由于sensor制造和chip的限制，我们需要⽤到图像插值（scaler/resize）技术，这种⽅法代价⼩，使⽤⽅便。
同时，该技术还可以放⼤⽤户希望看到的感兴趣区域。
图像缩放算法往往基于插值实现，常见的图像插值算法包括最近邻插值（Nearest-neighbor）、双线性插值（Bilinear）、双⽴⽅插值（bicubic）、lanczos插值、⽅向插值（Edge-directed interpolation）、example-based插值、深度学习等算法。
插值缩放的原理是基于⽬标分辨率中的点，将其按照缩放关系对应到源图像中，寻找源图像中的点（不⼀定是整像素点），然后通过源图像中的相关点插值得到⽬标点。
本篇⽂章，我们介绍Nearest-neighbor和Bilinear插值的原理及C实现。
插值算法原理如下：1. Nearest-neighbor最近邻插值，是指将⽬标图像中的点，对应到源图像中后，找到最相邻的整数点，作为插值后的输出。
如下图所⽰，P为⽬标图像对应到源图像中的点，Q11、Q12、Q21、Q22是P点周围4个整数点，Q12与P离的最近，因此P点的值等于Q12的值。
这⾥写图⽚描述由于图像中像素具有邻域相关性，因此，⽤这种拷贝的⽅法会产⽣明显的锯齿。
2. Bilinear双线性插值使⽤周围4个点插值得到输出，双线性插值，是指在xy⽅法上，都是基于线性距离来插值的。
如图1，⽬标图像中的⼀点对应到源图像中点P(x,y)，我们先在x⽅向插值：然后，进⾏y⽅向插值：可以验证，先进⾏y⽅向插值再进⾏x⽅向插值，结果也是⼀样的。
值得⼀提的是，双线性插值在单个⽅向上是线性的，但对整幅图像来说是⾮线性的。
3. C实现使⽤VS2010，⼯程包含三个⽂件，如下：main.cpp#include <string.h>#include <iostream>#include "resize.h"int main(){const char *input_file = "D:\\simuTest\\teststream\\00_YUV_data\\01_DIT_title\\data.yuv"; //absolute pathconst char *output_file = "D:\\simuTest\\teststream\\00_YUV_data\\01_DIT_title\\data_out2.yuv"; //absolute pathint src_width = 720;int src_height = 480;int dst_width = 1920;int dst_height = 1080;int resize_type = 1; //0:nearest, 1:bilinearresize(input_file, src_width, src_height, output_file, dst_width, dst_height, resize_type);return 0;}resize.cpp#include "resize.h"int clip3(int data, int min, int max){return (data > max) ? max : ((data < min) ? min : data);if(data > max)return max;else if(data > min)return data;elsereturn min;}//bilinear takes 4 pixels (2×2) into account/** 函数名： bilinearHorScaler* 说明：⽔平⽅向双线性插值* 参数：*/void bilinearHorScaler(int *src_image, int *dst_image, int src_width, int src_height, int dst_width, int dst_height){double resizeX = (double)dst_width / src_width;for(int ver = 0; ver < dst_height; ++ver){for(int hor = 0; hor < dst_width; ++hor){double srcCoorX = hor / resizeX;double weight1 = srcCoorX - (double)((int)srcCoorX);double weight2 = (double)((int)(srcCoorX + 1)) - srcCoorX;double dstValue = *(src_image + src_width * ver + clip3((int)srcCoorX, 0, src_width - 1)) * weight2 + *(src_image + src_width * ver + clip3((int)(srcCoorX + 1), 0, src_width - 1)) * weight1; *(dst_image + dst_width * ver + hor) = clip3((uint8)dstValue, 0, 255);}}}/** 函数名： bilinearVerScaler* 说明：垂直⽅向双线性插值* 参数：*/void bilinearVerScaler(int *src_image, int *dst_image, int src_width, int src_height, int dst_width, int dst_height){double resizeY = (double)dst_height / src_height;for(int ver = 0; ver < dst_height; ++ver){for(int hor = 0; hor < dst_width; ++hor){double srcCoorY = ver / resizeY;double weight1 = srcCoorY - (double)((int)srcCoorY);double weight2 = (double)((int)(srcCoorY + 1)) - srcCoorY;double dstValue = *(src_image + src_width * clip3((int)srcCoorY, 0, src_height - 1) + hor) * weight2 + *(src_image + src_width * clip3((int)(srcCoorY + 1), 0, src_height - 1) + hor) * weight1; *(dst_image + dst_width * ver + hor) = clip3((uint8)dstValue, 0, 255);}}}/** 函数名： yuv420p_NearestScaler* 说明：最近邻插值* 参数：*/void nearestScaler(int *src_image, int *dst_image, int src_width, int src_height, int dst_width, int dst_height){double resizeX = (double)dst_width /src_width; //⽔平缩放系数double resizeY = (double)dst_height / src_height; //垂直缩放系数int srcX = 0;int srcY = 0;for(int ver = 0; ver < dst_height; ++ver) {for(int hor = 0; hor < dst_width; ++hor) {srcX = clip3(int(hor/resizeX + 0.5), 0, src_width - 1);srcY = clip3(int(ver/resizeY + 0.5), 0, src_height - 1);*(dst_image + dst_width * ver + hor) = *(src_image + src_width * srcY + srcX);}}}void resize(const char *input_file, int src_width, int src_height, const char *output_file, int dst_width, int dst_height, int resize_type){//define and init src bufferint *src_y = new int[src_width * src_height];int *src_cb = new int[src_width * src_height / 4];int *src_cr = new int[src_width * src_height / 4];memset(src_y, 0, sizeof(int) * src_width * src_height);memset(src_cb, 0, sizeof(int) * src_width * src_height / 4);memset(src_cr, 0, sizeof(int) * src_width * src_height / 4);//define and init dst bufferint *dst_y = new int[dst_width * dst_height];int *dst_cb = new int[dst_width * dst_height / 4];int *dst_cr = new int[dst_width * dst_height / 4];memset(dst_y, 0, sizeof(int) * dst_width * dst_height);memset(dst_cb, 0, sizeof(int) * dst_width * dst_height / 4);memset(dst_cr, 0, sizeof(int) * dst_width * dst_height / 4);//define and init mid bufferint *mid_y = new int[dst_width * src_height];int *mid_cb = new int[dst_width * src_height / 4];int *mid_cr = new int[dst_width * src_height / 4];memset(mid_y, 0, sizeof(int) * dst_width * src_height);memset(mid_cb, 0, sizeof(int) * dst_width * src_height / 4);memset(mid_cr, 0, sizeof(int) * dst_width * src_height / 4);uint8 *data_in_8bit = new uint8[src_width * src_height * 3 / 2];memset(data_in_8bit, 0, sizeof(uint8) * src_width * src_height * 3 / 2);uint8 *data_out_8bit = new uint8[dst_width * dst_height * 3 / 2];memset(data_out_8bit, 0, sizeof(uint8) * dst_width * dst_height * 3 / 2);FILE *fp_in = fopen(input_file,"rb");if(NULL == fp_in){//exit(0);printf("open file failure");}FILE *fp_out = fopen(output_file, "wb+");//data readfread(data_in_8bit, sizeof(uint8), src_width * src_height * 3 / 2, fp_in);//Y componentfor(int ver = 0; ver < src_height; ver++){for(int hor =0; hor < src_width; hor++){src_y[ver * src_width + hor] = data_in_8bit[ver * src_width + hor];}}//c component YUV420Pfor(int ver = 0; ver < src_height / 2; ver++){for(int hor =0; hor < src_width / 2; hor++){src_cb[ver * (src_width / 2) + hor] = data_in_8bit[src_height * src_width + ver * src_width / 2 + hor];src_cr[ver * (src_width / 2) + hor] = data_in_8bit[src_height * src_width + src_height * src_width / 4 + ver * src_width / 2 + hor];}}//resizeif(0 == resize_type){nearestScaler(src_y, dst_y, src_width, src_height, dst_width, dst_height);nearestScaler(src_cb, dst_cb, src_width / 2, src_height / 2, dst_width / 2, dst_height / 2);nearestScaler(src_cr, dst_cr, src_width / 2, src_height / 2, dst_width / 2, dst_height / 2);}else if(1 == resize_type){bilinearHorScaler(src_y, mid_y, src_width, src_height, dst_width, src_height);bilinearHorScaler(src_cb, mid_cb, src_width / 2, src_height / 2, dst_width / 2, src_height / 2);bilinearHorScaler(src_cr, mid_cr, src_width / 2, src_height / 2, dst_width / 2, src_height / 2);bilinearVerScaler(mid_y, dst_y, dst_width, src_height, dst_width, dst_height);bilinearVerScaler(mid_cb, dst_cb, dst_width / 2, src_height / 2, dst_width / 2, dst_height / 2);bilinearVerScaler(mid_cr, dst_cr, dst_width / 2, src_height / 2, dst_width / 2, dst_height / 2);}else{nearestScaler(src_y, dst_y, src_width, src_height, dst_width, dst_height);nearestScaler(src_cb, dst_cb, src_width / 2, src_height / 2, dst_width / 2, dst_height / 2);nearestScaler(src_cr, dst_cr, src_width / 2, src_height / 2, dst_width / 2, dst_height / 2);}//data writefor(int ver = 0; ver < dst_height; ver++){for(int hor =0; hor < dst_width; hor++){data_out_8bit[ver * dst_width + hor] = clip3(dst_y[ver * dst_width + hor], 0, 255);}}for(int ver = 0; ver < dst_height / 2; ver++){for(int hor = 0; hor < dst_width / 2; hor++){data_out_8bit[dst_height * dst_width + ver * dst_width / 2 + hor] = clip3(dst_cb[ver * (dst_width / 2) + hor], 0, 255);data_out_8bit[dst_height * dst_width + dst_height * dst_width / 4 + ver * dst_width / 2 + hor] = clip3(dst_cr[ver * (dst_width / 2) + hor], 0, 255); }}fwrite(data_out_8bit, sizeof(uint8), dst_width * dst_height * 3 / 2, fp_out);delete [] src_y;delete [] src_cb;delete [] src_cr;delete [] dst_y;delete [] dst_cb;delete [] dst_cr;delete [] mid_y;delete [] mid_cb;delete [] mid_cr;delete [] data_in_8bit;delete [] data_out_8bit;fclose(fp_in);fclose(fp_out);}resize.h#ifndef RESIZE_H#define RESIZE_H#include <stdio.h>#include <string.h>typedef unsigned char uint8;typedef unsigned short uint16;int clip3(int data, int min, int max);void bilinearHorScaler(int *src_image, int *dst_image, int src_width, int src_height, int dst_width, int dst_height);void bilinearVerScaler(int *src_image, int *dst_image, int src_width, int src_height, int dst_width, int dst_height);void nearestScaler(int *src_image, int *dst_image, int src_width, int src_height, int dst_width, int dst_height);void resize(const char *input_file, int src_width, int src_height, const char *output_file, int dst_width, int dst_height, int resize_type);#endif效果⽐较将720x480分辨率图像放⼤到1080p，1：1截取局部画⾯如下，左边是最近邻放⼤的效果，右边是双线性效果，可以看到，双线性放⼤的锯齿要明显⽐最近邻⼩。
数值分析
这里h为常数，称为步长， , 和 分别为向前, 向后和中心差分算子.
以向前差分为例 利用一阶差分可定义二阶差分
2 f k f k 1 f k f k 2 2 f k 1 f k .
还可以定义m阶差分 m f k m 1 f k 1 m 1 f k ;
xi ƒ(xi) x0 x1 x2 x3 xn ƒ(x0) ƒ(x1) ƒ(x2) ƒ(x3) ƒ(xn) 一阶 均差 二阶均差 三阶均差 n阶均差
ƒ[x0, x1] ƒ[x1, x2] ƒ[x0, x1, x2] ƒ[x2, x3] ƒ[x1, x2, x3] ƒ[x0, x1, x2, x3] ƒ[xn-1, xn] ƒ[xn-2, xn-1, xn] ƒ[xn-3, xn-2, xn-1, xn] ƒ[x0, x1,…, xn]
0 ( x) 1 1 ( x ) x x0 2 ( x ) ( x x0 )( x x1 )
LL n ( x ) ( x x0 )( x x1 ) L ( x xn1 )
当增加一个节点xn+1时，只需加上基函数
n1 ( x xi ) 即可.
n n n j
f k ( I E ) f k (1)
1 2
n1
………… f [ x, x0 , ... , xn1 ] f [ x0 , ... , xn ] ( x xn ) f [ x, x0 , ... , xn ]
1 + (x x0) 2 + … … + (x x0)…(x xn1)
n1
f [ x, x0 ] f [ x0 , x1 ] ( x x1 ) f [ x, x0 , x1 ]
第6章2DEM内插方法与数据管理
逐点内插方法 以每一待定点为中心，定义一 个局部函数去拟合周围的数据 点。逐点内插法十分灵活，精 度较高，计算方法简单又不需 很大的计算机内存，但计算速 度可能比其它方法慢。
移动曲面拟合法
(l) 建立局部坐标
对DEM每一个格网点，将坐标原点 移至该DEM格网点P（Xp，Yp）
Xi Xi X Yi Yi Y p
整型量存贮
Zi INT (Zi Z0 ) 10 0.5
m
将高程数据减去一常数Z0
差分映射 相邻数据间的增量，数据范围较小， 可以利用一个字节存贮一个数据，使 数据压缩至原有存贮量的近四分之一
Z 0 Z 0
Z 0 Z1 Z 2 Z n 1 1 0 0
压缩编码
当根据各数出现的概率设计一定的编码， 用位数（bit）最短的码表示出现概率最 大的数，出现概率较小数用位数较长的 码表示，则每一数据所占的平均位数比 原来的固定位数（16或8）小 数据的平均最小 位数可用信息论 中熵的定义计算
H (d1d 2 d n ) pk log2 pk
A
误差方程式 若 A点是已知高程点，作为观测值， 以格网高程Zi,j…作为待定的未知数
vA (1 X )(1 Y )Zi, j X (1 Y )Zi 1, j
(1 X )YZi, j 1 XYZi 1, j 1 Z A
虚拟 1 n
pi Z i pi
i 1
多面函数法DEM内插
“任何一个圆滑的数学表面总是可以用一 系列有规则的数学表面的总和，以任意的 精度进行逼近。”也就是一个数学表面上 某点（X，Y）处高程Z的表达式为:
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
f(n+a)=(1-a)f(n)+af(n+1), 0<a<1
Note: when a=0.5, we simply have the average of two
x x0 x x1 L( x ) y0 y1 x0 x1 x1 x0
Numerical Examples
1D Zero-order (Replication)
f(n)
n
f(x) x
1D First-order Interpolation (Linear)
f(n)
n
f(x) x
Linear Interpolation Formula
Heuristic: the closer to a pixel, the higher weight is assigned Principle: line fitting to polynomial fitting (analytical formula) f(n) f(n+a) f(n+1)
Scenario I: Resolution Enhancement
Low-Res.
High-Res
High Resolution Imaging Science Experiment
High Resolution Imaging Science Experiment
/mro/mission/instruments/hirise/
%%%%%%%%%%%%%%%%%%%%%%%%%%% % How to implement image interpolation by yourself? % %%%%%%%%%%%%%%%%%%%%%%%%%%%% % zero-order interpolation (replication) [M,N]=size(x); z0=zeros(2*M,2*N); z0(1:2:2*M,1:2:2*N)=x; h=[1 1;1 1]; z0=filter2(h,z0); % display interpolated image figure,imshow(z0,[]);
% zero-order interpolation (replication) y0=interp2(x,'nearest');
% display interpolated image figure,imshow(y0,[]);
% first-order interpolation (bilinear interpolation) y1=interp2(x,'linear'); % display interpolated image figure,imshow(y1,[]); % check bilinear property x(1:2,1:2) y1(1:5,1:5) % bicubic interpolation y2=interp2(x,'cubic'); % display interpolated image figure,imshow(y2,[]);
% read the PGM image into MATLAB x=imread('fl_orig.pgm'); % display the image imshow(x); % convert image to double format x=double(x); whos x
% read the instruction of interp2 help interp2 (or imresize)
Pixel Replication
low-resolution image (100×100) high-resolution image (400×400)
Bilinear Interpolmage (100×100)
high-resolution image (400×400)
Scenario II: Image Inpainting
Non-damaged
Damaged
Image warping is the process of digitally manipulating an image such that any shapes portrayed in the image have been significantly distorted. Warping may be used for correcting image distortion as well as for creative purposes (e.g., morphing[1]).
f(n)=[0,120,180,120,0]
Interpolate at 1/2-pixel
f(x)=[0,60,120,150,180,150,120,60,0], x=n/2
Interpolate at 1/3-pixel
f(x)=[0,20,40,60,80,100,120,130,140,150,160,170,180,…], x=n/6
Graphical Interpretation of Interpolation at Half-pel
row
column
f(m,n)
g(m,n)
最近邻插值
Numerical Examples
a
zero-order a a a a c c c c b b d d b b d d
b
c
d first-order a (a+b)/2 (a+c)/2 (a+b+c+d)/4 c (c+d)/2 b (b+d)/2 d
真实值: 10.7238
Introduction
• What is image interpolation?
– An image f(x,y) tells us the intensity values at the integral lattice locations, i.e., when x and y are both integers – Image interpolation refers to the guess of intensity values at missing locations, i.e., x and y can be arbitrary – Note that it is just a guess (Note that all sensors have finite sampling distance)
– We want GOOD images
• If some block of an image gets damaged during the transmission, we want to repair it
– We want COOL images
• Manipulate images digitally can render fancy artistic effects as we often see in movies
已知函数表求满足:
L(x0)=y0, L(x1)=y1 的线性函数 L(x)
过两点直线方程
x f (x )
x0 y0
x1 y1
y1 y0 L( x ) y0 ( x x0 ) x1 x0
引子 求 115 的近似值(函数值: 10.7238)
11 10 115 10 (115 100) 10.7143 121 100
1D Third-order Interpolation (Cubic)*
f(n) n
f(x) x
1D Third-order Interpolation (Cubic)*
http://www.paulinternet.nl/?page=bicubic
压缩的概念:
观测的离散数据可以想象成现实中无穷多 信息的代表。通过给定数据求出插值函数意味 着用简单的规则代替无穷多信息。尽管期待这 种简单规则精确地反映实际情况是不现实的, 但是它可以充分接近实际。这一类压缩是有损 的压缩, 即它会产生误差。
From 1D to 2D
Engineers’ wisdom: divide and conquer 2D interpolation can be decomposed into two sequential 1D interpolations. The ordering does not matter (row-column = column-row) Such separable implementation is not optimal but enjoys low computational complexity “If you don’t know how to solve a problem, there must be a related but easier problem you know how to solve. See if you can reduce the problem to the easier one.” - rephrased from G. Polya’s “How to Solve It”
% check replication property x(1:4,1:4) z0(1:9,1:9)
% first-order interpolation (bilinear interpolation) [M,N]=size(x); z1=zeros(2*M,2*N); z1(1:2:2*M,1:2:2*N)=x; h=[1 2 1;2 4 2;1 2 1]/4; z1=filter2(h,z1); % display interpolated image figure,imshow(z1,[]);