Fourier Interpolation

Mathematics 2024. 3. 20. 21:37

$0\le x <L$ 구간에서 일정한 간격($\Delta = L/N$)으로 샘플링된 데이터 $\{y_k|k=0,...,N-1\}$이 주어졌을 때 이들의  DFT 계수 $\{Y_k | k=0, N-1\}$을 이용하면 보간함수를 만들 수 있다.

$$ y(x) = Y_0 + \sum_{0 <k <N/2}  \left( Y_k e^{+\frac {2\pi i}{L}  k x} + Y_{N-k} e^{ -\frac {2\pi i}{L}  k x} \right)  +  Y_{N/2} \cos\Big( \frac{\pi}{L}Nx\Big) $$

마지막 항(Nyquist term)은 $N$이 짝수일 때만 더해진다. 구체적인 유도과정은 https://kipl.tistory.com/406에 있다. 샘플링 데이터가 실수인 경우는 $Y_{N-k} = (Y_k)^*$가 성립하므로 표현이 더 간단해진다.

 

FFT를 이용한 영상의 미분

주어진 함수 $y(x)$를 구간 $0\le x < L$에서 일정한 간격으로 샘플링해서 얻은 $N$개의 data $\{y_n = y(nL/N)|n=0, 1,...,N-1\} $의 discrete Fourier transform(DFT) $Y_k$는 $$ Y_ k = \frac {1}{N} \sum _{n = 0}^{N-1} y_n e^{ -\frac {2\p

kipl.tistory.com

 

$x/L\to x$으로 하면 이 보간함수의 domain은 $[0,1]$로 제한된다.

// [x=0:1]
// Fourier coefficient (re[], im[]);
double FourierInterp(double x, std::vector<double>& re, std::vector<double>& im) {
    const int N = re.size();
    double fval = re[0];
    for (int k = N >> 1; k-->1;) {
        double ang = TWOPI * k * x;
        double cc = cos(ang);
        double ss = sin(ang);
        fval += re[k] * cc - im[k] * ss;     // 입력데이터가 실수인 경우:
        fval += re[N-k] * cc + im[N-k] * ss; // Y[N-k] = (Y[k])*
    }
    if (!(N&1))           //k = N/2 for even N;
        fval += re[N>>1] * cos(PI * N * x);
    return fval;
}

$f(x) = \sin^2(x)\cos(x)$을 $[0,3\pi]$ 구간에서 일정한 간격으로 얻은 30개의 데이터를 이용해서 보간함수를 만든다.

double f(double x) {
    return sin(x) * sin(x) * cos(x);
}
void Finterpolate() {
    const int nsamples = 30;
    std::vector<double> sample(nsamples);
    const double L = 3 * PI, dx = L / nsamples;
    for (int k = 0; k < nsamples; k++) sample[k] = f(k * dx);

    std::vector<double> re(sample.size(), 0), im(sample.size(), 0);
    // perform DFT;
    DFT(sample, re, im);

    FILE *fp = fopen("finp.txt", "w");
    const double N = 200;
    const double dy = dx * sample.size() / N;
    for (int k = 0; k < N; k++) { 
        double x = k * dy;
        fprintf(fp, "%f %f\n", x, FourierInterp(x / L, re, im));
    }
    fclose(fp);
}

 
 
 
 
 
 
 
 
 
 
 
 
728x90

'Mathematics' 카테고리의 다른 글

Generalized eigenvalues problem  (0) 2024.03.02
수치적으로 보다 정밀한 이차방정식의 해  (0) 2024.02.23
열방정식의 Green function  (0) 2024.02.12
지구의 나이는?  (0) 2024.02.11
n 차원 공의 부피  (2) 2024.02.07
Posted by helloktk
,

주어진 데이터를 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' 카테고리의 다른 글

Derivatives of Bezier Curves  (1) 2024.04.12
Least Squares Bezier Fit  (0) 2024.04.05
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
Posted by helloktk
,

일정한 간격 $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차 도함수가 연속적임을 요구하면

\begin{align} s = 1^\pm:~~~& 3A_1 + 2B_1 = 3A_2 + 2B_2+C_2 \end{align}

그리고 샘플링된 데이터가 모두 같은 경우 보간함수도 상수함수가 되는 것이 타당하므로

$$ g(x) = \sum_k K \left( \frac{x-x_k}{h}\right) = 1~~~\text{if} ~~\forall f_j = 1$$

을 만족시켜야 한다. $x_j <x<x_{j+1}$일 때 $x  =  x_j + sh,  ~(0< s<1)$로 쓸 수 있고, kernel이 반지름이 2인  support를 가지므로 

$$ g(x) =  K(s+1) + K(s) +  K(s-1) + K(s-2)=1$$

임을 알 수 있다. 위에서 주어진 $K(s)$을 대입해서 정리하면 다음과 같은 항등식을 얻는다.

$$ -1 + A_1 + 9 A_2 + B_1 + 5 B_2 + 3 C_2 + 2 D_1 +2 D_2 \\ +(-3 A_1 - 9 A_2 - 2 B_1 - 2 B_2) s + (3 A_1 + 9 A_2 + 2 B_1 + 2 B_2) s^2 =0$$

이 항등식의 계수가 0이 되어야 한다는 사실에서 2 개의 추가 조건을 얻으므로 총 8개 계수 중 2개가 미결정 free parameter로 남는다. 보통 이 두 계수는 $D_1=1-B/3$,  $D_2=4C + 4B/3$처럼 매개화한다. 이 경우 kernel 함수는

$$ K(s)  = \frac{1}{6} \left\{ \begin{matrix}   (12-9B-6C)|s|^3 +(-18+12B+6C) |s|^2 + (6-2B) & |s|< 1\\ (-B-6C)|s|^3 +(6B+30C)|s|^2 + (-12B-48C) |s| + (8B+24C) & 1\le |s|<2\\0 & \text{otherwise} \end{matrix} \right.$$

따라서 cubic spline kernel은 두 개의 파라미터 (B,C)에 의해서 정해진다. 또한 kernel 함수의 적분은 $B,C$에 상관없이 항상 1이어서 총 가중치의 합이 1임이 자동으로 보증된다.

$$\int_{-\infty}^\infty K(s)ds =1$$

이 중에는 이미지의 resampling에서 많이 사용되는 커널도 있는데, 잘 알려진 경우를 보면

$$ \begin{matrix} (B,C)=(0,1) & \text{Cardinal spline} \\ (B,C)=(0,1/2) & \text{Catmull-Rom spline } \\ (B,C)=(0,3/4) & \text{used in photoshop} \\ (B,C)=(1/3,1/3) & \text{ Mitchell-Netravali spline}  \\ (B,C)=(1,0) & \text{B-spline}\end{matrix}$$

$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;
}
 
 
 
 
 
 
 
 

 

 
 
 
 
 
 
 
 
 
 
728x90

'Image Recognition > Fundamental' 카테고리의 다른 글

Linear Least Square Fitting: perpendicular offsets  (0) 2024.03.22
Ellipse Fitting  (0) 2024.03.02
Bilateral Filter  (0) 2024.02.18
영상에서 Impulse Noise 넣기  (2) 2023.02.09
Canny Edge: Non-maximal suppression  (0) 2023.01.11
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' 카테고리의 다른 글

Least Squares Bezier Fit  (0) 2024.04.05
Why Cubic Splines?  (9) 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
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?  (9) 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
,