最近有人问起二维三次卷积插值算法（Cubic Convolution Interpolation）及其程序的问题，我想到这种插值算法也能用于图像的插值，所以也就有了兴趣，稍微了解了一下，并试着用Fortran实现了一下，同时也顺便复习了一下Fortran90的一下新特性。如果需要，敬请使用；发现Bug，烦请明示。

值得注意的是，我在程序中将原本为Nx×Ny的格点数据向正反方向都增加了节点，变为0:Nx+2×0:Ny+2。这样做的目的是为了计算每个节点x，y方向导数的方便，不再需要区分边界点与非边界点，统一使用中心差分公式即可。扩展节点数据的取法保证了边界点的导数由最近两点的差分决定。

注意：
如果需要计算任意位置的导数，需要直接利用插值公式的导数来计算，而不能先求出每个节点的导数再对导数进行插值。

参考文献

1.    韩复兴, 孙建国, 杨昊. 基于二维三次卷积插值算法的波前构建射线追踪. 吉林大学学报(地球科学版), 38(2), 2008.

2.  李清, 侯永军, 沈春林. 数字地形数据的二维三次卷积插值. 南京航空航天大学学报, 29(4), 1997.

2012-07-24 20:59:20 初稿
2012-07-26 09:39:47 更新注意
2012-11-26 16:18:26 修正代码错误，默认三阶精度边界条件

==============================J*e*r*k*w*i*n*@*g*m*a*i*l*.*c*o*m=============================
  1| Program CubConv
  2|     integer i, j, NgrdX, NgrdY, Nx, Ny
  3|     real*8  Xmin, Ymin, Xmax, Ymax, dXgrd, dYgrd, x, y, dx, dy, Rtmp, &
  4|     &       F, CubConvInterp
  5|     real*8, allocatable:: Zxy(:,:)
  6|
  7| ! 对最高二次的函数，插值应该给出精确值，故用作测试
  8|     F(x, y) = (x*x+y*y)
  9|
 10| ! 原始网格设定
 11|     NgrdX = 3
 12|     NgrdY = 3
 13|     dXgrd = 2.D0
 14|     dYgrd = 2.D0
 15|     Xmin = -1.D0
 16|     Ymin = -1.D0
 17|     Xmax = Xmin+dble(NgrdX-1)*dXgrd
 18|     Ymax = Ymin+dble(NgrdY-1)*dYgrd
 19|
 20|     allocate(Zxy(0:NgrdX+2, 0:NgrdY+2))
 21|
 22| ! 获取网格节点
 23|     Zxy = 0.D0
 24|     do i=1, NgrdX
 25|         x = Xmin+dble(i-1)*dXgrd
 26|         do j=1, NgrdY
 27|             y = Ymin+dble(j-1)*dYgrd
 28|             Zxy(i, j) = F(x, y)
 29|         end do
 30|     end do
 31|
 32| ! 处理边界点，三阶精度
 33|     Zxy(0, 1:NgrdY) = 3.D0*( Zxy(1, 1:NgrdY)-Zxy(2, 1:NgrdY) )+Zxy(3, 1:NgrdY)
 34|     Zxy(1:NgrdX, 0) = 3.D0*( Zxy(1:NgrdX, 1)-Zxy(1:NgrdX, 2) )+Zxy(1:NgrdX, 3)
 35|
 36|     Zxy(NgrdX+1, 1:NgrdY) = 3.D0*(Zxy(NgrdX, 1:NgrdY)-Zxy(NgrdX-1, 1:NgrdY)) &
 37|     &                      +Zxy(NgrdX-2, 1:NgrdY)
 38|     Zxy(1:NgrdX, NgrdY+1) = 3.D0*(Zxy(1:NgrdX, NgrdY)-Zxy(1:NgrdX, NgrdY-1)) &
 39|     &                      +Zxy(1:NgrdX, NgrdY-2)
 40|
 41|     Zxy(0, 0) = 3.D0*( Zxy(1, 0)-Zxy(2, 0) )+Zxy(3, 0)
 42|     Zxy(NgrdX+1, 0) = 3.D0*( Zxy(NgrdX, 0)-Zxy(NgrdX-1, 0) )+Zxy(NgrdX-2, 0)
 43|
 44|     Zxy(0, NgrdY+1) = 3.D0*( Zxy(0, NgrdY)-Zxy(0, NgrdY-1) )+Zxy(0, NgrdY-2)
 45|     Zxy(NgrdX+1, NgrdY+1) = 3.D0*(Zxy(NgrdX, NgrdY+1)-Zxy(NgrdX-1, NgrdY+1)) &
 46|     &                      +Zxy(NgrdX-2, NgrdY+1)
 47|
 48| ! 改变网格，测试插值，误差应为0
 49|     Nx = 50
 50|     Ny = 50
 51|     dx = (Xmax-Xmin)/dble(Nx-1)
 52|     dy = (Ymax-Ymin)/dble(Ny-1)
 53|     do i=1, Nx
 54|         x = Xmin+dble(i-1)*dx
 55|         do j=1, Ny
 56|             y = Ymin+dble(j-1)*dy
 57|             Rtmp = CubConvInterp(x, y, Xmin, Ymin, dXgrd, dYgrd, NgrdX, NgrdY, Zxy)
 58|             write(*, '(5F18.9)') x, y, F(x, y), Rtmp, Rtmp-F(x, y)
 59|         end do
 60|     end do
 61|
 62| End Program CubConv
 63| !
 64| !
 65| Real*8  Function CubConvInterp(Xpos, Ypos, Xmin, Ymin, dXgrd, dYgrd, Nx, Ny, Zxy)
 66|     integer i, j, Nx, Ny
 67|     real*8  Xpos, Ypos, Xmin, Ymin, dXgrd, dYgrd, u, v, Zxy(0:Nx+2, 0:Ny+2), &
 68|     &       X(4), Y(4), C(4, 4), D(4, 4), DX(4), DY(4), CDX(4)
 69|
 70|     data D(1, 1:4) / -1.D0,  2.D0, -1.D0, 0.D0 /
 71|     data D(2, 1:4) /  3.D0, -5.D0,  0.D0, 2.D0 /
 72|     data D(3, 1:4) / -3.D0,  4.D0,  1.D0, 0.D0 /
 73|     data D(4, 1:4) /  1.D0, -1.D0,  0.D0, 0.D0 /
 74|
 75|     i = int((Xpos-Xmin)/dXgrd)+1
 76|     j = int((Ypos-Ymin)/dYgrd)+1
 77|     u = (Xpos-(Xmin+dble(i-1)*dXgrd))/dXgrd
 78|     v = (Ypos-(Ymin+dble(j-1)*dYgrd))/dYgrd
 79|
 80|     C(1, 1:4) = [ Zxy(i-1, j-1), Zxy(i, j-1), Zxy(i+1, j-1), Zxy(i+2, j-1) ]
 81|     C(2, 1:4) = [ Zxy(i-1, j  ), Zxy(i, j  ), Zxy(i+1, j  ), Zxy(i+2, j  ) ]
 82|     C(3, 1:4) = [ Zxy(i-1, j+1), Zxy(i, j+1), Zxy(i+1, j+1), Zxy(i+2, j+1) ]
 83|     C(4, 1:4) = [ Zxy(i-1, j+2), Zxy(i, j+2), Zxy(i+1, j+2), Zxy(i+2, j+2) ]
 84|
 85|     X(1:4) = [u*u*u, u*u, u, 1.D0]
 86|     Y(1:4) = [v*v*v, v*v, v, 1.D0]
 87|
 88|     DX  = matmul(D, X)
 89|     DY  = matmul(D, Y)
 90|     CDX = matmul(C, DX)
 91|     CubConvInterp = 0.25D0*dot_product(DY, CDX)
 92| End Function CubConvInterp
============================================================================================

转载本文请联系原作者获取授权，同时请注明本文来自李继存科学网博客。
链接地址：https://blog.sciencenet.cn/blog-548663-595297.html

上一篇：伟大的公式
下一篇：Paul的教导