仓库源文，站点原文

layout: post title: 二维三次卷积插值算法及Fortran代码 categories:

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

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

代码

<table class="highlighttable"><th colspan="2">fortran</th><tr><td><div class="linenodiv" style="background-color: #f0f0f0; padding-right: 10px"><pre style="line-height: 125%"> 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</pre></div></td><td class="code"><div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span style="color: #AA22FF; font-weight: bold">Program </span>CubConv <span style="color: #00BB00; font-weight: bold">integer </span>i, j, NgrdX, NgrdY, Nx, Ny <span style="color: #00BB00; font-weight: bold">real</span><span style="color: #666666">*8</span> Xmin, Ymin, Xmax, Ymax, dXgrd, dYgrd, x, y, dx, dy, Rtmp, & & F, CubConvInterp <span style="color: #00BB00; font-weight: bold">real</span><span style="color: #666666">*8</span>, <span style="color: #AA22FF; font-weight: bold">allocatable::</span> Zxy(:,:) <span style="color: #008800; font-style: italic">! 对最高二次的函数，插值应该给出精确值，故用作测试</span> F(x, y) <span style="color: #666666">=</span> (x<span style="color: #666666">*</span>x<span style="color: #666666">+</span>y<span style="color: #666666">*</span>y) <span style="color: #008800; font-style: italic">! 原始网格设定</span> NgrdX <span style="color: #666666">=</span> <span style="color: #666666">3</span> NgrdY <span style="color: #666666">=</span> <span style="color: #666666">3</span> dXgrd <span style="color: #666666">=</span> <span style="color: #666666">2.</span>D0 dYgrd <span style="color: #666666">=</span> <span style="color: #666666">2.</span>D0 Xmin <span style="color: #666666">=</span> <span style="color: #666666">-1.</span>D0 Ymin <span style="color: #666666">=</span> <span style="color: #666666">-1.</span>D0 Xmax <span style="color: #666666">=</span> Xmin<span style="color: #666666">+</span><span style="color: #AA22FF">dble</span>(NgrdX<span style="color: #666666">-1</span>)<span style="color: #666666">*</span>dXgrd Ymax <span style="color: #666666">=</span> Ymin<span style="color: #666666">+</span><span style="color: #AA22FF">dble</span>(NgrdY<span style="color: #666666">-1</span>)<span style="color: #666666">*</span>dYgrd <span style="color: #AA22FF; font-weight: bold">allocate</span>(Zxy(<span style="color: #666666">0</span>:NgrdX<span style="color: #666666">+2</span>, <span style="color: #666666">0</span>:NgrdY<span style="color: #666666">+2</span>)) <span style="color: #008800; font-style: italic">! 获取网格节点</span> Zxy <span style="color: #666666">=</span> <span style="color: #666666">0.</span>D0 <span style="color: #AA22FF; font-weight: bold">do </span>i<span style="color: #666666">=1</span>, NgrdX x <span style="color: #666666">=</span> Xmin<span style="color: #666666">+</span><span style="color: #AA22FF">dble</span>(i<span style="color: #666666">-1</span>)<span style="color: #666666">*</span>dXgrd <span style="color: #AA22FF; font-weight: bold">do </span>j<span style="color: #666666">=1</span>, NgrdY y <span style="color: #666666">=</span> Ymin<span style="color: #666666">+</span><span style="color: #AA22FF">dble</span>(j<span style="color: #666666">-1</span>)<span style="color: #666666">*</span>dYgrd Zxy(i, j) <span style="color: #666666">=</span> F(x, y) <span style="color: #AA22FF; font-weight: bold">end do</span> <span style="color: #AA22FF; font-weight: bold"> end do</span> <span style="color: #008800; font-style: italic">! 处理边界点，三阶精度</span> Zxy(<span style="color: #666666">0</span>, <span style="color: #666666">1</span>:NgrdY) <span style="color: #666666">=</span> <span style="color: #666666">3.</span>D0<span style="color: #666666">*</span>( Zxy(<span style="color: #666666">1</span>, <span style="color: #666666">1</span>:NgrdY)<span style="color: #666666">-</span>Zxy(<span style="color: #666666">2</span>, <span style="color: #666666">1</span>:NgrdY) )<span style="color: #666666">+</span>Zxy(<span style="color: #666666">3</span>, <span style="color: #666666">1</span>:NgrdY) Zxy(<span style="color: #666666">1</span>:NgrdX, <span style="color: #666666">0</span>) <span style="color: #666666">=</span> <span style="color: #666666">3.</span>D0<span style="color: #666666">*</span>( Zxy(<span style="color: #666666">1</span>:NgrdX, <span style="color: #666666">1</span>)<span style="color: #666666">-</span>Zxy(<span style="color: #666666">1</span>:NgrdX, <span style="color: #666666">2</span>) )<span style="color: #666666">+</span>Zxy(<span style="color: #666666">1</span>:NgrdX, <span style="color: #666666">3</span>) Zxy(NgrdX<span style="color: #666666">+1</span>, <span style="color: #666666">1</span>:NgrdY) <span style="color: #666666">=</span> <span style="color: #666666">3.</span>D0<span style="color: #666666">*</span>(Zxy(NgrdX, <span style="color: #666666">1</span>:NgrdY)<span style="color: #666666">-</span>Zxy(NgrdX<span style="color: #666666">-1</span>, <span style="color: #666666">1</span>:NgrdY)) & & <span style="color: #666666">+</span>Zxy(NgrdX<span style="color: #666666">-2</span>, <span style="color: #666666">1</span>:NgrdY) Zxy(<span style="color: #666666">1</span>:NgrdX, NgrdY<span style="color: #666666">+1</span>) <span style="color: #666666">=</span> <span style="color: #666666">3.</span>D0<span style="color: #666666">*</span>(Zxy(<span style="color: #666666">1</span>:NgrdX, NgrdY)<span style="color: #666666">-</span>Zxy(<span style="color: #666666">1</span>:NgrdX, NgrdY<span style="color: #666666">-1</span>)) & & <span style="color: #666666">+</span>Zxy(<span style="color: #666666">1</span>:NgrdX, NgrdY<span style="color: #666666">-2</span>) Zxy(<span style="color: #666666">0</span>, <span style="color: #666666">0</span>) <span style="color: #666666">=</span> <span style="color: #666666">3.</span>D0<span style="color: #666666">*</span>( Zxy(<span style="color: #666666">1</span>, <span style="color: #666666">0</span>)<span style="color: #666666">-</span>Zxy(<span style="color: #666666">2</span>, <span style="color: #666666">0</span>) )<span style="color: #666666">+</span>Zxy(<span style="color: #666666">3</span>, <span style="color: #666666">0</span>) Zxy(NgrdX<span style="color: #666666">+1</span>, <span style="color: #666666">0</span>) <span style="color: #666666">=</span> <span style="color: #666666">3.</span>D0<span style="color: #666666">*</span>( Zxy(NgrdX, <span style="color: #666666">0</span>)<span style="color: #666666">-</span>Zxy(NgrdX<span style="color: #666666">-1</span>, <span style="color: #666666">0</span>) )<span style="color: #666666">+</span>Zxy(NgrdX<span style="color: #666666">-2</span>, <span style="color: #666666">0</span>) Zxy(<span style="color: #666666">0</span>, NgrdY<span style="color: #666666">+1</span>) <span style="color: #666666">=</span> <span style="color: #666666">3.</span>D0<span style="color: #666666">*</span>( Zxy(<span style="color: #666666">0</span>, NgrdY)<span style="color: #666666">-</span>Zxy(<span style="color: #666666">0</span>, NgrdY<span style="color: #666666">-1</span>) )<span style="color: #666666">+</span>Zxy(<span style="color: #666666">0</span>, NgrdY<span style="color: #666666">-2</span>) Zxy(NgrdX<span style="color: #666666">+1</span>, NgrdY<span style="color: #666666">+1</span>) <span style="color: #666666">=</span> <span style="color: #666666">3.</span>D0<span style="color: #666666">*</span>(Zxy(NgrdX, NgrdY<span style="color: #666666">+1</span>)<span style="color: #666666">-</span>Zxy(NgrdX<span style="color: #666666">-1</span>, NgrdY<span style="color: #666666">+1</span>)) & & <span style="color: #666666">+</span>Zxy(NgrdX<span style="color: #666666">-2</span>, NgrdY<span style="color: #666666">+1</span>) <span style="color: #008800; font-style: italic">! 改变网格，测试插值，误差应为0</span> Nx <span style="color: #666666">=</span> <span style="color: #666666">50</span> Ny <span style="color: #666666">=</span> <span style="color: #666666">50</span> dx <span style="color: #666666">=</span> (Xmax<span style="color: #666666">-</span>Xmin)<span style="color: #666666">/</span><span style="color: #AA22FF">dble</span>(Nx<span style="color: #666666">-1</span>) dy <span style="color: #666666">=</span> (Ymax<span style="color: #666666">-</span>Ymin)<span style="color: #666666">/</span><span style="color: #AA22FF">dble</span>(Ny<span style="color: #666666">-1</span>) <span style="color: #AA22FF; font-weight: bold">do </span>i<span style="color: #666666">=1</span>, Nx x <span style="color: #666666">=</span> Xmin<span style="color: #666666">+</span><span style="color: #AA22FF">dble</span>(i<span style="color: #666666">-1</span>)<span style="color: #666666">*</span>dx <span style="color: #AA22FF; font-weight: bold">do </span>j<span style="color: #666666">=1</span>, Ny y <span style="color: #666666">=</span> Ymin<span style="color: #666666">+</span><span style="color: #AA22FF">dble</span>(j<span style="color: #666666">-1</span>)<span style="color: #666666">*</span>dy Rtmp <span style="color: #666666">=</span> CubConvInterp(x, y, Xmin, Ymin, dXgrd, dYgrd, NgrdX, NgrdY, Zxy) <span style="color: #AA22FF; font-weight: bold">write</span>(<span style="color: #666666">*</span>, <span style="color: #BB4444">'</span>(<span style="color: #666666">5</span>F18.<span style="color: #666666">9</span>)<span style="color: #BB4444">'</span>) x, y, F(x, y), Rtmp, Rtmp<span style="color: #666666">-</span>F(x, y) <span style="color: #AA22FF; font-weight: bold">end do</span> <span style="color: #AA22FF; font-weight: bold"> end do</span> <span style="color: #AA22FF; font-weight: bold">End Program </span>CubConv <span style="color: #008800; font-style: italic">!</span> <span style="color: #008800; font-style: italic">!</span> <span style="color: #00BB00; font-weight: bold">Real</span><span style="color: #666666">*8</span> <span style="color: #AA22FF; font-weight: bold">Function </span>CubConvInterp(Xpos, Ypos, Xmin, Ymin, dXgrd, dYgrd, Nx, Ny, Zxy) <span style="color: #00BB00; font-weight: bold">integer </span>i, j, Nx, Ny <span style="color: #00BB00; font-weight: bold">real</span><span style="color: #666666">*8</span> Xpos, Ypos, Xmin, Ymin, dXgrd, dYgrd, u, v, Zxy(<span style="color: #666666">0</span>:Nx<span style="color: #666666">+2</span>, <span style="color: #666666">0</span>:Ny<span style="color: #666666">+2</span>), & & X(<span style="color: #666666">4</span>), Y(<span style="color: #666666">4</span>), C(<span style="color: #666666">4</span>, <span style="color: #666666">4</span>), D(<span style="color: #666666">4</span>, <span style="color: #666666">4</span>), DX(<span style="color: #666666">4</span>), DY(<span style="color: #666666">4</span>), CDX(<span style="color: #666666">4</span>) <span style="color: #AA22FF; font-weight: bold">data </span>D(<span style="color: #666666">1</span>, <span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">/</span> <span style="color: #666666">-1.</span>D0, <span style="color: #666666">2.</span>D0, <span style="color: #666666">-1.</span>D0, <span style="color: #666666">0.</span>D0 <span style="color: #666666">/</span> <span style="color: #AA22FF; font-weight: bold">data </span>D(<span style="color: #666666">2</span>, <span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">/</span> <span style="color: #666666">3.</span>D0, <span style="color: #666666">-5.</span>D0, <span style="color: #666666">0.</span>D0, <span style="color: #666666">2.</span>D0 <span style="color: #666666">/</span> <span style="color: #AA22FF; font-weight: bold">data </span>D(<span style="color: #666666">3</span>, <span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">/</span> <span style="color: #666666">-3.</span>D0, <span style="color: #666666">4.</span>D0, <span style="color: #666666">1.</span>D0, <span style="color: #666666">0.</span>D0 <span style="color: #666666">/</span> <span style="color: #AA22FF; font-weight: bold">data </span>D(<span style="color: #666666">4</span>, <span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">/</span> <span style="color: #666666">1.</span>D0, <span style="color: #666666">-1.</span>D0, <span style="color: #666666">0.</span>D0, <span style="color: #666666">0.</span>D0 <span style="color: #666666">/</span> i <span style="color: #666666">=</span> <span style="color: #AA22FF">int</span>((Xpos<span style="color: #666666">-</span>Xmin)<span style="color: #666666">/</span>dXgrd)<span style="color: #666666">+1</span> j <span style="color: #666666">=</span> <span style="color: #AA22FF">int</span>((Ypos<span style="color: #666666">-</span>Ymin)<span style="color: #666666">/</span>dYgrd)<span style="color: #666666">+1</span> u <span style="color: #666666">=</span> (Xpos<span style="color: #666666">-</span>(Xmin<span style="color: #666666">+</span><span style="color: #AA22FF">dble</span>(i<span style="color: #666666">-1</span>)<span style="color: #666666">*</span>dXgrd))<span style="color: #666666">/</span>dXgrd v <span style="color: #666666">=</span> (Ypos<span style="color: #666666">-</span>(Ymin<span style="color: #666666">+</span><span style="color: #AA22FF">dble</span>(j<span style="color: #666666">-1</span>)<span style="color: #666666">*</span>dYgrd))<span style="color: #666666">/</span>dYgrd C(<span style="color: #666666">1</span>, <span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">=</span> [ Zxy(i<span style="color: #666666">-1</span>, j<span style="color: #666666">-1</span>), Zxy(i, j<span style="color: #666666">-1</span>), Zxy(i<span style="color: #666666">+1</span>, j<span style="color: #666666">-1</span>), Zxy(i<span style="color: #666666">+2</span>, j<span style="color: #666666">-1</span>) ] C(<span style="color: #666666">2</span>, <span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">=</span> [ Zxy(i<span style="color: #666666">-1</span>, j ), Zxy(i, j ), Zxy(i<span style="color: #666666">+1</span>, j ), Zxy(i<span style="color: #666666">+2</span>, j ) ] C(<span style="color: #666666">3</span>, <span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">=</span> [ Zxy(i<span style="color: #666666">-1</span>, j<span style="color: #666666">+1</span>), Zxy(i, j<span style="color: #666666">+1</span>), Zxy(i<span style="color: #666666">+1</span>, j<span style="color: #666666">+1</span>), Zxy(i<span style="color: #666666">+2</span>, j<span style="color: #666666">+1</span>) ] C(<span style="color: #666666">4</span>, <span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">=</span> [ Zxy(i<span style="color: #666666">-1</span>, j<span style="color: #666666">+2</span>), Zxy(i, j<span style="color: #666666">+2</span>), Zxy(i<span style="color: #666666">+1</span>, j<span style="color: #666666">+2</span>), Zxy(i<span style="color: #666666">+2</span>, j<span style="color: #666666">+2</span>) ] X(<span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">=</span> [u<span style="color: #666666">*</span>u<span style="color: #666666">*</span>u, u<span style="color: #666666">*</span>u, u, <span style="color: #666666">1.</span>D0] Y(<span style="color: #666666">1</span>:<span style="color: #666666">4</span>) <span style="color: #666666">=</span> [v<span style="color: #666666">*</span>v<span style="color: #666666">*</span>v, v<span style="color: #666666">*</span>v, v, <span style="color: #666666">1.</span>D0] DX <span style="color: #666666">=</span> <span style="color: #AA22FF">matmul</span>(D, X) DY <span style="color: #666666">=</span> <span style="color: #AA22FF">matmul</span>(D, Y) CDX <span style="color: #666666">=</span> <span style="color: #AA22FF">matmul</span>(C, DX) CubConvInterp <span style="color: #666666">=</span> <span style="color: #666666">0.25</span>D0<span style="color: #666666">*</span><span style="color: #AA22FF">dot_product</span>(DY, CDX) <span style="color: #AA22FF; font-weight: bold">End Function </span>CubConvInterp </pre></div> </td></tr></table>

注意

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

参考

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