일정한 간격 $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' 카테고리의 다른 글

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
Statistical Region Merging  (0) 2022.06.11
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' 카테고리의 다른 글

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
Statistical Region Merging  (0) 2022.06.11
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
,

\begin{gather}BF[I]_{\bf p} = \frac{1}{W_\bf{p}}  \sum_{ {\bf q} \in S} G_{\sigma_s} (||{\bf p} - {\bf q}||) G_{\sigma_r}(|I_{\bf p} - I_{\bf q} |) I_{\bf q} \\  W_{\bf p} = \sum_{{\bf q}\in S} G_{\sigma_s} (||{\bf p}-{\bf q}||) G_{\sigma_r}(|I_{\bf p} - I_{\bf q} |) \\ G_\sigma ({\bf r}) = e^{ - ||\bf{r}||^2/2\sigma^2 }\end{gather}

 

Bilateral Filter
Gaussian Filter

 

smoothing based on the nonlinear heat eq

// sigmar controls the intensity range that is smoothed out. 
// Higher values will lead to larger regions being smoothed out. 
// The sigmar value should be selected with the dynamic range of the image pixel values in mind.
// sigmas controls smoothing factor. Higher values will lead to more smoothing.
// convolution through using lookup tables.
int BilateralFilter(BYTE *image, int width, int height, 
    double sigmas, double sigmar, int ksize, BYTE* out) {
    //const double sigmas = 1.7;
    //const double sigmar = 50.;
    double sigmas_sq = sigmas * sigmas;
    double sigmar_sq = sigmar * sigmar;
    //const int ksize = 7;
    const int hksz = ksize / 2;
    ksize = hksz * 2 + 1;
    std::vector<double> smooth(width * height, 0);
    // LUT for spatial gaussian;
    std::vector<double> spaceKer(ksize * ksize, 0);
    for (int j = -hksz, pos = 0; j <= hksz; j++) 
        for (int i = -hksz; i <= hksz; i++) 
            spaceKer[pos++] = exp(- 0.5 * double(i * i + j * j)/ sigmas_sq); 
    // LUT for image similarity gaussian;
    double pixelKer[256];
    for (int i = 0; i < 256; i++)
        pixelKer[i] = exp(- 0.5 * double(i * i) / sigmar_sq);

    for (int y = 0, imgpos = 0; y < height; y++) {
        int top = y - hksz;
        int bot = y + hksz;
        for (int x = 0; x < width; x++) {
            int left = x - hksz;
            int right = x + hksz;
            // convolution;
            double wsum = 0;
            double fsum = 0; 	
            int refVal = image[imgpos];
            for (int yy = top, kpos = 0; yy <= bot; yy++) {
                for (int xx = left; xx <= right; xx++) {
                    // check whether the kernel rect is inside the image;
                    if ((yy >= 0) && (yy < height) && (xx >= 0) && (xx < width)) {
                        int curVal = image[yy * width + xx];
                        int idiff = curVal - refVal;
                        double weight = spaceKer[kpos] * pixelKer[abs(idiff)];
                        wsum += weight;
                        fsum += weight * curVal;
                    }
                    kpos++;
                }
            }
            smooth[imgpos++] = fsum / wsum;
        }
    }

    for (int k = smooth.size(); k-- > 0;) {
        int a = int(smooth[k]);
        out[k] = a < 0 ? 0: a > 255 ? 255: a;
    }
    return 1;
}
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

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

Cubic Spline Kernel  (1) 2024.03.12
Ellipse Fitting  (0) 2024.03.02
영상에서 Impulse Noise 넣기  (2) 2023.02.09
Canny Edge: Non-maximal suppression  (0) 2023.01.11
Statistical Region Merging  (0) 2022.06.11
Posted by helloktk
,

열방정식은 매질에서 열이 전달되는 과정을 기술한다. 열은 온도차이가 있을 때 전달되는 에너지로 온도 분포 $u(\vec{r},t)$가 주어진 경우 주변으로 나가는 열에너지는 온도의 gradient $\nabla u$에 비례한다. 

$$ \text{heat flux: }\quad  \vec{j}= \kappa (u)\nabla u$$

$\kappa$는 매질에서 열이 얼마나 잘 전도는지를 표현하는 일반적으로 $u$에 의존할 수 있다. 균일하고 등방적인 매질인 경우는 $\kappa$가 일정한 상수에 해당한다. 국소적으로 단위부피당 단위시간당 빠져나가는 열에너지는 thermal flux의 divergence에 비례하고 이 값이 0이 아니면 그 지점의 온도 변화를 만든다. 따라서 온도는 다음과 같은 방정식을 만족한다.

$$ u_t = \nabla \cdot \vec{j}  = \nabla\cdot(\kappa \nabla u)$$

이 방정식은 물질의 확산 현상에도 적용이 가능한데 이 경우 $u$는 물질의 밀도분포에 해당할 것이다.

$\kappa$가 상수인 경우 이 방정식은 선형방정식으로 초기 온도분포 $u(x,t=0)=f(x)$가 주어지면 해는 사이즈가 $\sigma = \sqrt{2\kappa t}$인 가우시안 커널(heat kernel 또는 Greeen function)과 초기 온도분포의 convolution으로 표현할 수 있다. 1차원 선형 열방정식의 해는 구체적으로 다음과 같이 표현된다.

\begin{align} u(x,t) =  \int \frac{1}{\sqrt{4\pi \kappa t}} e^{- (x-y)^2 / 4\kappa t } f(y) dy = ( G_{\sigma} * f)(x,t)\\G_\sigma (x) = \frac{1}{\sqrt{2\pi\sigma} }  e^{- x^2/2\sigma^2} \end{align}

열방정식은 뜨거운 커피가 담긴 잔을 방에 놓았을 때 방 안의 온도가 어떻게 변할지도 알 수 있게 해 준다. 물론 경험적으로 커피의 온도는 내려가고 방안은 잔의 주변부터 점차로 데워져서 나중에서는 모든 부분이 일정한 온도를 가지는 평형상태에 도달한다. 이는  가우시안 커널의 크기가 $\sqrt{t}$에 비례하므로 시간이 흐를수록 모든 지점에서 온도는 일정한 값으로 수렴할 것임을 쉽게 예측할 수 있다.

이 열방정식을 이미지에 대해서도 적용할 수 있다. 이 경우 초기 온도분포는 주어진 이미지의 명암값에 해당한다. 열방정식의 우변을 이미지에 적용하면 Gaussian 블러링에 해당하고, 시간은 블러링의 순차적인 적용 횟수를 의미한다. 이미지가 블러링 되는 경우 영상이 담고 있는 정보는 점점 잃게 되는데 이를 방지하기 위해서는 영상의 에지 정보를 보존하는 방법으로 열방정식을 변형시켜 보자. 영상에서 에지영역에서는 gradient 값이 커지므로 이 지점에서는 커널사이즈를 작게 (열방정식 관점에서는 열전도도를 줄임) 만들면 블러링 효과가 줄어들게 된다. 이를 위해서 다음과 같은 열전도도 함수를 도입하자.

$$ \kappa (u) = \frac {1}{1+ |\nabla u|^2/\lambda^2 }$$

여기서 $\lambda$은 contrast 파라미터로 에지의 세기가 $|\nabla| \gg \lambda$ 영역에서는 블러링의 영향을 거의 받지 않게 되어 열방정식을 적용하더라도 에지에 대한 정보 손실은 거의 일어나지 않는다. 그러나 $\kappa $가 상수가 아닌 경우는 열방정식은 비선형이므로 해를 명시적으로 쓸 수 없고 오직 수치해석적으로만 구할 수 있다.

먼저 미분연산자(Sobel operator: $\lambda$의 설정은 미분연산자의 normalization을 고려해야 함)을 이용해서 이미지의 gradient을 구하면 flux $\vec{j}$을 구할 수 있고, 다시 $j_x$에 대해서 $x-$방향 미분, $j_y$에 대해서 $y-$방향 미분을 해서 $\vec{j}$의 divergence, 즉 이미지의 일반화된 Laplace 연산 결과를 얻을 수 있다. 이 결과에 시간간격을 곱한 양을 이전 이미지에 더하여 새로운 이미지를 얻는다. 이 과정을 반복적으로 수행하면 된다.

$$ u(x,y, t+\Delta t) = u(x, y, t) + \Delta t \times \nabla \cdot (\kappa (|\nabla u| ) \vec{j})$$

$\kappa$가 상수일 때 $t \to \kappa t$로 시간변수를 바꾸면 열방정식은 차원이 없는 형태로 되고, $\kappa <1$이므로 시간간격은 $\Delta t \ll 1$를 만족하도록 선택해야 한다.

 
 
 

 

$\lambda=50,~\Delta t=0.01$일 때 5 스텝 간격으로의 변화

 
double conductivity(double gmagsq, double lambdasq) {
    return 1.0/ (1.0 + gmagsq/lambdasq);
}
int LaplaceImage(double *img, int w, int h, double lambda, double *laplace) {
    const double lambdasq = lambda * lambda;
    // Sobel gradient 3x3
    const int nn[] = { -w - 1, -w, -w + 1, -1, 0, 1, w - 1, w, w + 1};
    const int sobelX[] = { -1, 0, 1, -2, 0, 2, -1, 0, 1};
    const int sobelY[] = { -1, -2, -1, 0, 0, 0, 1, 2, 1};
    std::vector<double> Gx(w * h, 0);
    std::vector<double> Gy(w * h, 0);
    std::vector<double> Gmag2(w * h, 0);
    const int xmax = w - 1, ymax = h - 1;
    // gradient-x, gradient-y, and mag of gradient.
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                double v = img[pos + nn[k]];
                sx += sobelX[k] * v;
                sy += sobelY[k] * v;
            }
            Gx[pos] = sx;
            Gy[pos] = sy;
            // gx^2 + gy^2;
            Gmag2[pos] = sx * sx + sy * sy;
        }
        pos++; // skip x=xmax;
    }
    // flux;
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double kappa = conductivity(Gmag2[pos], lambdasq);
            // multiply conductivity;
            Gx[pos] *= kappa;
            Gy[pos] *= kappa;
        }
        pos++; // skip x=xmax;
    }
    // divergence -> laplace;
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                sx += sobelX[k] * Gx[pos + nn[k]];
                sy += sobelY[k] * Gy[pos + nn[k]];
            }
            laplace[pos] = sx + sy;
        }
        pos++; // skip x=xmax;
    }
    return 1;
}
void EpSmooth(BYTE *image, int w, int h, BYTE *out) {
    const double lambda = 50.0;
    const double dt = 0.01; //time step;
    const int maxiter = 100;
    std::vector<double> fimage(w * h, 0);
    std::vector<double> laplace(w * h, 0);
    for (int k = fimage.size(); k-- > 0;) fimage[k] = image[k];
    for (int iter = maxiter; iter-- > 0;) {
        LaplaceImage(&fimage[0], w, h, lambda, &laplace[0]);
        // update image;
        for (int k = fimage.size(); k-- > 0; ) 
            fimage[k] += dt * laplace[k];   
    }
    // preparing output;
    for (int k = fimage.size(); k-- > 0; ) {
        int a = int(fimage[k]);
        out[k] = a > 255 ? 255: a < 0 ? 0: a;
    }
}

 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

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

Anisotropic Diffusion Filter  (0) 2024.02.23
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
,