- 2012-07-24 20:59:20 初稿
- 2012-07-26 09:39:47 更新注意
- 2012-11-26 16:18:26 修正代码错误, 默认三阶精度边界条件
最近有人问起二维三次卷积插值算法（Cubic Convolution Interpolation）及其程序的问题, 我想到这种插值算法也能用于图像的插值, 所以也就有了兴趣, 稍微了解了一下, 并试着用Fortran实现了一下, 同时也顺便复习了一下Fortran90的一下新特性. 如果需要, 敬请使用; 发现Bug, 烦请明示.
值得注意的是, 我在程序中将原本为Nx×Ny的格点数据向正反方向都增加了节点, 变为0:Nx+2×0:Ny+2. 这样做的目的是为了计算每个节点x, y方向导数的方便, 不再需要区分边界点与非边界点, 统一使用中心差分公式即可. 扩展节点数据的取法保证了边界点的导数由最近两点的差分决定.
代码
| fortran | |
|---|---|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 | Program CubConv
integer i, j, NgrdX, NgrdY, Nx, Ny
real*8 Xmin, Ymin, Xmax, Ymax, dXgrd, dYgrd, x, y, dx, dy, Rtmp, &
& F, CubConvInterp
real*8, allocatable:: Zxy(:,:)
! 对最高二次的函数，插值应该给出精确值，故用作测试
F(x, y) = (x*x+y*y)
! 原始网格设定
NgrdX = 3
NgrdY = 3
dXgrd = 2.D0
dYgrd = 2.D0
Xmin = -1.D0
Ymin = -1.D0
Xmax = Xmin+dble(NgrdX-1)*dXgrd
Ymax = Ymin+dble(NgrdY-1)*dYgrd
allocate(Zxy(0:NgrdX+2, 0:NgrdY+2))
! 获取网格节点
Zxy = 0.D0
do i=1, NgrdX
x = Xmin+dble(i-1)*dXgrd
do j=1, NgrdY
y = Ymin+dble(j-1)*dYgrd
Zxy(i, j) = F(x, y)
end do
end do
! 处理边界点，三阶精度
Zxy(0, 1:NgrdY) = 3.D0*( Zxy(1, 1:NgrdY)-Zxy(2, 1:NgrdY) )+Zxy(3, 1:NgrdY)
Zxy(1:NgrdX, 0) = 3.D0*( Zxy(1:NgrdX, 1)-Zxy(1:NgrdX, 2) )+Zxy(1:NgrdX, 3)
Zxy(NgrdX+1, 1:NgrdY) = 3.D0*(Zxy(NgrdX, 1:NgrdY)-Zxy(NgrdX-1, 1:NgrdY)) &
& +Zxy(NgrdX-2, 1:NgrdY)
Zxy(1:NgrdX, NgrdY+1) = 3.D0*(Zxy(1:NgrdX, NgrdY)-Zxy(1:NgrdX, NgrdY-1)) &
& +Zxy(1:NgrdX, NgrdY-2)
Zxy(0, 0) = 3.D0*( Zxy(1, 0)-Zxy(2, 0) )+Zxy(3, 0)
Zxy(NgrdX+1, 0) = 3.D0*( Zxy(NgrdX, 0)-Zxy(NgrdX-1, 0) )+Zxy(NgrdX-2, 0)
Zxy(0, NgrdY+1) = 3.D0*( Zxy(0, NgrdY)-Zxy(0, NgrdY-1) )+Zxy(0, NgrdY-2)
Zxy(NgrdX+1, NgrdY+1) = 3.D0*(Zxy(NgrdX, NgrdY+1)-Zxy(NgrdX-1, NgrdY+1)) &
& +Zxy(NgrdX-2, NgrdY+1)
! 改变网格，测试插值，误差应为0
Nx = 50
Ny = 50
dx = (Xmax-Xmin)/dble(Nx-1)
dy = (Ymax-Ymin)/dble(Ny-1)
do i=1, Nx
x = Xmin+dble(i-1)*dx
do j=1, Ny
y = Ymin+dble(j-1)*dy
Rtmp = CubConvInterp(x, y, Xmin, Ymin, dXgrd, dYgrd, NgrdX, NgrdY, Zxy)
write(*, '(5F18.9)') x, y, F(x, y), Rtmp, Rtmp-F(x, y)
end do
end do
End Program CubConv
!
!
Real*8 Function CubConvInterp(Xpos, Ypos, Xmin, Ymin, dXgrd, dYgrd, Nx, Ny, Zxy)
integer i, j, Nx, Ny
real*8 Xpos, Ypos, Xmin, Ymin, dXgrd, dYgrd, u, v, Zxy(0:Nx+2, 0:Ny+2), &
& X(4), Y(4), C(4, 4), D(4, 4), DX(4), DY(4), CDX(4)
data D(1, 1:4) / -1.D0, 2.D0, -1.D0, 0.D0 /
data D(2, 1:4) / 3.D0, -5.D0, 0.D0, 2.D0 /
data D(3, 1:4) / -3.D0, 4.D0, 1.D0, 0.D0 /
data D(4, 1:4) / 1.D0, -1.D0, 0.D0, 0.D0 /
i = int((Xpos-Xmin)/dXgrd)+1
j = int((Ypos-Ymin)/dYgrd)+1
u = (Xpos-(Xmin+dble(i-1)*dXgrd))/dXgrd
v = (Ypos-(Ymin+dble(j-1)*dYgrd))/dYgrd
C(1, 1:4) = [ Zxy(i-1, j-1), Zxy(i, j-1), Zxy(i+1, j-1), Zxy(i+2, j-1) ]
C(2, 1:4) = [ Zxy(i-1, j ), Zxy(i, j ), Zxy(i+1, j ), Zxy(i+2, j ) ]
C(3, 1:4) = [ Zxy(i-1, j+1), Zxy(i, j+1), Zxy(i+1, j+1), Zxy(i+2, j+1) ]
C(4, 1:4) = [ Zxy(i-1, j+2), Zxy(i, j+2), Zxy(i+1, j+2), Zxy(i+2, j+2) ]
X(1:4) = [u*u*u, u*u, u, 1.D0]
Y(1:4) = [v*v*v, v*v, v, 1.D0]
DX = matmul(D, X)
DY = matmul(D, Y)
CDX = matmul(C, DX)
CubConvInterp = 0.25D0*dot_product(DY, CDX)
End Function CubConvInterp
|
注意
如果需要计算任意位置的导数, 需要直接利用插值公式的导数来计算, 而不能先求出每个节点的导数再对导数进行插值.
参考
- 韩复兴, 孙建国, 杨昊. 基于二维三次卷积插值算法的波前构建射线追踪. 吉林大学学报(地球科学版), 38(2), 2008.
- 李清, 侯永军, 沈春林. 数字地形数据的二维三次卷积插值. 南京航空航天大学学报, 29(4), 1997.
- 二维三次卷积插值算法
