728x90

영상처리에서 영상에서 분리된 객체의 형상을 기술할 때 moment를 많이 사용한다. 영상의 모멘트는 픽셀 위치(의 정수 지수 함수)를 픽셀 값으로 가중치를 주어서 낸 평균이라고 간단하게 기술할 수 있고, 이미지가 $I(x, y)$로 주어진 경우(배경은 $I(x,y)=0$인 픽셀이다)

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

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

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

로 객체의 질량중심의 위치다. 0차 central moment는 원래의  moment와 같고, 1차 central moment는 질량중심이 원점이므로 0으로 주어진다. 이 central moment는 이미지에서 객체가 이동하더라도 동일한 값을 가진다.

 

그런데 같은 객체를 담은 이미지라도 확대하거나 축소하면 central moment의 값은 변하게 되므로 central moment를 계산해서 객체의 형상을 판별하는 작업에는 적당하지 않을 수 있다. 이미지를 전체적으로 축소하거나 확대하더라도 같은 값을 주는 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{align} \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 \end{align}

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

$$ (\eta_{20} - \eta_{02})^2 + 4 \eta_{11}^2 \quad (\text{up to overall scale})$$

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

$$ \eta_{20} + \eta_{02}$$

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

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

Image Moments  (0) 2021.12.04
Orientation 추정  (0) 2021.11.30
Poisson Image Editing  (0) 2021.08.03
Sampling Theorem  (0) 2021.05.12
Lanczos Resampling  (0) 2021.05.08
Interpolation Kernels  (0) 2021.05.05
Posted by helloktk

댓글을 달아 주세요

728x90

 

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

$$\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

 

Local Ridge Orientation

지문에서 ridge의 방향(orientation)은 gradient에 수직한 방향이다(그런데 벡터인 gradient와는 달리 ridge의 방향은 모호함이 있다. 시계방향 또는 반시계방향으로 90도 회전이 모두 동일한 ridge의 방향이

kipl.tistory.com

https://kipl.tistory.com/111

 

Ellipse Parameters

원뿔을 평면으로 잘랐을 때 나타나는 곡선인 conic section은 직교 좌표계에서 $(x, y)$에 대한 2차 형식으로 쓰인다. 이 conic section이 타원을 기술할 때 parameter {$a, b, c, d, e, f$}를 이용해서 타원의..

kipl.tistory.com

https://kipl.tistory.com/58

 

Object Orientation

영상에서 전경 물체가 어떤 방향으로 정렬이 되어있는가를 찾는 문제는 다양한 영상 인식 알고리즘에서 나타난다. 예를 들면, 영상에서 사람의 머리가 어떤 자세를 취하고 있는가를 묻는 것에

kipl.tistory.com

 

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

Image Moments  (0) 2021.12.04
Orientation 추정  (0) 2021.11.30
Poisson Image Editing  (0) 2021.08.03
Sampling Theorem  (0) 2021.05.12
Lanczos Resampling  (0) 2021.05.08
Interpolation Kernels  (0) 2021.05.05
Posted by helloktk

댓글을 달아 주세요

728x90

주어진 이미지(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)$로 쓰인다:

출처: Poisson Image Editing J. Matías Di Martino, Gabriele Facciolo, Enric Meinhardt-Llopis

 

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

Image Moments  (0) 2021.12.04
Orientation 추정  (0) 2021.11.30
Poisson Image Editing  (0) 2021.08.03
Sampling Theorem  (0) 2021.05.12
Lanczos Resampling  (0) 2021.05.08
Interpolation Kernels  (0) 2021.05.05
Posted by helloktk

댓글을 달아 주세요

728x90

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) .$$

 

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

Orientation 추정  (0) 2021.11.30
Poisson Image Editing  (0) 2021.08.03
Sampling Theorem  (0) 2021.05.12
Lanczos Resampling  (0) 2021.05.08
Interpolation Kernels  (0) 2021.05.05
Fowler Angle  (0) 2021.04.05
Posted by helloktk

댓글을 달아 주세요

728x90

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)Si(2\pi-3\omega)+(4\pi -3\omega)Si(4\pi-3\omega)\\+(2\pi +3\omega)Si(2\pi+3\omega)+(4\pi+3\omega)Si(4\pi+3\omega)\Big)$$

// 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;
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
        CLAMP(xx, 0, srcWidth - 1);
        s += weight * src_line[xx];
        weightSum += weight;
    }
    return s / weightSum; 
}

 

bicubic downsample(1/2): alias 발생

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

Poisson Image Editing  (0) 2021.08.03
Sampling Theorem  (0) 2021.05.12
Lanczos Resampling  (0) 2021.05.08
Interpolation Kernels  (0) 2021.05.05
Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Posted by helloktk

댓글을 달아 주세요

728x90

Interpolation은 이산적으로 주어진 샘플들 사이에서 존재하지 않는 값을 추정하는 작업이다.  한 지점에서 interpolation을 통해 값을 추정하기 위해서는 주변의 알려진 샘플 값들을 참조해야 하고 적절한 가중치를 주어야 한다. 또한 interpolation 함수는 충분히 부드럽게 주변 샘플 값들과 연결이 되어야 한다. 이런 관점에서 보면 interpolation은 주변 샘플에 적절한 가중치를 주는 kernel함수와 convolution으로 이해할 수 있고, 또한 샘플링 과정에서 잃어버린 정보를 샘플 데이터를 smoothing해서 복원하는 과정으로 해석할 수 있다.

kernel 함수가 $K(x)$고, 일정한 간격으로 주어진 1차원 샘플의 값이 $\{ (x_k, c_k)\}$일 때 interpolation 함수는 

$$ f(x)  =\sum_k c_k K \Big( \frac{x - x_k}{h}\Big)$$

의 형태로 표현할 수 있다. $h$는 샘플의 간격으로 이미지에서는 픽셀의 간격이므로 $h=1$이다.

interpolation함수는 샘플 위치를 통과해야 하므로, 

$$ f(x_j) = c_j = \sum_{k}  c_k K( x_j - x_k)  \quad \longrightarrow\quad K(x_j - x_k ) = \delta_{jk}$$

즉, $K(0)=1$이고, $K(i)=0~(i\ne 0)$임을 알 수 있다. 또한 kernel이 주는 가중치의 합이 1이어야 하므로

$$\int_{-\infty}^{\infty}  K(x) dx = 1$$이어야 한다.

영상처리에서 잘 알려진 interpolation 방법은 nearest-neighbor, linear, cubic interpolation이 많이 사용된다. 아래 그림은 kernel의 중심이 원점에 있는 경우를 그렸다.

Nearest neighbor:

$$ K(x) = \left\{ \begin{array}{ll} 1, & |x| < \frac{1}{2} \\ 0, & |x| \ge\frac{1}{2}\end{array} \right. ,  \quad\quad  H(\omega) =\int_{-\infty}^{\infty}  K(x) e^{-i\omega x}dx= \text{sinc} \Big( \frac{\omega}{2} \Big) $$

Linear interpolation: 

$$ K(x) = \left\{ \begin{array}{ll} 1 - |x| , & |x| < 1 \\ 0, & |x| \ge 1  \end{array} \right. ,  \quad\quad  H(\omega) = \text{sinc}^2 \Big(\frac{\omega}{2} \Big) $$

Cubic interpolation:

$$ K(x) = \left\{ \begin{array}{ll} \frac{1}{2}( 3|x|^3 - 5|x|^2 +2)  , & |x| < 1 \\ \frac{1}{2}(-|x|^3 +5 |x|^2 -8|x|+4 ) , & 1\le |x| <2   \\ 0 , & |x| \le |x| \end{array} \right. , \\ \quad\quad  H(\omega) = \frac{12}{\omega^2} \Big( \text{sinc}^2 \Big(\frac{\omega}{2}\Big) -\text{sinc}(\omega) \Big) -\frac{4}{ \omega^2} \Big( 2\text{sinc}^2 (\omega) -2 \text{sinc}(\omega) -\text{sinc}(2\omega) \Big) $$

Cubic interpolation kernel이 다른 nearest-neighbor, linear kernel 보다 low-pass 특성이 강해지므로 이미지가 더 잘 smoothing 되는 특성을 가질 것임을 알 수 있다.

double KernelCubic(double x) {
    double absx = fabs(x);
    if (absx < 1)  return 0.5 * (2 + absx * (absx * (3 * absx - 5)))
    if (absx < 2)  return 0.5 * (4 - absx * (absx * (absx + 5) - 8));
    return 0;
}

일반적인 cubic convolution kernel은 4개의 샘플데이터를 이용해 보간하므로 폭이 4(반지름 =2)이다. $x_k=-1,0,1,2$에서 4개의 샘플값 $c_0, c_1, c_2,c_3$이 주어진 경우, $0\le x <1$구간에서 cubic interpolation은  kernel의 중심을 원점에서 $x$로 평행이동하면,

$$f(x) = c_0 K(-1-x) + c_1 K(-x) + c_2 K(1-x) + c_3 K(2-x)$$

로 표현된다. 다음 예는 $0\le  x<1$ 구간에서 $(x-1)^2$(blue)을 cubic interpolation(red)을 한 결과이다.

 

 

일반적인 cubic  convolution kernel은 한 개의 모양을 조정하는 한 개의 조정 parameter을 가지는 연속이고 한 번 미분가능한 형태로 주어진다.

$$ K(x) = \left\{   \begin{array}{ll} (a+2) |x|^3 - (a+3) |x| ^2 +1, &   |x| < 1 \\ a|x|^3 - 5a |x|^2 + 8a |x| - 4a ,&  1 \le |x| <2 \\ 0 ,& |x| \ge 2 \end{array} \right. $$

실용적인 파라미터의 범위는 $[-3,0]$이고( $a<-3$이면 $|x|<1$ 구간에서 단조감소가 아니어서 0이 아닌 곳에서 1보다 큰 값을 갖게 된다), 보통 많이 사용하는 cubic interpolation은 $a=-0.5$에 해당한다. 일반적으로 $a$가 -3에 가까워지면 최소 위치가 더 깊어지므로 일종의 미분 연산자처럼 작용하게 되어 이미지에 적용할 때 sharpness가 증가하게 된다. $a$가 0에 가까워지면 최소 위치 깊이가 0에 가까워지므로 gaussian filter처럼 이미지를 blurring하는 효과가 발생한다.

Lanczos Interpolation: 

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

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

https://kipl.tistory.com/277

 

Fixed-Point Bicubic Interpolation

아날로그 신호로부터 디지털 데이터로 바꿀 때 보통 시간당(sound) 또는 공간적 거리당(image) 일정한 횟수로 데이터를 수집한다. 이 비율을 sampling rate이라 한다. 그런데 sampling rate을 바꾸어야 할

kipl.tistory.com

kipl.tistory.com/55

 

Bicubic Interpolation

이미지 처리에서 픽셀 좌표는 간격이 1인 2차원 그리드의 교차점으로 볼 수 있다. 이미지를 확대하거나 축소할 때, 픽셀과 픽셀 사이 구간에서 값이 필요한 경우가 생긴다. 간단히는 주변의 가장

kipl.tistory.com

 

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

Sampling Theorem  (0) 2021.05.12
Lanczos Resampling  (0) 2021.05.08
Interpolation Kernels  (0) 2021.05.05
Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Posted by helloktk

댓글을 달아 주세요

728x90

/*
** 두 점 A, B가 주어졌을 때 벡터 (B-A)가 x-축과 이루는 각을 [0,8)사이로 보내는  
** 간단한 연산의 단조증가함수를 만들어 실제의 각을 대신하도록 한다. 
** 실제 각을 다른 계산에 넣어 사용하지 않고 비교만 할 때 매우 유용하다.
** (0,  1,  2,   3,   4,   5,   6,   7,   8)은 각각 실제 각도
** (0, 45, 90, 135, 180, 225, 270, 315, 360)에 해당한다.
*/
double FowlerAngle(CPoint A, CPoint B) {
    return FowlerAngle(B.x - A.x, B.y - A.y);
}
double FowlerAngle(double dx, double dy) {
    double adx = dx < 0 ? -dx: dx;
    double ady = dy < 0 ? -dy: dy;
    int code = (adx < ady) ? 1: 0;
    if (dx < 0)  code += 2;
    if (dy < 0)  code += 4;
    switch (code) {
    case 0: if (dx == 0) return 0;   /*     dx = dy = 0     */
            else return ady / adx;   /*   0 <= angle <=  45 */
    case 1: return 2 - adx / ady;    /*  45 <  angle <=  90 */
    case 2: return 4 - ady / adx;    /* 135 <= angle <= 180 */
    case 3: return 2 + adx / ady;    /*  90 <  angle <  135 */
    case 4: return 8 - ady / adx;    /* 315 <= angle <  360 */
    case 5: return 6 + adx / ady;    /* 270 <= angle <  315 */
    case 6: return 4 + ady / adx;    /* 180 <  angle <= 225 */
    case 7: return 6 - adx / ady;    /* 225 <  angle <  270 */
    }
};

 

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

Lanczos Resampling  (0) 2021.05.08
Interpolation Kernels  (0) 2021.05.05
Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Posted by helloktk

댓글을 달아 주세요

728x90

각 전경 픽셀에 대해서 모든 배경 픽셀까지 거리를 구한 다음 최솟값을 할당하면 된다. 픽셀 수가 n = w * h 개면 time complexity는 O(n^2)이다 (이미지 폭이나 높이의 4승이므로 연산량이 상당하다.) linear time Euclidean distance transform은 kipl.tistory.com/268.

ListPlot3D[ Reverse@ImageData[pic], PlotRange -> Full]

mathematica를 이용한 3D plot

int brute_force_edt(double *img, int w, int h) {
    int n = w * h;
    int *list = new int [n];
    int fg_cnt = 0;
    int bg_ptr = n - 1;               //배경 픽셀 위치는 끝에서 역으로 채움;
    //전경과 배경 픽셀 위치를 분리;
    for (int i = 0; i < n; ++i)
        if (img[i]) {                 // foreground;
            list[fg_cnt++] = i;
            img[i] = DBL_MAX;
        } else                       // background;
            list[bg_ptr--] = i;

    for (int i = 0; i < fg_cnt; ++i) {     // 전경 픽셀;
        int xi = list[i] % w, yi = list[i] / w;
        for (int j = fg_cnt; j < n; ++j) { // 배경 픽셀까지 거리를 계산해서 최소값을 할당;
            int xj = list[j] % w, yj = list[j] / w;
            double dx  = xi - xj, dy = yi - yj;  //
            double dst = dx * dx + dy * dy;
            if (img[list[i]] > dst)
                img[list[i]] = dst;
        }
    }
    delete [] list;
    return fg_cnt;
}

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

Interpolation Kernels  (0) 2021.05.05
Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Posted by helloktk

댓글을 달아 주세요

728x90

평균이 $\lambda$인 Poisson 분포를 가지는 random number를 만드는 레시피는 [0,1]에서 추출된 균일한 random number 수열이 $u_1, u_2, u_3,...$로 주어질 때, $u_1  u_2 u_3....u_{k+1}\le  e^{-\lambda}$을 처음 만족하는 $k$를 찾으면 됨이 Knuth에 의해서 알려졌다.

int PoissonValue(double val) { 
    double L = exp(-val);
    int k = 0;
    double p = 1;
    do {
        k++;
        // generate uniform random number u in [0,1] and let p ← p × u.
        p *= double(rand()) / RAND_MAX;
    } while (p > L);
    return (k - 1);
}

이미지에서 각각의 픽셀 값은 영상신호를 받아들이는 ccd 센서에서 노이즈가 포함된 전기신호의 평균값이 구현된 것으로 볼 수 있다. 따라서 영상에 Poisson 노이즈를 추가하는 과정은 역으로 주어진 픽셀 값을 평균으로 가지는 가능한 분포를 찾으면 된다.(물론 이미지에 무관하게 노이즈를 추가할 수도 있다. 예를 들면 모든 픽셀에 대해서 같은 평균값을 갖는 포아송 노이즈를 주는 경우다)

void AddPoissonNoise(CRaster& raster, CRaster& noised) {
    if (raster.IsEmpty() || raster.GetBPP() != 8) return;
    CSize sz = raster.GetSize();
    noised = raster;
    srand(unsigned(time(0)));
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)noised.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            int a = PoissonValue(*p++);
            *q++ = a > 0xFF ? 0xFF: a;
        }
    }
}

 

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

Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Posted by helloktk

댓글을 달아 주세요

728x90

blue cell: obstacles, black cell = start( or end) 

$$\text{distance measure}:~\text{Manhattan distance}=|x_2 - x_1 | + |y_2 - y_1|$$

Grassfire 알고리즘을 써서 최단경로를 찾는 문제는 DFS로도 구현할 수 있으나 비효율적이고(왜 그런지는 쉽게 알 수 있다) BFS로 구현하는 것이 더 적합하다. 

void grassfire(CPoint q, int **map, int w, int h, int **dist) {
    for (int y = 0; y < h; y++) 
        for (int x = 0; x < w; x++) 
            dist[y][x] = INF;  //unvisited cells;

    std::queue<CPoint> Q;
    dist[q.y][q.x] = 0;     //start( or end) position: distance = 0;
    Q.push(q);
    while (!Q.empty()) {
        CPoint p = Q.front(); Q.pop();
        int distance = dist[p.y][p.x];
        if (p.x > 0     && map[p.y][p.x - 1] == 0 && dist[p.y][p.x - 1] == INF) {
            dist[p.y][p.x - 1] = distance + 1;
            Q.push(CPoint(p.x - 1, p.y)); 
        }
        if (p.x < w - 1 && map[p.y][p.x + 1] == 0 && dist[p.y][p.x + 1] == INF) {
            dist[p.y][p.x + 1] = distance + 1;
            Q.push(CPoint(p.x + 1, p.y));
        }
        if (p.y > 0     && map[p.y - 1][p.x] == 0 && dist[p.y - 1][p.x] == INF) {
            dist[p.y - 1][p.x] = distance + 1;
            Q.push(CPoint(p.x, p.y - 1));
        }
        if (p.y < h - 1 && map[p.y + 1][p.x] == 0 && dist[p.y + 1][p.x] == INF) {
            dist[p.y + 1][p.x] = distance + 1;
            Q.push(CPoint(p.x, p.y + 1));
        }
    }
}

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

Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Posted by helloktk

댓글을 달아 주세요

728x90
#define PIX_SORT(a,b) { if ((a) > (b)) PIX_SWAP((a), (b));} 
#define PIX_SWAP(a,b) { BYTE tmp = (a); (a) = (b); (b) = tmp;} 
BYTE opt_med9(BYTE p[9]) { 
    PIX_SORT(p[1], p[2]); PIX_SORT(p[4], p[5]); PIX_SORT(p[7], p[8]); 
    PIX_SORT(p[0], p[1]); PIX_SORT(p[3], p[4]); PIX_SORT(p[6], p[7]); 
    PIX_SORT(p[1], p[2]); PIX_SORT(p[4], p[5]); PIX_SORT(p[7], p[8]); 
    PIX_SORT(p[0], p[3]); PIX_SORT(p[5], p[8]); PIX_SORT(p[4], p[7]); 
    PIX_SORT(p[3], p[6]); PIX_SORT(p[1], p[4]); PIX_SORT(p[2], p[5]); 
    PIX_SORT(p[4], p[7]); PIX_SORT(p[4], p[2]); PIX_SORT(p[6], p[4]); 
    PIX_SORT(p[4], p[2]); return(p[4]); 
}
double ImageSharpness(BYTE *img, int w, int h) {
    const int npixs = w * h;
    const int xend = w - 1, yend = h - 1;
    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};
    //median filter;
    BYTE arr[9];
    BYTE *medimg = new BYTE [npixs];
    memset(medimg, 0, npixs);
    for (int y = 1, pos = y * w; y < yend; y++) {
        pos++;
        for (int x = 1; x < xend; x++, pos++) {
            for (int k = 0; k < 9; k++) arr[k] = img[pos + nn[k]];
            medimg[pos] = opt_med9(arr);
        }
        pos++;
    }
    // Tenenbaum gradient;
    double sharpness = 0;
    for (int y = 1, pos = y * w; y < yend; y++) {
        pos++;
        for (int x = 1; x < xend; x++, pos++) {
            double gx = 0, gy = 0;
            for (int k = 0; k < 9; k++) {
                int v = medimg[pos + nn[k]];
                gx += sobelX[k] * v;
                gy += sobelY[k] * v;
            }
            sharpness += gx * gx + gy * gy;
        }
        pos++;
    }
    delete [] medimg;
    return sharpness;
}

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

Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Insertion Sort  (0) 2021.02.24
Posted by helloktk

댓글을 달아 주세요

728x90

void selection_sort(BYTE data[], int n) {
    for (int i = 0; i < n - 1; i++) {
        int jmin = i;
        for (int j = i + 1; j < n; j++)
            if (data[j] < data[jmin]) jmin = j;
        if (jmin != i) SWAP(data[jmin], data[i]);
    }
}

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

Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Insertion Sort  (0) 2021.02.24
Optimized Median Search  (0) 2021.02.24
Posted by helloktk

댓글을 달아 주세요

728x90

#define SWAP(a, b) {BYTE tmp = (a); (a) = (b); (b) = tmp;}
void  bubble_sort(BYTE data[], int n) {
    for (int i = n - 1; i > 0; i--) {
        for (int j = 0; j < i; j++) {
            if (data[j] > data[j + 1]) 
                SWAP(data[j], data[j + 1]);
        }
    }
}

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

Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Insertion Sort  (0) 2021.02.24
Optimized Median Search  (0) 2021.02.24
Zhang-Suen Thinning Algorithm  (0) 2021.02.18
Posted by helloktk

댓글을 달아 주세요

728x90

void insertion_sort ( BYTE *data, int n ) {
    int j;
    for (int i = 1; i < n; i++ ) {
        BYTE remember = data[(j = i)];
        while ( --j >= 0 && remember < data[j] ) {
            data[j + 1] = data[j];
            data[j] = remember;
        }
    }
}

 

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

Selection Sort  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Insertion Sort  (0) 2021.02.24
Optimized Median Search  (0) 2021.02.24
Zhang-Suen Thinning Algorithm  (0) 2021.02.18
Is Power of 2  (0) 2021.02.12
Posted by helloktk

댓글을 달아 주세요

728x90

참고:  ndevilla.free.fr/median/median/src/optmed_bench.c

/*
 * Optimized median search on 9 values
 */
#define PIX_SORT(a,b) { if ((a) > (b)) PIX_SWAP((a), (b)); }
#define PIX_SWAP(a,b) { BYTE temp = (a); (a) = (b); (b) = temp; }
BYTE opt_med9(BYTE p[9]) {
    PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ; 
    PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[4]) ; PIX_SORT(p[6], p[7]) ; 
    PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ; 
    PIX_SORT(p[0], p[3]) ; PIX_SORT(p[5], p[8]) ; PIX_SORT(p[4], p[7]) ; 
    PIX_SORT(p[3], p[6]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[2], p[5]) ; 
    PIX_SORT(p[4], p[7]) ; PIX_SORT(p[4], p[2]) ; PIX_SORT(p[6], p[4]) ; 
    PIX_SORT(p[4], p[2]) ; return(p[4]) ;
}
BYTE opt_med25(BYTE p[]) { 
    PIX_SORT(p[0], p[1]) ;   PIX_SORT(p[3], p[4]) ;   PIX_SORT(p[2], p[4]) ;
    PIX_SORT(p[2], p[3]) ;   PIX_SORT(p[6], p[7]) ;   PIX_SORT(p[5], p[7]) ;
    PIX_SORT(p[5], p[6]) ;   PIX_SORT(p[9], p[10]) ;  PIX_SORT(p[8], p[10]) ;
    PIX_SORT(p[8], p[9]) ;   PIX_SORT(p[12], p[13]) ; PIX_SORT(p[11], p[13]) ;
    PIX_SORT(p[11], p[12]) ; PIX_SORT(p[15], p[16]) ; PIX_SORT(p[14], p[16]) ;
    PIX_SORT(p[14], p[15]) ; PIX_SORT(p[18], p[19]) ; PIX_SORT(p[17], p[19]) ;
    PIX_SORT(p[17], p[18]) ; PIX_SORT(p[21], p[22]) ; PIX_SORT(p[20], p[22]) ;
    PIX_SORT(p[20], p[21]) ; PIX_SORT(p[23], p[24]) ; PIX_SORT(p[2], p[5]) ;
    PIX_SORT(p[3], p[6]) ;   PIX_SORT(p[0], p[6]) ;   PIX_SORT(p[0], p[3]) ;
    PIX_SORT(p[4], p[7]) ;   PIX_SORT(p[1], p[7]) ;   PIX_SORT(p[1], p[4]) ;
    PIX_SORT(p[11], p[14]) ; PIX_SORT(p[8], p[14]) ;  PIX_SORT(p[8], p[11]) ;
    PIX_SORT(p[12], p[15]) ; PIX_SORT(p[9], p[15]) ;  PIX_SORT(p[9], p[12]) ;
    PIX_SORT(p[13], p[16]) ; PIX_SORT(p[10], p[16]) ; PIX_SORT(p[10], p[13]) ;
    PIX_SORT(p[20], p[23]) ; PIX_SORT(p[17], p[23]) ; PIX_SORT(p[17], p[20]) ;
    PIX_SORT(p[21], p[24]) ; PIX_SORT(p[18], p[24]) ; PIX_SORT(p[18], p[21]) ;
    PIX_SORT(p[19], p[22]) ; PIX_SORT(p[8], p[17]) ;  PIX_SORT(p[9], p[18]) ;
    PIX_SORT(p[0], p[18]) ;  PIX_SORT(p[0], p[9]) ;   PIX_SORT(p[10], p[19]) ;
    PIX_SORT(p[1], p[19]) ;  PIX_SORT(p[1], p[10]) ;  PIX_SORT(p[11], p[20]) ;
    PIX_SORT(p[2], p[20]) ;  PIX_SORT(p[2], p[11]) ;  PIX_SORT(p[12], p[21]) ;
    PIX_SORT(p[3], p[21]) ;  PIX_SORT(p[3], p[12]) ;  PIX_SORT(p[13], p[22]) ;
    PIX_SORT(p[4], p[22]) ;  PIX_SORT(p[4], p[13]) ;  PIX_SORT(p[14], p[23]) ;
    PIX_SORT(p[5], p[23]) ;  PIX_SORT(p[5], p[14]) ;  PIX_SORT(p[15], p[24]) ;
    PIX_SORT(p[6], p[24]) ;  PIX_SORT(p[6], p[15]) ;  PIX_SORT(p[7], p[16]) ;
    PIX_SORT(p[7], p[19]) ;  PIX_SORT(p[13], p[21]) ; PIX_SORT(p[15], p[23]) ;
    PIX_SORT(p[7], p[13]) ;  PIX_SORT(p[7], p[15]) ;  PIX_SORT(p[1], p[9]) ;
    PIX_SORT(p[3], p[11]) ;  PIX_SORT(p[5], p[17]) ;  PIX_SORT(p[11], p[17]) ;
    PIX_SORT(p[9], p[17]) ;  PIX_SORT(p[4], p[10]) ;  PIX_SORT(p[6], p[12]) ;
    PIX_SORT(p[7], p[14]) ;  PIX_SORT(p[4], p[6]) ;   PIX_SORT(p[4], p[7]) ;
    PIX_SORT(p[12], p[14]) ; PIX_SORT(p[10], p[14]) ; PIX_SORT(p[6], p[7]) ;
    PIX_SORT(p[10], p[12]) ; PIX_SORT(p[6], p[10]) ;  PIX_SORT(p[6], p[17]) ;
    PIX_SORT(p[12], p[17]) ; PIX_SORT(p[7], p[17]) ;  PIX_SORT(p[7], p[10]) ;
    PIX_SORT(p[12], p[18]) ; PIX_SORT(p[7], p[12]) ;  PIX_SORT(p[10], p[18]) ;
    PIX_SORT(p[12], p[20]) ; PIX_SORT(p[10], p[20]) ; PIX_SORT(p[10], p[12]) ;
    return (p[12]);
}
/*---------------------------------------------------------------------------
   Function :   kth_smallest()
   In       :   array of elements, # of elements in the array, rank k
   Out      :   one element
   Job      :   find the kth smallest element in the array
   Notice   :   use the median() macro defined below to get the median. 
                Reference:
                  Author: Wirth, Niklaus 
                   Title: Algorithms + data structures = programs 
               Publisher: Englewood Cliffs: Prentice-Hall, 1976 
    Physical description: 366 p. 
                  Series: Prentice-Hall Series in Automatic Computation 
 ---------------------------------------------------------------------------*/
#define SWAP(a, b) {BYTE t = (a); (a) = (b); (b) = t; }
BYTE kth_smallest(BYTE a[], int n, int k) {
    int l = 0, m = n - 1 ;
    while (l < m) {
        BYTE x = a[k] ;
        int i = l, j = m ;
        do {
            while (a[i] < x) i++ ;
            while (x < a[j]) j-- ;
            if (i <= j) {
                SWAP(a[i], a[j]) ;
                i++ ; j-- ;
            }
        } while (i <= j) ;
        if (j < k) l = i ;
        if (k < i) m = j ;
    }
    return a[k] ;
}
#define median(a, n) kth_smallest(a, n, (((n) & 1) ? ((n) / 2) : (((n) / 2) - 1)))

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

Bubble Sort  (0) 2021.02.24
Insertion Sort  (0) 2021.02.24
Optimized Median Search  (0) 2021.02.24
Zhang-Suen Thinning Algorithm  (0) 2021.02.18
Is Power of 2  (0) 2021.02.12
Flood-Fill Algorithm을 이용한 Connected Component Labeling  (1) 2021.02.10
Posted by helloktk

댓글을 달아 주세요

728x90

더보기
비교: Rosenfeld Algorithm(Graphics Gem IV)
int Thinning_2pass(BYTE *image, int w, int h) {
    const int xmax = w - 1, ymax = h - 1;
    const int nn[9] = {-w - 1,- w, -w + 1, 1, w + 1, w, w - 1, -1, 0};//clockwise;
    const BYTE FG = 255, BG = 0;
    bool *flag = new bool [w * h];
    int pass = 0, ok = 0;
    int nb[9];
    while (!ok) {
        ok = 1; pass = (pass + 1) % 2;
        for (int y = 1, pos = w; y < ymax; y++) {
            pos++; //x=0;
            for (int x = 1; x < xmax; x++, pos++) flag[pos] = false;
            pos++; //x=w-1;
        }
        for (int y = 1, pos = w; y < ymax; y++) {
            pos++;//x=0;
            for (int x = 1; x < xmax; x++, pos++) {
                if (image[pos] == FG) { //fg;
                    int count = 0;
                    for (int k = 0; k < 9; k++) 
                        if (image[pos + nn[k]] == FG) count++;
                    if (count > 2 && count < 8) {
                        for (int k = 0; k < 8; k++) nb[k] = image[pos + nn[k]];
                        nb[8] = nb[0]; //cyclic;
                        int trans = 0;
                        for (int k = 0; k < 8; k++)
                            if (nb[k] == BG && nb[k + 1] == FG) trans++;
                        if (trans == 1) {
                            if (pass == 0 && (nb[3] == BG || nb[5] == BG ||	
                                (nb[1] == BG && nb[7] == BG))) {
                                    flag[pos] = true; ok = 0;
                            } else {
                                if (pass == 1 && (nb[1] == BG || nb[7] == BG ||	
                                    (nb[3] == BG && nb[5] == BG))) {
                                        flag[pos] = true; ok = 0;
                                }
                            }
                        }
                    }//(2 <count<8);
                }
            }//for-x;
            pos++;//x = w - 1 skip;
        } //for-y;
        // remove flaged pixels;
        for (int y = 1, pos = w; y < ymax; y++) {
            pos++;//x = 0;
            for (int x = 1; x < xmax; x++, pos++) 
                if (flag[pos]) image[pos] = BG;
            pos++; //x=w-1;
        }
    }
    delete[] flag;
    return 1;
}

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

Insertion Sort  (0) 2021.02.24
Optimized Median Search  (0) 2021.02.24
Zhang-Suen Thinning Algorithm  (0) 2021.02.18
Is Power of 2  (0) 2021.02.12
Flood-Fill Algorithm을 이용한 Connected Component Labeling  (1) 2021.02.10
Gradient Covariance  (0) 2021.01.27
Posted by helloktk

댓글을 달아 주세요

728x90
bool IsPowerOf2(int n) {
    return n && !(n & (n - 1));
    // return n > 0 && (n & (n - 1)) == 0;
}
bool IsPowerOf2(int i) {
    if (i == 1) return true;
    if (i == 0 || i & 1) return false;
    return IsPowerOf2(i >> 1);
} // i == 0; 
bool IsPowerof2(int x){
    int i = 0;
    while ((1 << i) < x) i++;
    if (x == (1 << i)) return true;
    return false;
}
bool IsPowerOf2(int n) {
   if (n == 0) return false;
   return (ceil(log2(n)) == floor(log2(n)));
}
int NextPowerOf2(int n) { //32-bit;
    n--;
    n |= n >> 1;
    n |= n >> 2;
    n |= n >> 4;
    n |= n >> 8;
    n |= n >> 16;
    n++;
    return n;    
} // NextPowerOf2(5) -> 8; NextPowerOf2(8) -> 8;

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

Optimized Median Search  (0) 2021.02.24
Zhang-Suen Thinning Algorithm  (0) 2021.02.18
Is Power of 2  (0) 2021.02.12
Flood-Fill Algorithm을 이용한 Connected Component Labeling  (1) 2021.02.10
Gradient Covariance  (0) 2021.01.27
점증적인 cosine/sine 값 계산  (0) 2020.12.28
Posted by helloktk

댓글을 달아 주세요

728x90

 

// labeling using depth-first search; non-recursive version;
int GetConnectedComponents(BYTE *image, int w, int h, int *table) {
    int npixs = w * h;
    int label = 0;    // starting label = 1;
    int *stack = new int [npixs];
    // initialize the table;
    for (int pos = 0; pos < npixs; pos++) 
        table[pos] = image[pos] ? -1: 0;  // Foreground = -1; Background = 0;

    for (int pos = 0; pos < npixs; pos++) {
        if (table[pos] == -1) { // Foreground;
            ++label;            // assign next label;
            int top = -1;       // stack initialization;
            stack[++top] = pos;
            while (top >= 0) {
                int adj = stack[top--];
                int xx = adj % w;
                int yy = adj / w;
                if (table[adj] == -1) {// Foreground;
                    table[adj] = label;
                    // check 4-way connectivity;
                    if (xx + 1 < w) stack[++top] = adj + 1; //RIGHT;
                    if (yy + 1 < h) stack[++top] = adj + w; //BOTTOM;
                    if (yy > 0)     stack[++top] = adj - w; //TOP;
                    if (xx > 0)     stack[++top] = adj - 1; //LEFT;
                }
            }
        }
    }
    delete [] stack;
    return label;    // total # of CCs;
};

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

Zhang-Suen Thinning Algorithm  (0) 2021.02.18
Is Power of 2  (0) 2021.02.12
Flood-Fill Algorithm을 이용한 Connected Component Labeling  (1) 2021.02.10
Gradient Covariance  (0) 2021.01.27
점증적인 cosine/sine 값 계산  (0) 2020.12.28
Fast Float Sqrt  (0) 2020.12.27
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2021.02.10 22:30 신고  댓글주소  수정/삭제  댓글쓰기

    DSU랑 Flood Fill이랑 시간 같나요?

728x90

이미지 처리 과정에서 미분에 해당하는 그래디언트 필드(gradient field: $g_x$, $g_y$ )를 이용하면 이미지 상의 특징인 corner, edge, ridge 등의 정보를 쉽게 얻을 수 있다. 이미지의 한 지점이 이러한 특징을 가지는 특징점이 되기 위해서는 그래디언트 필드의 값이 그 점 주변에서 (3x3나 5x5정도 크기의 window) 일정한 패턴을 유지해야 한다. 이 패턴을 찾기 위해서 그래디언트 필드에 PCA를 적용해보자. 수평과 수직방향의 그래디언트 field인 $g_x$와 $g_y$ 사이의 covariance 행렬은 다음 식으로 정의된다:

$$ \Sigma = \left [ \begin {array}{cc} < g_x^2 > & < g_x g_y > \\    <g_x g_y> & <g_y^2 > \end {array}\right] =\left [ \begin {array}{cc} s_{xx} & s_{xy} \\ s_{xy} & s_{yy}\end {array}\right];$$

$<...> = \int_{W}(...) dxdy$는 픽셀 윈도에 대한 적분을 의미한다. $\Sigma$의 eigenvalue는 항상 음이 아닌 값을 갖게 되는데 (matrix 자체가 positive semi-definitive), 두 eigenvalue이 $λ_1$, $λ_2$면

$$λ_1 + λ_2 = s_{xx} + s_{yy} \ge 0, \quad \quad    λ_1  λ_2 = s_{xx} s_{yy} - s_{xy}^2 \ge0 $$

을 만족한다 (완전히 상수 이미지를 배제하면 0인 경우는 없다). eigenvalue $λ_1$, $λ_2$는 principal axis 방향으로 그래디언트 필드의 변동(분산)의 크기를 의미한다. edge나 ridge의 경우는 그 점 주변에서 잘 정의된 방향성을 가져야 하고, corner의 경우는 방향성이 없어야 한다. edge나 ridge처럼 일방향성의 그래디언트 특성을 갖거나 corner처럼 방향성이 없는 특성을 서로 구별할 수 있는 measure가 필요한데, $λ_1$과 $λ_2$를 이용하면 차원이 없는 measure을 만들 수 있다. 가장 간단한 차원이 없는 측도(dimensionless measure)는  eigenvalue의 기하평균과 산술평균의 비를 비교하는 것이다.

$$ Q = \frac { {λ_{1} λ_{2}} }{ \left( \frac {λ_{1}+λ_{2}}{2} \right)^2} = 4\frac { s_{xx} s_{yy} - s_{xy}^2}{(s_{xx} + s_{yy})^2};$$

기하평균은 산술평균보다도 항상 작으므로

$$ 0 \le Q \le 1 $$

의 범위를 갖는다. 그리고 $Q$의 complement로

$$P = 1-Q = \frac{(s_{xx}-s_{yy})^2 + 4 s_{xy}^2}{(s_{xx}+s_{yy})^2};$$를 정의할 수 있는 데 $0 \le P \le 1$이다. $Q$와 $P$의 의미는 무엇인가? 자세히 증명을 할 수 있지만 간단히 살펴보면 한 지점에서 $Q \rightarrow 1$이려면 $λ_{1} \approx λ_{2}$이어야 하고, 이는 두 주축이 동등하다는 의미이므로 그 점에서는 방향성이 없는 코너의 특성을 갖게 된다. 반대로 $Q \rightarrow 0$이면 강한 방향성을 갖게 되어서 edge(ridge) 특성을 갖게 된다.

 

실제적인 응용으로는 지문 인식에서 지문 영역을 알아내거나 (이 경우는 상당이 큰 윈도를 사용해야 한다) 또는 이미지 텍스쳐 특성을 파악하기 위해서는 이미지를 작은 블록으로 나누고 그 블록 내의 미분 연산자의 균일성을 파악할 필요가 있는데 이 차원이 없는 측도는 이미지의 상태에 상관없이 좋은 기준을 주게 된다.

 

참고 논문:

Image field categorization and edge/corner detection from gradient covariance
Ando, S.

Pattern Analysis and Machine Intelligence, IEEE Transactions on
Volume 22, Issue 2, Feb 2000 Page(s):179 - 190

 

** 네이버 블로그 이전;

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

Is Power of 2  (0) 2021.02.12
Flood-Fill Algorithm을 이용한 Connected Component Labeling  (1) 2021.02.10
Gradient Covariance  (0) 2021.01.27
점증적인 cosine/sine 값 계산  (0) 2020.12.28
Fast Float Sqrt  (0) 2020.12.27
2차원 Gaussian 분포 생성  (0) 2020.12.10
Posted by helloktk

댓글을 달아 주세요

728x90

이미지 처리에서 각도를 다루게 되면 삼각함수를 이용해야 한다. 일반적인 PC 상에서는 CPU의 성능이 크게 향상이 되어 있고, floating point 연산이 잘 지원이 되므로 큰 문제가 없이 삼각함수를 바로 이용할 수가 있다. 그러나 embeded 환경에서는 연산 속도가 그다지 빠르지 않고 floating point 연산은 에뮬레이션을 하는 경우가 대부분이다. 이런 조건에서 루프 내에서 연속적으로 삼각함수를 호출하는 것은 성능에 많은 영향을 주게 되므로 다른 방법을 고려해야 한다. 일반적으로 각도를 점진적으로 증가시키면 cosine이나 sine 값을 계산해야 할 필요가 있는 경우에는 아래와 같이 처리한다:

for (int i = 0; i < nsteps; ++i) {
    double angle = i * angle_increment;
    double cos_val = cos( angle ) ;
    double sin_val = sin( abgle ) ;
    // do something useful with cos_val, and sin_val;
    // .......................;
}

물론 미리 계산된 cosine/sine의 값을 가지고 처리할 수도 있다. 하지만 매우 작은 각도 단위로 계산하는 경우는 lookup table의 크기가 매우 커지게 되므로 직접 계산을 하는 것보다 이득이 없을 수 있다. 이런 경우에는 다음의  삼각함수의 관계식을 이용하면 매번 cosine이나 sine 함수의 호출 없이도 필요한 값을 얻을 수 있다;

$$\cos(\theta+\delta)-\cos(\theta)=- (2\sin^2(\delta/2) \cos(\theta) + \sin(\delta)\sin(\theta))= -(\alpha \cos(\theta) + \beta \sin (\theta));$$

$$\sin(\theta+\delta)-\sin(\theta)= -(2\sin^2(\delta/2) \sin(\theta) - \sin(\delta)\cos(\theta))= -(\alpha \sin(\theta) - \beta \cos (\theta));$$

여기서 $\alpha = 2 \sin^2(\delta/2)$, $\beta = \sin(\delta)$다. 처음에 한번 계산된 $\alpha$와 $\beta$ 값을 이용하면, 추가적인 계산이 없이 점진적으로 $\cos$과 $\sin$의 값을 계산할 수 있다.

    int i = 0 ;
    double cos_val = 1 ; // = cos(0);
    double sin_val = 0 ; // = sin(0);
    double beta = sin(angle_increment) ;
    double alpha = sin(angle_increment/2) ;
    alpha = 2 * alpha * alpha ; 
    do {
        // do something useful with cos_val and sin_val ;
        // ...........................;
        // calculate next values;

        double cos_inc = (alpha * cos_val + beta * sin_val) ;
        double sin_inc = (alpha * sin_val - beta * cos_val) ;
        cos_val -= cos_inc ;
        sin_val -= sin_inc ;
    } while (++i < nsteps);

이 방법은 이론상 정확한 방법이나 round-off error가 전파가 되어서 반복 회수가 커질수록 에러가 누적이 될 위험이 있다. 이를 방지하기 위해 중간에 카운터를 두어서 일정한 반복 이후에는 cosine/sine값 계산을 다시 하여 에러를 리셋시키면 된다.

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

Flood-Fill Algorithm을 이용한 Connected Component Labeling  (1) 2021.02.10
Gradient Covariance  (0) 2021.01.27
점증적인 cosine/sine 값 계산  (0) 2020.12.28
Fast Float Sqrt  (0) 2020.12.27
2차원 Gaussian 분포 생성  (0) 2020.12.10
PCA Line Fitting  (0) 2020.11.12
Posted by helloktk

댓글을 달아 주세요

728x90

 이미지 처리 연산 과정에서 제곱근 값을 계산해야 경우가 많이 생긴다. 이 경우 math-라이브러리의 sqrt() 함수를 사용하는 것보다 원하는 정밀도까지만 계산을 하는 근사적인 제곱근 함수를 만들어서 사용하는 것이 계산 비용 측면에서 유리할 수 있다. 

주어진 양수 $x$의 제곱근 $y = \sqrt{x}$를 구하는 것은 $F(y) = y^2 - x$의 근을 구하는 것과 같다. 프로그램적으로 근은 반복적인 방법에 의해서 구해지는데 많이 사용하는 알고리즘 중 하나가 Newton 방법이다. Newton 방법은 한 단계에서 근의 후보가 정해지면 그 값에서 F에 접하는 직선을 구한 후 그 직선이 가로축과 만나는 점을 다시 새로운 근으로 잡는 반복 과정을 일정한 에러 범위 내에 들어올 때까지 수행한다.

 

함수 $F(y)$의 근은 아래의 점화식에 의해서 구해지게 된다:

$$\tt y(n + 1) = y(n) - F(y(n)) / F'(y(n)), n = 0,1, 2,...;$$

제곱근의 경우에는 ($F(y) = y^2 - x$) 점화식은 

$$\tt y(n + 1) = (y(n) + x / y(n)) / 2;$$

의 형태로 주어진다. 그러나, 이 점화식에는 계산 부하가 상대적으로 큰 나누기 연산이 들어 있다. 그런데 $\sqrt{x} = x \times (1 / \sqrt{x})$로 표현할 수 있으므로 제곱근의 역수를 쉽게 계산하는 알고리즘을 찾으면 된다. $x$의 제곱근의 역수를 구하는 점화식은 (이 경우 $F(y) = 1 / y^2 - x$)

$$\tt y(n + 1) = y(n) * (3 - y(n) * y(n) * x) / 2;$$

이므로 전체적인 제곱근을 구하는 pseudo 코드는

float fast_sqrt(float x) {
    xhalf = 0.5F * x;
    y = initial_guess_of (1/sqrt(x));
    while (error is not small)
        // find approx inverse sqrt root of (x);
        y = y * 1.5 - y * y * y * xhalf;
    end while

    // finally return sqrt(x) = x * (1 / sqrt(x));
    return (x * y);
}

보통 iteration 횟수는 고정이 되어 있다. 남은 문제는 초기값 추정뿐이다. Newton 방법에서 잘 추정된 초기값은 몇 번의 반복으로도 원하는 정밀도를 얻을 수 있게 하지만 잘못된 추정은 수렴 속도를 더디게 한다. 또한 초기값의 추정 과정도 매우 간단해야 하는데, 좋은 추정을 하기 위해서는 IEEE에서 규정하는 double/float 포맷에 대한 정보와 수치해석 지식이 필요하다.

 

IEEE의 32비트 float 상수의 비트 구성은 $\tt 0 \sim 22$번째 비트까지는 소수점 아래의 유효숫자(Mantissa: $\tt M$)를 표기하고, $\tt 23$번째부터 $\tt 30$번째까지는 지수 부분(Exponent: $\tt E$)을 표시하는데 $\tt 127$만큼 더해진 값이다(+-지수를 표현하기 위한 것으로 실제 지수는 $\tt e = E - 127$). 마지막 $\tt 31$비트는 부호를 결정한다. 임의의 양의 float 상수 $\tt x$가 

$$\tt x = (1 + M) 2^{e} = (1 + M) 2^{E - 127};$$

로 쓰여지면, 정수로 변환된 비트 데이터는 $\tt ix = *(int~*)\&x = (E + M) 2^{23}$의 값을 갖는다. $\tt x$ 제곱근의 역수는

$$\tt \frac{1}{\sqrt{x}} = \frac{1}{\sqrt{1 + M}} 2^{- e / 2} ;$$

로 표현되므로 비트 데이터에서 지수부는 $\tt (- e / 2) + 127 = 190 - E / 2$로 표현된다. 

 

이제 필요한 것은 잘 선택된 $\tt 1/\sqrt{x}$의 초기값을 찾는 것으로, 첫 번째로 시도할 근사는 유효숫자 부분 $(\tt 1/\sqrt{1 + M})$을 제거하고, 단순히 지수 부분($\tt 2^{- e / 2}$)만 살리면 된다정수 표현으로 바꾸면 $\tt 190$이 $\tt 23$에서 $\tt 30$번째 비트를 채우도록 $\tt 23$을 시프트 하고, $\tt - E / 2$ 부분은 원래 float 상수의 정수 표현을 나누면 유효숫자 부분에 해당하는 정도만 에러가 발생한다:

   정수 표현 = (190 - E / 2) << 23;

                  = (190 << 23) - (x의 정수형 변환 & (0xFF << 23)) >> 1;

                  ~ (190 << 23) - (x의 정수형 변환) >> 1;    // 유효 숫자 / 2 만큼 에러 발생;

                  = 0x5F000000 - (*(int *)&x) >> 1;

유효숫자 부분은 좀 더 높은 초기 정밀도를 갖도록 근사를 할 수 있는데 결과만 쓰면 $\tt 0x5F000000$ 대신에 $\tt 0x5F3759DF$를 이용하면 된다.

//빠른 float 제곱근의 역수;
float fast_invsqrt(float x) {
    float xhalf = 0.5F * x;
    // initial guess of 1 / sqrt(x);
    int i = 0x5F3759DF - (*(int *)&x) >> 1;
    x = *(float *)&i;    
    //iterate
    x = x * (1.5F - x * x * xhalf);
    x = x * (1.5F - x * x * xhalf);
    // end of calculation of 1 / sqrt(x);
    return (x);
};
// 빠른 float 제곱근;
float fast_sqrt(float x) {
    float xhalf = 0.5F * x;
    //initial guess of 1 / sqrt(x);
    int i = 0x5F3759DF - (*(int *)&x) >> 1;
    float y = *(float *)&i;    
    //iterate
    y = y * (1.5F - y * y * xhalf);
    y = y * (1.5F - y * y * xhalf);
    //end of calculation of 1/sqrt(x);
    //return sqrt(x)=x * 1/sqrt(x);
    return (x * y)   
}; 

$ \log_2 (x) = e + \log_2 (1 + M)$이므로 x의 근삿값은 $\log_2 (1+M)$의 근삿값을 잘 선택하면 된다. 함수 $f(M)=\log_2 (1+M) - M$는 $M=1/\ln(2)-1$에서 최댓값 $\tt 0.0860713$을 가지므로

$$0\le \log_2 (1+M)-M \le 0.0860713$$

이다. 따라서 오차를 최댓값의 절반으로 잡아서 근사하면

$$\log_2 (1+M) \approx M+ \delta,\quad \delta = \frac{0.0860713}{2}=0.0430385$$ 

그리고 $\delta$는 $x$에 무관한 값이 된다. $x$의 비트 데이터를 정수로 표현하면

$$ i_x= 2^{23} (E+M) = 2^{23}(e + \delta +M) + 2^{23}(127-\delta) \approx 2^{23}\log_2(x) + 2^{23}(127-\delta).$$ 로 근사할 수 있다. 두 번째 항은 $x$에 무관한 값이다. 따라서, 제곱근의 역수에 해당하는 비트 데이터의 근사적인 정수 표현은

$$i_{\frac{1}{\sqrt{x}}} \approx 2^{23}\log_2 (1/\sqrt{x}) + 2^{23}(127-\delta)\approx \frac{3}{2} 2^{23} (127-\delta ) -\frac{1}{2} i_x$$

로 쓸 수 있고, 위에서 구한 $\delta$을 대입하면

$$\tt  i_{\frac{1}{\sqrt{x}}} \approx  0x5F37BCB6 - i_x>>1$$ 

을 얻을 수 있다. $\delta = 0.0450466$을 대입하면 위에 적힌 코드에서 사용한 $\tt 0x5F3759DF$가 됨을 알 수 있다. 그러나 어떤 근거로 이 값을 선택했는지는...

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

Gradient Covariance  (0) 2021.01.27
점증적인 cosine/sine 값 계산  (0) 2020.12.28
Fast Float Sqrt  (0) 2020.12.27
2차원 Gaussian 분포 생성  (0) 2020.12.10
PCA Line Fitting  (0) 2020.11.12
Histogram Equalization  (0) 2020.11.12
Posted by helloktk

댓글을 달아 주세요

728x90

$$g(x, y; \sigma)= \frac{1}{2\pi \sigma^2}\exp\left( - \frac{x^2 +y^2 }{2\sigma^2 }\right)$$

원점을 중심으로 분산이 1 인 정규분포를 갖는 평면상의 점을 발생시킬 필요가 있는 경우에 아래의 함수를 이용한다: Box-Muller method

내장 함수 rand()를 이용하면 $[0,1]$ 구간에서 균일한 분포를 만들 수 있다. $[0,1]\times[0,1]$ 영역에서 균일한 분포를 갖는 두 개의 1차원 랜덤 변수 $x_1$, $x_2$을 이용해서 확률 변수 $y_1$, $y_2$의 2차원 가우시안 분포를 주는 변환을 구하려면, 우선 확률 보존에 의해서 

$$ dx_1 dx_2 = \frac{1}{\sqrt{2\pi}} e^{-y_1^2 /2 } dy_1 \frac{1}{\sqrt{2\pi}} e^{-y_2^2 /2 } dy_2 $$

을 만족시켜야 한다. 극좌표 $\rho^2 = y_1^2 + y_2^2$, $\phi = \tan^{-1}(y_2/y_1)$을 도입하면 

$$ dx_1 dx_2 = e^{-\rho^2/2} \rho d\rho \frac{d\phi}{2\pi}$$로 쓸 수 있으므로 미소 변환은

$$ dx_1 = e^{-\rho^2/2} \frac{d\rho^2}{2},\quad dx_2 = \frac{d\phi}{2\pi}$$ 로 선택할 수 있다. 따라서

$$ x_1 = \exp[-\rho^2/2] = \exp[- (y_1^2 + y_2^2 )/2],$$

$$ x_2 = \frac{1}{2\pi} \tan^{-1} (y_2/y_1).$$

그리고, 역변환은

$$ y_1 = \sqrt{-2 \ln (x_1) } \cos (2\pi x_2),\quad y_2 = \sqrt{-2 \ln(x_1) } \sin (2\pi x_2).$$

삼각함수의 계산을 피하기 위해서 반지름이 1인 원 내부에서 정의되는 두 개의 균일한 랜덤 변수 $v_1, v_2$를 도입하면, 다음 변환에 의해서 반지름 1인 원 영역이 변의 길이가 1인 정사각형 영역으로 보내진다.

$$x_1 = R^2(=S) = v_1^2 + v_2^2, \quad 2\pi x_2 = \tan^{-1} (v_2/v_1).$$

따라서 반지름이 1인 원 내부의 랜덤 변수를 써서 2차원의 가우시안 분포를 만들어주는 변환은

$$ y_1 = \sqrt{-2\ln S} \frac{v_1}{\sqrt{S}}, \quad y_2 = \sqrt{-2\ln S} \frac{v_2}{\sqrt{S}}$$

// 참고;  Numerical Receipe;
#define drand() ((double)rand() / RAND_MAX)   // [0,1]
// mean = 0; sigma = 1;
void normal_2D(double *x, double *y) {
    while (1) {
        double V1 = -1. + (2. * drand());
        double V2 = -1. + (2. * drand());
        double S = V1 * V1 + V2 * V2;
        if (S < 1.) {
            double rad = sqrt(-2.*log(S) / S);
            *x = V1 * rad;
            *y = V2 * rad;
            /* if (mean, sigma) was given,
            *x = meanx + sigma * (*x);
            *y = meany + sigma * (*y);
            */
            return;
        }
    }
};

 

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

점증적인 cosine/sine 값 계산  (0) 2020.12.28
Fast Float Sqrt  (0) 2020.12.27
2차원 Gaussian 분포 생성  (0) 2020.12.10
PCA Line Fitting  (0) 2020.11.12
Histogram Equalization  (0) 2020.11.12
Least Square Fitting of Circle  (0) 2020.11.11
Posted by helloktk

댓글을 달아 주세요

728x90

평면 위에 점집합이 주어지고 이들을 잘 기술하는 직선의 방정식을 구해야 할 경우가 많이 발생한다. 이미지의 에지 정보를 이용해 선분을 찾는 경우에 hough transform과 같은 알고리즘을 이용하는 할 수도 있지만 수치해석적으로 직접 fitting을 할 수도 있다. 점집합의 데이터를 취합하는 과정은 항상 노이즈에 노출이 되므로 직선 위의 점뿐만 아니라 직선에서 (많이) 벗어난 outlier들이 많이 들어온다. 따라서 line-fitting은 이러한 outlier에 대해서 매우 robust 해야 한다. 데이터 fitting의 경우에 초기에 대략적인 fitting에서 초기 파라미터를 세팅하고, 이것을 이용하여서 점차로 정밀하게 세팅을 해나가는 반복적인 방법을 많이 이용한다. 입력 데이터가 $\{(x_i, y_i)| i=0,..., N-1\}$로 주어지는 경우에 많이 이용하는 최소자승법에서는 각 $x_i$에서 직선 상의 $y$ 값과 주어진 $y_i$의 차이(residual)의 제곱을 최소로 하는 직선의 기울기와 $y$ 절편을 찾는다. 그러나 데이터가 $y$축에 평행하게 분포하는 경우를 다루지 못하게 되며, 데이터 점에서 직선까지 거리를 비교하는 것이 아니라 $y$값의 차이만 비교하므로 outlier의 영향을 매우 심하게 받는다.

이러한 문제를 제거 또는 완화하기 위해서는 PCA(principal axis analysis)를 이용할 수 있다. 점들이 선분을 구성하는 경우, 선분 방향으로는 점 위치의 편차가 크지만 수직 방향으로는 편차가 상대적으로 작다. 따라서 평면에서 점 분포에 대한 공분산 행렬 $\Sigma$의 고윳값과 고유 벡터를 구하면, 큰 고윳값을 갖는 고유 벡터 방향이 선분의 방향이 될 것이다. 잘 피팅이 이루어지려면 두 고윳값의 차이가 커야 한다. 또한 outlier에 robust 한 피팅이 되기 위해서는 각 점에 가중치를 부여해서 공분산 행렬에 기여하는 가중치를 다르게 하는 알고리즘을 구성해야 한다. 처음 방향을 설정할 때는 모든 점에 동일한 가중치를 부여하여 선분의 방향을 구한 후 다음번 계산에서는 직선에서 먼 점이 공분산 행렬에 기여하는 weight를 줄여 주는 식으로 하면 된다. weight는 점과 직선과의 거리에 의존하나 그 형태는 항상 정해진 것이 아니다.

 

// PCA-방법에 의한 line-fitting;
int LineFit_PCA(POINT P[], int N, double weight[], double line[4]){
    // 초기화 시 weight[i] = 1.;
    double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0, sw = 0;
    for (int i = 0; i < N; i++) {
         int x = P[i].x, y = P[i].y ;
         sx  += weight[i] * x; 
         sy  += weight[i] * y;
         sxx += weight[i] * x * x; 
         sxy += weight[i] * x * y; 
         syy += weight[i] * y * y;
         sw  += weight[i];
    }
    // variances;
    double vxx = (sxx - sx * sx / sw) / sw;
    double vxy = (sxy - sx * sy / sw) / sw;
    double vyy = (syy - sy * sy / sw) / sw;
    // principal axis의 기울기;
    double theta = atan2(2 * vxy, vxx - vyy) / 2;
    line[0] = cos(theta);
    line[1] = sin(theta);
    // center of mass (xc, yc);
    line[2] = sx / sw;  
    line[3] = sy / sw;
    // line-eq:: sin(theta) * (x - xc) = cos(theta) * (y - yc);
    return (1);
};
void test_main(std::vector<POINT>& pts, double line_params[4]) {
    // initial weights = all equal weights;
    std::vector<double> weight(pts.size(), 1); 
    while (1) {
       LineFit_PCA(&pts[0], pts.size(), &weight[0], line_params) ;
       //(1) check goodness of line-fitting; if good enough, break loop;
       //(2) re-calculate weight, normalization not required.
     }
};

아래 그림은 weight를 구하는 함수로 $weight= 1 /\sqrt{1+dist\times dist}$를 이용하고, fitting 과정을 반복하여 얻은 결과다. 상당히 많은 outlier가 있음에도 영향을 덜 받는다. 파란 점이 outlier이고, 빨간 직선은 outlier가 없는 경우 fitting 결과고, 파란선은 outlier까지 포함한 fitting 결과다.

##: 네이버 블로그에서 이전;

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

Fast Float Sqrt  (0) 2020.12.27
2차원 Gaussian 분포 생성  (0) 2020.12.10
PCA Line Fitting  (0) 2020.11.12
Histogram Equalization  (0) 2020.11.12
Least Square Fitting of Circle  (0) 2020.11.11
Integer Sqrt  (0) 2020.11.11
Posted by helloktk

댓글을 달아 주세요

728x90

서로 다른 조건에서 생성이 된 두 개의 영상을 비교하기 위해서는 먼저 영상을 정규화시키어야 한다. 쉽게 생각할 수 있는 정규화의 방법은 두 개의 영상이 동일한 히스토그램을 같도록 만들면 된다. 어떻게 동일한 히스토그램을 만들 수 있는가? 간단히 생각할 수 있는 방법은 주어진 영상의 히스토그램을 전 픽셀에 걸쳐서 균일한 빈도를 보이는 히스토그램이 되도록 변환을 하는 것이다. 즉, 히스토그램을 평탄화시키는 것이다. 이 변환을 $F(x)$라고 하면, 입력 영상의 픽셀 값 $x$는 새로운 픽셀 값 $y = F(x)$로 바뀐다. 영상의 히스토그램은 픽셀 값의 빈도를 나타내므로 확률 밀도 함수로 해석할 수 있고 (픽셀 값을 연속 변수로 볼 때) 주어진 영상의 확률 밀도 함수를 $P(x)$라고 하자. 그러면 $P(x) dx$는 픽셀 값이 $[x, x+dx]$인 픽셀이 영상에서 나타날 확률을 의미한다. 변환 후에서는 임의의 픽셀 값이 균일해야 하므로 확률 밀도 함수는 상수 $a$이고, 픽셀 값의 구간이 $[0,255]$ 이므로 $a = 1/ 255$로 주어진다. 따라서, 변환이 확률을 보존하기 위해서는 $$P(x) dx = a dy$$을 만족해야 하므로 변환 후 픽셀 값 y와 변환 전 픽셀 값 x사이의 관계를 알 수 있다. $$\frac {dy}{ dx} = \frac {1}{a}  P(x) $$ 이 식을 적분하면, $$y = F(x) = \frac {1}{a} \int_0^x P(x') dx' =  255\times  \text {cumulative sum of } P(x)$$ 즉, 히스토그램을 평탄화시키는 변환은 원래 영상 히스토그램의 누적합을 전체 픽셀 수로 나눈 값, $$k \rightarrow F(k)=255 \times \frac {\sum_{i=0}^{k} H [i]}{ \sum_{i=0}^{255} H [i]}$$으로 결정된다. 영상의 히스토그램은 이산적이기 때문에 실제로 모든 픽셀 값에서 균일하게 나타나지 않지만, 구간 변환은 거의 일정하게 나타난다. 이 변환을 거친 영상의 누적 히스토그램(cumulative histogram)을 $y$-축 픽셀 값을 $x$-축으로 해서 그리면 거의 일직선으로 증가하는 모양을 나타낸다.

이 히스토그램의 평탄화는 픽셀 값이 좁은 영역에 몰려있는 영상을 전 영역에 분포하도록 한다. 인간은 영상의 밝기에 의해서 보다는 밝고 어두움의 변화의 크기에 의해서 인지도가 증가하게 되는데, 평탄화된 영상은 보다 쉽게 인식이 될 수 있는 구조를 가진다.

int histogrm_eq(BYTE *src, int width, int height, BYTE *dst) {
    int nn = width * height;
    int hist[256] = {0};
    int cum[256];
    for (int k = 0; k < nn; k++) hist[src[k]]++ ;
    cum[0] = hist[0];
    for (int k = 1; k < 256; k++) 
    	cum[k] = hist[k] + cum[k - 1];
    if (cum[255] == 0) return 0;  // null image;
    double a = 255. / cum[255];
    for (int k = 0; k < nn; k++) {
        int x = int(a * cum[src[k]] + 0.5);
        dst[k] = x > 255 ? 255: x < 0 ? 0 : x;
    }
    return 1;
}

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

2차원 Gaussian 분포 생성  (0) 2020.12.10
PCA Line Fitting  (0) 2020.11.12
Histogram Equalization  (0) 2020.11.12
Least Square Fitting of Circle  (0) 2020.11.11
Integer Sqrt  (0) 2020.11.11
Parabolic Interpolation in Peak Finding  (3) 2020.11.10
Posted by helloktk

댓글을 달아 주세요

728x90

점집합을 일반적인 2차 곡선으로 피팅하는 경우에 방정식은

$$ a x^2 + by^2 + cxy +d x + ey +f = 0$$

의 계수를 주어진 데이터를 이용하여서 구해야 한다. 실제 문제에서는 타원, 포물선 쌍곡 선등의 타입에 따라 몇 가지 제약 조건을 넣어 피팅을 한다. 원은 타원의 특별한 경우로 일반적으로 $a = b, c = 0$의 제약 조건이 필요하다. 그러나 보다 엄밀하게 제약을 하게 되면 $a = b = 1$의 추가 조건을 줄 수 있다. 이 경우는 점들이 모두 일직선에 있는 경우를 ($a = b = 0$) 취급할 수 없게 된다. 이 예외적인 경우를 제외하고는 최소자승법을 사용하면 계수를 매우 쉽게 구할 수 있기 때문에 많이 이용된다.

 

문제: 

$$x^2  + y^2 + A x + B  y + C = 0$$

에서 $A, B, C$를 최소자승법을 사용해서 구하라. 

 

주어진 점집합이 원 위의 점이면 우변이 0이 되어야 하나, 실제 데이터를 얻는 과정에서 여러 노이즈에 노출되므로 일반적으로 0이 되지 않는다. 최소자승법은 주어진 점들이 원에서 벗어나는 정도의 제곱 합이 최소가 되도록 하는 계수 $A, B, C$를 결정한다.  원과 점의 편차의 제곱합
$$ L=\sum_ i   \left |x_i^2 + y_i^2 + A x_i + B y_i + C \right|^2 , $$

의 극값을 찾기 위해서 $A, B,$ 그리고 $C$에 대해 미분을 하면

$$\frac{\partial L}{\partial A} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) x_i = 0, $$

$$\frac{\partial L}{\partial B} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) y_i = 0, $$

$$\frac{\partial L}{\partial C} = 2 \sum_i (x_i^2 + y_i^2 + A x_i + B y_i + C) = 0. $$

이 연립방정식을 풀면  $A, B, C$를 구할 수 있다. 우선 세번째 식에서 

$$ CN = -S_{x^2} - S_{y^2} - AS_x - BS_y ,$$

을 얻고, 이를 첫번째와 두번째 식에 각각 대입하면

$$A ( NS_{x^2} - S_x^2) + B ( N S_{xy} - S_x S_y ) =-N S_{x^3} - N S_{xy^2} + S_{x^2} S_x + S_{y^2} S_x, $$

$$A ( NS_{xy} - S_x S_y ) + B ( N S_{y^2} - S_y^2) = -N S_{y^3} - N S_{x^2 y}  +S_{x^2} S_y +S_{y^2} S_y, $$

을 얻을 수 있다. 다시 정리하면 두 개의 연립방정식

$$-a_1 A - a_2 B = 2c_1,$$

$$-b_1 A - b_2 B = 2c_2,$$

을 얻는다. $a_1, a_2 =  b_1, b_2, c_1, c_2$는 코드에서 정의되어 있다. 그리고

따라서, 추정된 원의 중심 $(c_x, c_y)$는 

$$ c_x = - \frac{A}{2} = \frac{c_1 b_2  - c_2 a_2}{ a_1 b_2 - a_2 b_1},$$

$$ c_y = - \frac{B}{2} = \frac{-c_1 b_1 + c_2 a_1}{a_1 b_2 -a_2 b_1},$$

로 주어지고, 반지름은 

$$r^2 =  c_x^2 +c_y^2 - C = c_x^2 + c_y^2 + \frac{1}{N}( S_{x^2}+S_{y^2}- 2c_x S_x - 2 c_y S_y)$$

로 주어진다.

/* 구현 코드 */
BOOL circleFit_LS(POINT Q[], int N, POINT *center, double *radius) {
    double sx  = 0.0,  sy = 0.0;
    double sx2 = 0.0, sy2 = 0.0, sxy  = 0.0;
    double sx3 = 0.0, sy3 = 0.0, sx2y = 0.0, sxy2 = 0.0;
    /* compute summations */
    for (int k = 0; k < N; k++) {
        double x = Q[k].x, xx = x * x;
        double y = Q[k].y, yy = y * y;
        sx   += x;       sy   += y;
        sx2  += xx;      sy2  += yy;      sxy  += x * y;
        sx3  += x * xx;  sy3  += y * yy;
        sx2y += xx * y;  sxy2 += yy * x;
    }
    /* compute a's,b's,c's */
    double a1 = 2.0 * (sx * sx - sx2 * N);
    double a2 = 2.0 * (sx * sy - sxy * N);
    double b1 = a2;
    double b2 = 2.0 * (sy * sy - sy2 * N);
    double c1 = (sx2 + sy2) * sx - (sx3 * N + sxy2) * N;
    double c2 = (sx2 + sy2) * sy - (sy3 * N + sx2y) * N;
    double det = a1 * b2 - a2 * b1;
    if (fabs(det) < 1.e-10)                /*collinear한 경우임;*/
        return FALSE;
    /* floating point center */
    double cx = (c1 * b2 - c2 * b1) / det;
    double cy = (a1 * c2 - a2 * c1) / det;
    /* compute radius squared */
    double radsq = cx * cx + cy * cy + (sx2 + sy2 - 2 * (sx * cx  + sy * cy)) / N;
    *radius = sqrt(radsq);
    /* integer center */
    center->x = int(cx + .5) ;
    center->y = int(cy + .5) ;
    return TRUE;
}

 

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

PCA Line Fitting  (0) 2020.11.12
Histogram Equalization  (0) 2020.11.12
Least Square Fitting of Circle  (0) 2020.11.11
Integer Sqrt  (0) 2020.11.11
Parabolic Interpolation in Peak Finding  (3) 2020.11.10
Histogram Matching  (0) 2012.11.03
Posted by helloktk

댓글을 달아 주세요

728x90

정수의 제곱근을 정수 연산만 가지고 계산해야 하는 경우가 많이 있는데 (특히, 임베디드 환경에서 프로그래밍을 하는 경우) 이를 구현하는 방법에 대해서 알아보자. 32비트 정수의 제곱근은 16비트 이하이므로 가장 최상위의 비트 값이 세팅이 된 정수를 가지고 제곱근 조건을 만족시키는지 확인하고, 만족시키면 다음의 하위 비트 값이 세팅이 된 정수 값을 기존의 값에 더해서 다시 제곱근 조건을 만족시키는지 확인하는 과정을  마지막 비트까지 반복하면 된다. 그런데 가능한 가장 최상위의 비트 값을 갖는 정수는 0x8000 (= 1<<15) 이므로 양의 정수 x의 제곱근을 구하는 과정은 아래와 같다.

           unsigned s = 0x8000;  //= (1 << 15);
           unsigned t = 0;
           do {                          //새로 세팅된 값이 제곱근 조건을 만족하는가?
                if ((t + s) * (t + s) <= x ) t = t + s;  //만족하면 하위 비트값을 더함;
                s = s >> 1;                                //다음 하위 비트로 내려감;
           } while (s != 0) ;
           return (t);

이 알고리즘은 16번의 반복으로 항상 끝이 나지만, 중간에 곱셈 연산이 들어가 있다. 이 곱셈 연산을 없애는 방법을 알아보자. 먼저           

 (t + s) * (t + s) = t * t + 2 * t * s + s * s;

이므로 if 문 ( ) 내부 식 양변에서 t^2를 빼도록 하자. 이 경우 t의 값이 갱신이 되면

 x - t^2 => x - (told + sold)^2 => x - told * told - (2 * told * sold + sold * sold) ;

이므로, x 값에서도 역시 2 * t * s + s * s의 값을 제거해야 한다. 따라서

                  unsigned nr = 2 * t * s + s * s;
                  if (nr <= x) {
                         x = x - nr;
                         t = t + s;
                  }

로 고쳐야 한다. 이제, s * s 와 2 * t * s을 별도의 변수를 도입하여 컨트롤하면,

          unsigned m = s * s;    // m = 0x40000000;
          unsigned r = 2 * t * s ;// => 2 * (told + sold) * sold
                                        // = 2 * told * sold + 2sold * sold 
                                        // = r + m + m = nr + m;
          do {
                  nr = m + r ;
                  if (nr <= x) {
                         x = x - nr ;
                         r = nr + m;
                         t = t + s;
                  } ;
                  s = s >> 1;
                  m = m >> 2;
                  r = r >> 1;
          } while (s != 0);

로 사용이 가능하다. 그런데 여기서 r = 2 * t * s인데, 마지막으로 루프를 빠져나올 때 s가 1인 상태에서 s / 2를 하므로 (실수 연산을 고려하면) 이때 r 값은 t 값과 같다. 따라서 t 변수를 사용할 필요가 없고, 마찬가지로 s 변수도 사용할 필요가 없어진다 (m을 컨트롤 변수로 이용). 

//16-step integer_sqrt;
unsigned isqrt( unsigned x ) {
       unsigned m = 0x40000000 ;   //=(1 << 30)
       unsigned r = 0, nr;
       do {
            nr = m + r ;
            if (nr <= x) {
                   x -= nr ;
                   r = nr + m;
            }
            m >>= 2;
            r >>= 1;
       } while (m != 0);
       return (r) ;
};

결과는 항상 <= sqrt(x) 이다. 실제로 빠를까?

 

한번 더 다듬도록 하자. 나중에 r을 시프트 하지 말고 먼저 하면, r = 2 * t * s => r / 2 = t * s 이므로, if 문의 r은 r =  r + m이 된다. 정리하면

unsigned isqrt( unsigned x ) {
       unsigned m = 0 x 40000000 ;   //= (1 << 30)
       unsigned r = 0, nr;
       do {
            nr = m + r ;
            r >>= 1;
            if (nr <= x) {
                   x -= nr ;
                   r += m;
            }
            m >>= 2;
       } while (m != 0);
       return (r) ;
};

이 된다. 이는 잘 알려진 integer sqrt의 값을 구하는 함수이다.

 

다시 한번 더 형태를 바꾸면, r = r + m 항을 보면, r은 m을 누적하는 변수이고, m은 하나의 비트만 세팅이 된 변수이므로 덧셈을 비트 or 연산으로 바꿀수 있다.

             r += m   ==> r |= m;

그리고 매크로를 도입하여 루프를 풀어쓰기 하면, 아래의 코드를 얻을 수 있다.

unsigned isqrt( unsigned x ) {
       unsigned m = 0x40000000 ;   //= (1 << 30)
       unsigned r = 0, nr;
#define STEP(k)             \
           nr = m + r;      \
           r <<= 1;         \
           if (nr <= x) {   \
                x -= nr;    \
                r |= m;     \
           }                 
       STEP(15) ;  m>>=2;
       STEP(14) ;  m>>=2;
       STEP(13) ;  m>>=2;
       STEP(12) ;  m>>=2;
       STEP(11) ;  m>>=2;
       STEP(10) ;  m>>=2;
       STEP(9)  ;  m>>=2;
       STEP(8)  ;  m>>=2;
       STEP(7)  ;  m>>=2;
       STEP(6)  ;  m>>=2;      
       STEP(5)  ;  m>>=2;
       STEP(4)  ;  m>>=2;
       STEP(3)  ;  m>>=2;
       STEP(2)  ;  m>>=2;
       STEP(1)  ;  m>>=2;
       STEP(0)  ;
       return (r) ;
};

이 코드 또한 인터넷 상에 ARM 개발자 사이트나 게임 개발자 사이트에 널리 알려진 코드이다.

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

Histogram Equalization  (0) 2020.11.12
Least Square Fitting of Circle  (0) 2020.11.11
Integer Sqrt  (0) 2020.11.11
Parabolic Interpolation in Peak Finding  (3) 2020.11.10
Histogram Matching  (0) 2012.11.03
삼각형의 외접원: 외접원의 중심 2  (2) 2012.10.20
Posted by helloktk

댓글을 달아 주세요

728x90

히스토그램처럼 균일한 간격으로 주어진 데이터에서 피크 값의 위치를 찾기 위해서는 먼저 노이즈 제거를 위해 몇 번의 smoothing 과정을 적용해야 한다 (구체적인 방법은 필요에 따라 다양하다). smoothing 과정을 거친 히스토그램에서 피크 위치를 찾는 것은 쉬운 작업이다. 그런데 경우에 따라 그 위치를 sub-order까지 구해야 필요가 생긴다. 예를 들면, 실수 값 데이터 시퀀스 대한 히스토그램을 만들려면 먼저 정수로 양자화시켜야 가능하다. 이 양자화된 정보를 담고 있는 히스토그램에서 피크의 위치를 찾으려 할 때 양자화 과정에서 잃어버린 정밀도를 복원하려면 interpolation을 써야 한다. 간단하게 parabolic 근사로 피크의 위치를 찾는 경우를 생각하자. 이 방법은 수학적으로 단순하지만 영상처리 알고리즘에서 많이 이용이 되고 있다. 데이터 시퀀스가 균일한 간격에서 주어졌으므로 계산은 $-1, 0, 1$의 세 군데 위치에서 중앙 근방에 피크가 나타나는지를 고려하면 된다. 세 점에서 히스토그램 값이 각각 $h_m, h_0, h_p$일 때 이들을 지나는 이차 곡선의 피크 위치를 찾자. 주어진 이차 곡선을

$$y = a  (x - c)^2 + b$$

꼴로 쓰면 $c$가 0에서 벗어난 정도를 나타낸다.

$$(-1, h_m) \quad \rightarrow \quad h_m = a(-1-c)^2 +b; \\(0, h_0) \quad \rightarrow \quad h_0 = a c^2 +b; \\ (+1, h_p) \quad \rightarrow \quad h_m = a(+1-c)^2 +b; $$ 이므로 $$h_m - h_p = 4ac, \\ h_m + h_p -2 h_0=2 a,\\ \therefore ~c = \frac { h_m - h_p}{2(h_m + h_p -2 h_0 )}$$

아래 코드는 피크의 위치가 중앙점에서 벗어난 정도를 준다.

bool parabolicInterpolate(double hm, double h0, double hp, double *c) {
      double a = hm + hp - 2 * h0;
      if (a >= 0) return false; // not a parabola(a==0), not convex up(a>0);
      *c = 0.5 * (hm - hp) / a;
      if (*c < -0.5 || *c > 0.5) return false; //  too far;
      return true;
}

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

Least Square Fitting of Circle  (0) 2020.11.11
Integer Sqrt  (0) 2020.11.11
Parabolic Interpolation in Peak Finding  (3) 2020.11.10
Histogram Matching  (0) 2012.11.03
삼각형의 외접원: 외접원의 중심 2  (2) 2012.10.20
삼각형의 외접원: 외접원의 중심  (0) 2012.10.19
Posted by helloktk

댓글을 달아 주세요

  1. 2021.03.24 22:29  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

  2. 2021.03.25 06:47  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

728x90

동일한 물체나 풍경이라도 촬영에 사용된 센서의 종류나 조명 조건 때문에 사진이 더 어둡거나 반대로 더 밝게 표현될 수 있다. 이런 경우 사진을 잘 찍힌 사진과 비교하여 보정하려면 어떤 기준을 가지고 작업을 해야 할까? 동일한 내용의 영상이라면 단순히 히스토그램을 늘리거나 줄이는 과정을 (histogram stretching) 사용할 수 있다. 그러나 같은 물체나 풍경을 담지 않은 기준 영상의 히스토그램과 같은 내용의 히스토그램을 되도록 사진을 조작하고 싶으면 histogram stretching으로는 안된다. 

 

입력 영상의 주어진 픽셀 값이 기준 영상의 어떤 픽셀 값에 해당하게 변환을 해야 하는가를 결정해야 한다. 변환 $f$는 두 영상의 픽셀 값의 확률분포를 보존하도록 만들어진다. 입력 영상의 히스토그램 구간 $(i-1, i]$이 기준 영상의 히스토그램 구간 $(f(i-1), f(i)]$으로 변환될 때 각각 구간에서 픽셀이 나타날 확률이 동일해야 한다. 영상에서 픽셀 값은 이산적이므로 확률보다는 확률의 누적합을 이용하는 것이 더 편리하다. 입력 영상의 히스토그램 구간 $[0, i]$에서 픽셀이 나타날 누적 확률과 같은 누적 확률을 주는 기준 영상의 히스토그램 구간이 $[0, k=f(i)]$일 때 

$$i(영상)\rightarrow k = f(i) (기준영상)$$

으로 변환하면 된다. 이 과정은 히스토그램 stretching에도 그대로 적용이 된다. 누적 히스토그램(cumulative histogram)은 주어진 픽셀 값까지 확률 합을 쉽게 구할 수 있게 하는 일종의 lookup table이다. cumulative histogram의 $i$-번째 값은 픽셀 값이 $0$부터 $i$까지 히스토그램의 누적합이기 때문이다. 따라서 히스토그램 매칭(histogram matching, histogram specification)은 같은 cumulative histogram의 값을 갖는 픽셀 값 사이의 변환을 의미한다. 두  영상의 픽셀 수가 같지 않은 경우에는 cumulative histogram를 쓸 수 없다. 이 경우에는 전체 픽셀 수로 정규화시킨 cumulative distribution function(cdf)을 사용하면 된다.

 

 

히스토그램 평탄화(histogram equalization)은 각 그레이 값이 나타날 확률이 균일한 기준 영상을 (즉, 1차 cumulative histogram이 직선인 영상) 사용한 히스토그램 매칭에 해당한다. 따라서 기준 영상이 필요하지 않는다. 

 

 

//map을 lookup table로 이용해서 영상을 변환할 수 있다.
void matchHistogram(int srcHist[256], int refHist[256], int map[256]) {
    double cumSrc[256], cumRef[256];
    cumSrc[0] = srcHist[0];
    cumRef[0] = refHist[0];
    for (int i = 1; i < 256; i++) {
        cumSrc[i] = cumSrc[i - 1] + srcHist[i];
        cumRef[i] = cumRef[i - 1] + refHist[i];
    }
    // normalize ==> make CDF; 두 비교 영상이 같은 크기가 아닐 수 있으므로 필요하다;
    for (int i = 0; i < 256; i++) {
        cumSrc[i] /= cumSrc[255];
        cumRef[i] /= cumRef[255];
    }
    // 두 비교영상의 cumulant histogram이 가장 가까운(같은) 픽셀끼리 변환한다;
    for (int i = 0; i < 256; i++) {
        int k = 255;  
        while (cumRef[k] > cumSrc[i]) k--;
        map[i] = k < 0 ? 0 : k;
        // 좀 더 smooth한 interpolation이 필요함;
    }
};

실제 사용은 grey 이미지인 경우:

더보기
int histogramSpec(BYTE *src, int ws, int hs, BYTE *ref, int wr, int hr, BYTE *out) {
    int srcHist[256], int refHist[256], int map[256];
    makeHistogram(src, ws, hs, srcHist);
    makeHistogram(ref, wr, hr, refHist);
    matchHistogram(srcHist, refHist, map);
    for (int k = ws * hs - 1; k >= 0; k--)
        out[k] = map[src[k]];
    return 1;
}

 

Posted by helloktk

댓글을 달아 주세요

728x90

삼각형의 외접원을 기하학적인 방법이 아닌 대수적인 방법을 이용해서 구해보자. 중심이 $(x_c, y_c)$이고 반지름이 $R$인 원의 방정식은

$$(x-x_c)^2 + (y- y_c)^2 = R^2$$

으로 주어진다. 이 식에서 세 개의 변수 $x_c$, $y_c$, $R$를 결정해야는데, 각 변수에 대해서 2차 방정식의 형태로 주어져 있으므로 좀 더 쉬운 1차식 형태가 되도록 변형을 해보자. 방정식을 전개하면

$$(2x) x_c + (2y) y_c + R^2 - x_c^2 - y_c^2 = x^2 + y^2  \\ \longrightarrow ~(2x) x_c + (2y) y_c +d = x^2 +y^2,\quad d=R^2 -x_c^2 -y_c^2.$$

이제 원의 방정식은 $x_c$, $y_c$, $d$에 대한 1차식으로 정리되었다. 삼각형의 외접원의 방정식은 세 꼭짓점을 통과하므로 세 꼭지 짐 $A(x_1, y_1)$, $B(x_2, y_2)$, $C(x_3, y_3)$을 위의 식에 넣으면 3개의 방정식을 얻는다. 이 방정식을 행렬로 표현하면

 

따라서, Cramer의 규칙을 적용하면 다음과 같이 외접원의 중심 (또한 $d$을 구해서 외접원의 반지름을 구할 수 있다)을 얻을 수 있다.

 

determinant를 정리하면,

 

을 얻을 수 있다. 이 식의 분모는 세 꼭짓점에 동일 직선 위에 있지 않으면 0이 아니다. 한 꼭짓점과 다른 꼭짓점 간 좌표의 차이 4개 $(x_2-x_1, y_2-y_1, x_3-x_1, y_3-y_1)$와 좌표의 합 4개 $(x_1+x_2, y_1+y_2, x_1+x_3, y_1+y_3)$를 계산하여 중심점의 좌표를 얻을 수 있다. 

//

구현된 코드는 http://kipl.tistory.com/113; 또는

int circumcenter2(CPoint P, CPoint Q, CPoint R, double *xc, double *yc) {
    double A = Q.x - P.x, B = Q.y - P.y;
    double C = R.x - P.x, D = R.y - P.y;
    double E = A * (P.x + Q.x) + B * (P.y + Q.y);
    double F = C * (P.x + R.x) + D * (P.y + R.y);
    double G = 2. * ( A * D - B * C);
    if (G != 0.) { // 세 점이 일직선에 놓이지 않는 경우; 이 경우만 외접원이 정의된다;
        *xc = (D * E - B * F) / G;
        *yc = (A * F - C * E) / G;
        return 1;
    } else return 0;
}

분모 $G$의 계산을 $G = 2 ( A  (R.y - Q.y) - B  (R.x - Q.x))$ 로 쓰는 경우가 많다. 위의 행렬식을 보면 같은 결과지만 식을 좀 더 대칭적으로 보이게 한다: $G = 2.(A  (D - B) - B  (C - A)) = 2 (A D - B C)$.

Posted by helloktk

댓글을 달아 주세요

  1. 지나가던이 2013.09.12 14:09  댓글주소  수정/삭제  댓글쓰기

    d가 있어서 비선형방정식인데 determinant에 관한 법칙들과 crammer's rule을 적용하는게 가능한가요?

    • helloktk 2013.09.14 01:37 신고  댓글주소  수정/삭제

      미지수를 xc, yc,R 로 보면 비선형 방정식이지만, xc, yc, d를 변수로 보면 선형방정식 형태입니다 (d를 구하면 R은 구할 수 있지요.)

728x90

평면 위 세 점 $A$, $B$, $C$가 삼각형을 만들 때 외접원 중심을 구해보자. 점 $C$을 원점으로 하면, $A$, $B$는 각각

의 두 벡터로 표현이 된다 (원점을 $C$점으로 평행 이동했다고 생각하면 된다). 이 두 벡터에 각각 수직이고 삼각형과 같은 평면에 놓인 ($\vec a$, $\vec b$를 각각 시계 방향으로 회전시켜 만든) 두 벡터 $\vec m$, $\vec n$를 다음과 같이 정의하자.

그러면, 변 $CA$의 중점을 지나면서 수직인 직선의 방정식은

또, $CB$의 중점을 지나면서 수직인 직선의 방정식은 

이 두 직선의 교점이 외접원의 중심이 된다. 매개변수 $s$, $t$를 구하기 위해서 두 식을 빼면

 매개변수 $s$는 $\vec b$와 $\vec n$이 수직인 조건을 이용하면

따라서, 외접원의 중심은  

여기서, 3개 벡터의 외적의 성질을 이용해서 (평면에서 외적은 숫자이다)

가 항등식으로 성립하고, 또한 $\vec a$, $\vec m$이 서로 수직이면서 길이가 같고, $\vec m$ 도 $\vec b$에 수직이면 길이가 같으므로

따라서, 

여기서, 윗 식의 분모를 보면

이므로, 이 값이 0이 아니려면 세 점이 일직선에 있지 않으면 된다.

그런데 이 벡터는 점 $C$를 원점으로 하여서 계산을 한 것이므로, 원래의 좌표계에 대한 식으로 바꾸려면 $(C_x, C_y)$를 더해 주어야 한다.

또는, 성분으로 표현하면

int circumcenter ( CPoint A, CPoint B, CPoint C, double *xc, double *yc) {
    double ax = A.x - C.x, ay = A.y - C.y ;
    double bx = B.x - C.x, by = B.y - C.y ;
    double asq = ax * ax + ay * ay;
    double bsq = bx * bx + by * by;
    double ccw = ax * by - ay * bx;
    if ( ccw != 0. ) { //세 점임 일직선 위에 있지 않는 경우; 이 경우만 외접원이 정의됨;
        *xc = C.x + ( by * asq - ay * bsq ) / ( 2 * ccw ) ;
        *yc = C.y + ( -bx * asq + ax * bsq ) / ( 2 * ccw ) ;
        return 1;
    } else return 0;
}

 

Posted by helloktk

댓글을 달아 주세요