주어진 데이터를 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$이다. 이제 

$$ \kappa [f] = \int_a^b (h''(x) - S''(x))^2 = \kappa [h]  + \kappa [S]  - 2 \int_a^b h''(x) S''(s) dx $$

그런데, cubic spline의 경우 각 node에서 2차 미분이 연속이므로  

\begin{align} \int_a^b h''(x) S''(x) dx &= h'(x) S''(x)\Big|_a^b - \int _a^b h'(x) S'''(x)dx \\ &=h'(b)S''(b)-h'(a)S''(b) - \sum_i \int_{x_i}^{x_{i+1}} h'(x) S_i '''(x) dx \\  &=h'(b)S''(b)-h'(a)S''(a)- \sum_i \int_{x_i}^{x_{i+1} } h'(x) (6 a_i) dx \\ &= h'(b)S''(b)-h'(a)S''(a)- 6\sum_i a_i \left[  h(x_{i+1})- h(x_i) \right] \\&= h'(b)S''(b)-h'(a)S''(a) \end{align}  $$ \therefore  \quad \kappa [f] = \kappa[h]+\kappa[S]+ h'(b)S''(b)-h'(a) S''(a)$$

임의의 곡선 $h$에 대해서 $\kappa[h]\ge 0$이므로, 양끝에서 2차 미분이 $S''(a)=S''(b)=0$으로 고정된 natural cubic spline은  

$$ \kappa[f] \ge \kappa [S]$$

이므로 주어진 샘플링 데이터를 곡률척도가 가장 작게 즉, 가장 부드럽게 보간하는 곡선임을 알 수 있다.

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

'Computational Geometry' 카테고리의 다른 글

Natural Cubic Spline  (0) 2024.03.11
Approximate Distance Between Ellipse and Point  (0) 2024.03.08
Distance from a Point to an Ellipse  (0) 2024.03.06
Data Fitting with B-Spline Curves  (0) 2021.04.30
Closest Pair of Points  (0) 2021.04.27
Posted by helloktk
,
 

구간 $[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에서는 샘플링된 데이터값을 가져야 한다. 

\begin{align}     (a) ~~& S(t_j) = f_j \\ (b)~~& S_{j+1}(t_{j+1}) = S_j( t_{j+1}) \\ (c)~~& S'_{j+1}(t_{j+1}) = S'_j(t_{j+1}) \\ (d)~~& S''_{j+1}(t_{j+1})  = S''_{j}(t_{j+1}) \end{align}

$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인 조건만 부여된다.

 

$j-$번째 구간의 삼차함수를 

$$S_j(t) = a_j + b_j (t - t_j) + c_j (t-t_j)^2 + d_j (t - t_j)^3$$

로 쓰면 (a)에서 

$$ S_j (t_j) = a_j = f_i,~~ j=0,1,..n-2$$

(b)에서 

$$ a_{j+1} =a_j+  b_j + c_j +d_j,~~j=0,1,...,n-2$$

(c)에서 

$$b_{j+1} = b_j + 2c_j+ 3d_j,~~j=0,1,...,n-2 $$

(d)에서 

$$ c_{j+1} = c_j+ 3d_j,~~j=0,1,...,n-2$$

이므로 (b)와 (c)에서 $d_j$을 소거하면

$$  c_j = 3(a_{j+1}-a_j) -b_{j+1} - 2b_j$$

그리고 (a)에서 $$d_j = -2(a_{j+1}-a_j) + b_j  + b_{j+1}$$ 이므로 $b_j$에 대해 정리하면 다음과 같은 점화식을 얻을 수 있다.

$$ b_{j+1} + 4b_j + b_{j-1}= 3(a_{j+1}- a_{j-1})=3(f_{j+1}-f_{j}),~~j=1,2,...,n-2$$

물론 $j=0$일 때는 (note, $S''_0(t_0) = 0 \to c_0=0$) $a_1 =a_0+b_0+d_0, b_1 = b_0 + 3d_0$이므로

$$ b_1 + 2b_0 = 3(f_1 - f_0) $$

$j=n-1$일 때 계수는 $S_{n-1}(t_{n-1}) = f_{n-1}$이고 $S''_{n-1}(t_{n-1} )=c_{n-1}=0$이므로 $$ 2b_{n-1} + b_{n-2} = 3(f_{n-1}-f_{n-2})$$

임을 알 수 있다. 따라서 계수를 구하는 과정은 $b_j~(j=0,...,n-1)$을 구하는 것으로 결정된다. 이를 행렬로 표현하면

$$ \begin{bmatrix} 2&1& \\1&4&1\\ &1&4&1 \\ & & & \cdots \\ &&&&1&4&1  \\ && & & & 1 &2 \end{bmatrix} \begin{bmatrix} b_0\\b_1\\ \vdots \\   b_{n-2} \\b_{n-1}\end{bmatrix}=\begin{bmatrix} 3(f_1-f_0) \\ 3(f_2-f_0) \\  \vdots \\  3(f_{n-1}- f_{n-3}) \\3(f_{n-1} - f_{n-2}) \end{bmatrix}$$

와 같다. band 행렬은 upper triangle로 변환한 후 역치환과정을 거치면 쉽게 해를 구할 수 있다.

 

평면에서 주어진 점들을 보간하는 곡선은 이들 점을 표현하는 곡선의 매개변수를 일정한 간격으로 나누어서 샘플링된 결과로 취급하면, $x, y$ 성분에 대해서 각각 natural cubic spline를 구하여 얻을 수 있다. 

struct Cubic {
    double a,b,c,d;  /* a + b*t + c*t^2 +d*t^3 */
    Cubic(double a_, double b_, double c_, double d_) 
        : a(a_), b(b_), c(c_), d(d_) { }
    /* evaluate;*/
    double eval(double t) {
        return (((d*t) + c)*t + b)*t + a;
    }
};
void calcNaturalCubic(std::vector<double>& x, std::vector<Cubic>& Spline) {
    std::vector<double> gamma(x.size());
    std::vector<double> delta(x.size());
    std::vector<double> D(x.size());
    /* solve the banded equation:
    [2 1       ] [  D[0]]   [3(x[1]   - x[0])  ]
    |1 4 1     | |  D[1]|   |3(x[2]   - x[0])  |
    |  1 4 1   | | .    | = |      .           |
    |    ..... | | .    |   |      .           |
    |     1 4 1| | .    |   |3(x[N-1] - x[N-3])|
    [       1 2] [D[N-1]]   [3(x[N-1] - x[N-2])]
    ** make the banded matrix to an upper triangle;
    ** and then back sustitution. D[i] are the derivatives at the nodes.
    */
    const int n = x.size() - 1;  // note n != x.size()=N;
    // gamma;
    gamma[0] = 0.5;
    for (int i = 1; i < n; i++)
        gamma[i] = 1/(4-gamma[i-1]);
    gamma[n] = 1/(2-gamma[n-1]);
    // delta;
    delta[0] = 3*(x[1]-x[0])*gamma[0];
    for (int i = 1; i < n; i++) 
        delta[i] = (3*(x[i+1]-x[i-1])-delta[i-1])*gamma[i];
    delta[n] = (3*(x[n]-x[n-1])-delta[n-1])*gamma[n];
    // D;
    D[n] = delta[n];
    for (int i = n; i-->0;)
        D[i] = delta[i] - gamma[i]*D[i+1];

    /* compute the coefficients;*/
    Spline.clear(); Spline.reserve(n);
    for (int i = 0; i < n; i++)
        Spline.push_back(Cubic(x[i], D[i], 3*(x[i+1]-x[i])-2*D[i]-D[i+1],
            2*(x[i]-x[i+1]) + D[i] + D[i+1])) ;
}
void NaturalCubicSpline(std::vector<CPoint>& points, 
                        std::vector<CPoint>& curve) {
    curve.clear();
    if (points.size() < 2) return;
    std::vector<double> xp(points.size()), yp(points.size());
    for (int i = points.size(); i-->0;)
        xp[i] = points[i].x, yp[i] = points[i].y;

    std::vector<Cubic> splineX, splineY;
    calcNaturalCubic(xp, splineX);
    calcNaturalCubic(yp, splineY);
#define STEPS 12
    curve.reserve(splineX.size() * STEPS + 1);
    curve.push_back(CPoint(int(splineX[0].eval(0) + 0.5), int(splineY[0].eval(0) + 0.5)));
    for (int i = 0; i < splineX.size(); i++) {
        for (int j = 1; j <= STEPS; j++) {
            double t = double(j) / STEPS;
            curve.push_back(CPoint(int(splineX[i].eval(t) + 0.5), int(splineY[i].eval(t) + 0.5)));
        }
    }
}

 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

'Computational Geometry' 카테고리의 다른 글

Why Cubic Splines?  (2) 2024.03.16
Approximate Distance Between Ellipse and Point  (0) 2024.03.08
Distance from a Point to an Ellipse  (0) 2024.03.06
Data Fitting with B-Spline Curves  (0) 2021.04.30
Closest Pair of Points  (0) 2021.04.27
Posted by helloktk
,
 

타원은 베지어 곡선을 이용해서 근사적으로 표현할 수 있음을 잘 알고 있다(https://kipl.tistory.com/313). 표준형 타원의 경우 1사분면에서 표현만 알면 대칭성에 의해서 나머지 사분면에서는 바로 구할 수 있으므로 거리를 구하는 점이 1사분면에 있는 경우만 알아보면 충분하다. 장축 반지름이 $a$, 단축 반지름이 $b$인 타원은 다음과 같이 베이지 곡선으로 표현된다.

 

Bezier Curve Approximation of an Ellipse

Bezier 곡선을 이용한 원의 근사처럼 타원을 근사해보도록 하자. 원점을 중심으로 하고 장축 반지름이 $a$, 단축 반지름이 $b$인 타원을 1사분에서 3차 Bezier curve을 이용해서 근사하려면 4개의 control

kipl.tistory.com

$$\mathbf {B}(t) = (1-t^3) \mathbf {P}_0 + 3t(1-t)^2 \mathbf {P}_1 + 3t^2 (1-t) \mathbf {P}_2 + t^3 \mathbf {P}_3 = \left(\begin {array}{c} a x(t) \\ b y(t) \end {array}\right) $$

$$ x(t) = 3k (1-t)^2 t + 3 (1-t) t^2 + t^3 \\ y(t) = x(1-t) ,  \quad 0 \le t \le 1$$

$k$ 값은 $t=1/2$일 때 $(a/\sqrt {2}, b/\sqrt {2})$을 통과하는 조건을 부여하여 정하면

$$ k = \frac {4}{3}(\sqrt {2}-1)= 0.5522847498...$$임을 알 수 있다. 

 

한 점 ${\bf P}=(p,q)$에서 베지어 곡선 위의 한 점 ${\bf B}(t)$까지의 거리가 최단거리가 되기 위해서는 다음 직교 조건

$$f(t)= ({\bf B}(t) - {\bf P})\cdot {\bf B}'(t) = 0 \\ f(t)=a^2x(t)x'(t)-b^2y(t)y'(t)-ap x'(t)-bqy'(t)=0$$

을 만족하는 근을 찾으면 된다. 풀어야 할 방정식이 5차이므로 직접적으로 근을 찾는 것은 불가능하므로 Newton-Raphson 방법을 사용하도록하자. 앞선 정확한 거리 계산에서와는 달리 초기 $t$ 값 설정에 민감하게 의존하지 않는다.

//cubic bezier_x: P0(0, 1), P1(k, 1),P2(1, k), P3(1, 0);
double B(double t) { //t in [0:1]
    const double k = 4.*(sqrt(2.)-1)/3;
    return t*(3*k + t*(3 - 6*k + (-2 + 3*k)*t));
}
// derivative of B(t);
double DB(double t) {
    const double k = 4.*(sqrt(2.)-1)/3;
    return 3*k + t*(6 - 12*k + (-6 + 9*k)*t);
}
// derivative of DB(t);
double D2B(double t) {
    const double k = 4.*(sqrt(2.)-1)/3;
    return 6 - 12*k + (-12 + 18*k)*t;
}
// ellipse radii=(a, b);
double dist2EllipseBezier3(double p, double q, double a, double b,
                           double& xt, double& yt) {
    if (a == b) return dist2Circle(p, q, a, xt, yt);
    double x = fabs(p), y = fabs(q);
    const double eps = 0.001 / max(a, b);
    double t = 0.5;  // mid
    while (1) {
        // Newton-Raphson;
        double f = a*a*B(t)*DB(t)-b*b*B(1-t)*DB(1-t)-a*x*DB(t)+b*y*DB(1-t);
        if (f == 0) break;
        double df = a*a*(DB(t)*DB(t)+B(t)*D2B(t))+b*b*(DB(1-t)*DB(1-t)+B(1-t)*D2B(1-t))
                    -a*x*D2B(t)-b*y*D2B(1-t);
        double dt = f / df;
        t -= dt;
        t = max(0, min(1, t));
        if (abs(dt) < eps ) break;
    }
 
    xt = a * B(t); yt = b * B(1-t);
    xt = p >= 0? xt: -xt;
    yt = q >= 0? yt: -yt;
    return hypot(p - xt, q - yt);
}
 
 
 
 
 
 
 
 
 

이심률이 큰 경우와 작은 경우

 

 
 
 
 
 
 
 
 
 
728x90

'Computational Geometry' 카테고리의 다른 글

Why Cubic Splines?  (2) 2024.03.16
Natural Cubic Spline  (0) 2024.03.11
Distance from a Point to an Ellipse  (0) 2024.03.06
Data Fitting with B-Spline Curves  (0) 2021.04.30
Closest Pair of Points  (0) 2021.04.27
Posted by helloktk
,

 

 

타원  외부의 한 지점에서 타원까지의 최단거리를 구하자. 타원 피팅 문제에서 피팅 에러를 예측하는데 유용하게 사용할 수 있다. 일반적인 타원은 타원의 중심과 회전각을 구해서 원점으로 평행이동 후 다시 반대로 회전시켜서 표준형으로 만든 후 구하면 되기 때문에 표준형에 대해서만 고려해도 충분하다.(외부점도 동일하게 평행이동 시키고 회전시켜야 한다). 또한 표준형 타원의 경우도  대칭성에 의해서 외부점이 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$이어야 한다.

$$ \left(-\frac{b^2 x}{a^2 y} \right)  \frac{y-q}{x-p} = -1 \quad \to\quad (a^2 -b^2 )xy + b^2 qx - a^2 y p = 0$$

접점 $(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{c} =\left( \frac{ a^2 - b^2 }{a} \cos ^3 t , \frac{b^2 - a^2}{b} \sin^3 t\right)  $$

만약 타원의 한 지점에서 접촉원 중심과 외부점을 연결하는 직선이 접선에 수직이 아닌 경우에는 그 지점은 최단점이 아니다. 이때는 이 직선과 타원의 교점까지의 거리가 좀 더 짧아진다. 다시 이 교점에서 접촉원을 구한 후 외부점과 중심을 연결하는 선분이 교차하는 새로운 교점을 찾는 과정을 반복하면 점점 교점은 최단점에 접근함을 알 수 있다. 이를 이용하면 좀 더 효율적인 점과 타원사이의 최단거리를 구하는 알고리즘을 만들 수 있다. 타원 위의  한 점 $\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$로 내부하는 지점이다. 따라서 새로운 교점의 위치는

$$ \vec{r}' \approx \vec{c} + \frac{r}{q} \vec{q} $$

이다. 물론 이 지점이 정확히 타원 위에 있지는 않다. 타원 위에 있는 점은 각 축의 반지름으로 각성분을 나누면 단위원 위에 있어야 한다는 사실을 이용하면 $\vec{r}'$에 가장 가까이 있는 타원 상의 한 지점을 얻을 수 있다.  즉, 

$$  \cos t = \frac{x'/a}{\sqrt{(x'/a)^2+(y'/b)^2}},\quad    \sin t = \frac{y'/b}{\sqrt{(x'/a)^2 + (y'/b)^2 }}$$

인 $t$을 구하면 새로운 타원 위의 지점은

$$ x''= a \cos (t), \quad y'' = b \sin t$$

로 된다.

double dist2Ellipse2(double p, double q, double a, double b, double& xt, double& yt) {
    if (a==b) return dist2Circle(p, q, a, xt, yt);//원의 경우에는 별도로;
    const double px = abs(p), double py = abs(q);
    const double eps = 0.01/max(a,b);
    const double cc = a*a - b*b;
    double cost, sint; 
    //근이 여러개 존재하므로 초기값 설정에 주의해야 한다;
    if ((px*px/a/a + py*py/b/b) < 1) //타원에 포함되는 경우;
    	cost = sint = 1/sqrt(2.); //  pi/4에서 출발;
    else { //타원바깥에 있는 경우;
        cost = cos(atan2(py/b, px/a));
    	sint = sin(atan2(py/b, px/a));
    }
    double t = 1, told = t;
    while (1) {
        // 곡률중심;evoute
        double cx = cc * pow(cost, 3.) / a;
        double cy = - cc * pow(sint, 3.) / b;
        // 곡률중심에 대한 타원 위의 점(교점)
        double x = a*cost - cx, y = b*sint - cy;
        // 곡률중심에 대한 외부점 변위
        double qx = px - cx, qy = py - cy;
        //거리;
        double r = hypot(x, y), q = hypot(qx, qy);
        cost = min(1, max(0, (qx * r / q + cx) / a));
        sint = min(1, max(0, (qy * r / q + cy) / b));
        double t = hypot(cost, sint); //1이 되면 좋지만 안되는 경우;
        cost /= t; sint /= t;
        if (abs(t-told) < eps) break;
        told = t;
    }
    xt = abs(a *cost), yt = abs(b *sint);
    xt = p >= 0 ? xt: -xt;
    yt = q >= 0 ? yt: -yt;
    return hypot(p - xt, q - yt);
}

이심률이 큰 경우

728x90

'Computational Geometry' 카테고리의 다른 글

Natural Cubic Spline  (0) 2024.03.11
Approximate Distance Between Ellipse and Point  (0) 2024.03.08
Data Fitting with B-Spline Curves  (0) 2021.04.30
Closest Pair of Points  (0) 2021.04.27
DDA Algorithm  (0) 2021.04.25
Posted by helloktk
,

$m$개의 control point $\{\mathbf{Q}_i\}$가 주어지면  $p$차의 B-Spline curve는 basis 함수를 $N_{i, p}(t)$을 써서

$$ \mathbf{B}(t) = \sum_{i=0}^{m-1} N_{i, p}(t) \mathbf{Q}_i $$

로 표현할 수 있다. 이를 이용해서 일정한 순서로 샘플링된 평면상의 $N$개의 입력점 $\{ \mathbf{P}_i \}$을 찾아보자. B-spline 곡선을 이용하면 이 문제는 control 점을 찾는 문제가 된다. 곡선이 입력점을 잘 표현하기 위해서는 곡선과 입력점과의 차이를 최소로 하는 control 점을 찾아야 한다:

$$    \mathbf{ Q}^* = \text{argmin}(L), \quad\quad L:= \sum_{k = 0}^{N-1} | \mathbf{B}(t_k) - \mathbf{P}_k|^2 $$

여기서 $\{ t_k| k=0,1,...,N-1\}$는 입력점이 얻어지는 sample time으로 $0= t_0\le t_1\le...\le t_{N-1}= 1$로 rescale 할 수 있다. 

행렬을 이용해서 식을 좀 더 간결하게 표현할 수 있다. 

$$L = \sum_{k = 0}^{N-1} | \mathbf{B}(t_k) - \mathbf{P}_k|^2 = \sum_{k = 0}^{N-1} \left| \sum_{i=0}^{m-1} N_{i, p}(t_k) \mathbf{Q}_i -  \mathbf{P}_k \right|^2 = \sum_{k = 0}^{N-1} \left| \sum_{i=0}^{m-1} A_{ki} \mathbf{Q}_i - \mathbf{P}_k \right|^2, \\ \quad A_{ki} = N_{i, p}(t_k) $$

로 쓸 수 있으므로, $\hat {Q}= (\mathbf{Q}_0, \mathbf{Q}_1,..., \mathbf{Q}_{m-1})^t$, $\hat {P} = (\mathbf{P}_0, \mathbf{P}_1,..., \mathbf{P}_{N-1})^t$인 벡터, 그리고 $ \hat {A} = (A_{ki})$인 행렬로 보면

$$ L= \left| \hat{A} \cdot \hat {Q} - \hat {P} \right|^2 =  (\hat{Q}^t \cdot \hat {A}^t -\hat{P}^t ) \cdot (\hat{A}\cdot \hat{Q} - \hat{P}).$$

위 식의 값을 최소로 하는 최소자승해는 $\hat {Q}^t$에 대한 미분 값이 0인 벡터를 찾으면 된다; $$ \frac {\partial L}{\partial \hat {Q}^t } = \hat {A}^t \cdot \hat {A} \cdot \hat {Q} - \hat {A}^{t} \cdot \hat {P} = 0.$$ 이 행렬 방정식의 해는

$$ \hat {Q}^* = ( \hat {A} ^t \cdot \hat {A})^{-1} \cdot ( \hat {A}^t \cdot \hat {P})$$ 로 표현된다. $\hat {A} ^t \cdot \hat {A}$가 (banded) real symmetric ($m\times m$) 행렬이므로 Cholesky decomposion을 사용하면 쉽게 해를 구할 수 있다. $\hat{A}$가 banded matrix 형태를 가지는 이유는 basis가 local support에서만 0이 아닌 값을 취하기 때문이다. 

open b-spline(cubic)

int BSplineFit_LS(std::vector<CPoint>& data,
                  int degree,             // cubic(3); 
                  int nc,                 // num of control points;
                  std::vector<double>& X,
                  std::vector<double>& Y) // estimated control points;
{
    // open b-spline;
    std::vector<double> knot((nc - 1) + degree + 2);
    for (int i = 0; i <= nc + degree; i++) knot[i] = i;
    
    std::vector<double> t(data.size());       // parameter;
    double scale = (knot[nc] - knot[degree]) / (data.size()-1);
    for (int i = data.size(); i-->0;) 
        t[i] = knot[degree] + scale * i;

    std::vector<double> A(ndata * nc);
    for (int i = data.size(); i-->0;)
        for (int j = 0; j < nc; j++)
            A[i * nc + j] = Basis(j, degree, &knot[0], t[i]); //A(i,j)=N_j(t_i)

    // S = A^t * A; real-symmetric matrix;
    std::vector<double> Sx(nc * nc);
    std::vector<double> Sy(nc * nc);
    for (int i = 0; i < nc; i++) {
        for (int j = 0; j < nc; j++) {
            double s = 0;
            for (int k = data.size(); k-->0;)
                s += A[k * nc + i] * A[k * nc + j];
            Sx[i * nc + j] = s;
        }
    }
    //copy;
    for (int i = nc*nc; i-->0;) Sy[i] = Sx[i];
    // X = A^t * P.x;  Y = A^t * P.y
    for (int i = 0; i < nc; i++) {
        double sx = 0, sy = 0;
        for (int k = data.size(); k-->0;) {
            sx += A[k * nc + i] * data[k].x;
            sy += A[k * nc + i] * data[k].y;
        };
        X[i] = sx; Y[i] = sy;
    };
    // solve real symmetric linear system; S * x = X, S * y = Y;
    // solvps(S, X) destories the inputs;
    // ccmath-2.2.1 version;
    int res1 = solvps(&Sx[0], &X[0], nc);
    int res2 = solvps(&Sy[0], &Y[0], nc);
    return res1 == 0 && res2 == 0;
};

**네이버 블로그 이전;

 
728x90

'Computational Geometry' 카테고리의 다른 글

Approximate Distance Between Ellipse and Point  (0) 2024.03.08
Distance from a Point to an Ellipse  (0) 2024.03.06
Closest Pair of Points  (0) 2021.04.27
DDA Algorithm  (0) 2021.04.25
B-Spline  (0) 2021.04.25
Posted by helloktk
,