Fortranで3次元の線形補完ルーチンを書く機会があったのでメモ。1次元のケースから確認し、2次元、3次元に拡張する流れで進める。

1次元の線形補間

１次元グリッド $x_1, x_2, ... x_n$ 上で値 $v_1, v_2, ... v_n$ がそれぞれ定義されているとする.

ここでは, $x_1$ と $x_2$ の間にある $x_i$ における値 $v_i$ を求めてみる.

まず距離として以下の量を定義する.

$\begin{align} l_1 &= x_i - x_1 \\ l_2 &= x_2 - x_i \\ l &= l_1 + l_2 \end{align}$

$(x_i,v_i)$ が $(x_1,v_1)$ と $(x_2,v_2)$ を結んだ直線上に位置しているので, $v_i$ は $v_1$ と $v_2$ の距離に逆比例する重み付き平均で求められる. 式で書けば

$\begin{align} v_i &= v_1 \frac{l_2}{l} + v_2 \frac{l_1}{l} \\ &= v_1 \frac{x_2-x_i}{x_2 - x_1} + v_2 \frac{x_i - x_1}{x_2-x_1} \end{align}$

となる. このイメージを持ったまま多次元に拡張すると分かりやすい.

２次元の線形補間

２次元に拡張した場合を考える.
まずx軸方向で１次元線形補間を行い(図中の灰色点), 次にy軸方向で１次元線形補間を行う(図中の赤色点)と考えれば良い. 式で書くと

$\begin{align} v_{ii} &= \frac{y_2 - y_i}{y_2 - y_1}\left( v_{11} \frac{x_2-x_i}{x_2 - x_1} + v_{21} \frac{x_i - x_1}{x_2-x_1} \right) + \frac{y_i - y_1}{y_2 - y_1}\left( v_{12} \frac{x_2-x_i}{x_2 - x_1} + v_{22} \frac{x_i - x_1}{x_2-x_1} \right) \\ &= \frac{1}{S} \{ (x_2-x_i)( y_2 - y_i )v_{11} + (x_i - x_1)( y_2 - y_i )v_{21} \\ &\quad + (x_2-x_i)(y_i - y_1)v_{12} + (x_i - x_1)(y_i - y_1)v_{22} \} \end{align}$

ただし $S=(x_2-x_1)(y_2-y_1)$ としている. 式を眺めると添字の関係性に規則性があることがわかる.

３次元の線形補間

３次元の図を書くのが大変なのであとは規則性に従って公式を導く.

$\begin{align} v_{ii} = \frac{1}{V} \{ (x_2-x_i)(y_2 - y_i )(z_2 - z_i )v_{111} + (x_2-x_i)(y_i - y_1 )(z_2 - z_i )v_{121}\\ + (x_i-x_1)(y_i-y_1)(z_2-z_i)v_{221} + (x_i-x_1)(y_2-y_i)(z_2-z_i)v_{211}\\ + (x_i-x_1)(y_2-y_i)(z_i-z_1)v_{212} + (x_i-x_1)(y_i-y_1)(z_i-z_1)v_{222}\\ + (x_2-x_i)(y_i-y_1)(z_i-z_1)v_{122} + (x_2-x_i)(y_2-y_i)(z_i-z_1)v_{112} \} \end{align}$

ただし $V=(x_2-x_1)(y_2-y_1)(z_2-z_1)$ としている.

コード

Fortranによる３次元の線形補間アルゴリズムを載せます.

module interpol

    implicit none
    
    type InterpolCoeffs3D
        integer :: i = 0
        integer :: j = 0
        integer :: k = 0
        real(8) :: b111 = 0.d0
        real(8) :: b121 = 0.d0
        real(8) :: b112 = 0.d0
        real(8) :: b122 = 0.d0
        real(8) :: b211 = 0.d0
        real(8) :: b212 = 0.d0
        real(8) :: b221 = 0.d0
        real(8) :: b222 = 0.d0
    end type InterpolCoeffs3D
    
    contains
    
    subroutine xyz_interpol3D_coeff(xmin,dx,nx,ymin,dy,ny,zmin,dz,nz,xout,yout,zout,c,err)
        real(8), intent(in)           :: xmin
        real(8), intent(in)           :: dx
        integer, intent(in)           :: nx
        real(8), intent(in)           :: ymin
        real(8), intent(in)           :: dy
        integer, intent(in)           :: ny
        real(8), intent(in)           :: zmin
        real(8), intent(in)           :: dz
        integer, intent(in)           :: nz
        real(8), intent(in)           :: xout
        real(8), intent(in)           :: yout
        real(8), intent(in)           :: zout
        type(InterpolCoeffs3D), intent(out) :: c
        integer, intent(out), optional      :: err
    
        real(8) :: x1, x2, y1, y2, z1, z2, xp, yp, zp, dV
        integer :: i, j, k, err_status
    
        err_status = 1
        xp = max(xout,xmin)
        yp = max(yout,ymin)
        zp = max(zout,zmin)
        i = floor((xp-xmin)/dx)+1
        j = floor((yp-ymin)/dy)+1
        k = floor((zp-zmin)/dz)+1
    
        if ((((i.gt.0).and.(i.le.(nx-1))).and.((j.gt.0).and.(j.le.(ny-1)))).and.((k.gt.0).and.(k.le.(nz-1)))) then
            x1 = xmin + (i-1)*dx
            x2 = x1 + dx
            y1 = ymin + (j-1)*dy
            y2 = y1 + dy
            z1 = zmin + (k-1)*dz
            z2 = z1 + dz
            dV = dx*dy*dz
    
            c%b111 = ((x2 - xp) * (y2 - yp) * (z2 - zp))/dV
            c%b121 = ((x2 - xp) * (yp - y1) * (z2 - zp))/dV
            c%b221 = ((xp - x1) * (yp - y1) * (z2 - zp))/dV
            c%b211 = ((xp - x1) * (y2 - yp) * (z2 - zp))/dV
            c%b212 = ((xp - x1) * (y2 - yp) * (zp - z1))/dV
            c%b222 = ((xp - x1) * (yp - y1) * (zp - z1))/dV
            c%b122 = ((x2 - xp) * (yp - y1) * (zp - z1))/dV
            c%b112 = ((x2 - xp) * (y2 - yp) * (zp - z1))/dV
            c%i = i
            c%j = j
            c%k = k
            err_status = 0
        endif
    
        if(present(err)) err = err_status
    
    end subroutine xyz_interpol3D_coeff
    
    end module interpol
    
    program test
        use interpol
        implicit none

        real(8) :: xmin,dx,ymin,dy,zmin,dz,xout,yout,zout,rnd,vout
        integer :: nx,ny,nz,err,ix,iy,iz,i,j,k
        type(InterpolCoeffs3D) :: c
        real(8),allocatable :: v(:,:,:)

        xmin = 1.0d0
        dx   = 0.5d0
        ymin = 2.0d0
        dy   = 0.4d0
        zmin = 3.0d0
        dz   = 0.6d0
        nx   = 10
        ny   = 15
        nz   = 5

        allocate(v(nx,ny,nz))

        do ix = 1,nx
            do iy = 1,ny
                do iz = 1,nz
                    call random_number(rnd)
                    v(ix,iy,iz) = rnd
                end do
            end do
        end do

        call xyz_interpol3D_coeff(xmin,dx,nx,ymin,dy,ny,zmin,dz,nz,xout,yout,zout,c,err)

        i = c%i
        j = c%j
        k = c%k
        vout = c%b111*v(i,j,k) + c%b121*v(i,j+1,k) + c%b221*v(i+1,j+1,k) + c%b211*v(i+1,j,k) + &
               c%b212*v(i+1,j,k+1) + c%b222*v(i+1,j+1,k+1) + c%b122*v(i,j+1,k+1) + c%b112*v(i,j,k+1)

        write(*,*) vout

    end program test

参考

リンク

どすえのブログ

ソフトウェア開発ブログ

Fortranによる３次元線形補間 | 3D Linear Interpolation

1次元の線形補間

２次元の線形補間

３次元の線形補間

コード

参考