주어진 데이터를 interpolation을 할 때 3차 spline을 많이 사용한다. 그럼 왜 3차 일까? 데이터를 연결하는 spline은 되도록이면 부드럽게 연결되어야 한다. 곡선이 부드럽게 그려지기 위해서는 급격한 꺾임이 없어야 하는데 얼마나 급격히 꺾이는가는 곡률에 비례하고, 곡률은 그 지점에서 함수의 두 번 미분값에 비례한다.(곡선의 곡률: https://kipl.tistory.com/105) 따라서 어떤 함수 $f(x)$가 구간 $[a,b]$에서 얼마나 부드럽게 연결되는가에 대한 척도는 곡선의 각부분에서 곡률의 크기의 제곱을 를 더한 다음 (음수가 아닌) 적분이 제공할 수 있다.
$$ \kappa[f] = \int_a^ b( f'' (x))^2 dx $$
$a \le x \le b$ 구간에서 균일하게 샘플링된 데이터가 있고 이들을 보간하는 3차 spline이 $S(x)=\{S_i(x) = a_i x^3 + b_i x^2 + c_i x +d_i\} $라고 하자. 또 두 번 이상 미분가능한 임의의 함수 $f(x)$도 주어진 샘플링 데이터를 지나가는 보간함수라고 하자. 이 경우 natural cubic spline이 주어진 sampling 데이터를 지나가면서 가장 부드럽게 이어지는 곡선임을 보일 수 있다.
$$ \kappa [f] \ge \kappa [S] \qquad \forall~ f$$
Spline $S$와 함수 $f$의 차이를 $h(x)= f(x) - S(x)$라 하면 각 node에서 $h(x_i)=0$이다. 이제
일정한 간격 $h$마다 샘플링된 데이터 $\{ (x_k, f_k) \}$를 이용해서 이들 데이터를 표현하는 spline를 구해보자. spline은 주어진 샘플링 데이터을 통과할 필요는 없으므로 일반적으로 interpolation 함수는 아니다. 이 spline은 샘플링 데이터와 kernel이라고 불리는 함수의 convolution 형태로 표현할 수 있다.
$$ g(x) = \sum_k f_k K \left( \frac{x-x_k}{h}\right)$$
이미지의 resampling 과정에서 spline를 이용하는데 이때 사용 가능한 kernel의 형태와 그 효과를 간단히 알아보자.
3차 spline kernel은 중심을 기준으로 반지름이 2인 영역 $(-2,1),(-1,0), (0,1), (1,2)$에서만 0이 아닌 piecewise 삼차함수다. 그리고 이 함수는 우함수의 특성을 갖는다. 따라서 가능한 형태는
$$ K(s) = \left\{ \begin{matrix} A_1|s|^3 + B_1 |s|^2 +C_1 |s| + D_1 & |s| <1 \\ A_2 |s|^3 + B_2 |s|^2 + C_2 |s| + D_2 & 1 \le |s|<2 \\ 0 & \text{otherwise} \end{matrix} \right. $$처럼 쓸 수 있다. 계수를 완전히 결정하기 위해서는 8개의 조건이 필요한다. 우선 우함수이므로 원점에서 미분값이 제대로 정의되려면 $$C_1=0$$도 만족해야 한다, 그리고 각 node에서 연속성을 요구하면
\begin{align} s=1^\pm:~~~& A_1+B_1 +D_1 = A_2 +B_2+C_2 +D_2 \\ s=2^\pm:~~~& 8A_2 +4B_2 +2C_2 +D_2 =0 \end{align} 임을 알 수 있다. 또한 각 node에서 부드럽게 연결되기 위해서 1차 도함수가 연속적임을 요구하면
$B=0$인 경우는 $s=0$일 때 1이고, $|s|=1,2$일 0이므로 interpolation kernel($K(i-j)=\delta_{ij}$)에 해당한다. 그리고 $B=0, C=1/2$인 경우인 Catmul-Rom spline은 node에서 2차 도함수까지도 연속이므로 샘플링 데이터를 생성한 원 아날로그 함수에 $O(h^3)$이내에서 가장 유사하게 근사함을 보일 수도 있다.
// Mitchell Netravali Reconstruction Filter
// B = 0 C = 0 - Hermite B-Spline interpolator
// B = 0, C = 1/2 - Catmull-Rom spline
// B = 1/3, C = 1/3 - Mitchell Netravali spline
// B = 1, C = 0 - cubic B-spline
double MitchellNetravali(double x, double B, double C) {
x = fabs(x);
if (x >= 2) return 0;
double xx = x*x;
if (x >= 1) return ((-B - 6*C)*xx*x
+ (6*B + 30*C)*xx + (-12*B - 48*C)*x
+ (8*B + 24*C))/6;
if (x < 1) return ((12 - 9*B - 6*C)*xx*x +
(-18 + 12*B + 6*C) * xx + (6 - 2*B))/6;
}
구간 $[a=t_0, ...t_{n-1}=b]$에서 일정한 간격(꼭 일정한 간격일 필요는 없지만 여기서는 1로 잡음)으로 샘플링된 데이터 $\{ f_0, ...f_{n-1} \}$이 있다. $n-1$ 개의 각 구간에서 데이터를 보간하는 삼차다항식 함수 모임 $$S(t)= \{S_j (t)| t_j \le t <t_{j+1} , ~~j=0,1,2...,n-2 \} $$을 찾아보자. 전체 구간에서 연속적이고 충분히 부드럽게 연결되기 위해서는 우선 각 node에서 연속이어야 하고, 또한 1차 도함수와 2차 도함수도 연속이 되어야 한다. 물론 각 node에서는 샘플링된 데이터값을 가져야 한다.
$n-1$개 구간에서 각각 정의된 3차 함수를 결정하려면 $4 \times (n-1)$개의 조건이 필요하다. (a)에서 $n$개, (b), (c), (d)에서 $3\times (n-2)$개의 조건이 나오므로 총 $n+3\times(n-2)=4n-6$개의 조건만 생긴다. 삼차식을 완전히 결정하기 위해서는 2개의 추가적인 조건을 부여해야 하는데 보통 양끝에서 2차 도함수값을 0으로 하거나(natural boundary) 또는 양끝에서 도함수 값을 특정한 값으로 고정시킨다(clamped boundary). 여기서는 양끝의 2차 도함수를 0으로 한 natural cubic spline만 고려한다. 그리고 $S_j(t)$가 $n-2$개의 구간에 대해서만 정의되어 있지만, 끝점을 포하는 구간에서도 정의하자. 이 경우 $S_{n-1}(t), t\ge t_{n-1}$에 부여된 조건은 $t_{n-1}$에서 $S_{n-2}(t)$와 연속, 미분연속 그리고 2차 도함수가 0인 조건만 부여된다.
타원은 베지어 곡선을 이용해서 근사적으로 표현할 수 있음을 잘 알고 있다(https://kipl.tistory.com/313). 표준형 타원의 경우 1사분면에서 표현만 알면 대칭성에 의해서 나머지 사분면에서는 바로 구할 수 있으므로 거리를 구하는 점이 1사분면에 있는 경우만 알아보면 충분하다. 장축 반지름이 $a$, 단축 반지름이 $b$인 타원은 다음과 같이 베이지 곡선으로 표현된다.
타원 외부의 한 지점에서 타원까지의 최단거리를 구하자. 타원 피팅 문제에서 피팅 에러를 예측하는데 유용하게 사용할 수 있다. 일반적인 타원은 타원의 중심과 회전각을 구해서 원점으로 평행이동 후 다시 반대로 회전시켜서 표준형으로 만든 후 구하면 되기 때문에 표준형에 대해서만 고려해도 충분하다.(외부점도 동일하게 평행이동 시키고 회전시켜야 한다). 또한 표준형 타원의 경우도 대칭성에 의해서 외부점이 1 사분면에 있는 경우만 살펴보면 충분하다.
외부점 $(p, q)$에 타원 상의 한 점 $(x, y)$까지의 거리가 최단 거리가 되기 위해서는 두 점을 잇는 선분과 타원에서 접선이 서로 직교해야 한다는 사실에서 $(x,y)$ 위치를 정할 수 있다.
우선 타원 상의 한 점 $(x,y)$에서 접선의 기울기는
$$ \frac{dy}{dx} = - \frac{b^2 x}{a^2 y}$$
이고, $(p,q)$까지 선분의 기울기가 $$ m =\frac{ y - q}{x-p}$$이다. 서로 직교하기 위해서는 이 둘의 곱이 $-1$이어야 한다.
접점 $(x, y)$는 타원 위에 있다는 점을 고려하면 $x$나 $y$에 대한 4차 방정식을 풀면 답을 analytic하게 쓸 수 있다. 하지만 직접적으로 이용하기에는 너무 복잡한 형태이므로 대부분의 경우는 Newton-Raphson 방법과 같은 반복적인 수치해석 기법을 이용해서 답을 얻는다. 이를 위해 타원을
$$ x= a \cos t, \quad y = b \sin t$$
로 매개화하면 최단거리를 주는 점을 찾는 문제는 다음 $f(t)$의 근을 찾는 문제로 환원이 된다.
$$ f(t) = (a^2 -b^2) \cos t \sin t - ap \sin t + bq \cos t = 0, \quad 0 \le t \le\frac{\pi}{2}$$
$z= \tan(t/2)$로 치환하면 이 식의 근을 구하는 문제는 $z$에 대한 4차 방정식의 근을 구하는 문제와 동일함을 보일 수 있다.
/*(a, b)=radii of ellipse; 타원에 포함되는 점에 대해서는 초기조건 설정에 주의;*/
double dist2Ellipse(double p, double q, double a, double b, double& xt, double& yt) {
if (a == b) return dist2Circle(p, q, a, xt, yt); //원의 경우 별도 처리;
if (x == 0 && y == 0) //원점의 경우 별도처리;
if (a > b) { xt = 0, yt = b; return b;}
else { xt = a, yt = 0; return a;}
double x = fabs(p), y = fabs(q);
const double eps = 0.01 / max(a, b);
const int maxiters = 10;
double t;
if ((x*x/a/a + y*y/b/b) < 1) {// 타원 내부의 경우 초기조건의 조정이 필요함;
double t1 = atan2(y/b, x/a), t2 = atan(1.);
if (a > b) t = max(t1, t2);
else t = min(t1, t2);
} else t = atan2(y/b, x/a);
for (int iter = 0; iter < maxiters; iter++) {
double c = cos(t), s = sin(t);
double f = (a*a - b*b)*c*s - a*x*s + b*y*c;
if (f == 0) break;
double Df = (a*a - b*b)*(c*c - s*s) - a*x*c - b*y*s;
double dt = -f/Df;
t += dt;
t = max(0, min(2*atan(1.), t)); //[0:pi/2];
if (abs(dt) < eps ) break;
}
xt = a * cos(t); yt = b * sin(t);
xt = p >= 0? xt: -xt;
yt = q >= 0? yt: -yt;
return hypot(p-xt, q-yt);
}
위 알고리즘은 Newton-Raphson iteration이 수렴하지 않을 수도 있다. 좀 더 좋은 방법은 외부점에서 타원에 내린 수선의 발에서 접촉원(osculating circle)을 구하면 수선의 연장선은 접촉원의 중심을 지난다.
만약 타원의 한 지점에서 접촉원 중심과 외부점을 연결하는 직선이 접선에 수직이 아닌 경우에는 그 지점은 최단점이 아니다. 이때는 이 직선과 타원의 교점까지의 거리가 좀 더 짧아진다. 다시 이 교점에서 접촉원을 구한 후 외부점과 중심을 연결하는 선분이 교차하는 새로운 교점을 찾는 과정을 반복하면 점점 교점은 최단점에 접근함을 알 수 있다. 이를 이용하면 좀 더 효율적인 점과 타원사이의 최단거리를 구하는 알고리즘을 만들 수 있다. 타원 위의 한 점 $\vec{e}$에서 접촉원의 중심을 구하면 중심 $\vec{c}$와 외부점$\vec{p}$을 연결하는 직선과 타원의 교점은 대략적으로 $\vec c$와 $\vec p$를 $|\vec{r}=\vec{e}-\vec{c}|: |\vec{q}=\vec{p}-\vec{c}| - |\vec{e}-\vec{c}|=r:q-r$로 내부하는 지점이다. 따라서 새로운 교점의 위치는