Geometry & Recognition

열방정식에서 매질에서 열의 확산은 특별히 선호된 방향이 없이 모든 방향에 대해서 동등하다. 물론 매질의 열전도도가 다르면 국속적으로 열확산이 균일하지 않을 수는 있어도 방향을 따지는 않는다. 그러면 방향에 따라 열의 확산을 정도를 다르게 설정하려면 어떻게 해야 할까? 이는 열방정식을 이용해서 이미지의 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;
};

영상처리에서 $2\times2$ 행렬의 고유값과 고유벡터를 수치적으로 구해야 하는 상황이 많이 생긴다. 이 경우 2차인 특성방정식을 풀어서 고유값을 구하게 되는데, 두 고유값의 차이가 매우 큰 경우 단순한 근의 공식을 적용한 결과는 수치적으로 부정확한 값을 줄 수 있다. 이를 피하는 방법을 간단히 살펴보자.

이차방정식 

$$x^2 - b x +c =0$$

의 두 해는 

\begin{align}x_1 = \frac{b + d}{2},\quad x_2 = \frac{b - d}{2},\qquad d= \sqrt{b^2 - 4c}\end{align}

로 표현된다.  고유값의 특성방정식으로 보면 계수 $b$는 두 고유값의 합, 판별식 $d$는 두 고유값의 차이다. 그런데 $b ^2 \gg c $ 판별식 $d$는 $b$와 매우 근접한 수치이므로 거의 비슷한 두 수의 뺄셈으로 표현된 $x_2$의 수치계산 결과는 매우 부정확해질 수 있다. 예로 

$$ x^2 - 10^8 x +1=0$$

의 해를 프로그램으로 구하면 

	double b = 1.0e8, c = 1.;
	double d = sqrt(b * b - 4. * c);
	double x1 = (b + d)/2., x2 = (b - d)/2.;
	TRACE("x1=%e, x2=%e", x1, x2);

프로그램을 돌리면 얻게 되는 결과는

$$ \tt  x1=1.000000e+008, \qquad x2=7.450581e-009$$

그런데 테일러 전개를 사용하면 $x_2$는 

$$x_2 = \frac{b}{2}\left( 1 - \sqrt{1-\frac{4c}{b^2}} \right)  \approx \frac{b}{2}\times \frac{2c}{b^2} = \frac{c}{b}=10^{-8}$$ 

이어야 한다. 따라서 $x_2$를 바로 근의 공식으로 구할 경우 25%정도의 오차가 생긴다. 이런 경우에는 근의 공식에 의한 직접적인 계산보다는 이차방정식의 특성을 이용하는 것이 수치적으로 더 정교한 답을 얻을 수 있다. 두 근의 곱이 $x_1 x_2 =c$이므로 $x_2$는 

$$x_2 = c / x_1\approx 10^{-8}$$

을 통해서 구하는 것이 더 수치적으로 정교하다.  일반적으로 두 근의 합  $b$가 양수인 경우는  $x_2$를, 음수인 경우는 $x_1$을 두 근의 곱을 이용해서 계산하는 것이 더 정교한 수치값을 얻을 수 있다.

// Find two roots of a quadratic eq: x * x - b * x + c = 0;
void Quadratic2(double b, double c, double& x1, double& x2) {
    double d = sqrt(b * b - 4. * c);
    x1 = (b + d)/ 2., x2 = (b - d)/ 2.;
    TRACE("x1=%e, x2=%e\r\n", x1, x2);
    if (b > 0) {
        if (x1 != 0) x2 = c / x1;
        TRACE("Improved: x1=%e, x2=%e\r\n", x1, x2);
    } else {
        if (x2 != 0) x1 = c / x2;
        TRACE("Improved: x1=%e, x2=%e\r\n", x1, x2);
    }
}

이 경우 결과는 훨씬 더 정교한 값을 보여준다.

$$\tt Improved: x1=1.000000e+008, \quad x2=1.000000e-008$$