Geometry & Recognition

영상처리에서 영상에서 분리된 객체의 형상을 기술할 때 moment를 많이 사용한다. 영상은 각 픽셀 위치에서 컬러 값이 할당된 일종의 2 변수 함수로 생각할 수 있다. 영상의 모멘트는 픽셀 위치(의 단항식)를 컬러값으로 가중치를 주어서 낸 평균이라고 간단하게 기술할 수 있다.  이미지가 $I(x, y)$로 주어진 경우 $p+q$-차 moment는

$$m_{pq} = \iint x^p y^q I(x, y) dxdy, \quad p, q=0,1,2,3...$$

로 정의한다. 그러나 이 정의는 똑같은 형상이라도 원점에서 얼마나 떨어진 위치에 있는가에 따라 결과가 달라지므로 (원점의 선택에 의존되는 정의다. 객체를 구성하는 픽셀 개수를 알려주는 $m_{00}$는 불변이다.) 형상 기술에 적합한 형태가 아니다. 이 문제는 객체의 질량중심에 대한 상대 위치에 대한 weighed mean인 central moment를 feature 값으로 이용하면 된다:

$$ \mu_{pq} = \iint (x - \overline{x})^p ( y - \overline{y})^q I(x, y) dx dy$$

여기서 

$$\overline{x} = m_{10}/m_{00}, \quad \overline{y} = m_{01}/m_{00}$$

로 객체의 질량중심(center of mass)의 위치다. 0차 central moment는 원래의  moment와 같고, 1차 central moment는 질량중심이 원점이므로 0으로 주어진다: $\mu_{10}=\mu_{01}=0$. 이 central moment는 이미지에서 객체가 이동하더라도 동일한 값을 가진다.

그런데 같은 객체를 담은 이미지라도 확대하거나 축소하면 central moment의 값은 변하게 되므로 영상의 scale이 달라지는 경우는 central moment를 feature로 한 형상 판별하는 작업에는 적당하지 않을 수 있다. 이미지를 전체적으로 축소하거나 확대하더라도 같은 값을 주는 scale invariant인 moment 정의가 필요한데, 이는 central moment를 정규화시켜 사용하면 된다.

우선 이미지를 $\lambda$ 만큼 확대하거나 축소하면 변환된 이미지는 $I'(x', y') = I(x'/\lambda, y'/\lambda)$로 표현된다. 이 이미지에서 central moment를 계산하면

\begin{align} \mu'_{pq} &= \iint (x'- \overline{x'})^p ( y' - \overline{y'})^q I(x'/\lambda , y'/\lambda ) dx'dy'\\ &= \lambda^{p+q+2} \iint (x- \overline{x} )^p(y -\overline{y})^q I(x, y) dx dy \\ &= \lambda^{p+q+2} \mu_{pq} \end{align}

즉, 이미지를 scaling하면 $\mu_{pq}$는 scaling factor $\lambda$의 ${p+q+2}$ 지수승만큼 변하게 된다. 이 변화를 없애기 위해서 양변을 제일 단순한 central moment인 $\mu_{00}^\gamma$로 나누어서 전체적인 scale factor가 사라지게 만들자. 그런데 $$\mu'_{00} = \lambda ^2 \mu_{00}$$

이므로

$$ \frac{ \mu'_{pq} }{ (\mu'_{00})^\gamma} = \lambda^{p+q+2-2\gamma} \frac{\mu_{pq}}{ (\mu_{00})^\gamma } $$

$\gamma$를

$$ \gamma = \frac{p+q}{2}+1$$

로 선택하면 다음의 normalized central moment는 scaling에 대해서 불변인 성질을 가진다.

$$ \eta_{pq} = \frac{\mu_{pq}}{\mu_{00}^\gamma}$$

이진 이미지에서 지름 $2a$인 원과 한 변이 $2a$인 정사각형에 대해서 2차 central moment을 구해보면

$$\text{circle: }~  \mu_{20}= \mu_{02}=\frac{\pi}{4} a^4 = 0.7854 a^4 ,~~\mu_{11}=0,$$

$$\text{rectangle: }~\mu_{20}=\mu_{02}=\frac{4}{3} a^4 = 1.3333 a^4, ~~\mu_{11}=0.$$

그리고 $\mu_{00}= \pi a^2~(\text{circle})$, $\mu_{00}= 4 a^2 ~(\text{rectangle})$이므로 2차 normalized central moment는

$$\text{circle: }~ \eta_{20}=\eta_{02} = \frac{1}{4\pi} = 0.07956,\quad \eta_{11}=0, $$

$$\text{rectangle:}~\eta_{20}=\eta_{02}=\frac{1}{12}=0.08333,\quad \eta_{11} = 0.$$

normalized central moment는 크기에 의존하지는 않지만 원과 사각형에서 분명한 차이를 보이므로 이미지에서 분리된 객체를 구별하는 기준으로 삼을 수 있을 것이다.

물론 normalzed central moment는 객체를 회전시키면 그 값이 다시 변하게 된다. 하지만 이들을 잘 조합하면 회전 불변인 성질까지 추가되는 invariant moment를 구성할 수 있다. 이차 central moment의 경우 질량중심에 대한 회전변환을 고려할 때 서로 다르게 변환하는 부분이 섞여 있는데 이를 분리하면

\begin{gather} \mu_{20}= \iint \left[ \frac{1}{2}(x^2-y^2) + \frac{1}{2}(x^2+y^2) \right] f(x, y) dxdy =T_{11} +\delta_{11} S ,\\ \mu_{11}= \iint xy I(x, y) dxdy = T_{12}+\delta_{12}S, \\ \mu_{02} = \iint \left[ -\frac{1}{2}(x^2 - y^2) +\frac{1}{2}(x^2+y^2) \right] dx dy= T_{22} + \delta_{22} S ,\\ T_{ij} = \iint \left[x_i x_j -\frac{1}{2}(x^2+y^2) \delta_{ij} \right] I(x, y) dxdy=\frac{1}{2} \left(\begin{array}{cc}\mu_{20} - \mu_{02} & 2\mu_{11} \\ 2\mu_{11}& -(\mu_{20}-\mu_{02}) \\\end{array} \right),\\ S= \iint (x^2 + y^2) I(x, y)dx dy =\mu_{20}+\mu_{02}   \end{gather}

로 쓸 수 있는데, $T_{ij}$는 회전 변환에 대해서 rank가 2인 tensor의 성분을 구성하고, $S$는 $rank=0$인 scalar가 된다.  traceless인 $rank=2$인 텐서를 가지고 만들 수 있는 회전에 대해 invariant 한 양은 tensor의 determinant ($\det(T)$)인데, 이를 normalized central moment로 표현하면 (scaling까지 고려하기 위해서)

$$ \text{det}(T) \propto (\eta_{20} - \eta_{02})^2 + 4 \eta_{11}^2 \quad (\text{up to a numerical factor})$$

로 쓸 수 있고, $rank = 0$ 인 부분($S$)은 그 자체로 회전에 불변이다. (물리적으로는 질량중심축에 대한 rotational inertial에 대응하는 양이다). 

$$ \eta_{20} + \eta_{02}= \frac{1}{\mu_{00}^2} \iint (x^2 + y^2 ) I(x, y) dxdy$$

다른 회전 invariant 한 양은 order가 3 이상인 central moment를 이용해서 만든 rank=3 이상인 tensor의 invariant을 이용해서 얻을 수 있다.

평면에서 주어진 벡터장의 orientation을 찾는 문제는 영상처리 알고리즘에서 자주 접하게 된다. 벡터의 방향은 두 성분의 부호와 상대적인 크기에 따라 달라지지만, 기준선에 대해 상대적으로 기울어진 정도를 나타내는 orientation은 정반대 방향의 두 벡터 $(v_x, v_y)$와 $(-v_x, -v_y)$에 같은 값이 부여되고, 그 값은 벡터의 $x$ 성분과 $y$ 성분의 비의 arctangent 값

$$\theta=\tan^{-1} \left( \frac {v_y}{ v_x }\right)$$

로 계산할 수 있다. 

영상처리에서는 영상에 내재하는 잡음에 의한 영향을 줄이기 위해 한 지점에서 orientation을 추정할 때 보통 그 지점 주변의 벡터 성분의 평균을 이용한다. 주변에서 정반대 방향의 두 벡터가 있는 경우 이 두 벡터는 기하학적으로 같은 orientation을 주지만 더하는 경우 서로 상쇄되어 평균에는 기여가 없으므로 위 식을 사용하면 잘못된 예측을 줄 수 있다. 따라서 잡음을 고려한 상황에서 좀 더 robust 하게 orientation을 추정할 수 있는 방법이 있어야 한다. 벡터 성분의 상대적인 부호만 고려하는 식으로 바꾸기 위해서 $\tan \theta$ 대신에 $\tan (2\theta)$를 고려하자.

\[ \tan (2\theta) = \frac {2\tan(\theta)}{1-\tan^2(\theta)}=\frac {2v_x v_y}{v_x^2 - v_y^2}. \]

분모에서는 각 성분의 제곱, 분자는 두 성분의 곱으로 표현되므로 성분 사이의 상대부호가 같은 경우에는 우측식은 같은 값을 주므로 분모, 분자를 주변 평균값 $v_x v_y ~\longrightarrow ~<v_x v_y>$, $v_x^2 - v_y^2 ~\longrightarrow ~<v_x^2> - < v_y^2>$으로 대체하여도 올바른 orientation을 주게 된다. orientation 각도는

\[ \theta = \frac {1}{2} \tan^{-1}\left( \frac {2 <v_x v_y>}{ <v_x^2> - <v_y^2>}  \right)  \]

으로 주어진다. 실제 계산은 인자가 singular해지는 경우를 피하기 위해서 $\text {atan2}()$ 함수를 사용한다.

https://kipl.tistory.com/293

https://kipl.tistory.com/111

https://kipl.tistory.com/58

주어진 이미지(destination image)의 일부분($\Omega$)을 다른 이미지(source image: $g(x, y)$)로 자연스럽게 대체하는 seamless blending의 한 기법인 Possion Image Editing을 살펴보자. 수학적으로는 destination image의 한 영역 $\Omega$에서 source image를 보간(interpolation)하는 함수 $f(x,y)$를 찾는 작업이다. $f$는 보간 영역의 경계($\partial\Omega$)에서는 destination image와 같은 컬러 값($f^{*}$)을 가져야 한다.

그럼 보간의 기준은 무엇이 되어야 하는가? blending이 잘되기 위해서는 보간 영역에서는 source image와 같은 형상으로 보여야 한다. 이미지에서 형상의 인식은 컬러값 자체가 아니라 컬러의 변화인 gradient에 의해서 결정이 된다. 따라서 보간 함수는 source image의 gradient를 최대한 유지하도록 선택되어야 한다. 이제 seamless blending을 주는 보간 함수를 찾는 문제는 source image가 blending 영역의 경계에서는 destination 이미지의 color값을 가지면서 영역 내부에서는 source image의 gradient 값을 최대한 유지하는 함수를 찾는 변분 문제로 귀결된다: source image의 gradient를 $\mathbf{v}(x, y)=\nabla g= (\partial g/\partial x, \partial g / \partial y)$라 하면,

$$ \underset{f}{\text{argmin}} \iint_{\Omega}  | \nabla f - \mathbf{v}|^2 dxdy, \quad\text {with} \quad f = f^*~\text {on}~\partial \Omega.$$ 

위 식을 $f$에 대해서 variation을 취하면$$\iint 2(\nabla \delta f)\cdot(\nabla f - \mathbf{v}) dxdy = - 2 \iint \delta f ( \nabla^2 f - \nabla \cdot \mathbf{v}) dxdy + 2 \iint  \nabla \cdot \left[ \delta f (\nabla f - \mathbf{v})\right] dxdy    $$을 얻을 수 있는데 마지막 항은 total derivative이므로 기여가 없다. 따라서 보간함수 $f$는 다음 방정식을 만족해야 한다:

$$ \nabla^2 f = \nabla \cdot \mathbf{v}, \quad \text {or} \quad \frac {\partial^2 f}{\partial x^2} + \frac {\partial ^2 f}{\partial y^2 } = \frac {\partial v_x }{\partial x}+ \frac {\partial v_y}{\partial y}.$$

보간함수를 찾는 과정은 domain $\Omega$에서 Dirichlet boundary condition $f=f^* ~\text {on} ~\partial\Omega$이 부여된 Poisson 방정식의 해를 찾는 것과 같다는 것을 알려준다.

$\nabla \cdot\mathbf{v}= \nabla^2 g$이므로 Poisson 방정식의 해를 $f(x,y)= g(x,y) + h(x,y)$로 쓰면 풀어야 할 문제는 다음의 Laplace 방정식으로 해를 구하는 문제로 바뀐다.

$$ \nabla^2 h = 0, \quad ~~ h|_{\partial\Omega} = (f^* - g)|_{\partial\Omega}$$

** 1차원 예: $f^*(x)=1-\frac{x}{\pi}$, $g(x)=\frac{1}{2} \cos (4x)$ 인 경우, $f(x)=g(x) + h(x)$로 쓰면, $h(x)$는 $h(x=0) = f^*(0) - g(0) = 0.5$, $h(\pi)=f^*(\pi)-g(\pi) = -\frac{1}{2}$의 경계조건을 만족하는 $d^2h/dx^2 = 0$의 해로 주어지는데, $h(x)=\frac{1}{2}- \frac{x}{\pi}$임을 쉽게 알 수 있다. 따라서 $f(x)= \frac{1}{2}-\frac{x}{\pi} +\frac{1}{2} \cos(4x)$로 쓰인다:

comb 함수: 일정한  간격($T$)으로 주어진 message를 샘플링하는 함수.

$$\text{comb}_T(t) := \sum_{n=-\infty}^{\infty} \delta (t- nT)$$

주기가 $T$인 함수다: 

$$ \text{comb}_T(t) = \text{comb}_T(t+T)$$

따라서 Fouries series 전개가 가능하다.

$$\text{comb}_T(t) = \sum_{n=-\infty}^{ \infty}   c_n e^{i 2\pi n t /T}$$

계수 $c_n$은?

\begin{align} c_n :=&\frac{1}{T} \int_{-T/2}^{T/2} \text{comb}_T(t) e^{-i 2\pi n t /T}dt \\ =&  \frac{1}{T}\int_{-T/2}^{T/2} \delta (t) e^{ -i 2\pi n t/T} dt  \\ =&\frac{1}{T} e^{-i 2\pi n (0) /T } = \frac{1}{T}.\end{align}

따라서, $\text{comb}_T(t)$의 Fourier series 전개는 

$$\text{comb}_T(t) = \frac{1}{T} \sum_{n = -\infty}^{\infty} e^{i2\pi  n t/T}.$$

frequency domain에서 delta 함수의 역 Fourier transform은 정의에 의해서

$$ {\cal F}^{-1} [\delta (f-f_0)] = \int \delta (f-f_0) e^{i 2\pi t f }df = e^{i 2\pi t f_0}$$

그럼 $\text{comb}_T(t)$의 Fourier transform은 어떻게 표현되는가?

\begin{align} {\cal F}[\text{comb}_T(t)]=& \frac{1}{T} \sum {\cal F}[ e^{i2\pi n t/T}] \\  =& \frac{1}{T} \sum {\cal F}[ {\cal F}^{-1} [ \delta (f - n/T)]] \\ =& \frac{1}{T} \sum_{n = -\infty}^{\infty} \delta(f - n/T).\end{align} $\text{comb}$ 함수의 Fourier 변환은 frequency domain에서 $\text{comb}$ 함수이고 (up to constant factor),  time domain에서 주기가 $T$일 때 frequency domain에서는 $ 1/T$의 주기를 가진다.

주어진 message $m(t)$에서 일정한 간격 $T$로 샘플링된 message $m_s(t)$는 $\text{comb}$ 함수를 이용하면

$$m_s (t) : = m(t) \text{comb}_T(t) = \sum m(nT) \delta (t- nT)$$

로 표현된다.

양변에 Fourier transform을 적용하면,

\begin{align} M_s(f) = {\cal F} [m_s ] =& {\cal F}[m(t) \text{comb}_T(t)]= {\cal F}[m] * {\cal F}[\text{comb}_T]  \\  =& \frac{1}{T} \sum \int \delta(f' - n /T) M(f- f') df' \\ =& \frac{1}{T} \sum_{n=- \infty}^{\infty}  M(f - n/T) . \end{align}

따라서 message의 spectrum이 band-limited이고, $\text{band-width} \le \frac{1}{2} f_s = \frac{1}{2T}$인 조건을 만족하면 샘플링된 데이터를 이용해서 원 신호를 복원할 수 있다.

이 경우에 low-pass filter 

$$ H(f/f_s) := T \cdot \text{rect}(f/f_s)=\left\{ \begin{array}{ll} T ,& |f/f_s| < 1/2 \\ 0 , & |f/f_s| >1/2,\end{array}\right. $$

을 sampled massage의 Fourier transform에 곱해주면, 원 message의 Fourier transform을 얻는다:

$$ M(f) = H(f) M_s(f).$$

그런데 $M_s(f)$는 frequency domain에서 주기가 $f_s = 1/T$인 주기함수이므로 Fourier series로 표현할 수 있다:($\text{comb}_T$ 함수와 같은 방식으로 하면 계수를 쉽게 찾을 수 있다: Poisson summation formula)

$$M_s(f) = \frac{1}{T}\sum M(f- n f_s ) = \sum_{-\infty}^\infty   m(nT) e^{-i 2\pi nf T}$$

따라서, 

\begin{align} M(f) = & H(f/f_s)   M_s(f) \\ =& H(f/f_s)  \sum  m(nT) e^{-i 2\pi nfT} \\ =& \sum m(nT) \Big(\text{rect}(f/f_s )T  e^{-i 2\pi nTf} \Big) \\=& \sum m(nT) {\cal F} \left[  \text{sinc}  \left( \pi  \frac{t-nT}{T} \right) \right]\end{align}

이므로 양변에 역 Fourier transform을 하면 sampled 된 message $\{m(nT)|n\in Z\}$를 이용해서 원 message를 복원할 수 있는 식을 얻을 수 있다(Whittaker-Shannon interpolation):

$$m(t) = \sum_{n=-\infty}^{\infty}  m(nT) \,\, \text{sinc}\left(  \pi \frac{t-nT}{T} \right) .$$

Lanczos kernel:  a sinc function windowed by the central lobe of a second, longer, sinc function(wikipedia)

$$K(x) = \left\{ \begin {array}{ll}   \text {sinc}(\pi x) \text {sinc}\Big( \frac {\pi x}{a}\Big), & |x| < a \\ 0 ,& |x| \ge a\end {array} \right. = \text{sinc}(\pi x) \text{sinc}\Big( \frac{\pi x}{a}\Big)\text{Box}(|x|<a) $$

$$H(\omega)=\frac{1}{2\pi^2\sqrt{2\pi}}\Big( -(2\pi-3\omega)\text{Si}(2\pi-3\omega)+(4\pi -3\omega)\text{Si}(4\pi-3\omega)\\+(2\pi +3\omega)\text{Si}(2\pi+3\omega)+(4\pi+3\omega)\text{Si}(4\pi+3\omega)\Big).$$

여기서 sine intergral은 다음과 같이 정의된다:

$$\text{Si}(x)=\int_0^x \frac{\sin(t)}{t} dt.$$

$x=0$ 근방에서는 Tayor 전개식을 이용하는 것이 더 효율적인다.

Lanczos kernel을 사용해서 일정한 간격으로 sampling 된 데이터 $\{f_i \}$를 보간하려고 할 때 보간함수 $F(x)$는 Lanczos kernel과 데이터의 convolution으로 표현할 수 있다.

$$ F(x) = \sum_{i \in \text{window}} f_i K(x - i) $$

Interpolation 함수는 영상의 resampling 과정에서 사용할 수 있다. 주어진 이미지를 확대하거나 축소하려면 기존 이미지의 픽셀과 픽셀 사이 구간에서 픽셀 값을 알아내야 하는데, 이때 interpolation 함수를 이용하여 필요한 데이터를 얻는 과정인 resampling 한다.

다음 코드는 Lanczos3 kernel을 사용해서 이미지를 확대하거나 축소할 때 행에 대해서 필요한 resampling을 수행한다. 여기서는 kernel의 중심이 소스 이미지의 픽셀과 픽셀의 가운데에 놓이도록 설정하였다($-0.5$의 의미). 이는 작은 영상을 크게 확대할 때 가장자리에서 왜곡이 되는 것을 방지하기 위해서 인데 큰 영상을 확대/축소할 때는 영향이 없다. 이 interpolation은 separable이므로 2차원의 경우 행방향으로 진행한 후 그 결과를 다시 열 방향으로 진행하는 2-pass 방식으로 사용해도  된다.

// windowed sinc(x); window=sinc(x/3);
static double Lanczos3(double x){
    if (x < 0) x = -x; // symmetric;
    x *= PI;
    if (x < 0.01)    return 1. + (- 5./ 27. + 
        (14. / 1215. - 17. * x * x / 45927.) * x * x) * x * x; 
    else if (x < 3.) return 3. * sin(x) * sin(x / 3.) / x / x;
    else     	     return 0;
};
// interpolation in the horizonal direction: single channel or gray image;
// x = pixel position in a destination image;
double Lanczos3_line(BYTE *src_line, int srcWidth, int x, double xscale) {
    double halfWidth;
    if (xscale > 1.) halfWidth = 3.;
    else             halfWidth = 3. / xscale;

    double centerx = double(x) / xscale - 0.5;  //center loc of filter in the src_line;
    int left  = (int)floor(centerx - halfWidth);
    int right = (int)ceil(centerx + halfWidth);
    if (xscale > 1) xscale = 1;
    double s = 0;
    double weightSum = 0;
    for (int ix = left; ix <= right; ix++) {   
        double weight = Lanczos3((centerx - ix) * xscale);
        int xx = ix;         // for ix<0 || ix>=srcWidth: repeat boundary pixels
        xx = min(max(xx, 0), srcWidth - 1));
        s += weight * src_line[xx];
        weightSum += weight;
    }
    return s / weightSum; 
}

이차원 이미지의 경우 수평방향의  적용한 후 다시 수직방향으로 적용하면 된다.