二维三次卷积插值算法及Fortran代码

类别:    标签: fortran 算法   阅读次数:   版权: (CC) BY-NC-SA

最近有人问起二维三次卷积插值算法(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

注意

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

参考

  1. 韩复兴, 孙建国, 杨昊. 基于二维三次卷积插值算法的波前构建射线追踪. 吉林大学学报(地球科学版), 38(2), 2008.
  2. 李清, 侯永军, 沈春林. 数字地形数据的二维三次卷积插值. 南京航空航天大学学报, 29(4), 1997.
  3. 二维三次卷积插值算法
随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: Paul的教导
后一篇: 我的博客文本编辑器——EmEditor

访问人次(2015年7月 9日起): | 最后更新: 2017-06-23 03:12:17 UTC | 版权所有 © 2008 - 2017 Jerkwin