有理样条

2024-10-20

有理样条（共4篇）

有理样条篇1

经典样条理论多是在多项式空间进行的,文献[1]以实例说明了有理函数逼近效果更佳,因此,有理函数逼近作为多项式逼近的推广多年来愈来愈引起人们的关注,近年来有理逼近已成为数值逼近,函数逼近论,CAGD等领域的研究热点,并已在无限维动力系统中得以应用[2]。

有理样条函数是多项式样条的一种自然推广,特别是对带极点的函数,现给出的带极点有理样条函数逼近是合适的。文献[3]介绍了1980年以来,有理样条的一些结果,但这些结果一般还将导致非线性方程组,计算复杂。为了避免有理样条插值的这个困难,基于四维带极点有理函数空间U4建立一种有理样条函数。

1 带极点有理函数空间

定义1[4] 设[a,b]是一个实紧区间,给定 $p_{1}, p_{2}, \dots, p_{d} \in \bar{R} [a, b]$ 和一组正整数n1,n2,…,nd,定义函数组μi(1≤i≤n)如下(其中 $\bar{R} = R \cup {\infty}$ , $\sum_{i = 1}^{d}$ ni=n):

$μ_{k j} (x) = {\begin{matrix} x^{j - 1} ‚ & (1 \leq j \leq n_{1}, k = 1) \\ \frac{1}{(x - p_{k})^{j}} ‚ & (1 \leq j \leq n_{k}, k = 2, \dots, d) \end{matrix} (1)$

记线性包Rpn=span{μkj} $_{k = 1, j = 1}^{d, n_{k}}$ ,则称Rpn为带极点的有理函数空间。称广义数列 $\bar{p} = (\overset{n_{1}}{\overset{⌢}{p_{1}, \dots, p_{1}}}, \overset{n_{2}}{\overset{⌢}{p_{2}, \dots, p_{2}}}, \dots, \overset{n_{d}}{\overset{⌢}{p_{d}, \dots, p_{d}}})$ 为极点列(其中p1=∞)。

易见μi(1≤i≤n)线性无关,因此Rpn为n维向量函数空间。

文献[5]中给出了空间Rpn上满足

$λ_{s} (p) = λ_{s} (f) ‚ s = 1, 2, \dots, n (2)$

的有理函数p,其中广义泛函

$\begin{array}{l} λ_{s} = λ_{i j} = \frac{d^{j - 1} (\cdot)}{d x^{j - 1}} |_{x = x_{i}} \\ (i = 1, 2, \dots, k . j = 1, 2, \dots, m_{k}) (3) \end{array}$

这里s= $\sum_{l = 1}^{i - 1}$ ml+j(m0=0)。

该文献还证明了问题的存在唯一性,进一步推导了误差表示。

2 四阶带极点有理样条函数

由于高次函数的不稳定性,实际应用中常常利用低阶函数空间来进行插值。下面在具有两个极点的四维空间U4=span $1, x, \frac{1}{x - a + ε}, \frac{1}{b + ε - x}$ 上研究有理样条问题。

2.1 C1类分片低次有理插值函数

定理1 对[a,b]做分割a=x0<x1<…<xn=b,每个区间 $[x_{i - 1}, x_{i}] (i = 1, 2, \dots, n)$ 上相应的C1类低次Hermite有理插值表达式为

$\begin{array}{l} p (t) = (f_{0} (\frac{t - x_{i - 1}}{h_{i}}), f_{1} (\frac{t - x_{i - 1}}{h_{i}}), \\ h_{i} g_{0} (\frac{t - x_{i - 1}}{h_{i}}), h_{i} g_{1} (\frac{t - x_{i - 1}}{h_{i}})) (\begin{matrix} f (x_{i - 1}) \\ f (x_{i}) \\ f^{'} (x_{i - 1}) \\ f^{'} (x_{i}) \end{matrix}) (4) \end{array}$

基函数分别为

$\begin{array}{l} f_{0} (x) = ε^{2} + ε + 1 + (- 2 ε^{2} - 2 ε - 1) x - \\ \frac{(ε + 1)^{2} ε^{2}}{x + ε} + \frac{(ε + 1)^{2} ε^{2}}{1 + ε - x} ‚ \\ f_{1} (x) = - ε (ε + 1) + (2 ε^{2} + 2 ε + 1) x + \\ \frac{(ε + 1)^{2} ε^{2}}{x + ε} - \frac{(ε + 1)^{2} ε^{2}}{1 + ε - x} ‚ \\ g_{0} (x) = ε (ε + 1) - ε (ε + 1) x - \\ \frac{(ε + 1)^{3} ε^{2}}{(2 ε + 1) (x + ε)} + \frac{(ε + 1)^{2} ε^{3}}{(2 ε + 1) (1 + ε - x)} ‚ \\ g_{1} (x) = - ε (ε + 1) x - \frac{(ε + 1)^{2} ε^{3}}{(2 ε + 1) (x + ε)} + \\ \frac{(ε + 1)^{3} ε^{2}}{(2 ε + 1) (1 + ε - x)} 。 \end{array}$

容易验证这组规范基满足

(1) f0(x)+f1(x)=1;

(2) fi(j)=gi′(j)=δij, (i,j=0,1);

fi′(j)=gi(j)=0, (i,j=0,1)。

定义2 函数组(f0(x),f1(x),g0(x),g1(x))称为空间U4上的规范低次有理Hermite基。

对分片低次Hermite有理插值p(t)有下面的误差定理。

定理2 设M(t)=(t-xi-1+εhi)((1+ε)hi-t+xi-1),w(x)=(x-xi-1)2(x-xi)2,则区间[xi-1,xi]上满足条件p(xi) = yi,p′(xi) = yi′ (i=0,1,2,…,n)的C1类分片低次Hermite有理插值p(t)的插值余项为

$f (t) - p (t) = \frac{w (t)}{4! Μ (t)} (Μ (ξ) f (ξ))^{(4)} ‚ ξ \in (x_{i - 1}, x_{i}) (5)$

且 |R(t)|=|f(t)-p(t)

$\begin{array}{l} \leq \frac{h_{i}^{2}}{384 ε (1 + ε)} \times \\ \max_{x_{i - 1} \leq ξ \leq x_{i}} | (Μ (ξ) f (ξ))^{(4)} | 。 \end{array}$

2.2 C2类有理插值样条函数

定义3 给定一个分划Δ:a=x0<x1<…<xn=b,函数s(x)满足

ⅰ)在每个区间 $[x_{i - 1}, x_{i}] (i = 1, 2, \dots, n)$ 上,

$\begin{array}{l} s (t) = a_{i 0} + a_{i 1} \frac{t - x_{i - 1}}{h_{i}} + a_{i 2} \frac{1}{(\frac{t - x_{i - 1}}{h_{i}}) + ε} + \\ a_{i 3} \frac{1}{1 + ε - (\frac{t - x_{i - 1}}{h_{i}})} (6) \end{array}$

ⅱ)s(t)∈C2[a,b]。

则称s(t)为区间[a,b]上关于分划Δ的四阶带极点有理样条函数。若还满足在节点xi(i=1,2,…,n)处有s(xi)=yi,则称s(t)为[a,b]上的关于Δ的四阶带极点有理插值样条函数。该样条空间称为四阶带极点有理样条函数空间,记为R3,2 (2)(Δ; U4 )。

下面给出常见边界条件下样条存在唯一性,同时给出样条的计算方法。

定理3 给定常见的几种边界条件:

$(a) {\begin{cases} s^{'} (x_{0}) = y_{0^{'}} \\ s^{'} (x_{n}) = y_{n^{'}} \end{cases} ‚ (b) {\begin{cases} s^{˝} (x_{0}) = y_{0^{″}} \\ s^{˝} (x_{n}) = y_{n^{″}} \end{cases} ‚$

$\begin{array}{l} (c) {\begin{cases} s^{'} (x_{0}) = y_{0^{'}} \\ s^{˝} (x_{n}) = y_{n^{″}} \end{cases} 或 {\begin{cases} s^{˝} (x_{0}) = y_{0^{″}} \\ s^{'} (x_{n}) = y_{n^{'}} \end{cases}, \\ (d) s^{˝} (x_{- 1}) = s^{˝} (x_{n + 1}) = 0 。 \end{array}$

则在上述边界条件下四阶带极点有理插值样条函数存在且唯一。

证明由于定理1中样条函数满足内节点处一阶连续,即s(t)∈C1[0,1],由内节点处二阶连续 s″(xi-0)=s″(xi+0),得

$\begin{array}{l} \frac{h_{i + 1}}{h_{i} + h_{i + 1}} m_{i - 1} + \frac{2 ε^{2} + 2 ε + 1}{ε (ε + 1)} m_{i} + \frac{h_{i}}{h_{i} + h_{i + 1}} m_{i + 1} = \\ \frac{3 ε^{2} + 3 ε + 1}{ε (ε + 1)} [\frac{h_{i + 1}}{h_{i} + h_{i + 1}} \frac{y_{i} - y_{i - 1}}{h_{i}} + \frac{h_{i}}{h_{i} + h_{i + 1}} \frac{y_{i + 1} - y_{i}}{h_{i + 1}}] (7) \end{array}$

设

${\begin{cases} λ_{i} = \frac{h_{i + 1}}{h_{i} + h_{i + 1}}, μ_{i} = 1 - λ_{i} \\ l_{i} = 2 + ε^{- 1} (ε + 1)^{- 1} (i = 1, \dots, n - 1) \\ c_{i} = (3 + ε^{- 1} (ε + 1)^{- 1}) \times \\ [λ_{i} \frac{y_{i} - y_{i - 1}}{h_{i}} + μ_{i} \frac{y_{i + 1} - y_{i}}{h_{i + 1}}] (8) \end{cases}$

则方程(7)可化为

$λ_{i} m_{i - 1} + l_{i} m_{i} + μ_{i} m_{i + 1} = c_{i} (i = 1, 2, \dots, n - 1) (9)$

n-1个方程n+1个未知数,边界条件给出2个附加方程。由边界(a)(b)(c)(d)得到的附加方程

$\begin{array}{l} l_{0} m_{0} + μ_{0} m_{1} = c_{0}, \\ λ_{n} m_{n - 1} + l_{n} m_{n} = c_{n} 。 \end{array}$

因此二阶有理样条插值函数的求解归结为下面方程组的求解。

$\begin{array}{l} (\begin{matrix} l_{0} & μ_{0} \\ λ_{1} & 2 + ε^{- 1} (1 + ε)^{- 1} & μ_{1} \\ ⋮ & ⋮ & ⋱ \\ λ_{n - 1} & 2 + ε^{- 1} (1 + ε)^{- 1} & μ_{n - 1} \\ λ_{n} & l_{n} \end{matrix}) \times \\ [m_{0} m_{1} \dots m_{n - 1} m_{n}]^{Τ} = [c_{0} c_{1} \dots c_{n - 1} c_{n}]^{Τ} (10) \end{array}$

(10)式中λi,li,μi,ci(i=1,2,…,n-1)l0,ln见式(8), μ0,λn,c0,cn由各边界条件可得。方程组(10)的系数矩阵是严格对角占优的。由基尔希高林(Gershgorin)定理知上述有理样条函数若规定了型值yi(i=0,1,…,n),则它的解存在且唯一,同时解法是稳定的。

可以证明上述定理所得R3,2 (2)(Δ ;U4 )中有理样条函数的误差估计可达到 $ο (h^{\frac{3}{2}})$ 。

摘要：主要研究带极点有理样条函数空间——R3,2(2)(Δ;U4),不仅证明了样条函数的存在唯一性,而且还给出了其计算方法。该方法利用Gershgorin定理,由追赶法求解,解法稳定。

关键词：带极点有理空间,C1类分片有理插值,带极点有理样条,边界条件

参考文献

[1]Langer R E.On numerical approximation.Madison:The University of Wisconsin Press,1959:25—43

[2]Angel R S,Guillermo L L.Approximation of transfer of infinite di-mensional dynamical systems by Tational interpolants with prescribed poles.Journal of Mathematical Analysis and Application,2000;244(1):147—168

[3]王仁宏,朱功勤.有理函数逼近及其应用.北京:科学出版社,2004

[4]Muehlbach G.Interpolation by Cauchy-Vandermonde systems and ap-plications.Journal of Computational and Applied Mathematics,2000;122(2):203—222

[5]谭高山.Cauchy-Vandermonde空间上插值问题的研究.南京晓庄学院学报,2007;3:14—17

二次有理B样条G2连续插值篇2

关键词：二次有理B样条,型值点,控制顶点,相对曲率,权因子,插值,过渡曲线

0引言

曲线造型是计算机辅助几何设计和图形学的基础, 其典型代表是参数曲线, 如Bézier曲线和B样条曲线[1,2]等, 尤其是NURBS已作为统一的形状数学描述方法纳入工业产品几何定义的STEP国际标准。在CAD/CAM中用有理多项式函数表示曲线、曲面已经越来越广泛, 但有理插值由于涉及权因子和控制顶点的反求, 至今仍是有理函数曲线曲面表示领域中没有完善解决的问题。

由于几何连续性与参数的选取以及具体的参数化无关, 它着眼于形状内在几何特征的描述, 只要求较弱的限制条件, 为形状定义与形状控制提供了额外的自由度, 也为人机交互进行形状设计提供了更有效的工具, 因此几何连续性已逐步成为光顺性研究的重要依据, 成为曲线曲面造型的一个发展趋势。文献[3]给出了基于曲率参数的NURBS曲线插值;文献[4] 给出了三次B样条曲线间G2连续条件;文献[5] 给出了三次有理Bézier样条曲线G2光滑拼接条件;文献[6] 构造了G2连续的三次有理Bézier插值曲线;文献[7] 给出了有理二次Bézier曲线的G2连续条件, 构造了G2连续的二次有理Bézier插值曲线。上述文献讨论的是B样条、NURBS曲线以及Bézier样条曲线的连续条件。二次有理Bézier曲线是二次有理B样条曲线的特例, 二次有理B样条曲线G2连续拼接的条件是否与二次有理Bézier类似?对于曲线的调整能否通过改变某段曲线的首末端曲率来实现, 而不是改变控制顶点或权因子?

本文讨论了二次有理B样条G2连续条件并给出了一种二次有理B样条G2连续插值曲线的构造方法:通过给出某段曲线的首端的曲率κPk (0) 和该段曲线的首端的切矢量的方向角θ以及该段曲线的权因子ωk, 利用几何连续条件来得到其余各段曲线的控制顶点及权因子, 从而生成整条插值曲线。并且给出了过渡曲线的构造方法, 可以得到闭的G2连续插值曲线。该方法构造的插值曲线可以通过改变κPk (0) 、θ和ωk的值对曲线进行修改, 不必做反求运算, 方便实际应用。

1二次有理B样条曲线及其导数

二次有理B样条曲线为[8,9]:

undefined (1)

其中P0、P1、P2为控制顶点;ω0、ω1、ω2为控制顶点的权因子。通常情况下, ωi>0。设ω0=ω2=1, ω1=ω, Bi, 2 (u) 是二次均匀B样条基函数:

undefined

由式 (1) 计算P (u) 的曲率[10]和对u的各阶导数在端点的值分别为:undefined (2)

undefined (3)

undefined (4)

undefined (5)

undefined (6)

2二次有理B样条G2光滑拼接条件

任意给定空间中的6个点P0、P1、P2、Q0、Q1、Q2, 分别以P0、P1、P2和Q0、Q1、Q2为控制顶点, ω1、ω2为权因子, 可构造两条二次有理B样条曲线:

C1:undefined

C2:undefined

现考虑它们的拼接条件。

定理1 两条曲线C1 (u) 、C2 (u) 的连续拼接的充要条件是:undefined。

证明两条有理B样条曲线连续拼接的充要条件是两曲线具有公共的连接点[11], 即P1 (1) =P2 (0) , 由二次有理B样条基函数的性质, 该式即:

undefined (7)

定理得证。

定理2 两条曲线C1 (u) 、C2 (u) 满足如下条件:存在α>0使得undefined, undefined, 则两条曲线C1 (u) 、C2 (u) 在公共连接点实现G1光滑拼接。

证明两条有理二次B样条G1光滑拼接的充要条件是:

(1) 连续拼接充要条件;

(2) 两条曲线在公共连接点处具有公共的切矢方向[11]。

条件 (1) 是前面讨论过的undefined, 而条件 (2) 则为:

P′1 (1) =αP′2 (0) α>0 (8)

根据式 (3) 、式 (5) 可以计算出:

undefined

因而式 (8) 变为:

undefined (9)

式 (7) 与式 (9) 构成了两条二次有理B样条曲线G1光滑拼接的充要条件, 故定理2得证。

定理3 若两曲线C1 (u) 、C2 (u) 满足如下条件:

undefined

且存在α满足:

undefined (10)

则两条曲线C1 (u) 、C2 (u) 在公共连接点处实现G2光滑拼接。

证明两曲线C1 (u) 、C2 (u) 的G2光滑拼接的充要条件是:

(1) 在拼接点处满足G1连续条件;

(2) 拼接点处曲率相等。

由条件 (1) 知:

P1 (1) =P2 (0) P′1 (1) =αP′2 (0) α>0

由条件 (2) 可得 :

undefined (11)

将式 (3) -式 (6) 代入式 (11) , 可得:

undefined

结合G1连续的条件, 式 (7) 、式 (9) 、式 (10) 即为两曲线C1 (u) 、C2 (u) 的G2光滑拼接的充要条件, 从而定理3得证。

3G2连续的二次有理B样条插值曲线

问题:给定平面上一组型值点b1, b2, …, bn, 现要在每两个型值点bi, bi+1;bi+1, bi+2之间各构造一段二次有理B样条曲线Pi, Pi+1, ωi, ωi+1分别为其权因子, 且Pi分别以bi, bi+1为首末端点, Pi+1分别以bi+1, bi+2为首末端点, Pi, Pi+1满足G2连续, 且要求第k (1≤k≤n-1) 段曲线满足曲率约束。

令曲线段Pi (u) =Pi (x (u) , y (u) ) (i=1, …, n-1) 的控制顶点为P3i-2, P3i-1, P3i, 由二次有理B样条曲线的端点性质可知:

undefined (12)

undefined (13)

bk的坐标记为 (xbk, ybk) , 且满足xbk≤xbk+1;Pi的坐标记为 (xPi, yPi) 。设第k段曲线的相对曲率为κPk (u) 可表示为:

undefined

则第k段曲线的首端的相对曲率κPk (0) 可表示为:

undefined (14)

将式 (3) 、式 (4) 代入并结合式 (12) 、式 (13) 整理可得:

undefined (15)

如果考虑的是第k段曲线的曲率约束, 为简化计算过程, 可取bk为原点。由于曲率与曲线的平移旋转以及比例缩放无关。因此简化后的结果不会影响结果的一般性。记第k段曲线的首端的切矢的方向角为θ, 并规定逆时针方向为正。则式 (15) 可化简为:

undefined (16)

又:

undefined (17)

undefined (18)

则由式 (16) -式 (18) , 得 P3k-1的坐标可以由κPk (0) 和θ以及ωk表示:

xP3k-1=1+ωkωkybk+1cosθ-xbk+1sinθ8κPk (0) cosθ

yP3k-1=1+ωkωkybk+1cosθ-xbk+1sinθ8κPk (0) sinθ (19)

由式 (12) 、式 (13) 可得:

P3k-2= (1+ωk) bk-ωkP3k-1 (20)

P3k= (1+ωk) bk+1-ωkP3k-1 (21)

此时第k段插值曲线的控制顶点以及权因子已全部求得, 即可确定第k段曲线。

要使得第k+1段曲线与第k段曲线在公共点达到G2连续, 只要满足:

(1) Pk (1) =Pk+1 (0)

(2) Pk′ (1) =αk+1Pk+1′ (0)

(3) κPk (1) =κPk+1 (0)

将式 (3) -式 (6) 、式 (12) 、式 (13) 带入上述条件中, 整理可得:

undefined (22)

其中:

a1= (xP3k-1-xP3k-2) (yPk (1) -yPk (0) )

b1= (yP3k-1-yP3k-2) (xPk (1) -xPk (0) )

c1= (xP3k-1-xP3k-2) (yPk+1 (1) -yPk+1 (0) )

d1= (yP3k-1-yP3k-2) (xPk+1 (1) -xPk+1 (0) )

由条件 (2) 可得:

undefined (23)

由二次有理B样条端点性质可得:

P3k+1= (1+ωk+1) bk+1-ωk+1P3k+2 (24)

P3k= (1+ωk+1) bk+2-ωk+1P3k+2 (25)

要使得第k-1段曲线与第k段曲线在公共点达到G2连续, 只要满足:

(1) Pk-1 (1) =Pk (0)

(2) Pk-1′ (1) =αkPk′ (0)

(3) κPk-1 (1) =κPk (0)

同理可得:

undefined (26)

其中:

a2= (xP3k-1-xP3k-2) (yPk-1 (1) -yPk-1 (0) )

b2= (yP3k-1-yP3k-2) (xPk-1 (1) -xPk-1 (0) )

c2= (xP3k-1-xP3k-2) (yPk (1) -yPk (0) )

d2= (yP3k-1-yP3k-2) (xPk (1) -xPk (0) )

undefined (27)

P3k-5= (1+ωk-1) bk-1-ωk-1P3k-4 (28)

P3k-3= (1+ωk-1) bk-ωk-1P3k-4 (29)

此时, 各段插值曲线的控制顶点以及权因子已全部求得, 则由式 (1) 可依次确定各段插值于bi、bi+1 (i=1, 2, …, n-1) 的曲线, 由式 (1) 可得到整条G2连续的二次有理B样条插值曲线。

综上所述, 绘制过给定型值点的G2连续的二次有理B样条的步骤为:

Step1 输入型值点的坐标b1, ..., bn;

Step2 输入第k段曲线的首端的端点的曲率κPk (0) 和该段曲线首末端点连线与切矢所形成的角度θ以及权因子ωk;

Step3 计算由式 (19) , 式 (20) , 式 (21) 决定的第k段插值曲线的控制顶点的坐标P3k-2, P3k-1, P3k;

Step4 计算由式 (23) 至式 (29) 决定的其他控制顶点坐标P3j-5, P3j-4, P3j-3, P3j+1, P3j+2, P3j+3以及权因子ωj-1, ωj+1 (j=1, …, k-1, k+1, n-1) ;

Step5 按式 (1) 依次绘制第k (k=1, ..., n-1) 段曲线;

Step6 若所得的曲线不合要求, 则返回Step2, 修改κPk (0) 和θ以及ωk, 重复以上步骤, 直到满意为止。

4插值曲线的调节

由式 (19) 可知, 若κPk (0) 或θ或者ωk改变时, 第k段曲线的控制顶点会随之改变, 那么由G2连续条件可得所有的控制顶点均会发生变化, 即可通过改变κPk (0) 或θ或者ωk的取值对曲线进行调整。

下面给出实例来说明相对曲率, 夹角以及权因子对插值曲线的调节。

若给平面上的一列型值点, 分别为b1=[-10, 10], b2=[0, 0], b3=[5,7], b4=[13, 2], b5=[19, 6], b6=[24, 8], b7=[30, 12], b8=[35, 16], 如图1所示, 实线表示第二段曲线的首端曲率κP2 (0) =0.8, 该段曲线的首端的切矢量的方向角θ=-0.8, 权因子ω=[1,2,3,5]的插值曲线, 虚线表示第二段曲线的首端曲率κP2 (0) =0.9, 该段曲线的首端的切矢量的方向角θ=-0.8, 权因子ω=[1,2,3,5]的插值曲线。如图2所示, 实线表示第二段曲线的首端曲率κP2 (0) =0.8, 该段曲线的首端的切矢量的方向角θ=-0.8, 权因子ω=[1,2,3,5]的插值曲线, 虚线表示第二段曲线的首端曲率κP2 (0) =0.8, 该段曲线的首端的切矢量的方向角θ=-2, 权因子ω=[1,2,3,5]的插值曲线。

上述都是以改变某一段插值曲线的首端的曲率或者切矢角来改变整条插值曲线的形状。若改变某一段的权因子ωi, 那么曲线段Pi以及后面的插值曲线段将会发生变化, 如图3所示。实线表示的插值曲线的第二段曲线的首端曲率κP2 (0) =0.8, 该段曲线的首端的切矢量的方向角θ=-0.8, 权因子ω=[1,2,3,5], 虚线表示的是插值曲线的第二段曲线的首端曲率κP2 (0) =0.8, 该段曲线的首端的切矢量的方向角θ=-0.8, 权因子ω=[1,2,3,5]。

5闭曲线的构造

给定两条有理二次B样条曲线C1、C3, 其控制顶点分别为P0、P1、P2、V0、V1、V2;权因子分别为ω1、ω3。现欲在这两条曲线段之间构造一条曲线段C2, 与C1、C3均为G2连续拼接。设C2的控制顶点为Q0、Q1、Q2;权因子为ω2。

由二次有理B样条的表达式可得:

C1:undefined

C2:undefined

C3:undefined

由于C1, C2;C2, C3均满足G2连续条件, 我们可以得到:

(1) P1 (1) =P2 (0) ; (2) P1′ (1) =αP2′ (0) ;

(3) κP1 (1) =κP2 (0) ; (4) P2 (1) =P3 (0) ;

(5) P2′ (1) =βP3′ (0) ; (6) κP2 (1) =κP3 (0) 。

b1、b2、b3分别表示C1、C2、C3的首端端点;b4表示C3的末端端点。

将式 (3) -式 (6) 代入到上述条件中, 可得到:

undefined

由条件 (2) 、 (5) 可得, Q1在P1P2与V0V1的交点上, 可以通过求解P1P2与V0V1的直线方程组得到Q1的坐标。再根据上述条件, 可以求得:

undefined

Q0= (1+ω2) b2-ω2Q1Q2= (1+ω2) b3-ω2Q1

这时过渡曲线C2的权因子与控制顶点已全部得到, 也就可以得到整条过渡曲线。

下面给出闭曲线的例子。给定一列点, b1=[0, 0], b2=[2,6], b3=[8,9], b4=[5, -5], 第一段插值曲线的首端端点的曲率为κP1 (0) =-1, 该点的切矢方向角为θ=2.2, 权因子为ω=[3,5]。首先可以构造一条过型值点b1、b2、b3、b4的插值曲线, 然后构造一条过b1、b4的过渡曲线, 并且使得过渡曲线与插值曲线达到G2连续, 如图4所示。

6结束语

本文给出了二次有理B样条G2连续拼接的充要条件, 给出了G2连续的二次有理B样条插值曲线的构造方法, 并给出了过渡曲线的构造方法, 从而可以得到G2连续闭插值曲线。用这种方法得到的插值曲线具有如下优点:

(1) 如果给定某一段曲线的首端的端点的相对曲率和该段曲线的首端切矢量的方向角, 由G2连续条件可得到插值曲线的其余控制顶点, 它不需要做近似计算或者求解大型的方程组, 从而减少了计算量;

(2) 可以通过调节某一段曲线首端端点曲率以及该点的切矢方向角对曲线的整体形状进行修改, 如果希望局部修改插值曲线, 可以通过修改曲线的权因子, 那么曲线段Pi以及后面的曲线的形状将会发生改变, 有效而且直观;

(3) 由过渡曲线的构造方法可以得到闭的插值曲线, 而不需要大量的计算。

参考文献

[1]Tiller W.Rational B-splines for curves and surface representation[J].IEEE Computer Graph Application, 1983, 3 (1) :61-69.

[2]Imre J.Cubic parametric curves of given tangent and curvature[J].CAD, 1998, 30 (1) :1-16.

[3]程永福, 朱功勤.三次有理Bézier样条曲线G2连续拼接条件[J].合肥工业大学学报, 2007, 30 (3) :397-399.

[4]杨莉, 晁翠华, 贾晓.G2连续的三次有理Bézier样条插值曲线[J].机械科学与技术, 2000, 19 (3) :512-513.

[5]陈宝平, 尹志凌.基于有理二次Bézier曲线的G2连续的插值曲线[J].内蒙古大学学报, 2004, 35 (4) :464-466.

[6]苏晓红, 李东, 王宇颖.基于曲率参数的NURBS曲线插值[J].哈尔滨工业大学学报, 2001, 33 (1) :108-111.

[7]冯仁忠, 王仁宏.三次B样条曲线间G2连续条件[J].大连理工大学学报, 2003, 43 (4) :407-411.

[8]朱心雄.自由曲线曲面造型技术[M].北京:科学出版社, 2002.

[9]施法中.计算机辅助几何设计与非均匀有理B样条 (CAGD&NURBS) [M].北京:高等教育出版社, 2001.

[10]彭家贵, 陈卿.微分几何[M].北京:高等教育出版社.2002.

函数的分段有理二次B样条插值篇3

科学和工程计算中,函数的近似表示一直是一个重要课题。近似方法一般可归结为插值、逼近和拟合三种基本类型,经历长期发展,函数逼近方法[1,2,3]十分丰富。诸如Taylor公式、Lagrange-Newton插值、多项式样条方法、最佳平方逼近、最小二乘拟合、连分式插值、有理样条表示等等,它们在函数近似表示方面均发挥着不同程度的作用。

有理B样条曲线[4,5,6,7]是计算机辅助几何设计中常用的一元函数的参数形式,它具有良好的几何性质以及比较成熟的算法。文献[8]在讨论代数曲线的分段逼近时,采用“单权因子”的标准形式的Bézier,提出了代数曲线的“肩点”到有理Bézier曲线的“肩点”的最小距离逼近算法,以解决代数曲线的参数化问题。考虑到传统的插值方法较少考虑函数的几何性质,同时也很少见到以参数方程给出的插值函数,这里将函数的有理B样条形式引入函数近似表示领域,给出一种基于有理二次B样条曲线的插值算法。

根据函数曲线的凹凸性和单调性,对函数进行分段,定义分段函数的“凸包”,并由“凸包”来确定控制多边形;重点介绍基于“肩点”等“关键点”信息的“双权因子”插值算法;相对于“单权因子”逼近算法[8],误差有明显改善;最后以数值实验说明了算法的有效性。

1 函数的分段

1.1 函数的几何特征分段

已知函数:

$\tilde{F}$ :y=f(x) x∈[a0,an] (1)

设ai∈[a0,an],i=1,2,…n-1为函数拐点的横坐标,将定义域[a0,an]分割为:

[a0,an]=[a0,a1]∪[a1,a2]∪…

∪[an-2,an-1]∪[an-1,an] (2)

在区间:

[ai,ai+1] i=0,1,…,n-1 (3)

上,第i段函数表示为:

$\tilde{F}_{i}$ :y=f(x) x∈[ai,ai+1] (4)

设分段函数式(4)在区间式(3)上的极值点为ami,将区间式(3)进一步分割(若已经单调,则不再分割)为如下形式:

[ai,ami]∪[ami,ai+1],i=0,1,…n-1 (5)

则分段函数式(4)进一步分割为如下两段:

y=f(x) x∈[ai,ami] x∈[ami,ai+1] (6)

这样,所给函数在每个分段区间上具有固定的凹凸性和单调性。

1.2 分段函数的三角形凸包

假定分割后的一段曲线为:

F:y=f(x) x∈[a,b] (7)

设它的两个端点为V0(xa,ya)、V2(xb,yb),求出两端点处切线,记为:

T0:a0x+b0y+c0=0 T2:a2x+b2y+c2=0 (8)

容易证明,切线T0、T2一定相交于一点,记交点为V1(x1,y1)。这样,由三点V0、V1、V2形成了分段函数曲线F的三角形凸包(如图1所示)。

例如,给定函数(如图2所示)。

C:y=x3+x2ex-sin(2πx)+1 x∈[-0.5,0.5] (9)

它有一个拐点(-0.0077,1.0484)和两个极值点(-0.2538,2.0333)、(0.2296,0.0866),将函数分割为四段,各段及其三角形凸包如图3所示。

2 函数分段二次插值与逼近

2.1 “单权因子”标准有理二次Bézier曲线逼近

设控制顶点为Qi(xi,yi),i=0,1,2,权因子为wi(>0),i=0,1,2,则有理二次Bézier曲线P(u),u∈[0,1]可以表示为:

${\begin{cases} x (u) = \frac{w_{0} x_{0} (1 - u)^{2} + 2 w_{1} x_{1} u (1 - u) + w_{2} x_{2} u^{2}}{w_{0} (1 - u)^{2} + 2 w_{1} u (1 - u) + w_{2} u^{2}} \\ y (u) = \frac{w_{0} y_{0} (1 - u)^{2} + 2 w_{1} y_{1} u (1 - u) + w_{2} y_{2} u^{2}}{w_{0} (1 - u)^{2} + 2 w_{1} u (1 - u) + w_{2} u^{2}} \end{cases}$

按照Farin和Worsey提出的Bézier曲线的参数化标准形式[9],取w0=w2=1,为确定式(2)中余下的权因子w1,文献[8]中采用“肩点”到“肩点”距离最小逼近方法。具体方案是:已知凸包 $V_{i} (\tilde{x}_{i}, \tilde{y}_{i}), i = 0, 1, 2$ ;设F的“肩点”为SF(xF,yF),P(u)的“肩点”为:

$\begin{array}{l} S_{Ρ} = Ρ (1 / 2) = S (w_{1}) = (x_{1 / 2}, y_{1 / 2}) \\ = (\frac{x_{0} + 2 w_{1} x_{1} + x_{2}}{2 (1 + w_{1})}, \frac{y_{0} + 2 w_{1} y_{1} + y_{2}}{2 (1 + w_{1})}) \end{array}$

记d2(SF,SP)=(xF-x1/2)2+(yF-y1/2)2,由 $\frac{\partial d^{2}}{\partial w_{1}} = 0$ 来确定权因子w1,结果为:

$w_{1} = \frac{1}{2} \cdot \frac{(x_{0} + x_{2} - 2 x_{F}) + α (y_{0} + y_{2} - 2 y_{F})}{(x_{F} - x_{1}) + α (y_{F} - y_{1})}$

其中, $x_{i} = \tilde{x}_{i}, y_{i} = \tilde{y}_{i}, i = 0, 1, 2$ ; $α = \frac{y_{0} + y_{2} - 2 y_{1}}{x_{0} + x_{2} - 2 x_{1}}$ 。

2.2 “双权因子”有理二次均匀B样条插值

2.2.1 有理二次均匀B样条曲线

设控制顶点为Qi(xi,yi),权因子为wi(>0),i=0,1,2;则有理二次均匀B样条曲线[4,5,6,7]P(u),u∈[0,1]可以表示为:

考虑到计算复杂性等因素,这里选取:

w0=1 (11)

并称由此对应的B样条曲线为“准标准形式”。

2.2.2 控制顶点的确定

已知F的三角形凸包V0V1V2,设 $V_{i} = (\tilde{x}_{i}, \tilde{y}_{i}), i = 0, 1, 2$ ;取控制顶点Q1=V1,则:

$x_{1} = \tilde{x}_{1} y_{1} = \tilde{y}_{1}$ (12)

由P(0)=V0,P(1)=V2得 $\frac{Q_{0} + w_{1} Q_{1}}{1 + w_{1}} = V_{0}, \frac{w_{1} Q_{1} + w_{2} Q_{2}}{w_{1} + w_{2}} = V_{2}$ ,即有:

x0=(1+w1)x～0-w1x～1

y0=(1+w1)y～0-w1y～1 (13)

x2=(w1+w2)x～2-w1x～1w2

y2=(w1+w2)y～2-w1y～1w2 (14)

这样,控制顶点将借助于下文的权因子w1、w2的确定而确定(如图4所示),曲线P(u)中的w1、w2的自由变化对应曲线段A的插值曲线族(如图5所示)。

2.3 有理二次B样条插值

选取F上的“关键点”,如“肩点”[9]SF(xF,yF)(如图6所示),它指的是F上到线段V0V1具有最大距离的点,“中点” $Μ (\frac{(a + b)}{2}, f (\frac{(a + b)}{2}))$ ,“n等分点”等,同时选取B样条曲线P(u)上的“关键点”,如“肩点”SP(xP,yP)、参数域中点P(1/2)(记为S1/2)、参数域的“n等分点”等。这时,插值条件可以在“关键点”间的对等关系中选取一种,如SP=SF;S1/2=M等等,这里我们选取。

S1/2=SF (15)

考虑到式(10)、(11)两式有:

$\frac{x_{0} + 6 x_{1} w_{1} + x_{2} w_{2}}{1 + 6 w_{1} + w_{2}} = x_{F}$ ; $\frac{y_{0} + 6 y_{1} w_{1} + y_{2} w_{2}}{1 + 6 w_{1} + w_{2}} = y_{F}$ (16)

将式(13)、式(14)代入式(16)可以得到关于w1、w2的方程组,求出w1、w2,并由式(12)、式(13) 、式(14)确定控制顶点为Qi(xi,yi),i=0,1,2,最后由式(10)确定P(u)。

2.4 插值函数的性质

1) 保持原有函数的凹凸性不变;

2) 保持原有函数的单调性不变;

3) 在函数的分段点处保持G1连续性[6];

4) 插值曲线通过原始曲线的某个“关键点”。

3 分段插值算法

3.1 插值误差

在几何学中,曲线段C和它的近似表示Ca之间的误差经常使用Hausdorff距离[9]来描述:

$e (C, C_{a}) = d i s (C, C_{a}) = \max_{Ρ \in C} \min_{Q \in C_{a}} d (Ρ, Q)$

这种距离便于理论分析,但不便于计算。

设F为分段待插函数,P(u)为分段插值函数,下面给出一种易于计算的误差概念。

定义1 在参数域[0,1]上的任意取一点 $u_{i} = \frac{i}{n}, 0 \leq i \leq n, n \in Ν$ ,记P(u)上的对应点为(x(ui),y(ui)),F上的对应点为(x(ui),yi),定义插值误差为:

$e (F, Ρ (u)) = \max_{i} | y_{i} - y (u_{i}) |$ (17)

3.2 算法描述

输入原始函数 $\tilde{F}$ 以及误差限δ,输出误差小于δ的分段有理二次B样条曲线P(u)。

1) 由式(1)至式(6)对原始函数进行分段;

2) 计算分段函数曲线F两端切线T0,T2及其交点V1,确定三角形凸包V0V1V2;

3) 计算F的“肩点”SF以及P(u)上的点S1/2;

4) 根据式(12)、式(13)、式(14)、式(16)计算权因子w1,w2和控制顶点为Qi(xi,yi),i=0,1,2;

5) 依据式(10)计算插值曲线P(u);

6) 若插值误差e(F,P(u))<δ,则计算过程中止;否则,以“肩点” SF为界,将F一分为二,递归调用上述插值算法,直至满足给定的误差要求。

3.3 算法的收敛性

定理1 在算法的实施过程中,F和P(u)之间的Hausdorff距离收敛于0。即如果经过k次分割后的Hausdorff距离为Ek(F,P(u)),则:

$\lim_{k \to + \infty} E_{k} (F, Ρ (u)) = 0$ (18)

证明如图7所示,由函数的分段方法及其凸包形成过程可知,ΔV0V1V2是确定的。设∠V1V0V2=α,∠V1V2V0=β,0<α,β≤π/2,记边V0V2的高为h,则h=|V0V1|sinα,h=|V2V1|sinβ,记“肩点”SF到边V0V2的距离为d。设经过SF的切线分别和V0V1、V2V1相交于P1、P2,自P1、P2分别向V0SF、V2SF引垂线,记P1H1=d11,P2H2=d12。

又设∠P1V0SF=α1,∠P2V2SF=β1,0<α1,β1<π/2。显然E0(F,P(u))<d<h,经过在“肩点”处一次分割后:

d11=|P1V0|sinα1=a1|V1V0|sinα1 0<a1<1

d12=|P2V2|sinβ1=b1|V1V2|sinβ1 0<b1<1

经过k次分割后:

$d_{k 1} = a_{1} a_{2} \dots a_{k} | V_{1} V_{0} | \sin α_{k} 0 < a_{i} < 1 i = 1, 2 ‚ \dots ‚ k$

$d_{k 2} = b_{1} b_{2} \dots b_{k} | V_{1} V_{2} | \sin β_{k} 0 < b_{j} < 1 j = 1, 2 ‚ \dots ‚ k$

设 $a = \max_{i} {a_{i}}$ , $b = \max_{j} {b_{j}}$ ,显然0<a,b<1,于是:

$\lim_{k \to + \infty} d_{k 1} < \lim_{k \to + \infty} a^{k} | V_{1} V_{0} | \sin α_{k} = 0$

$\lim_{k \to + \infty} d_{k 2} < \lim_{k \to + \infty} b^{k} | V_{1} V_{2} | \sin β_{k} = 0$

从而, $\lim_{k \to + \infty} E_{k} (F, Ρ (u)) = 0$ 。

4 数值实验

图8为函数C分四段(从左至右)“双权因子”的一步插值效果,图9为其误差函数。一步插值误差 $e_{1} (F, Ρ (u)) = \max_{i} | y_{i} - y (u_{i}) | < 6 \times 10^{- 3}$ 。

图10为函数C分四段(从左至右)“单权因子”的一步逼近效果,图11为其误差函数。一步插值误差 $e_{1} (F, Ρ (u)) = \max_{i} | y_{i} - y (u_{i}) | < 15 \times 10^{- 3}$ 。

显然,插值算法的误差不到逼近算法的一半,同时误差分布比较均匀,近似效果明显优于逼近算法。

5 结语

基于函数的单调性和凹凸性分段的函数“双权因子”有理二次B样条插值算法,充实了插值方法的理论内涵,具有一定的实用价值。算法可以以任意精度近似表示已知函数。误差估计方法易于理解,计算简单。“双权因子”插值误差远小于“单权因子”逼近误差,而且误差分布比较均匀。但误差缩小的理论原因有待进一步探讨。

参考文献

[1] Wang Nengchao.Concise tutorial of computational approach[M].High Educational Press,2004.

[2] Wang Renhong.Numerical approximation[M].High Educational Press,1999.

[3] Yang Fengxiang,Zhai Ruicai, Sun Jing. Numerical Analysis [M]. Tinajin University Press, 1996.

[4] Zhong Xinxiong. Modeling Technology of Free Formed Curves and Surfaces[M]. Beijing, Science Press,2000.

[5] Shi Fazhong.Computer Aided Geometric Design an Non Uniform Rational Basic Spline[M]. Beijing,Higher Educational Press,2001.

[6] Les Piegl,Wayne Tiller. The NURBS Book [M]. New York: Springer,1997.

[7] Gerald Farin. Curves and Surfaces for Computer Aided Geometric Design [M].New York:ACADEMIC PRESS, 2000.

[8] Gao Xiao-Shan, Li Ming. Rational quadratic approximation to real algebraic curves[J].Computer Aided Geometric Design,2004,21:805-828.

有理样条篇4

数字图像融合是把多个不同模式的图像传感器获得的同一场景的多幅图像或同一传感器在不同时刻获得的同一场景的多幅图像合成为一幅图像的过程。目前图像融合技术现已广泛地应用于军事、遥感、计算机视觉、医学图像处理等领域中[1]。

图像融合的方法主要分为三种:简单的图像融合方法、基于塔式分解的图像融合方法和基于小波变换的图像融合方法。Nichols在文献[2]中提出一种将三幅不同的红外图像通过线性变换映射到RGB色度空间中。这种图像融合具有实现简单、融合速度快的特点。基于塔式分解的融合方法是首先对参加融合的各源图像分别进行塔式分解,对塔式分解后的图像按照一定的规则进行融合,再进行逆塔式变换得到融合图像。常用的图像塔式分解有Burt提出的拉普拉斯塔式分解[3]、Toet提出的低通比率塔式分解[4]、Akerman提出的Gaussian塔式分解[5],这些都是图像的多尺度、多分辨率分解。基于上述塔式分解的图像融合方法都能取得良好的融合效果。但是,这些塔式分解都是图像的冗余分解,分解后各层数据有冗余和相关性。因此,图像的塔式分解后的数据总量均比原来被分解的图像源有相应的增加,而且图像的上述塔式分解都是无方向性的。

文献[6,7,8]提出了基于小波变换的图像融合方法,即先对图像进行小波分解,形成各自的多尺度描述;然后在各个小波分解图像的相应子图像上按照一定的融合规则进行融合,形成融合后的图像多尺度描述;最后进行小波逆变换重构融合后的图像。因为小波是非冗余的,小波分解是有方向的,所以基于小波变换的图像融合方法可以获得效果更好的融合图像。王宏提出了一种基于图像块分割的多聚集图像融合方法[6],对两幅输入图像,若对比度大时就选择对比度大的块,若对比度相近时,就对该块对应像素的灰度值进行加权组合形成融合图像。这种方法关键在于如何选择权函数。

根据融合在处理阶段的不同,图像融合通常分为像素级图像融合、特征级图像融合和决策级图像融合。像素级图像融合是最低层次的图像融合,该层次的融合准确性最高,能够提供其他层次上的融合处理所不具有的细节信息,是当前的研究热点。

本文提出的图像融合属于像素级图像融合,利用有理三次样条函数构造权函数,对图像进行融合,给出相应的实验结果并同其他的方法进行比较。

1 多聚焦图像融合算法

多聚焦图像分为清晰和模糊两部分, 并认为输入原图像的清晰部分和对应的模糊部分成像的差异不大,大小基本相同。大多数情况下输入图像成像的差异是不能忽略的,因此将多聚焦图像分为清晰区域、模糊区域和两者交界的边界区域。为了便于辨别出图像的上述3个区域,将多聚焦图像分解成若干个块区域,然后对每个块区域进行判断,以确定其所属区域,再进行相应的处理。

1.1 图像块区域的局部对比度的判定与划分

在多聚焦图像融合中,如果一幅输入图像的某个区域是聚焦区域,则这个区域使用聚焦度量函数MMF(Multi-focus Measure Function)得到的值就要比对另一幅输入图像的相同区域使用聚焦度量函数得到的值要大。因此聚焦度量函数的的定义非常关键。传统的聚焦度量函数有方差、空间频率、梯度能量、拉普拉斯能量等[1],考虑到计算的复杂度和有效性,下面定义一种新多聚焦度量函数:

‖·‖为一种范数,通常取为‖·‖∞,‖·‖2,‖·‖1,本文采用2-范数来度量,N是窗口的尺寸。若输入图像是RGB空间的真彩图像,则 $\overset{⇀}{f} (x, y)$ 表示(x,y)处颜色的向量值函数,若是灰度值图像,则 $\overset{⇀}{f} (x, y)$ 退化为标量函数,表示(x,y)处的灰度值。

对图像分块后,利用聚焦度量函数计算分块图像的能量值MMFk,记X图像的第k块能量值MMF $_{k}^{X}$ ,若MMF $_{k}^{A}$ >MMF $_{k}^{B}$ ,则A的第k块为清晰区域,B的相应第k块为模糊区域,反之MMF $_{k}^{A}$ <MMF $_{k}^{B}$ ,则B的第k块为清晰区域,A的相应第k块为模糊区域。

将输入的图像划分为成清晰区域和模糊区域,为保证融合后的效果,将清晰区域与模糊区域相邻的所有块区域划为边界区域, 以此得到图像的3个不同的区域划分。

1.2 图像块区域的清晰度的定义

图像清晰度(Image Clarity)指图像上各种细节及其边界的清晰程度。现定义图像的清晰度如下:

其中i,j为整数,M=2m+1,N=2n+1,2m,2n为邻近(x,y)的块区域长宽。通常为了便于计算取m=n=1或2,‖·‖的计算同式(1)。

图像区域块Ak的清晰度定义为:

其中,k是以(x,y)为中心的窗口半径。

1.3 权函数的构造

Burt在文献[3]中提出了一种加权平均法,是多幅图像融合的简单有效方法,其优点是简单省时,适合实时处理。当输入的两幅图像很相似时,合成的图像就取两幅图像的平均值;当两幅图像差距很大时,就选择最显著的那幅图像。受以上思想的启发,本文将利用有理三次样条函数构造权函数。

权函数ω(s)的选择是本文提出的图像融合算法的关键,应当具有如下性质:

其中m(x,y)=IC(Ak)-IC(Bk),表示以(x,y)为中心,半径为k的A,B图像k区域块的清晰度差。

现构造有理三次样条函数如下:

设a=x0<x1<…<xn=b是区间[a,b]的一个划分,fi是函数f(x)在xi(i=1,2,…,n)处的函数值。记:

$h_{i} = x_{i + 1} - x_{i} Δ_{i} = (f_{i + 1} - f_{i}) / h_{i} θ = \frac{(x - x_{i})}{h_{i}} ‚ i = 1, \dots, n - 1 (7)$

有理三次样条函数构造为:

其中:

Pi(θ)=rifi(1-θ)3+{(2ri+ti)fi+rihidi}θ(1-θ)2+

{(ri+2ti)fi+1-tihidi+1}θ2(1-θ)+tifi+1θ3

Qi(θ)=(1-ri)θ+tiθθ∈[0,1]

ti>ri>0是三次有理样条函数的形状参数。在区间[xi,xi+1],i=1,...,n-1上满足:

R(xi)=fiR(xi+1)=fi+1

R′(xi)=diR′(xi+1)=di+1

di,i=1,...,n一般是未知的,计算方法如下:

$\begin{array}{l} d_{i} = {\begin{cases} 0 f_{i + 1} - f_{i - 1} = 0 \\ Δ_{i} Δ_{i - 1} (x_{i + 1} - x_{i - 1}) / (f_{i + 1} - f_{i - 1}) o t h e r w i s e \end{cases} i = 2, 3, \dots, n - 1 \\ d_{1} = {\begin{cases} 0 f_{3} - f_{1} = 0 \\ (Δ_{1})^{2} (x_{3} - x_{1}) / (f_{3} - f_{1}) o t h e r w i s e \end{cases} \end{array}$

$d_{n} = {\begin{cases} 0 f_{n} - f_{n - 2} = 0 \\ (Δ_{n - 1})^{2} (x_{n} - x_{n - 2}) / (f_{n} - f_{n - 2}) o t h e r w i s e \end{cases}$

当区间[a,b]被等距划分时,di与f′i的误差是O(h2)阶的。若权函数ω(s)采用式(8)过点(-T,0),(0,0.5),(T,1)来构造,显然满足权函数的两条性质,但是为了提高权函数的精度,首先过(-T,0),(0,0.5),(T,1)三点构造一个分段多项式:

用式(9)预估出权函数ω(s)在s=±T/3,±2T/3处的值0.7777,0.2222,0.9444,0.0556,然后通过七点:

(-T,0),(-2T/3,0.0556),(-T/3,0.2222),(0,0.5),(T/3,0.7778),(2T/3,0.9444)(T,1)

构造一个分段有理三次样条函数 $R (θ) = R (\frac{s - s_{i}}{s_{i + 1} - s_{i}}) ‚ i = 0, \dots, 5$ 。图1显示权函数图像示例。

当形状参数ri,ti满足 $\frac{r_{i}}{t_{i}} \geq \frac{(d_{i + 1} - Δ_{i})}{(Δ_{i} - d_{i})}$ ,i=1,2,…,n-1时,所构造的有理三次样条函数式(8)是保单调的,因此在利用式(8)构造权函数选择形状参数ri,ti时注意约束条件。

1.4 多聚焦图像融合处理的算法

当两幅多聚焦输入图像是彩色图像时,(x,y)处的像素值用向量值函数表示,即:

$\overset{⇀}{f} (x, y) = (r (x, y), g (x, y), b (x, y))$

若是灰度图像, $\overset{⇀}{f} (x, y)$ 退化为标量函数,表示(x,y)处的灰度值。彩色图像的融合在RGB空间完成,避免了颜色空间的转换运算。

Step1 将图像分块,利用式2计算块区域的能量,将整个图像划分成清晰区域和模糊区域,清晰区域与模糊区域相邻的所有块区域划为边界区域,以此得到图像的3个不同的区域划分。

Step2 由于输入图像A的清晰区域对应图像B的模糊区域,图像A的模糊区域对应图像B的清晰区域,因此在进行融合处理时, 直接选取清晰块区域作为融合后的相应块区域。

Step3 对于边界区域,计算两幅待融合图像对应块的清晰度,用IC(Ak)和IC(Bk)分别表示输入图像A,B在第k块处的清晰度,并计算清晰度差m(x,y)=IC(Ak)-IC(Bk)。

Step4 若|m(x,y)|>T,则:

其中 $\overset{⇀}{f}_{F} (x, y)$ 为融合后图像在(x,y)处的像素值, $\overset{⇀}{f}_{A} (x, y)$ , $\overset{⇀}{f}_{B} (x, y)$ 分别为输入图像A,B在(x,y)处的像素值。T为清晰度阈值。

若|m(x,y)|≤T

则

$\overset{⇀}{f}_{F} (x, y) = ω \overset{⇀}{f}_{A} (x, y) + (1 - ω) \overset{⇀}{f}_{B} (x, y) (11)$

其中ω,1-ω分别是输入图像A,B在(x,y)处的权值,是关于清晰度差m(x,y)的函数,用所构造的有理样条计算权值ω。

2 数字图像融合实验及实验结果分析

为了验证本文提出的算法有效性,将本文算法与其他算法[6,7,8]的实验比较。

若取定有理样条的形状参数ri,ti,i=1,2,…,n-1,对于不同的阈值T,融合的图像效果也各不相同,如图2所示。

图3(c)-(h)分别为平均值加权、FSD-Pyramid、Laplacian-Pyramiad、主成分分析、Gradient-Pyramid和Ratio-Pyramid 算法与本文所提出的算法融合源图像(a),(b)的结果。从融合的图像可以看出,本文提出算法融合的效果相对较好。图4为参考图像。图5(a)为参考图像的前聚焦彩色图像,(b)为参考图像的后聚焦彩色图像,(d)-(f)分别为离散小波,加权平均和FSD pyramid方法融合的结果,(c)为本文提出的方法融合的结果。从视觉效果看,(c)的纹理细节和图像边缘保持较好。

为客观地评价图像融合的实验效果和质量,引入一些性能评价参数。下文中F表示融合后图像的像素值,M和N分别表示图像的行数和列数。

(1) 熵

熵值的大小表示图像所含的平均信息量的多少,图像的熵定义为:

$Η = - \sum_{i = 1}^{Μ} Ρ_{i} \log_{2} Ρ_{i}$

其中H表示图像的熵;M表示图像总的灰度级数。熵值越大,图像所含信息量越大,信息越丰富。

(2) 空间频率

图像的行频率定义为:

$R F = \sqrt{\frac{1}{Μ \times Ν} \sum_{i = 1}^{Μ} \sum_{j = 2}^{Ν} ∥ \overset{⇀}{F} (i, j) - \overset{⇀}{F} (i, j - 1) ∥_{2}^{2}}$

图像的列频率定义为:

$C F = \sqrt{\frac{1}{Μ \times Ν} \sum_{j = 1}^{Ν} \sum_{i = 2}^{Μ} ∥ \overset{⇀}{F} (i, j) - \overset{⇀}{F} (i - 1, j) ∥_{2}^{2}}$

于是图像的空间频率定义为:

$S F = \sqrt{R F^{2} + C F^{2}}$

空间频率反映了一幅图像空间域的总体活跃程度,空间频率越大,图像的融合效果就越好。

(3) 平均梯度

图像的平均梯度定义为:

$\begin{array}{l} A G = \sum_{i = 1}^{Μ - 1} \sum_{j = 1}^{Ν - 1} [∥ \overset{⇀}{F} (i, j) - \overset{⇀}{F} (i + 1, j) ∥_{2}^{2} + ∥ \overset{⇀}{F} (i, j) - \\ \overset{⇀}{F} (i, j + 1) ∥_{2}^{2}) / 2]^{\frac{1}{2}} / ((Μ - 1) (Ν - 1)) \end{array}$

平均梯度是敏感反映图像对微小细节反差和纹理变化特征表达的能力,同时也反映了图像的清晰度。平均梯度越大,图像的层次越多,融合的图像就越清晰。

从表1可以看出,在信息熵、空间频率和平均梯度这几个反映图像融合效果的评价参数方面,本文提出的图像融合算法相比其他算法效果较好。

与图4相比,各种方法融合后的峰值信噪比PSNR对比如2所示。

由表2可以看出本文提出的方法对于彩色图像融合效果相对较好。

3 结语

多媒体和互联网的迅猛发展以及各类图像传感器的广泛使用,传统的图像处理方法己经不能适应日益剧增的图像源和人们对图像越来越高的精度要求。图像融合技术能够有效提高图像质量、提取更多的信息,特别是彩色图像融合。本文针对具有不同聚焦点的图像的特点,首先根据图像块区域的局部对比度对图像进行划分,其次给出图像块区域的清晰度的定义,在融合过程中对于那些属于清晰区域和模糊区域的直接选取源图像中清晰度高的对应部分,对于边界区域采用有理三次样条加权平均得到对应的图像融合部分。最后通过灰度图像和彩色图像实验证明这种方法是有效的。

摘要：针对多聚焦图像融合,从人眼视觉对比度以及多聚焦图像清晰与模糊像素之间的边界点出发,提出一种基于有理三次样条的图象融合算法。将多聚焦图像划成三个区域:清晰区域、模糊区域和边界区域。在边界区域,利用有理三次样条函数构造权函数,并以边界区域像素为中心,用清晰区域的像素和模糊区域像素做加权组合,计算边界区域像素的灰度值。经实验验证,提出的图像融合方法具有较强的灵活性和可靠性,具有较好的融合效果。

关键词：图像融合,多聚焦图像,有理三次样条,权函数,清晰度

参考文献

[1]覃征,鲍复民,李爱国,等.数字图像融合[M],西安:西安交通大学出版社,2004:1-60.

[2]Nichols L W,Mar J L.Conversion of infrared images to visible in color[J].Applied Optics,2008,7(9):1757-1762.

[3]Burt P J,Adelson E H.The Laplacian Pyramid as a Compact ImageCode[J].IEEE,Transcations on Communications,1983,31(4):532-540.

[4]Toet A,Ruyven L J.Merging and Visual Images by a Contract Pyramid[J].IEEE,Proceeding Visual Image Signal Process,1994,141(3):220-232.

[5]Akerman A.Pyramidal techniques for multi-sensor fusion[J].Senor Fu-sion,1992:124-131.

[6]王宏,敬忠良,李建勋.一种基于图像块分割的多聚焦图像融合方法[J].上海交通大学学报,2003,37(11):1743-1750.

[7]李树涛,王耀南,张昌凡.基于视觉特性的多聚焦图像融合[J].电子学报,2001,29(12):1699-1701.

[8]蒲恬,方庆赭,倪国强.基于对比度的多分辨图像融合[J].电子学报,2000,28(12):116-118.

【有理样条】推荐阅读：

有理函数05-12