파이썬 2차 보간법 완벽 예제!! - 면 보간, 면 피팅, interp2d, bisplrep

파이썬 3차 보간 애니메이션

파이썬 3차 보간 대표이미지

 
 

개요

    아래와 같이 9개의 점을 데이터로 측정했다고 치자. 그리고 보간법(Interpolation)을 통해 이 점들을 포함하는 하나의 곡면을 찾으려고 한다. 

XYZ 공간에 그려진 곡면에 대한 9개의 샘플링 포인트

   
 
    나는 10개 측정치에 대한 함수를 전혀 모르고, 이것들이 어떠한 곡면에 대한 데이터라는 것만 알기 때문에 이 점들의 개형을 최대한 잘 표현할 수 있는 곡면을 보간법으로 구해보려고 한다. 적용 알고리즘은 B-Spline Interpolation 이다.
 
 

B-Spline Interpolation

  B-Spline Interpolation에 관해서는 위키피디아(https://en.wikipedia.org/wiki/B-spline에 잘 정의되어 있으니 참고 바란다. 간단히 설명하자면, 어떤 복잡한 개형의 데이터를 한 번에 피팅?하지 않고 이를 조각조각 내어서 그 조각을 n차 함수로 각각 피팅한 것이라고 보면 된다. 

B-Spline 보간법, (출처: 위키피디아https://en.wikipedia.org/wiki/B-spline)

    위와 같이 B-Spline은 여러 개의 점들을 보간할 때 중간에 "매듭(빨강과 파랑 사이사이)"을 만들고 이 매듭은 n차 조각 함수들을 매끄럽게 이어주는 역할을 한다. 더 자세한 설명을 하면 좋을텐데 나는 수학은 전문가가 아니므로 일단 여기까지만 대충? 설명하도록 하고 코드로 넘어가보겠다!
 
 

코드구현

    우선 측정을 통해 얻은 샘플 포인트는 아래와 같다. 총 10개의 점이고 9개는 원점 대칭, 1개는 비대칭 점이다. 

import plotly.express as px 

xf = np.array([ 0.        ,  1.        ,  0.70710678,  0.        , -1.        ,
               -0.70710678, -0.70710678, -0.        ,  0.70710678,  0.7       ])
yf = np.array([ 0.        ,  0.        ,  0.70710678,  0.70710678,  0.        ,
                0.70710678, -0.70710678, -0.70710678, -0.70710678,  0.        ])
zf = np.array([ 0.    , -0.002 , -0.002 , -0.001 , -0.002 , -0.002 , -0.002 ,
               -0.001 , -0.002 , -0.0005])
               
fig = px.scatter_3d(x=xf,y=yf,z=zf, size=[2 for i in xf])
fig.show()

    위 코드를 실행한 결과는 아래와 같다. 주피터 노트북에서 실행 후 이리저리 돌려보자. 

B-spline Interpolation을 위한 10개의 샘플 포인트

 
    이제 이 점들을 부드럽게 이어줄 면을 구현해 볼 차례다. 
먼저 scipy 라이브러리의 interpolate 모듈에서 bisplrep과 bisplev(https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.bisplev.html)를 임포트 한다. 
bisplrep은 Bivariate B-Spline Interpolation Representation 함수이고, bisplev는 bisplrep으로 생성한 함수의 값을 evaluation 하는 함수이다.  
    즉, bisplrep으로 먼저 해당 포인트에 대한 tck[tx,ty, c, kx,ky]를 생성한 뒤 tck값을 bisplev함수에 좌표값과 함께 전달하면 보간된 함수의 z 값이 생성된다. 

반응형

from scipy.interpolate import bisplrep, bisplev

# 보간 함수에 xyz 값 입력, 
# kx, ky : 보간 차수 
# s : Smoothing Factor
tck = bisplrep(x=xf,y=yf,z=zf, kx=2, ky=2, s=0.1)

# 면 형성을 위해 새로운 x,y값 meshgrid 생성 
xx = np.linspace(-1,1,500)
yy = np.linspace(-1,1,500)
xxi, yyi = np.meshgrid(xx,yy, indexing='ij')

# bisplev에 새로운 x,y 값 입력을 위해 flatten 
xxf = xxi.flatten()
yyf = yyi.flatten()

# bisplev 함수에 새로운 x,y 값과 보간 함수 정보를 담은 tck값 전달 
zz = bisplev(xx, yy, tck)
zzf = zz.flatten()

# 3d scatter로 보간 상태 확인 
fig = px.scatter_3d(x=xf,y=yf,z=zf)
fig.update_traces(marker=dict(size=3, color='black'))
fig.show()

    이제 보간이 제대로 되었는지 plotly로 확인해보자. 

fig = go.Figure(data=[go.Surface(x=xxi, y=yyi, z=zz)])
fig.add_scatter3d(x=xf, y=yf, z=zf, mode='markers', marker=dict(size=2, color='black'))
fig.show()

면 보간 1차 결과

    일단 1차 결과는 원하는 면이 제대로 잡히지 않았고 뭔가 면이 부자연스럽다. (1) 실제 보간을 원하는 범위를 제대로 선택하고, (2) 보간 면의 차수(kx, ky)도 조절해보고자 한다. 차수를 올리려면 점이 더 필요하다. 따라서 대칭성을 이용하여 동일 위치에 같은 점을 몇 개 더 추가하였다. (약간의 편법이다.)

# 보간 차수(kx,ky)를 늘리기 위해 추가 점 생성 
xf = np.array([ 0.        ,  1.        ,  0.70710678,  0.        , -1.        ,
               -0.70710678, -0.        , -1.        , -0.70710678, -0.        ,
                1.        ,  0.70710678,  0.        ,  0.7       ])
yf = np.array([ 0.        ,  0.        ,  0.70710678,  0.70710678,  0.        ,
                0.70710678,  0.70710678,  0.        , -0.70710678, -0.70710678,
                0.        , -0.70710678, -0.70710678,  0.        ])
zf = np.array([ 0.    , -0.002 , -0.002 , -0.001 , -0.002 , -0.002 , -0.001 ,
               -0.002 , -0.002 , -0.001 , -0.002 , -0.002 , -0.001 , -0.0005])

    그리고 부자연스러운 면을 제거하기 위해 아래와 같이 실제 면의 바운더리를 설정하는 함수를 추가하였다. 바운더리는 아래 그림과 같이 지름이 1인 원형 + Y방향의 길이가 0.707..인 직사각형으로 된 면의 AND 조건이다.  

def IsInBoundary(x,y):
    r = 1
    if x**2+y**2 <= r**2:
        if -1 <= x <= 1 and -0.70718 <= y <= 0.70718:
            return True
    else:
        return False

   

면의 바운더리 설정 (matplotlib) : 직사각형 + 원형 AND 조건

    자, 이제 이 바운더리 내부의 점만 골라내고, B-Spline의 차수를 kx=3, ky=2, s=0.1로 하여 다시 보간을 해본다. 

# 보간 면 차수 kx=3, ky=2 재설정 
tck = bisplrep(x=xf,y=yf,z=zf, kx=3, ky=2, s=0.1)

# 신규 점에대해 meshgrid 생성 
xx = np.linspace(-1,1,200)
yy = np.linspace(-1,1,200)
xxi, yyi = np.meshgrid(xx,yy, indexing='ij')

xxf = xxi.flatten()
yyf = yyi.flatten()

# 보간된 새로운 좌표값 생성 
zz = bisplev(xx, yy, tck)
zzf = zz.flatten()

# 보간점들 중 바운더리 내부의 점만 유효값으로 인정 
for i, xt in enumerate(xx):
    for j, yt in enumerate(yy):
        if not IsInBoundary(xt,yt):
            xxi[i,j] = np.nan
            yyi[i,j] = np.nan
            zz[i,j] = 0

# 차트 그리기 
fig = go.Figure(data=[go.Surface(x=xxi, y=yyi, z=zz)])
fig.add_scatter3d(x=xf, y=yf, z=zf, mode='markers', marker=dict(size=2, color='black'))
fig.show()

    위 코드의 결과는 아래 그림과 같다. B-Spline이기 때문에 실제 측정치와 보간면 사이 값에 에러가 있을 수 있으나 대체로 잘 된 것 같아 보인다. 에러값을 평가하고 싶으면 bisplrep 함수의 인수에 full_output=1로 전달하고 tck, fp, ier, msg 네 개의 값 중 fp 값을 보면 된다. fp는 보간면의 weighted sum of squared residual(잔차에 대한 스퀘어 루트 가중평균) 값이다. 설명이 길어지기 때문에 여기서 별도로 다루지는 않겠다. 

 
 
 
 
    도움되셨다면 하트(♥) 부탁드리고, 더 궁금한 사항은 댓글로 남겨주세요 :)