타원은 베지어 곡선을 이용해서 근사적으로 표현할 수 있음을 잘 알고 있다(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
,

 

 

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

일반적인 conic section 피팅은 주어진 데이터 $\{ (x_i, y_i)\}$를 가장 잘 기술하는 이차식

$$F(x, y) = ax^2 + bxy +cy^2 + dx +ey +f=0 $$

의 계수 ${\bf u^T}= (a,b,c,d,e,f)$을 찾는 문제이다. 이 conic section이 타원이기 위해서는 2차항의 계수 사이에 다음과 같은 조건을 만족해야 한다.

$$\text{ellipse constraint:}~~ ac - b^2/4 >0$$

그리고 얼마나 잘 피팅되었난가에 척도가 필요한데 여기서는 주어진 데이터의 대수적 거리 $F(x,y)$을 이용하자. 주어진 점이 타원 위의 점이면 이 값은 정확히 0이 된다. 물론 주어진 점에서 타원까지의 거리를 사용할 수도 있으나 이는 훨씬 복잡한 문제가 된다.  따라서 해결해야 하는 문제는

\begin{gather}L = \sum _{i}  \left( ax_i^2 + bx_i y_i + cy_i^2 +dx_i + e y_i +f\right)^2 - \lambda( 4ac-b^2-1) \\= \left|\begin{pmatrix}x_0^2& x_0y_0 & y_0^2 & x_0 & y_0 & 1\\ x_1^2 & x_1 y_1& y_1^2 & x_1 & y_1 & 1 \\ x_2^2 & x_2y_2& y_2^2 & x_2& y_2 & 1\\ &&\vdots \\\end{pmatrix}\begin{pmatrix}a\\b\\c\\d\\e\\f \end{pmatrix}  \right|^2 -\lambda \left({\bf  u^T} \begin{pmatrix} 0& 0& 2&0&0&0\\ 0 &-1&0 &0 &0 &0\\ 2&0&0&0&0&0\\0&0&0&0&0&0 \\0&0&0&0&0&0\\0&0&0&0&0&0&  \end{pmatrix} \bf u -1\right) \\ =\bf u^T D^TD u -\lambda (u^T C u -1)\\ = \bf u^T S u -\lambda (u^T C u-1)\end{gather}

을 최소화시키는 계수 벡터 $\bf u$를 찾는 것이다. 여기서 제한조건으로 $4ac - b^2 =1= \bf u^T C u$로 설정했다. 

$\bf u^T$에 대해서 미분을 하면 

$$ \frac{\partial L}{\partial \bf u^T} =  \bf S u -\lambda C u=0$$

즉, 주어진 제한조건 $4ac - b^2=1$하에서 대수적 거리를 최소화시키는 타원방정식의 계수 $\bf u$를 구하는 문제는 scattering matrix $\bf S=D^T D$에 대한 일반화된 고유값 문제로 환원이 된다.

$$  \bf S u =\lambda C u \\ u^T C u =1$$

이 문제의 풀이는 직전의 포스팅에서 다른 바 있는데 $\bf S$의 제곱근 행렬 $\bf Q=S^{1/2}$를 이용하면 된다. 주어진 고유값 $\lambda$와 고유벡터 $\bf u$가 구해지면 대수적 거리는 $$\bf u^T S u = \lambda$$

이므로 이를 최소화시키기 위해서는 양의 값을 갖는 고유값 중에 최소에 해당하는 고유벡터를 고르면 된다. 그런데 고유값 $\lambda$의 부호별 개수는 $\bf C$의 고유값 부호별 개수와 동일함을 보일 수 있는데 (Sylverster's law of inertia),  $\bf C$의 고유값이 $\{-2,-1,2,0,0,0\}$이므로 $\lambda>0$인 고유값은 1개 뿐임을 알 수 있다. 따라서 $\bf S u = \lambda C u$를 풀어서 얻은  유일한 양의 고유값에 해당하는 고유벡터가 원하는 답이 된다.

https://kipl.tistory.com/370

 

Least Squares Fitting of Ellipses

일반적인 이차곡선은 다음의 이차식으로 표현이 된다: $$ F(x, y)=ax^2 + bxy + cy^2 +d x + ey + f=0$$ 6개의 계수는 모두 독립적이지 않고 어떤 종류의 이차곡선인가에 따라 제약조건이 들어온다. 주어진

kipl.tistory.com

https://kipl.tistory.com/565

 

Generalized eigenvalues problem

$\bf S$가 positive definite 행렬이고, $\bf C$는 대칭행렬일 때 아래의 일반화된 eigenvalue 문제를 푸는 방법을 알아보자. $$\bf S u = \lambda C u$$ 타원을 피팅하는 문제에서 이런 형식의 고유값 문제에 부딛

kipl.tistory.com

 Ref: https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/ellipse-pami.pdf

 

 
double FitEllipse(std::vector<CPoint>& points, double einfo[6] ) {     
    if ( points.size() < 6 ) return -1;
    double eigvals[6];
    std::vector<double> D(6 * points.size());
    double S[36];/*  S = ~D * D  */
    double C[36];
    double EIGV[36];/* R^T; transposed orthogonal matrix;*/

    double offx = 0, offy = 0;
    /* shift all points to zero */
    for(int i = points.size(); i--> 0; ) {	
        offx += points[i].x;
        offy += points[i].y;        	
    }
    offx /= points.size(); 
    offy /= points.size();

    /* for the sake of numerical stability, scale down to [-1:1];*/
    double scale = 0; 
    for (int i = points.size(); i-->0; ) {
        if (points[i].x > scale) scale = points[i].x;
        if (points[i].y > scale) scale = points[i].y;
    }
    double invscale = 1 / scale;
    /* ax^2 + bxy + cy^2 + dx + ey + f = 0*/
    /* fill D matrix rows as (x*x, x*y, y*y, x, y, 1 ) */
    for(int i = points.size(); i--> 0; ) {	
        double x = points[i].x - offx; x *= invscale; 
        double y = points[i].y - offy; y *= invscale;
        D[i*6 + 0] = x*x; D[i*6 + 1] = x*y;
        D[i*6 + 2] = y*y; D[i*6 + 3] = x;
        D[i*6 + 4] = y;   D[i*6 + 5] = 1;		
    }			

    /* scattering matrix: S = ~D * D (6x6)*/
    for (int i = 0; i < 6; i++) 
        for (int j = i; j < 6; j++) { /*upper triangle;*/
            double s = 0;
            for (int k = points.size(); k-- > 0; ) 
                s += D[k*6 + i] * D[k*6 + j];
            S[i*6 + j] = s;
        }
    for (int i = 1; i < 6; i++) /*lower triangle;*/
        for (int j = 0; j < i; j++) 	
            S[i*6 + j] = S[j*6 + i] ;
    
    /* fill constraint matrix C */
    for (int i = 0; i < 36 ; i++ ) C[i] = 0;
    C[12] =  2 ;//2x0 
    C[2 ] =  2 ;//0x2 
    C[7 ] = -1 ;//1x1

    /* find eigenvalues/vectors of scattering matrix; */
    double RT[36];	/* each row contains eigenvector; */
    JacobiEigens ( S, RT, eigvals, 6, 0 );
    /* create R and INVQ;*/
    double R[36];
    for (int i = 0; i < 6 ; i++) {
        eigvals[i] = sqrt(eigvals[i]);
        for ( int k = 0; k < 6; k++ ) {
            R[k*6 + i] = RT[i*6 + k];  /* R = orthogonal mat = transpose(RT);*/
            RT[i*6 + k] /= eigvals[i]; /* RT /= sqrt(eigenvalue) row-wise)*/
        }
    }
    /* create INVQ=R*(1/sqrt(eigenval))*RT;*/
    double INVQ[36];
    _MatrixMul(R, RT, 6, INVQ);

    /* create matrix INVQ*C*INVQ */
    double TMP1[36], TMP2[36];
    _MatrixMul(INVQ, C, 6, TMP1 );
    _MatrixMul(TMP1, INVQ, 6, TMP2 );
    
    /* find eigenvalues and vectors of INVQ*C*INVQ:*/
    JacobiEigens ( TMP2, EIGV, eigvals, 6, 0 );
    /* eigvals stores eigenvalues in descending order of abs(eigvals);*/
    /* search for a unique positive eigenvalue;*/
    int index = -1, count = 0;
    for (int i = 0 ; i < 3; i++ ) {
        if (eigvals[i] > 0) {
            index = i; // break;
            count++;
        }
    }
    /* only 3 eigenvalues must be non-zero 
    ** and only one of them must be positive;*/
    if ((count != 1) || (index == -1)) 
        return -1;
     
    /* eigenvector what we want: u = INVQ * v */
    double u[6]; 
    double *vec = &EIGV[index*6];
    for (int i = 0; i < 6 ; i++) {
        double s = 0;
        for (int k = 0; k < 6; k++) s += INVQ[i*6 + k] * vec[k];
        u[i] = s;
    }
    /* extract shape infos;*/
    PoseEllipse(u, einfo);
    /* recover original scale; center(0,1) and radii(2,3)*/
    for (int i = 0; i < 4; i++) einfo[i] *= scale;
    /* recover center */
    einfo[0] += offx; 
    einfo[1] += offy;
    return FitError(points, offx, offy, scale, u);
};
 
728x90

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

Linear Least Square Fitting: perpendicular offsets  (0) 2024.03.22
Cubic Spline Kernel  (1) 2024.03.12
Bilateral Filter  (0) 2024.02.18
영상에서 Impulse Noise 넣기  (2) 2023.02.09
Canny Edge: Non-maximal suppression  (0) 2023.01.11
Posted by helloktk
,

$\bf S$가 positive definite 행렬이고, $\bf C$는 대칭행렬일 때 아래의 일반화된 eigenvalue 문제를 푸는 방법을 알아보자. 
$$\bf  S u = \lambda C u$$

타원을 피팅하는 문제에서 이런 형식의 고유값 문제에 부딛히게 된다. 이 경우 $\bf S$는 scattering matrix이고, $\bf C$는 타원피팅에 걸리는 제한조건때문에 나온다. $\bf S$가 positive definite이므로  $\bf S$의 eigenvalues $\{\sigma_i > 0 \}$는 모두 0보다 크고, eigenvector를 이용하면 

$$ \bf S = R \Lambda R^T, ~~~~\Lambda=\text{diag}(\sigma_1, ...,\sigma_n)$$

처럼 분해할 수 있다. $\bf R$은 $\bf S$의 eigenvector를 열로 가지는 행렬로 orthogonal 행렬이다: $\bf R^{-1}=R^T$.

 

이제 $\bf S$의 제곱근 행렬을 $\bf Q, ~Q^2 =S$라면 

$$ \bf Q= R \Lambda^{1/2} R^T,~~~~\Lambda^{1/2} = \text{diag}( \sqrt{\sigma_1},...,\sqrt{\sigma_n})$$

임을 쉽게 확인할 수 있다. $\bf Q$을 이용하면 구하려는 고유값 문제는 

$$ \bf QQ u = \lambda C ~\to ~  Qu = \lambda Q^{-1} C Q Q^{-1}~~\to ~~ v = \lambda Q^{-1} C Q^{-1} v$$

이므로 $\bf Q^{-1} C Q^{-1}$의 고유값 문제 $(1/\lambda, \bf Qu)$로 단순화됨을 알 수 있다. $\bf Q$의 역행렬이 $\bf Q^{-1} = R \Lambda^{-1/2}R^T$임을 쉽게 체크할 수 있으므로 직접적으로 역행렬을 계산할 필요가 없어진다.

$$\bf \Lambda^{-1/2} = \text{diag}( 1/\sqrt{\sigma_1},...,1/\sqrt{\sigma_n})$$

 
 
 
728x90

'Mathematics' 카테고리의 다른 글

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

열방정식에서 매질에서 열의 확산은 특별히 선호된 방향이 없이 모든 방향에 대해서 동등하다. 물론 매질의 열전도도가 다르면 국속적으로 열확산이 균일하지 않을 수는 있어도 방향을 따지는 않는다. 그러면 방향에 따라 열의 확산을 정도를 다르게 설정하려면 어떻게 해야 할까? 이는 열방정식을 이용해서 이미지의 smoothing을 하는 필터에서도 동일한 상황이 생길 수 있다. 특히 에지를 보존하면서 smoothing을 할 때 에지에 수직방향으로는 확산이 잘 안일어나게 억제해야 할 필요가 있다. 이를 위해 2차원 열방정식의 우변을 살펴보자. 우변은 온도 $u(x,y;t)$의 Laplacian이 들어오는데 보다 일반적인 Hessian matrix $H$을 써서 표현하면 그 의미가 더 명확해진다.

$$ H = \begin{pmatrix} u_{xx} & u_{xy}\\ u_{xy} & u_{yy} \end{pmatrix}$$

$$ u_{xx} + u_{yy} = (1, 0).H. \begin{pmatrix} 1 \\0\end{pmatrix}+ (0,1) . H. \left( \begin{matrix} 0 \\1 \end{matrix} \right) = \text{Tr}\left ((e_1 \otimes  e_1^T + e_2 e \otimes _2^T).H\right)$$

즉, Laplacian은 평면에서 두 단위 벡터  $e_1=(1,0)^T$과 $e_2 =(0,1)^T$ 각각에 Hessian을 작용한 결과를 다시 자신과 내적을 한 값을 동일하게 더하여 나온다. 따라서 $x$방향이나 $y$방향에 무관하게 등방적으로 작용하게 된다. 그러나 $x$방향과 $y$방향에 다른 가중치 $\alpha, \beta\ne \alpha$를 둔다면 등방성은 사라지게 될 것이고 가중치가 큰 방향으로 열의 확산이 더 잘 일어나게 된다.

$$ \alpha u_{xx} + \beta u_{yy}= \alpha \times (1, 0).H. \left(\begin{matrix} 1 \\0\end{matrix} \right)+ \beta(0,1) . H. \left( \begin{matrix} 0 \\1 \end{matrix} \right) = \text{Tr}\left((\alpha e_1 \otimes e_1^T + \beta e_2 \otimes e_2^T).H \right)$$  

일반적으로 $x,y$ 방향의 단위벡터 대신에 서로 직각이 단위벡터를 사용해서 비등방성을 줄 수 있다. 예를 들면 에지 방향에 수직한 방향으로 확산을 억제하고 싶은 경우에는 그래디언트의 공변행렬 $\left(\begin{matrix} u_x^2 & u_x u_y\\ u_x u_y& u_y^2 \end{matrix}\right)$의 두 고유벡터 $v_1$, $v_2$를 이용할 수 있다. 이 경우 고유값이 큰 쪽에 해당하는 고유벡터의 방향이 에지에 수직방향이고, 작은 쪽이 에지 방향을 가리킨다. 

$$    \text{Tr}\left( (c_1 v_1 \otimes v_1^T + c_2 v_2 \otimes v_2^T ).H\right) = \text{Tr}(T.H) \\T= c_1 v_1 \otimes v_1^T + c_2 v_2 \otimes v_2^T $$

이 경우 비등방적 확산은 텐서필드 $T$에 의해서 제어가 되고 열(확산)방정식은 다음과 같이 간단한 형식으로 표현할 수 있다.

$$ u_t = \text{Tr} (T.H)$$

텐서필드는 고유벡터의 텐서곱을 선형결합해서 만드는데 계수 $c_1$과 $c_2$를 어떻게 선택하는가에 따라 확산의 유형이 달라진다. 에지를 보존하도록 컨트롤하기 위해서는 고유값이 큰 방향으로 가중치가 작게 되도록 선택을 해야 한다. 많이 사용되는 예로는 고유값이 각각 $\lambda_1 < \lambda_2$일 때

$$ c_ 1 = (1 + (\lambda_1 + \lambda_2)/s)^{-p_1},  ~~c_2 = (1+ (\lambda_1 + \lambda_2)/s)^{-p_2}, ~~0<p_1 < p_2$$ 

로 잡는다.

 
 
 
 
 
 
 
 
 
 

void AnisoDiffusion(double **image, int width, int height) {
    // preparation;
    // ....

    // main iteration loop
    for (int iter = 0; iter < max_iter; iter++) {
        // compute gradients;
        Gradient(image, width, height, Gx, Gy);
        // compute the tensor field T, used to drive the diffusion
        for (int x = width; x-- > 0;) {
            for (int y = height; y-- > 0;) {			
                // covariance matrix;
                double a = Gx[x][y] * Gx[x][y]; 
                double b = Gx[x][y] * Gy[x][y];
                double d = Gy[x][y] * Gy[x][y];
                // eigenvalues of symm matrix [a,b;b,d]
                double lambda1, lambda2;
                Eigenvalues(a, b, d, lambda1, lambda2); // 0 <= lambda1 <= lambda2;
                // eigenvector for lambda1; ev1=(u,v);
                double u, v;
                Eigenvector(a, b, d, lambda1, u, v);
                // eigenvector for lambda2: ev2=(-v,u);

                double val1 = lambda1 / drange2; 
                double val2 = lambda2 / drange2;
                double f1 = pow(1 + val1 + val2, -a1);
                double f2 = pow(1 + val1 + val2, -a2);
                // T = f1 * (ev1.ev1^T) + f2* (ev2.ev2^T)
                T[0][x][y] = f1 * u * u + f2 * v * v; //(0,0)
                T[1][x][y] = (f1 - f2) * u * v;       //(0,1)=(1,0)
                T[2][x][y] = f1 * v * v + f2 * u * u; //(1,1)
            }
        }
        double xdt = 0.0;
            
        // compute the velocity and update the iterated image
        for (int x = width; x-- > 0;) {
            int px = max(x - 1, 0);
            int nx = min(x + 1, width - 1);
            for (int y = height; y-- > 0;) {
                int py = max(y - 1, 0);
                int ny = min(y + 1, height);
                double Ipp = image[px][py]; 
                double Ipc = image[px][y] ; 
                double Ipn = image[px][ny]; 
                double Icp = image[x] [py]; 
                double Icc = image[x] [y] ; 
                double Icn = image[x] [ny]; 
                double Inp = image[nx][py]; 
                double Inc = image[nx][y] ;
                double Inn = image[nx][ny];
                // Hessian matrix of image; H=[Ixx,Ixy;Ixy,Iyy]
                double Hxx = Inc + Ipc - 2 * Icc;
                double Hyy = Icn + Icp - 2 * Icc;
                double Hxy = (Ipp + Inn - Ipn - Inp) / 4;
                // veloc = trace(T.H)
                veloc[x][y] = T[0][x][y] * Hxx + 2 * T[1][x][y] * Hxy + T[2][x][y] * Hyy;
            }
        }

        // find xdt coefficient
        if (dt > 0) {
            double max1 = veloc[0][0], min1 = veloc[0][0];
            for (int x = width; x-- > 0;) {
                for (int y = height; y-- > 0;) {
                    if (veloc[x][y] > max1) max1 = veloc[x][y];
                    if (veloc[x][y] < min1) min1 = veloc[x][y];
                }
            }
            xdt = dt / max(fabs(max1), fabs(min1)) * drange0;
        } else xdt = -dt;

        // update image
        for (int x = width; x-- > 0;) {
            for (int y = height; y-- > 0;) {
                image[x][y] += veloc[x][y] * xdt;
                // normalize image to the original range
                if (image[x][y] < initial_min) image[x][y] = initial_min;
                if (image[x][y] > initial_max) image[x][y] = initial_max;
            }
        }
    }
    // clean memories;
};
 
 
 
 
 
 
 
 
 
728x90

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

Edge Preserving Smoothing  (0) 2024.02.14
Watershed Segmentation  (0) 2021.02.27
Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (3) 2021.02.15
Fixed-Point Bicubic Interpolation  (1) 2021.01.19
Posted by helloktk
,