Geometry & Recognition

Interpolation은 이산적으로 주어진 샘플들 사이에서 존재하지 않는 값을 추정하는 작업이다.  한 지점에서 interpolation을 통해 값을 추정하기 위해서는 주변의 알려진 샘플 값들을 참조해야 하고 적절한 가중치를 주어야 한다. 또한 interpolation 함수는 충분히 부드럽게 주변 샘플 값들과 연결이 되어야 한다. 이런 관점에서 보면 interpolation은 주변 샘플에 적절한 가중치를 주는 kernel 함수와 convolution으로 이해할 수 있고, 또한 샘플링 과정에서 잃어버린 정보를 샘플 데이터를 smoothing해서 복원하는 과정으로 해석할 수 있다.

kernel 함수가 $K(x)$고, 일정한 간격으로 주어진 1차원 샘플의 값이 $\{(x_k,c_k)\}$일 때 interpolation 함수는 

$$f(x)=\sum _k{c}_{k}K(\frac{x-x_k}{h})$$

의 형태로 표현할 수 있다. $h$는 샘플의 간격으로 이미지에서는 픽셀의 간격이므로 $h=1$이다.

interpolation함수는 샘플 위치를 통과해야 하므로, 

$$f(x_j)=c_j=\sum _k{c}_{k}K(x_j-x_k)⟶K(x_j-x_k)={\delta }_{jk}$$

즉, $K(0)=1$이고, $K(i)=0(i\ne 0)$임을 알 수 있다. 또한 kernel이 주는 가중치의 합이 1이어야 하므로

$${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}K(x)dx=1$$이어야 한다.

영상처리에서 잘 알려진 interpolation 방법은 nearest-neighbor, linear, cubic interpolation이 많이 사용된다. 아래 그림은 kernel의 중심이 원점에 있는 경우를 그렸다.

Nearest neighbor:

$$K(x)=\{\begin{array}{ll}1,& |x|<\frac{1}{2}\\ 0,& |x|\ge \frac{1}{2}\end{array},H(\omega )={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}K(x)e^{-i\omega x}dx=\text{sinc}(\frac{\omega }{2})$$

Linear interpolation: 

$$K(x)=\{\begin{array}{ll}1-|x|,& |x|<1\\ 0,& |x|\ge 1\end{array},H(\omega )={\text{sinc}}^2(\frac{\omega }{2})$$

Cubic interpolation:

$$K(x)=\{\begin{array}{ll}\frac{1}{2}(3|x{|}^3-5|x{|}^2+2),& |x|<1\\ \frac{1}{2}(-|x{|}^3+5|x{|}^2-8|x|+4),& 1\le |x|<2\\ 0,& |x|\ge 2\end{array},$$

$$H(\omega )=\frac{12}{{\omega }^2}({\text{sinc}}^2(\frac{\omega }{2})-\text{sinc}(\omega ))$$$$-\frac{4}{{\omega }^2}(2{\text{sinc}}^2(\omega )-2\text{sinc}(\omega )-\text{sinc}(2\omega ))$$

Cubic interpolation kernel이 다른 nearest-neighbor, linear kernel 보다 low-pass 특성이 강해지므로 이미지가 더 잘 smoothing 되는 특성을 가질 것임을 알 수 있다.

double KernelCubic(double x) {
    double absx = fabs(x);
    if (absx < 1)  return 0.5 * (2 + absx * (absx * (3 * absx - 5)))
    if (absx < 2)  return 0.5 * (4 - absx * (absx * (absx + 5) - 8));
    return 0;
}

일반적인 cubic convolution kernel은 4개의 샘플데이터를 이용해 보간하므로 폭이 4(반지름 =2)이다. $x_k=-1,0,1,2$에서 4개의 샘플값 $c_0,c_1,c_2,c_3$이 주어진 경우, $0\le x<1$구간에서 cubic interpolation은  kernel의 중심을 원점에서 $x$로 평행이동하면,

$$f(x)=c_{0}K(-1-x)+c_{1}K(-x)+c_{2}K(1-x)+c_{3}K(2-x)$$

로 표현된다. 다음 예는 $0\le x<1$ 구간에서 $(x-1{)}^2$(blue)을 cubic interpolation(red)을 한 결과이다.

일반적인 cubic  convolution kernel은 한 개의 모양을 조정하는 한 개의 조정 parameter을 가지는 연속이고 한 번 미분가능한 형태로 주어진다.

$$K(x)=\{\begin{array}{ll}(a+2)|x{|}^3-(a+3)|x{|}^2+1,& |x|<1\\ a|x{|}^3-5a|x{|}^2+8a|x|-4a,& 1|x|<2\\ 0,& |x|\ge 2\end{array}$$

실용적인 파라미터의 범위는 $[-3,0]$이고( $a<-3$이면 $|x|<1$ 구간에서 단조감소가 아니어서 0이 아닌 곳에서 1보다 큰 값을 갖게 된다), 보통 많이 사용하는 cubic interpolation은 $a=-0.5$에 해당한다. 일반적으로 $a$가 -3에 가까워지면 최소 위치가 더 깊어지므로 일종의 미분 연산자처럼 작용하게 되어 이미지에 적용할 때 sharpness가 증가하게 된다. $a$가 0에 가까워지면 최소 위치 깊이가 0에 가까워지므로 gaussian filter처럼 이미지를 blurring하는 효과가 발생한다.

Lanczos Interpolation: 

$$K(x)=\{\begin{array}{ll}\text{sinc}(\pi x)\text{sinc}(\pi x/a),& |x|<a\\ 0,& |x|\ge a\end{array}$$

$$H(\omega )=\frac{1}{2{\pi }^2\sqrt{2\pi }}(-(2\pi -3\omega )\text{Si}(2\pi -3\omega )+(4\pi -3\omega )\text{Si}(4\pi -3\omega )$$ $$+(2\pi +3\omega )\text{Si}(2\pi +3\omega )+(4\pi +3\omega )\text{Si}(4\pi +3\omega ))$$

https://kipl.tistory.com/277

kipl.tistory.com/55

아날로그 신호로부터 디지털 데이터로 바꿀 때 보통 시간당(sound) 또는 공간적 거리당(image) 일정한 횟수로  데이터를 수집한다. 이 비율을 sampling rate이라 한다. 그런데 sampling rate을 바꾸어야 할 경우가 있다. 1초마다 한 번씩 데이터를 모았는데 실제로 0.5초마다 수집된 데이터가 필요가 한 경우에는 기존의 수집된 데이터와 데이터 사이에서 원 아날로그 신호 값을 알아내야 한다. sampling rate을 바꿀 때 기존의 데이터를 이용해서 원 아날로그 신호에서 수집한 것처럼 새로운 데이터를 만들어 내는 과정이 resampling이다. 이미지에서는 원본을  확대하는 up sampling이나 축소하는 down sampling 등이 resampling의 일종이다. resampling을 하기 위해서는 기존 데이터를 이용해서 데이터와 데이터 사이에서의 값을 추정하는 interpolation이 필요하게 된다. 많이 사용하는 interpolation 방법으로는 nearest neighbor interpolation, linear interpolation, cubic interpolation, spline interpolation 등이 있다.  여기서는 cubic interpolation을 이용한 이미지 resampling을 살펴본다.

$x$축 위의 4 점 $x=-1$, $x=0$, $x=1$, $x=2$에서 값이 $p_0$, $p_1$, $p_2$, $p_3$으로 주어진 경우 $0\le x\le 1$ 구간에서 cubic interpolation 함수는 

$$f(x)=a{x}^3+b{x}^2+cx+d$$

라고 쓰면, 우선 $x=0$에서 $p_1$, $x=1$에서 $p_2$를 지나므로

$$f(0)=d=p_1$$

$$f(1)=a+b+c+d=p_2$$

를 만족해야 하고, 양끝에서 도함수 값은 주변 입력점의 평균변화율로 주면, 다음식을 만족한다.

$$f^{\prime }(0)=c=\frac{p_2-p_0}{2}$$

$$f^{\prime }(1)=3a+2b+c=\frac{p_3-p_1}{2}$$

이를 이용해서 $a,b,c,d$를 구하면,

$$a=\frac{1}{2}(-p_0+3p_1-3p_2+p_3)$$

$$b=\frac{1}{2}(2p_0-5p_1+4p_2-p_3)$$

$$c=\frac{1}{2}(-p_0+p_2)$$

$$d=p_1.$$

따라서 cubic interpolation 함수는

$$\begin{array}{rl}f(x)& =\frac{1}{2}(-x^3+2x^2-x)p_0+\frac{1}{2}(3x^3-5x^2+2)p_1\\ & +\frac{1}{2}(-3x^3+4x^2+x)p_2+\frac{1}{2}(x^3-x^2)p_3\end{array}$$

로 쓰인다. 이는 다음처럼 정의된 kernel 함수를(Catmull-Rom filter) convolution한 형태로도 표현할 수 있다.

$$K(x)=\{\begin{array}{ll}\frac{1}{2}(3|x{|}^3-5|x{|}^2+2)& |x|<1\\ \frac{1}{2}(-|x{|}^3+5|x{|}^2-8|x|+4)& 1\le |x|<2\\ 0& \text{otherwise}\end{array}$$

$$f(x)=p_{0}K(-1-x)+p_{1}K(-x)+p_{2}K(1-x)+p_{3}K(2-x),(0\le x<1)$$

여기서는 좀 더 간단히 표현하기 위해서 $$f(x)=m_0{p}_0+m_1{p}_1+m_2{p}_2+m_3{p}_3$$ 로 표현한다. 물론 $m_i$는 $0\le x<1$의 함수이다. $f(x)$는 그림에서 보는 것처럼 $p_1$과 $p_2$를 통과한다.

구간 끝에서 도함수를 $$f^{\prime }(0)=\frac{p_2-p_0}{2},f^{\prime }(1)=\frac{p_3-p_1}{2}$$로 설정했기 때문에 인접 구간(eg. $[p_2,p_3]$)에서의 interpolation과 경계에서 같은  도함수를 가지게 되어 곡선이 꺽임없이 부드럽게 연결이 되는 특성을 가진다(예를 들면, $p_2$와 $p_3$을 보간하는 cubic interpolation함수는 $p_2$에서 미분값이 $(p_3-p_1)/2$인데, 이는 $p_1,p_2$구간의 interpolation  함수가 $p_2$에서 가지는 미분값과 같다).

이를 2차원 공간으로 확장하면, 평면 위의 격자점 $\{(x,y)|x,y=-1,0,1,2\}$에서 16개의 값 $\{p_{ij}|i,j=0,1,2,3\}$가 주어진 경우 정사각형 영역 $0\le x,y\le 1$에서 bicubic interpolation은 먼저 4개의 행 $y=-1,0,1,2$에서 각각 $x$ 방향의 cubic interpolation 된 값을 구한다:

$$\begin{array}{r}{q}_0=f(x,y=-1)=m_0{p}_{00}+m_1{p}_{10}+m_2{p}_{20}+m_3{p}_{30}\\ q_1=f(x,y=0)=m_0{p}_{01}+m_1{p}_{11}+m_2{p}_{21}+m_3{p}_{31}\\ q_2=f(x,y=1)=m_0{p}_{02}+m_1{p}_{12}+m_2{p}_{22}+m_3{p}_{32}\\ q_3=f(x,y=2)=m_0{p}_{03}+m_1{p}_{13}+m_2{p}_{23}+m_3{p}_{33}\end{array}$$

4개 행에서 구해서 값을 이용해서 다시 $y$ 방향으로 cubic interpolation 한 값을 구하면 bicubic interpolation이 완성된다: $$q=f(x,y)=q_0{n}_0+q_1{n}_1+q_2{n}_2+q_3{n}_3$$ 여기서 $n_i$는 $m_i$에서 $x$를 $y$로 치환한 값이다.

원본 이미지를 확대/축소 또는 변형 등의 변환 과정을 거쳐서 출력 이미지를 만들 때 원본 이미지의 픽셀 값을 resampling 하는 과정이 들어온다. 이때 원본 이미지의 픽셀과 픽셀 사이의 값이 필요한 경우가 발생하여 적절한 interpolation이 필요한데 많이 쓰이는 interpolation 중의 하나가 bicubic interpolation이다. 아래는 bicubic interpolation을 이용한 이미지 resampling을 수행하는 코드다.

interpolation에서 실수(float) 연산을 하지 않고 정수 연산만 사용하면서도 일정한 정밀도로 소수점 아래 부분의 정보를 유지하기 위해 나누기 전에 미리 256을 곱한 후에 나눈다(256 대신에 적당히 큰 수를 선택해도 되는데 2의 지수승을 잡으면 곱셈/나눗셈을 shift 연산으로 대체할 수 있는 장점이 있다). 이렇게 하면 나눈 몫에 $\mathtt{0}\mathtt{x}\mathtt{F}\mathtt{F}$ 비트 마스크를 적용할 때 남는 부분이 소수점 아랫부분을 표현한다. 정밀도는 $1/256$ 까지다. 중간 과정에서 소수점 이하 부분끼리 곱셈을 할 때는 항상 $256$으로 나누어서 $255$ 이하가 되도록 만들어야 한다. 최종 결과에서도 다시 $256$으로 나누기를 해주어야 된다. bicubic인 경우는 $x/y$ 양방향으로 적용하므로 $256\times 256$으로 나누어야 하고 cubic interpolation 계수에서 $1/2$이 들어오므로 추가로 $4$만큼 더 나누어 주어야 한다(코드의 마지막 결과에서 shift 연산 "$\mathtt{>>}\mathtt{18}$"이 들어온 이유다).

bicubic interpolation을 적용할 때 $\mathtt{y}\mathtt{=}\mathtt{0}$이나 $\mathtt{x}\mathtt{=}\mathtt{0}$에서는 이전 행이나 열이 없으므로 자신을 반복하는 방식으로 처리해 주어야 하고, 또 $\mathtt{y}\mathtt{=}\mathtt{r}\mathtt{o}\mathtt{w}\mathtt{s}\mathtt{-}\mathtt{2}\mathtt{,}\mathtt{r}\mathtt{o}\mathtt{w}\mathtt{s}\mathtt{-}\mathtt{1}$이나 $\mathtt{x}\mathtt{=}\mathtt{c}\mathtt{o}\mathtt{l}\mathtt{s}\mathtt{-}\mathtt{2}\mathtt{,}\mathtt{c}\mathtt{o}\mathtt{l}\mathtt{s}\mathtt{-}\mathtt{1}$일 때도 비슷한 처리가 있어야 한다.

int resample_bicubic ( BYTE *src, int cols, int rows,
                       BYTE *des, int newcols, int newrows ) {
    if (cols < 4 || rows < 4)
        return resample_bilinear(src, cols, rows, des, newcols, newrows);

    int ixn = cols - 4;       
    BYTE *pa, *pb, *pc, *pd;
    for (int j = 0; j < newrows; j++) {
        int yy = ( ( j * rows ) << 8 ) / newrows;
        int yp = yy >> 8;                        // src pixel y-position;
        int dy = yy & 0xFF;
        int dy2 = (dy * dy) >> 8;
        int dy3 = (dy2 * dy) >> 8;
        int n0 = -dy3 + 2 * dy2 - dy;
        int n1 = 3 * dy3 - 5 * dy2 + 2 * 256;
        int n2 = -3 * dy3 + 4 * dy2 + dy;
        int n3 = dy3 - dy2;
        //
        pb = src + yp * cols;                  //current line;
        if (yp == 0) pa = pb;
        else         pa = pb - cols;           //previous line;
        
        if (yy < rows - 2) {
            pc = pb + cols;                    //next line;
            pd = pc + cols;                    //next-next line;
        } else if (yp < rows - 1) {
            pc = pb + cols;        
            pd = pc;
        } else 
            pd = pc = pb;

        for (int i = 0; i < newcols; i++) {
            int xx = ( ( i * cols ) << 8 ) / newcols;
            int xp = xx >> 8;                    // src pixel x-position;
            int dx = xx & 0xFF;
            int dx2 = ( dx * dx ) >> 8;
            int dx3 = ( dx2 * dx ) >> 8;
            int m0 = -dx3 + 2 * dx2 - dx;
            int m1 = 3 * dx3 - 5 * dx2 + 2 * 256;
            int m2 = -3 * dx3 + 4 * dx2 + dx;
            int m3 = dx3 - dx2;
            int p = (xp == 0) ? 0 : (xp < ixn) ? xp - 1: ixn;    // p+3 <= ixn+3=cols-1;
            int a = ((m0 * pa[p] + m1 * pa[p + 1] + m2 * pa[p + 2] + m3 * pa[p + 3]) * n0 +
                     (m0 * pb[p] + m1 * pb[p + 1] + m2 * pb[p + 2] + m3 * pb[p + 3]) * n1 +
                     (m0 * pc[p] + m1 * pc[p + 1] + m2 * pc[p + 2] + m3 * pc[p + 3]) * n2 +
                     (m0 * pd[p] + m1 * pd[p + 1] + m2 * pd[p + 2] + m3 * pd[p + 3]) * n3) >> 18;
            *des++ = (a > 0xFF) ? 0xFF: (a < 0) ? 0: a;
        }
    }
	return 1;
}

bicubic interpolation을 하기 위해서는 4점이 필요하므로 폭이나 높이가 이보다 작은 경우는 bilinear interpolation을 사용한다. 다음 코드는 fixed-point bilinear interpolation을 구현한 코드다.

int resample_bilinear(BYTE *src, int cols, int rows,
                      BYTE *des, int newcols, int newrows ) {
    for (int j = 0; j < newrows; j++ ) {
        int yy = ((j * rows ) << 8) / newrows;
        int yp = yy >> 8;   // integer part; src y-position;
        int dy = yy & 0xFF; // fractional part;
        BYTE *curr = src + yp * cols;
        BYTE *next = curr + cols;
        for (int i = 0; i < newcols; i++) {
            int xx = ((i * cols ) << 8) / newcols;
            int xp = xx >> 8;       //src x-position;
            int dx = xx & 0xFF;
            int p00 = curr[xp];
            int p10 = curr[xp + 1];
            int p01 = next[xp];
            int p11 = next[xp + 1];
            int val = ((p11 * dx + p01 * (256 - dx)) * dy
                    + (p10 * dx + p00 * (256 - dx)) * (256 - dy)) >> 16;
            *des++ = val > 255 ? 0xFF: val < 0 ? 0 : val;
        }
    }
    return 1;
}

kipl.tistory.com/55

두 개의 점이 주어지는 경우 이들을 잇는 가장 단순한 곡선은 직선이다. 또 직선은 두 점을 잇는 가장 짧은 거리를 갖는 곡선이기도 하다. 그러면 $N$개의 점이 주어지는 경우는 이들 모두 지나는 곡선을 어떻게 찾을까? 구간별로 직선으로 연결을 시켜서 하나의 곡선을 구성할 수 있으나 부드럽게 이어지는 형태의 곡선은 아니다. 곡선이 smooth 함은 그 곡선의 곡률이 잘 정의된다는 의미이다(곡률은 곡선의 이차 미분에 의존한다). 주어진 점들을 연결하는 충분히 짧고 smooth 한 곡선을 찾는 문제는 곡률에 의한 에너지를 최소화시키는 해를 구하는 문제로 환원할 수 있다. 이 문제의 해는 잘 알려진 cubic spline이다. 그러나 cubic spline은 지나는 점(컨트롤 점)의 변동에 대해서 안정적이 아니다. 컨트롤 점 중에서 한 점만 변해도 곡선이 전역적(global)으로 변하는 특성을 나타낸다. 컨트롤 점의 변동이 있을 때 그 점 주위에서만 변화하는 국소 특성을 갖는 곡선을 찾기 위해서는 smoothness 조건을 완화해야 한다. 그렇지만 충분히 부드러운 곡선이어야 하므로 적어도 일차의 미분 값의 연속성 정도는 보장되어야 한다. 이런 특성을 갖는 삼차 곡선이 Catmull-Rom spline이다. 이름은 Edwin Catmull와 Raphie Rom의 이름에서 유래한다.

특징:

• 주어진 점들 위에서 정의됨(interpolation)
• 국소적인 변동 특성(노이즈에 안정적)
• 일차 미분의 연속성

두 점 $P_0$와 $P_1$을 지나는 일반적인 3차의 곡선을 적으면

$$P(t)=a_0+a_{1}t+a_2{t}^2+a_3{t}^3$$

이 곡선이 $t=0$에서 $P_0$을 지나고, $t=1$에서 $P_1$을 지난다고 하더라도, 완전히 결정되지 않는다. 따라서 추가적인 조건을 주어야 하는데, 여기서는 $P_0$와 $P_1$에서의 미분 값이 주어진다고 가정한다. 즉,

$$P(0)=P_0,P(1)=P_1,P^{\prime }(0)=P_0^{\prime },P^{\prime }(1)=P_1^{\prime }$$

이 조건에서 계수들은

$$\begin{gathered} a_0=P_0,a_1=P_0^{\prime }, \\ {a}_2=3(P_1-P_0)-2P_0^{\prime }-P_1^{\prime }, \\ {a}_3=2(P_0-P_1)+P_0^{\prime }+P_1^{\prime } \end{gathered}$$

로 주어진다;

$$P(t)=\left[\begin{array}{cccc}1& t& t^2& t^3\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ -3& 3& -2& -1\\ 2& -2& 1& 1\end{array}\right]\left[\begin{array}{c}{P}_0\\ P_1\\ P_0^{\prime }\\ P_1^{\prime }\end{array}\right]$$

그러나  미분 값을 직접 정해서 넣는 것은 문제가 있다.

4 이상의 점들이 주어지는 경우에는 컨트롤상에서의 미분 값을 그 점 주위의 두 점을 있는 직선의 기울기 값으로 근사 시킬 수 있다:

$$P_i^{\prime }⟶(P_{i+1}-P_{i-1})/2$$

이 미분 값을 사용하면 네 개의 컨트롤 점 $\{P_{i-1},P_i,P_{i+1},P_{i+2}\}$이 주어진 경우 $P_i(\leftarrow t=0)$와 $P_{i+1}(\leftarrow t=1)$ 사이를 보간하는 곡선은

$$\begin{array}{rl}P(t)& =\left[\begin{array}{cccc}1& t& t^2& t^3\end{array}\right]\left[\begin{array}{rrrr}1& 0& 0& 0\\ 0& 0& 1& 0\\ -3& 3& -2& -1\\ 2& -2& 1& 1\end{array}\right]\left[\begin{array}{c}{P}_i\\ P_{i+1}\\ (P_{i+1}-P_{i-1})/2\\ (P_{i+2}-P_i)/2\end{array}\right]\\ & =\frac{1}{2}\left[\begin{array}{cccc}1& t& t^2& t^3\end{array}\right]\left[\begin{array}{rrrr}0& 2& 0& 0\\ -1& 0& 1& 0\\ 2& -5& 4& -1\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{P}_{i-1}\\ P_i\\ P_{i+1}\\ P_{i+2}\end{array}\right]\end{array}$$

로 표현된다. 다시 정리하면

$$P(t)=\frac{1}{2}[2P_i+(-P_{i-1}+P_{i+1})t+(2P_{i-1}-5P_i+4P_{i+1}-P_{i+2})t^2$$

$$+(-P_{i-1}+3P_i-3P_{i+1}+P_{i+2})t^3]$$

로 쓸 수 있는데, 이는 각각의 입력점 사이 구간에서 매개변수를 이용해서 표현하기 좋은 꼴이다. 만약 매개변수를 knots으로 주어지는 경우는 다음과 같은 형식으로 바꾸어서 사용하는 것이 편리하다.

$$P(t)=b_0(t)P_{i-1}+b_1(t)P_i+b_2(t)P_{i+1}+b_3(t)P_{i+2}$$

$$\begin{array}{rl}{b}_0(t)& =\frac{1}{2}(-3t^3+2t^2-t)\\ b_1(t)& =\frac{1}{2}(3t^3-5t^2+2)\\ b_2(t)& =\frac{1}{2}(-3t^3+4t^2+t)\\ b_3(t)& =\frac{1}{2}(t^3-3t^2+3t-1)\end{array}$$

열린 곡선으로 보간을 하는 경우 처음과 끝 구간에는 미분 값을 구할 수 없으므로 정의가 안되지만, 나머지 구간에서는 주어진 점들을 연결하는 충분히 부드러운 곡선이 된다. 양 끝 구간의 보간은 동일점을 반복함으로써 전체적인 구간에서 잘 정의되게 만들 수 있다(끝점에서 기울기를 그 점과 인접점과의 기울기/2로 잡는 것과 같은 효과임). 닫힌 곡선인 경우에는 모든 구간에서 잘 정의가 된다. Catmull-Rom spline은 Bezier나 B-Spline 곡선의 경우와 다르게 컨트롤 점의 convex hull 내부에서만 정의되는 특성을 따르지 않는다.

CPoint CRSpline(double t, CPoint p1, CPoint p2, CPoint p3, CPoint p4) {
    double tt = t * t ;
    double ttt = tt * t ;
    double x = 0.5 * ((-p1.x + 3 * p2.x - 3 * p3.x + p4.x) * ttt
        + (2 * p1.x - 5 * p2.x + 4 * p3.x - p4.x) * tt
        + (-p1.x + p3.x) * t
        + 2 * p2.x);
    double y = 0.5 * ((-p1.y + 3 * p2.y - 3 * p3.y + p4.y) * ttt
        + (2 * p1.y - 5 * p2.y + 4 * p3.y - p4.y) * tt
        + (-p1.y + p3.y) * t
        + 2 * p2.y);
    return CPoint(int(x + .5), int(y + .5)) ;
}

//open spline의 drawing;
void DrawCatmullRom(std::vector<CPoint>& Q, CDC* pDC) {  
#define STEPS (20)
    if (Q.size() < 4) return ;
    CPen red(PS_SOLID, 1, RGB(0xFF, 0, 0));
    CPen *pOld = pDC->SelectObject(&red);
    const int n = Q.size();
    for (int i = 0; i < n - 1; i++) {
        pDC->MoveTo(Q[i]);
        // i = 0 인 경우에는 처음 점을 반복, i = n - 2인 경우에는 마지막 점을 반복..
        int ip = max(i - 1, 0);
        int inn = min(i + 2, n - 1);
        for (int it = 1; it <= STEPS; it++)
            pDC->LineTo(CRSpline(double(it)/STEPS, Q[ip], Q[i], Q[i + 1], Q[inn]));
    };
    pDC->SelectObject(pOld);
}

**네이버 블로그 이전;