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' 카테고리의 다른 글

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

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' 카테고리의 다른 글

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

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

segmentation result

1. 경계영역 처리를 보완해야 한다.

2. distance-map을 충분히 smooth하게 만들어야 한다: mean-filter 반복 적용.

3. 참고: imagej.net/Distance_Transform_Watershed의 구현보다 단순하지만 성능은 떨어진다.

#define insideROI(x, y) ((x) >= 0 && (x) < w && (y) >= 0 && (y) < h)
#define BACKGND 0xFF
int watershed(float *distMap, int w, int h, int *rgnMap) {
    const int xend = w - 1, yend = h - 1;
    const int nx[] =  {-1, 0, 1, -1, 0, 1, -1, 0, 1};
    const int ny[] =  {-1, -1, -1, 0, 0, 0, 1, 1, 1};
    // processing image;
    float **procImg = new float * [h];
    procImg[0] = new float [h * w];
    for (int y = 1; y < h; y++) procImg[y] = procImg[y - 1] + w;
    // presmooth; 10 iterations;
    float **pp = new float* [h];
    for (int y = 0; y < h; y++) pp[y] = &distMap[y * w];
    mean_filter(pp, w, h, 15, procImg); 

    // 입력 영상에서 경계는 global maximum으로 설정;
    for (int y = 0; y < h; y++ ) procImg[y][0] = procImg[y][xend] = BACKGND;
    for (int x = 0; x < w; x++ ) procImg[0][x] = procImg[yend][x] = BACKGND;
    // find local minima for each pixel;
    int minCnt = 1;
    for (int y = 0, pos = 0; y < h; y++) {
        for (int x = 0; x < w; x++, pos++) {
            float minVal = procImg[y][x];
            if (minVal == BACKGND) {
                rgnMap[pos] = 1;	// background(white);
                continue;
            }
            int minIdx = 4; //current position;
            for (int k = 0; k < 9; k++) {
                int xx = x + nx[k], yy = y + ny[k];
                if (insideROI(xx, yy) && procImg[yy][xx] <= minVal) {
                    minVal = procImg[yy][xx];
                    minIdx = k;
                }
            }
            if (minIdx == 4) rgnMap[pos] = ++minCnt; // center(x, y) is a new local minimum;
            else             rgnMap[pos] = -minIdx;  // direction of local-min =  "-(dir)"
        }
    }
    // follow gradient downhill for each pixel;
    // reset rgnMap to have the id's of local minimum connected with current pixel;
    for (int y = 1, pos = y * w; y < yend; y++ ) {
        pos++;
        for (int x = 1; x < xend; x++, pos++ ) {
            int minpos = pos;
            while ( rgnMap[minpos] <= 0 ) { //it is not a local minimum.
                switch ( rgnMap[minpos] ) {
                case  0: minpos += -w - 1; break; // top-lef = downhill direction;
                case -1: minpos += -w;	   break; // top
                case -2: minpos += -w + 1; break; // top-right;
                case -3: minpos--;         break; // left;
                case -5: minpos++;         break; // right;
                case -6: minpos += w - 1;  break; // bot-left;
                case -7: minpos += w;      break; // bot;
                case -8: minpos += w + 1;  break; // bot-right;
                }
            }
            // speed-up: use a stack to record downhill path.
            rgnMap[pos] = rgnMap[minpos]; //assign the id of a local minimum connected with current pixel;
        }
        pos++;
    }
    // rgnMap는 각 픽셀과 연결되는 국소점의 id을 알려주는 lookup table이다;
    delete [] procImg[0]; delete [] procImg;
    delete [] pp;
    return ( minCnt );
}

 

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

Watershed Segmentation  (0) 2021.02.27
Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (0) 2021.02.15
Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
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; }
pixelvalue 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]) ;
}
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

지문에서 ridge의 방향(orientation)은 gradient에 수직한 방향이다(그런데 벡터인 gradient와는 달리 ridge의 방향은 모호함이 있다. 시계방향 또는 반시계방향으로 90도 회전이 모두 동일한 ridge의 방향이다).  gradient의 방향각이 $\alpha$일 때 ridge의 방향각은 $\theta =\frac{\pi}{2} + \alpha$로 정의하자. 방향각을 주어진 gradient 성분 $\nabla I = (g_x, g_y)$로 표현하기 위해서

$$ \sin 2\alpha = 2 \sin \alpha \cos \alpha = 2 \frac{g_x}{\sqrt{g_x^2 + g_y^2}}\frac{g_y}{\sqrt{g_x^2 + g_y^2}}$$

$$\cos 2\alpha = \cos^2 \alpha - \sin^2\alpha = \frac{g_x^2}{g_x^2 + g_y^2} -\frac{g_y^2}{g_x^2 +g_y^2}$$ 

임을 이용하자. $$\tan 2\alpha = \frac{ 2g_x g_y}{ g_x^2 - g_y^2}$$로 쓸 수 있으므로 ridge의 방향각은

$$ \theta =\frac{\pi}{2} + \frac{1}{2} \arctan \frac{2g_x g_y}{g_x^2 - g_y^2}$$

임을 알 수 있다. $(g_x, g_y) \rightarrow (-g_x, -g_y)$이더라도 ridge의 방향은 동일하므로 $\theta$의 범위는 $0\le \theta < \pi$로 제한된다.

 

이미지의 한 지점에서 ridge의 방향은 그 점을 중심으로 하는 적당한 크기의 윈도 평균값으로 추정한다. gradient 성분은 Sobel 연산자나 Prewitt 연산자를 convolution해서 얻을 수 있다. 따라서 지문 이미지 상의 픽셀 좌표 $(i,j)$에서 ridge의 방향은

$$\theta_{ij} = \frac{\pi}{2} +\frac{1}{2}\arctan \frac{2 G_{xy}}{G_{xx} - G_{yy}}$$

$$ G_{xy} = \sum_{W} \text{gradX}_{ij} \times \text{gradY}_{ij}$$

$$ G_{xx} = \sum_{W} \text{gradX}_{ij} \times \text{gradX}_{ij}$$

$$ G_{yy} = \sum_{W} \text{gradY}_{ij} \times \text{gradY}_{ij}$$

아래 결과는 5x5 크기의 Gaussian filter를 적용한 지문 이미지에 17x17 크기의 윈도를 겹치지 않게 잡아서 평균을 계산한 결과다. $G_{xx}+G_{yy}$가 최대값의 $20\%$보다 작은 영역은 방향을 표시하지 않았다 (red:0~45, yellow:45~90, magenta: 90~135, cyan: 135~180).

int RidgeOrientaion(BYTE *img, int w, int h, int bsz/*17*/, int *omap) {
    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};
    const int xmax = w - 1, ymax = h - 1;
    int *buffer = new int [2 * w * h];
    memset(buffer, 0, sizeof(int) * 2 * w * h);
    int *gradx = &buffer[0];
    int *grady = gradx + (w * h);
    for (int y = 1, pos = w; y < ymax; y++) {
        pos++; //x = 0;
        for (int x = 1; x < xmax; x++, pos++) {
            int sx = 0, sy = 0; 
            for (int k = 0; k < 9; k++) {
                int v = img[pos + nn[k]];
                sx += sobelX[k] * v;
                sy += sobelY[k] * v;
            }
            gradx[pos] = sx;
            grady[pos] = sy;
        }
        pos++; // x = xmax;
    }
    const int ny = h / bsz;
    const int nx = w / bsz;
    for (int iy = 0; iy < ny; iy++) {
        for (int ix = 0; ix < nx; ix++) {
            int sxy = 0, sxx = 0, syy = 0; 
            int pos = (iy * w + ix) * bsz; //starting position of block(ix, iy)
            int *dx = &gradx[pos];
            int *dy = &grady[pos];
            for (int jj = 0; jj < bsz; jj++) {
                for (int ii = 0; ii < bsz; ii++) {
                    int gx = dx[jj * w + ii];
                    int gy = dy[jj * w + ii];
                    sxy += gx * gy;
                    sxx += gx * gx;
                    syy += gy * gy;
                }
            }
            omap[iy * nx + ix] = 0.5 * (PI + atan2(double(2.0 * sxy), double(sxx - syy)));
        }
    }
    //...find fingerprint regions and draw orientational flow;
    delete[] buffer;
    return 1;
}

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

Watershed Segmentation  (0) 2021.02.27
Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (0) 2021.02.15
Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
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

/* "Graphics Gems IV" 의 code를 수정함;*/
int CLAHE ( BYTE* image, int width, int height,
            BYTE gMin/*0*/, BYTE gMax/*255*/, int rgn_nx/*16*/, int rgn_ny/*16*/,
            int bins/*256*/, double fcliplimit )
{
    int* map_array	= new int [ rgn_nx * rgn_ny * bins ];
    int rgn_xsz		= width / rgn_nx;
    int rgn_ysz		= height / rgn_ny; 
    int rgn_area	= rgn_xsz * rgn_ysz;
    int clipLimit	= int( fcliplimit * ( rgn_xsz * rgn_ysz ) / bins );
    clipLimit = ( clipLimit < 1 ) ? 1 : clipLimit;
    //  
    /* calculate greylevel mappings for each contextual region */
    for (int iy = 0; iy < rgn_ny; iy++) {
        for (int ix = 0; ix < rgn_nx; ix++) {
            int *hist = &map_array[bins * ( iy * rgn_nx + ix )];
            int start = (iy * rgn_ysz) * width + ix * rgn_xsz;
            RgnHistogram (&image[start], width, rgn_xsz, rgn_ysz, hist, bins);
            ClipHistogram ( hist, bins, clipLimit );
            //convert hist to cuimulative histogram normalized to [gMin=0, gMax=255]
            MapHistogram ( hist, gMin, gMax, bins, rgn_area );
        }
    }

    /* bilinearInterp greylevel mappings to get CLAHE image */
    int szx, szy, ixl, ixr, iyt, iyb;
    BYTE *pimg = image;
    for (int iy = 0; iy <= rgn_ny; iy++ ) {
        if ( iy == 0 ) {					 
            szy = rgn_ysz / 2; iyt = iyb = 0;
        } else if ( iy == rgn_ny ) {				 
            szy = rgn_ysz / 2; iyt = rgn_ny - 1; iyb = iyt;
        } else {
            szy = rgn_ysz; iyt = iy - 1; iyb = iyt + 1;
        }
        for (int ix = 0; ix <= rgn_nx; ix++ ) {
            if ( ix == 0 ) {				 
                szx = rgn_xsz / 2; ixl = ixr = 0;
            } else if ( ix == rgn_nx ) {			 
                szx = rgn_xsz / 2; ixl = rgn_nx - 1; ixr = ixl;
            } else {
                szx = rgn_xsz; ixl = ix - 1; ixr = ixl + 1;
            }
            // cumulative histogram data;
            int *LU = &map_array[bins * ( iyt * rgn_nx + ixl )];
            int *RU = &map_array[bins * ( iyt * rgn_nx + ixr )];
            int *LB = &map_array[bins * ( iyb * rgn_nx + ixl )];
            int *RB = &map_array[bins * ( iyb * rgn_nx + ixr )];
            BilinearInterp (pimg, width, LU, RU, LB, RB, szx, szy );
            pimg += szx;			  
        }
        pimg += ( szy - 1 ) * width;
    }
    delete [] map_array;				
    return 0;						 
}
더보기
void ClipHistogram ( int* hist, int gLevels, int clipLimit ) {
    int excess = 0;
    for (int i = 0; i < gLevels; i++ ) { /* calculate total number of excess pixels */
        int diff =  hist[i] - clipLimit;
        if ( diff > 0 ) excess += diff;	 /* excess in current bin */
    };

    /* clip histogram and redistribute excess pixels in each bin */
    int incr = excess / gLevels;		/* average binincrement */
    int upper =  clipLimit - incr;		/* bins larger than upper set to cliplimit */
    for (int i = 0; i < gLevels; i++ ) {
        if ( hist[i] > clipLimit ) hist[i] = clipLimit; /* clip bin */
        else {
            if ( hist[i] > upper ) {			/* high bin count */
                excess -= hist[i] - upper;
                hist[i] = clipLimit;
            } else {							/* low bin count */
                excess -= incr;
                hist[i] += incr;
            }
        }
    }

    while ( excess ) { /* redistribute remaining excess  */
        int start = 0;
        while ( excess && start < gLevels) {
            int step = gLevels / excess;
            if ( step < 1 ) step = 1;		 /* step size at least 1 */
            for (int i = start; i < gLevels && excess; i += step) 
                if (hist[i] < clipLimit) { 
                    hist[i]++; excess--;
                }
            start++;
        }
    }
}

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

Watershed Segmentation  (0) 2021.02.27
Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (0) 2021.02.15
Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
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

아날로그 신호로부터 디지털 데이터로 바꿀 때 보통 시간당(sound) 또는 공간적 거리당(image) 일정한 횟수로  데이터를 수집한다. 이 비율을 sampling rate이라 한다. 그런데 sampling rate을 바꾸어야 할 경우가 있다. 1초마다 한 번씩 데이터를 모았는데 실제로 0.5초마다 수집된 데이터가 필요가 한 경우에는 기존의 수집된 데이터와 데이터 사이에서 원 아날로그 신호 값을 알아내야 한다. sampling rate을 바꿀 때 기존의 데이터를 이용해서 원 아날로그 신호에서 수집한 것처럼 새로운 데이터를 만들어 내는 과정이 resampling이다. 이미지에서는 원본을  확대하는 up sampling이나 축소하는 down sampling 등이 resampling의 일종이다. resampling을 하기 위해서는 기존 데이터를 이용해서 데이터와 데이터 사이에서의 값을 추정하는 interpolation이 필요하게 된다. 많이 사용하는 interpolation 방법으로는 nearest neighbor interpolation, linear interpolation, cubic interpolation, spline interpolation 등이 있다.  여기서는 cubic interpolation을 이용한 이미지 resampling을 살펴본다.

 

$x$축 위의 4 점 $x=-1$, $x=0$, $x=1$, $x=2$에서 값이 $p_0$, $p_1$, $p_2$, $p_3$으로 주어진 경우 $0 \le x \le1$ 구간에서 cubic interpolation 함수는 

$$ f(x)=ax^3 +b x^2 + cx +d$$

라고 쓰면, 우선 $x=0$에서 $p_1$, $x=1$에서 $p_2$를 지나므로

$$f(0) = d = p_1, \quad\quad f(1) = a + b+ c+d = p_2$$

를 만족해야 하고, 양끝에서 도함수 값은 주변 입력점의 평균변화율로 주면, 다음식을 만족한다.

$$  f'(0) = c = \frac{p_2 - p_0}{2}, \quad \quad f'(1) = 3a + 2b + c =\frac{p_3 - p_1}{2}$$

이를 이용해서 $a,b,c,d$를 구하면,

$$a= \frac{1}{2}( -p_0 +3 p_1 -3p_2 +p_3), ~b = \frac{1}{2}(2p_0 -5 p_1 + 4p_2 -p_3),~ c=\frac{1}{2} (-p_0 +p_2), ~d=p_1. $$

따라서 cubic interpolation 함수는

\begin{align} f(x) &=\frac {1}{2}\big(- x^3 + 2x^2 -x\big) p_0+\frac {1}{2}\big(3x^3 - 5x^2 +2 \big) p_1\\&+\frac {1}{2}\big( -3x^3 + 4 x^2 + x\big) p_2 +\frac {1}{2}\big( x^3 - x^2 \big) p_3 \end{align}

로 쓰인다. 이는 다음처럼 정의된 kernel 함수를(Catmull-Rom filter) convolution한 형태로도 표현할 수 있다.

$$ 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 & \text{otherwise} \end{array} \right. $$

$$  f(x) = p_0 K(-1-x) + p_1 K(-x) + p_2 K(1-x) + p_3 K(2-x),\quad (0\le x <1)$$

여기서는 좀 더 간단히 표현하기 위해서 $$ f(x) = m_0 p_0 + m_1 p_1 + m_2 p_2 + m_3 p_3 $$ 로 표현한다. 물론 $m_i$는 $0\le x <1$의 함수이다. $f(x)$는 그림에서 보는 것처럼 $p_1$과 $p_2$를 통과한다.

구간 끝에서 도함수를 $$f'(0) = \frac{p_2 - p_0}{2}, ~f'(1) = \frac{p_3 - p_1}{2}$$로 설정했기 때문에 인접 구간(eg. $[p_2,p_3]$)에서의 interpolation과 경계에서 같은  도함수를 가지게 되어 곡선이 꺽임없이 부드럽게 연결이 되는 특성을 가진다(예를 들면, $p_2$와 $p_3$을 보간하는 cubic interpolation함수는 $p_2$에서 미분값이 $(p_3 - p_1)/2$인데, 이는 $p_1,p_2$구간의 interpolation  함수가 $p_2$에서 가지는 미분값과 같다).

 

이를 2차원 공간으로 확장하면, 평면 위의 격자점 $\{(x, y)| x, y=-1,0,1, 2\}$에서 16개의 값 $\{ p_{ij} |i, j=0,1,2, 3\}$가 주어진 경우 정사각형 영역 $ 0\le x, y \le 1$에서 bicubic interpolation은 먼저 4개의 행 $y=-1,0,1,2$에서 각각 $x$ 방향의 cubic interpolation 된 값을 구한다:

\begin{align} q_0 =f(x, y=-1)= m_0 p_{00} + m_1 p_{10} + m_2 p_{20} + m_3 p_{30} \\ q_1 =f(x, y=0)= m_0 p_{01} + m_1 p_{11} + m_2 p_{21} + m_3 p_{31} \\ q_2 =f(x, y=1)= m_0 p_{02} + m_1 p_{12} + m_2 p_{22} + m_3 p_{32} \\ q_3 =f(x, y=2)= m_0 p_{03} + m_1 p_{13} + m_2 p_{23} + m_3 p_{33} \end{align}

4개 행에서 구해서 값을 이용해서 다시 $y$ 방향으로 cubic interpolation 한 값을 구하면 bicubic interpolation이 완성된다: $$ q =f(x,y)= q_0 n_0 + q_1 n_1 + q_2 n_2 + q_3 n_3$$ 여기서 $n_i$는 $m_i$에서 $x$를 $y$로 치환한 값이다.

 

원본 이미지를 확대/축소 또는 변형 등의 변환 과정을 거쳐서 출력 이미지를 만들 때 원본 이미지의 픽셀 값을 resampling 하는 과정이 들어온다. 이때 원본 이미지의 픽셀과 픽셀 사이의 값이 필요한 경우가 발생하여 적절한 interpolation이 필요한데 많이 쓰이는 interpolation 중의 하나가 bicubic interpolation이다. 아래는 bicubic interpolation을 이용한 이미지 resampling을 수행하는 코드다.

 

interpolation에서 실수(float) 연산을 하지 않고 정수 연산만 사용하면서도 일정한 정밀도로 소수점 아래 부분의 정보를 유지하기 위해 나누기 전에 미리 256을 곱한 후에 나눈다(256 대신에 적당히 큰 수를 선택해도 되는데 2의 지수승을 잡으면 곱셈/나눗셈을 shift 연산으로 대체할 수 있는 장점이 있다). 이렇게 하면 나눈 몫에 $\tt 0xFF$ 비트 마스크를 적용할 때 남는 부분이 소수점 아랫부분을 표현한다. 정밀도는 $1/256$ 까지다. 중간 과정에서 소수점 이하 부분끼리 곱셈을 할 때는 항상 $256$으로 나누어서 $255$ 이하가 되도록 만들어야 한다. 최종 결과에서도 다시 $256$으로 나누기를 해주어야 된다. bicubic인 경우는 $x/y$ 양방향으로 적용하므로 $256\times 256$으로 나누어야 하고 cubic interpolation 계수에서 $1/2$이 들어오므로 추가로 $4$만큼 더 나누어 주어야 한다(코드의 마지막 결과에서 shift 연산 "$\tt >> 18$"이 들어온 이유다).

 

bicubic interpolation을 적용할 때 $\tt y=0$이나 $\tt x=0$에서는 이전 행이나 열이 없으므로 자신을 반복하는 방식으로 처리해 주어야 하고, 또 $\tt y=rows-2, rows-1$이나 $\tt x=cols-2, cols-1$일 때도 비슷한 처리가 있어야 한다.

int resample_bicubic ( BYTE *src, int cols, int rows,
                       BYTE *des, int newcols, int newrows ) {
    if (cols < 4 || rows < 4)
        return resample_bilinear(src, cols, rows, des, newcols, newrows);

    int ixn = cols - 4;       
    BYTE *pa, *pb, *pc, *pd;
    for (int j = 0; j < newrows; j++) {
        int yy = ( ( j * rows ) << 8 ) / newrows;
        int yp = yy >> 8;                        // src pixel y-position;
        int dy = yy & 0xFF;
        int dy2 = (dy * dy) >> 8;
        int dy3 = (dy2 * dy) >> 8;
        int n0 = -dy3 + 2 * dy2 - dy;
        int n1 = 3 * dy3 - 5 * dy2 + 2 * 256;
        int n2 = -3 * dy3 + 4 * dy2 + dy;
        int n3 = dy3 - dy2;
        //
        if (yp == 0) pa = src + yp * cols;
        else         pa = src + (yp - 1) * cols;
        pb = src + yp * cols;
        if (yy < rows - 2) {
            pc = src + (yp + 1) * cols;
            pd = src + (yp + 2) * cols;
        } else if (yp < rows - 1) {
            pc = src + (yp + 1) * cols;
            pd = pc;
        } else {
            pc = src + yp * cols;
            pd = pc;
        }
        for (int i = 0; i < newcols; i++) {
            int xx = ( ( i * cols ) << 8 ) / newcols;
            int xp = xx >> 8;                    // src pixel x-position;
            int dx = xx & 0xFF;
            int dx2 = ( dx * dx ) >> 8;
            int dx3 = ( dx2 * dx ) >> 8;
            int m0 = -dx3 + 2 * dx2 - dx;
            int m1 = 3 * dx3 - 5 * dx2 + 2 * 256;
            int m2 = -3 * dx3 + 4 * dx2 + dx;
            int m3 = dx3 - dx2;
            int p = (xp == 0) ? 0 : (xp < ixn) ? xp - 1: ixn;    // p+3 <= ixn+3=cols-1;
            int a = ((m0 * pa[p] + m1 * pa[p + 1] + m2 * pa[p + 2] + m3 * pa[p + 3]) * n0 +
                     (m0 * pb[p] + m1 * pb[p + 1] + m2 * pb[p + 2] + m3 * pb[p + 3]) * n1 +
                     (m0 * pc[p] + m1 * pc[p + 1] + m2 * pc[p + 2] + m3 * pc[p + 3]) * n2 +
                     (m0 * pd[p] + m1 * pd[p + 1] + m2 * pd[p + 2] + m3 * pd[p + 3]) * n3) >> 18;
            *des++ = (a > 0xFF) ? 0xFF: (a < 0) ? 0: a;
        }
    }
	return 1;
}

bicubic interpolation을 하기 위해서는 4점이 필요하므로 폭이나 높이가 이보다 작은 경우는 bilinear interpolation을 사용한다. 다음 코드는 fixed-point bilinear interpolation을 구현한 코드다.

더보기
int resample_bilinear(BYTE *src, int cols, int rows,
                      BYTE *des, int newcols, int newrows ) {
    for (int j = 0; j < newrows; j++ ) {
        int yy = ((j * rows ) << 8) / newrows;
        int yp = yy >> 8;   // integer part;
        int dy = yy & 0xFF; // fractional part;
        for (int i = 0; i < newcols; i++) {
            int xx = ((i * cols ) << 8) / newcols;
            int xp = xx >> 8;
            int dx = xx & 0xFF;
            int pos = yp * cols + xp; //src position;
            int p00 = *(src + pos);
            int p10 = *(src + pos + 1);
            int p01 = *(src + pos + cols);
            int p11 = *(src + pos + cols + 1);
            int val = ((p11 * dx + p01 * (256 - dx)) * dy
                    + (p10 * dx + p00 * (256 - dx)) * (256 - dy)) >> 16;
            *des++ = val > 255 ? 0xFF: val < 0 ? 0 : val;
        }
    }
    return 1;
}

원본 컬러 이미지
2.5배 확대 영상

kipl.tistory.com/55

 

Bicubic Interpolation

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

kipl.tistory.com

 

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

Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (0) 2021.02.15
Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Posted by helloktk

댓글을 달아 주세요

728x90

Distance Transform은 이미지 상의 각 픽셀에서 가장 가까이 놓인 물체(전경)의 경계까지 거리를 구하는 연산이다. 보통의 경우에 이진 영상에만 적용이 된다. 구해진 거리의 정보를 이용해서 영상에 놓인 전경 물체의 skeleton을 구하거나 (OCR), 로봇의 팔 등을 자동 제어하는 데 있어서 움직이는 경로에 놓인 방해물의 위치를 예측할 수 있으며, 또한 물체까지의 경로 찾기(path finding) 등에 이용이 된다. distance transform은 사용되는 거리 측도(measure)에 따라 그 특성이 달라진다. 우선 가장 단순한 Euclidean 측도를 이용하도록 하자.

 

distance transform은 전경의 경계선을 모두 찾아서 배경의 각 픽셀에서 거리의 최솟값을 할당하면 된다. 그러나 경계가 많은 픽셀로 구성이 되는 이미지에서는 연산 비용이 매우 커지게 된다. 여기서는 이미지의 크기에 선형으로 비례하는 알고리즘에 대해서 살펴보도록 하겠다. 픽셀 수에 선형이기 위해서는 알고리즘이 raster 스캔 방식을 따라야 한다.

 

Euclidean 측도를 사용하는 distance transform은 한 점 $(x, y)$에서 전경 픽셀 $(i, j)$까지 거리의 제곱이 $z = (x-i)^2 + (y-j)^2$로 표현된다. 이미지 상의 각 전경 픽셀 꼭짓점으로 하는 paraboloid를 생각하면 distance transform은 한 픽셀에서 가장 밑에 위치한 paraboloid를 찾는 문제로 환원된다. 그러나 이차원 영역에서 이를 찾기는 쉽지 않은 작업이므로 이를 수평과 수직 방향으로 분리가 가능한지를 찾아보아야 한다. 우선 $x$-방향으로 스캔하면서 $x$-축으로 사영된 거리를 저장하면서 다시 $y$-축 방향으로 스캔하면서 저장된 $x$-방향의 거리의 제곱의 정보를 이용하면 된다. 이것은  $x$-방향으로 스캔할 때 거리의 제곱 값을 이미지의 그레이 값으로 저장하면 보다 쉽게 해결이 된다. 그러면 $y$-방향으로 스캔할 때는 꼭짓점이 이전에 $x$-방향 스캔 시 저장된 거리의 제곱만큼 이동이 된 이차곡선들 중에서 맨 아래에 놓인 것을 찾으면 된다.  따라서 2차원 distance transform은 1차원의 distance transform을 행방향 그리고 열방향으로 2번 시행을 하면 얻을 수 있다.

 

1차원의 distance transform은 일직선에 놓인 각 픽셀 위치를 꼭짓점으로 하는 이차 곡선들 중에서 가장 아래에 놓인 곡선을 선택하여 그 값을 취하면 된다. 거리의 최솟값을 구할 때 배경에 해당되는 픽셀은 선택되지 않기 위해서 배경 위치에서 이차곡선에 매우 큰 (수학적으로는 무한대) offset을 준다:

$$\text {1-dim distance transform}(p)=\underset {q\in\text {pixels}}{\text {min}}  \left( (q-p)^2+ F(q) \right)$$

초기 offset $F(q)$ 값은 전경 픽셀에서는 0, 배경 픽셀에서는 +INFINITY로 설정한다.  전경의 경계에서 형성이 된 2차 곡선 가장 낮게 놓이므로 그것에 의해서 거리가 결정이 된다.

일반적인 F(q)의 값이 주어지는 경우에는 문제는 각 픽셀점에서 보는 가장 낮게 놓인 이차곡선이 어느 것이고, 그것이 성립이 되는 영역을 찾는 문제로 귀결이 된다. 위치 $p$, $q$에 꼭짓점이 있는 이차곡선의 교차점이

$$s=\frac { (F(p) + p^2)-(F(q) + q^2) }{ 2(p-q)}$$ 표현되므로 쉽게 가장 낮은 2차 곡선의 영역을 찾을 수 있다.

 

아래 그림은 일차원 이미지의 distance transform의 결과이다: 붉은 선은 이미지(0:전경, 255:배경)이고, 녹색 선은 거리 값을 나타내는데 직관적으로 맞음을 확인할 수 있다.

네이버 블로그 이전

이 일차원 distance transform을 각각의 행에 대해서 한 후에, 그 결과를 다시 각각의 열에 또다시 적용하면 원하는 결과를 얻는다. 아래의 그림은 글씨에 대해서 skeleton을 구하기 위해서 distance transform을 적용하였다(skeleton을 얻기 위해서 반대로 함: 흰색=배경, 검은색=전경)

** 참고: en.wikipedia.org/wiki/Distance_transform

** cs.brown.edu/people/pfelzens/dt/에서 관련 논문과 원 코드를 찾을 수 있다;

#define INF     (1E+20)
#define SQR(a)	((a) * (a))
/* distanceTransform of 1d function using the Euclidean distance */
struct para_envelop {
    int x;		   // x-coord of the vertex of parabolar;
    double start;  // 현재 parabola가 가장 아래 높인 구간의 시작;
};
double para_intersection(int q, int p, double *dblimg) {
    return ((dblimg[q] + SQR(q)) - (dblimg[p] + SQR(p))) / (2. * q - 2. * p);
}
void DistanceTransform_1D(double *dblimg, int n, double* dist) {
    para_envelop *env = new para_envelop [n + 1] ;
    int k = 0;				// index of right_most parabola lower envelop;
    env[k].x = 0; env[k].start = -INF;
    env[k + 1].start = +INF;
    for (int q = 1; q < n; ++q) {
        double s = para_intersection(q, env[k].x, dblimg);
        // q에 꼭지점이 있는 parabola가 가장 밑에 있는 영역의 시작 위치(start)를 찾음;
        while (s <= env[k].start) s = para_intersection(q, env[--k].x, dblimg);
        env[++k].x = q ;    // 가장 아래에 있는 parabola 갱신;
        env[k].start = s ;
        env[k + 1].start = +INF;
    }
    k = 0;
    for (int q = 0; q < n; ++q) {
        while (env[k + 1].start < q) k++;//현재 위치에서 가장 낮은 parabola를 찾음;
        // Euclidean 거리 정보를 기록;
        dist[q] = SQR(q - env[k].x) + dblimg[env[k].x];
    }
    delete[] env;
}

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

Contrast Limited Adaptive Histogram Equalization (CLAHE)  (0) 2021.02.15
Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Posted by helloktk

댓글을 달아 주세요

728x90

FFT는 입력 데이터의 개수$(N)$가 2의 지수승으로 주어질 때  $O(N\log N)$의 연산만으로 빠르게 DFT을 수행하는 알고리즘이다. 주어진 $N$개의 data $\{F_0, F_1, ..., F_{N-1}\}$의  DFT $\{a_0, a_1,..., a_{N-1}\}$는

$$ a_k = \sum_{n = 0}^{N-1} F_n \exp\Big( -i \frac{2\pi nk }{N} \Big)=\sum_{n = 0}^{N-1} F_n (W^n )^k , \quad W = \exp\Big( -i \frac{2\pi}{N}\Big)$$ 로 표현된다. $N$이 2의 지수승으로 주어지면 이 합은 $n$이 짝수인 항과 홀수인 항으로 아래처럼 분리할 수 있다:

\begin{align} a_k &= \sum _{n=0}^{N/2-1} F_{2n} (W^{2n})^k+W^{k} \sum_{n=0}^{N/2-1} F_{2n+1}(W^{2n})^k\\ &=\text{(N/2개 even 항에 대한 DFT)} + W^k \times \text{(N/2개 odd 항에 대한 DFT)}\\ &= e_k + W^k o_k \end{align}

이는 원래의 문제가 처음 입력 데이터의 절반만 가지고 처리하는 두 개의 동일한 문제로 줄일 수 있음을 보여준다. 전형적인 divide and conquer 기법을 적용할 수 있는 구조이므로 매우 빠르게 연산을 수행할 수 있다(Cooley-Tukey algorithm). 더불어 주기성 때문에 $a_k$의 계산도 절반만 하면 된다: 

$$a_{k+ N/2} = e_k - W^k o_k, \quad k = 0,1,..., N/2-1.$$

역변환의 경우에는 twiddle factor $W$를 $W^*$로 바꾸면 되므로(전체 normalization factor $\frac{1}{N}$가 덧붙여진다) 코드에서 쉽게 수정할 수 있다: $\tt W_{im}  \rightarrow  - { W}_{im}$. 아래는 FFT를 재귀적으로 구현한 코드이다. 비재귀적 구현은 kipl.tistory.com/22에서 찾을 수 있다.

void split(int N, double* data) {
    double *tmp = new double [N / 2];
    for (int i = 0; i < N / 2; i++)   //copy odd elements to tmp buffers;
        tmp[i] = data[2 * i + 1];
    for (int i = 0; i < N / 2; i++)   // move even elements to lower half;
        data[i] = data[2 * i];            
    for (int i = 0; i < N / 2; i++)   // move odd elements to upper half;
        data[i + N / 2] = tmp[i];     
    delete [] tmp;
}
void fft(int N, double* re, double* im) {
    if (N < 2) return;
    else {
        split(N, re);                             // divide;
        split(N, im);        
        fft(N / 2, &re[    0], &im[    0]);       // conquer;
        fft(N / 2, &re[N / 2], &im[N / 2]);
        for (int k = 0; k < N / 2; k++) {
            double Ere = re[k];
            double Eim = im[k];
            double Ore  = re[k + N / 2];
            double Oim  = im[k + N / 2];
            double Wre =  cos(2. * PI * k / N);   // real part of twiddle factor:W
            double Wim =  -sin(2. * PI * k / N);  // imag part of twiddle factor:W
            double WOre = Wre * Ore - Wim * Oim;  // WO = W * O
            double WOim = Wre * Oim + Wim * Ore;
            re[k] = Ere + WOre;
            im[k] = Eim + WOim;
            re[k + N / 2] = Ere - WOre;
            im[k + N / 2] = Eim - WOim;
        }
    }
    return;
}

 6개의 주파수$\{f=2,5,9,11,21,29\text{Hz}\}$가 섞인 신호에서 [0,1]초 사이에 일정한 간격으로 sampling한  64개의 data을 FFT한 결과 (Nyquist frequency=32Hz);

$$s(t)=\sum_{i=0}^{5} (i+1) \cos(2\pi f_i t) $$

int fft_test_main() {
    const double PI = 4. * atan(1.);
    const int samples = 64;
    double re[samples], im[samples];
    const int nfreqs = 6;
    double freq[nfreqs] = {2, 5, 9, 11, 21, 29};
    for ( int i = 0; i < samples; i++ ) {
        double signal = 0;
        for ( int k = 0; k < nfreqs; k++ ) 
            signal += (k + 1) * cos ( 2.0 * PI * freq[k] * i / samples );
        re[i] = signal;
        im[i] = 0.;
    }
    fft ( samples, &re[0], &im[0] );
    for ( int i = 0; i < samples; i ++ ) TRACE ( "%d, %f, %f\n", i, re[i], im[i] );
    return 0;
}

FFT결과:

real part (imaginary part = 0)

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

Fixed-Point Bicubic Interpolation  (0) 2021.01.19
Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Posted by helloktk

댓글을 달아 주세요

728x90

Edge를 보존하는 low-pass filter로 edge의 방향으로 픽셀 평균값으로 대체한다. 방향은 8방향 (45도) 기준이다. 에지의 방향은 그림처럼 중심 픽셀을 기준으로 주변의 8개를 픽셀을 이용해서 45도씩 회전하면서  차이(평균 변화율)를 계산하여 가장 작은 값을 주는 방향이 된다. 중심 픽셀이 line[1][1]일 때

     0도 에지: abs(line[1][0] - line[1][2])가 최소

   45도 에지: abs(line[0][0] - line[2][2])가 최소

   90도 에지: abs(|ine[0][1] - line[2][1])가 최소

 135도 에지: abs(line[0][2] - line[2][0])가 최소 

edge 방향이 정해졌으면 에지 방향의 3 픽셀의 평균값으로 중심 픽셀 값을 대체하면 에지는 보존하면서 평균 필터를 적용한 효과를 얻을 수 있다.

 

void smooth_ep_3x3(BYTE *src, int width, int height, BYTE* dst) {
    const int xend = width - 2, yend = height - 2;
    BYTE *line[3];
    line[0] = src; line[1] = line[0] + width; line[2] = line[1] + width;
    BYTE *dst_line = dst + width;        // starting dst row;
    for (int y = 0; y < yend; ++y) {
        for (int x = 0; x < xend; ++x) {
            int diff1 = line[1][x] - line[1][x + 2];
            if (diff1 < 0) diff1 = - diff1;
            int diff2 = line[0][x] - line[2][x + 2];
            if (diff2 < 0) diff2 = - diff2;
            int diff3 = line[0][x + 1] - line[2][x + 1];
            if (diff3 < 0) diff3 = - diff3;
            int diff4 = line[0][x + 2] - line[2][x];
            if (diff4 < 0) diff4 = - diff4;
            
            if (diff1 <= diff2 && diff1 <= diff3 && diff1 <= diff4)         //0-도
                dst_line[x + 1] = (line[1][x] + line[1][x + 1] + line[1][x + 2]) / 3;
            else if (diff2 <= diff3 && diff2 <= diff4)                      //45-도
                dst_line[x + 1] = (line[0][x] + line[1][x + 1] + line[2][x + 2]) / 3;
            else if (diff3 <= diff4)                                        //90-도;
                dst_line[x + 1] = (line[0][x + 1] + line[1][x + 1] + line[2][x + 1]) / 3;
            else                                                            //135-도;
                dst_line[x + 1] = (line[0][x + 2] + line[1][x + 1] + line[2][x]) / 3;
                
            dst_line += width;                              //move to next dst line;
        }
        // increases src line ptr;
        BYTE *tptr = line[2] + width;
        line[0] = line[1]; line[1] = line[2]; line[2] = tptr;
    }
};

9번 반복 적용 결과: 지문처럼 규칙적인 패턴이 나오는 경우는 Gabor 필터를 이용하면 보다 좋은 결과를 얻을 수 있다.

네이버 블로그 이전

iteration에 따른 correlation의 계산:

iter=0, corr w.r.t original=0.925144, corr w.r.t previous=0.925144
iter=1, corr w.r.t original=0.903661, corr w.r.t previous=0.985120
iter=2, corr w.r.t original=0.888224, corr w.r.t previous=0.994821
iter=3, corr w.r.t original=0.876582, corr w.r.t previous=0.997403
iter=4, corr w.r.t original=0.866254, corr w.r.t previous=0.998183
iter=5, corr w.r.t original=0.856620, corr w.r.t previous=0.998596
iter=6, corr w.r.t original=0.847365, corr w.r.t previous=0.998841
iter=7, corr w.r.t original=0.838400, corr w.r.t previous=0.998981
iter=8, corr w.r.t original=0.829703, corr w.r.t previous=0.999088

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

Distance Transform  (0) 2021.01.16
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Posted by helloktk

댓글을 달아 주세요

728x90

컬러를 양자화하는 방법으로 median-cut알고리즘 이외에도 Octree을 이용한 알고리즘도 많이 쓰인다. (페인트샵 프로를 보면 이 두 가지 방법과 더불어 간단한 websafe 양자화 방법을 이용한다). Octree는 8개의 서브 노드를 가지는 데이터 구조이다. 이것은 R, G, B의 비트 평면에서의 비트 값(각각의 레벨에서 8가지가 나옴)을 서브 노드에 할당하기 편리한 구조이다. 풀컬러를 이 Octree를 이용해서 표현을 하면 깊이가 9인 트리가 생성이 된다 (root + 각각의 비트 평면 레벨 = 8^8 = num of true colors). Octree  quantization은 이 트리의 깊이를 조절하여서 leaf의 수가 원하는 컬러수가 되도록 조절하는 알고리즘이다. 전체 이미지를 스캔하면서 트리를 동적으로 구성하고 난 후에 트리의 leaf의 수가 원하는 컬러수보다 작으면 각각의 leaf가 표현하는 컬러를 팔레트로 만들면 된다. 그러나 컬러를 삽입하면서 만들어지는 leaf의 수가 원하는 컬러의 수보다 많아지는 경우에는 가장 깊은 leaf 노드를 부모 노드에 병합시키는 작업을 하여서 (이러고 나면 그 부모 노드가 이젠 leaf가 된다) leaf 노드의 수를 줄이는 과정을 시행한다. 이 병합 과정은 원래의 leaf노드들의 컬러 값의 평균을 취하여 부모 노드에 전달한다. 따라서 Octree 알고리즘은 RGB 컬러에서 상위 비트가 컬러에 대표성을 갖도록 하위 비트 값을 병합하여서 컬러를 줄이는 방법을 이용한 것이다.
       
특정 RGB 비트 평면에서 R, G, B비트 값을 이용하여서 child-노드의 인덱스를 만드는 방법은

    child-node index = (R-비트<<2)|(G-비트<<1)|(B-비트) ;

 

     (R, G, B)=(109,204,170)의 경우;

      node |
      level | 0 1 2 3 4 5 6 7
     -------------------------
           R | 0 1 1 0 1 1 0 1
           G | 1 1 0 0 1 1 0 0
           B | 1 0 1 0 1 0 1 0
     ------|------------------
      child | 3 6 5 0 7 6 1 4
     index |
  
Octree 알고리즘은 컬러 이미지를 스캔하면서 트리를 동적으로 구성할 수 있으므로 median-cut 알고리즘보다도 더 메모리 사용에 효과적인 방법을 제공한다. 그리고 트리 깊이를 제한하여서 일정 깊이 이상이 안되도록 만들 수 있어서 빠른 알고리즘을 구현할 수 있다. 최종적으로  팔레트를 구성하는 컬러 값은 leaf노드에 축적된 컬러의 평균값이 된다 (일반적으로 leaf-노드가 나타내는 비트 값과는 불일치가 생긴다). 이 알고리즘은 재귀적인 방법을 이용하면 쉽게 구현이 가능하고, 또한 웹상에서 구현된 소스도 구할 수 있다.

 

참고 논문: "A Simple Method for Color Quantization: Octree Quantization." by M. Gervautz and W. Purgathofer 1988.

256 컬러 이미지(트리 깊이 = 5): RGB 값이 주어지면 전체 트리에서 해당 leaf을 찾아서 팔레트 인덱스를 얻는다. 이와는 다른 방법으로는 팔레트가 주어졌으므로 주어진  RGB 값과 가장 가까운 유클리디안 거리를 주는 팔레트의 인덱스를 할당하는 방법도 있다. 이 방법을 이용하면 아래의 결과와 약간 차이가 생긴다.

양자화된 컬러 이미지
L2 error;

** 네이버 블로그에서 이전;

source code(C++): web.archive.org/web/20050306011057/www.drmay.net/octree/

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

FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Posted by helloktk

댓글을 달아 주세요

728x90

컬러 영상을 처리할 때 가장 흔히 사용하는 컬러 표현은 RGB 컬러이다. 이것은 R, G, B에 각각 8-비트를 할당하여 256-단계를 표현할 수 있게 하여, 전체적으로 256x256x256=16777216가지의 컬러를 표현할 수 있게 하는 것이다. 그러나 인간의 눈은 이렇게 많은 컬러를 다 구별할 수 없으므로 24-비트 RGB 컬러를 사용하는 경우는 대부분의 경우에 메모리의 낭비와 연산에 오버헤드를 가져오는 경우가 많이 생긴다. RGB 컬러 영상을 R, G, B를 각각 한축으로 하는 3차원의 컬러 공간에서의 벡터(점)로 표현이 가능하다. 컬러 영상의 픽셀들이 RGB삼차원의 공간에 골고루 펴져 있는 경우는 매우 드물고, 많은 경우에는 이 컬러 공간에서 군집(groups)을 이루면서 분포하게 된다. 하나의 군(group)은 유사한 RGB 값을 갖는 픽셀들로 구성이 되므로, 이 군에 포함이 되는 픽셀들에게 (군의 중앙에 해당하는) 대표적인 컬러 값을 대체하면 그 군에 포함이 된 픽셀은 이젠 RGB공간에서 한 점으로 표현이 되고, RGB공간상에는 픽셀 수만큼의 점이 있는 것이 아니라, 대표 RGB 값에 해당하는 점만이 존재하게 된다. 따라서 적당한 Lookup테이블(colormap)을 이용하면, 적은 메모리 공간만을 가지고도 원본의 컬러와 유사한 컬러를 구현할 수 있다.

 

이제 문제는 원래의 컬러 공간을 어떻게 군집화 하는가에 대한 것으로 바뀌었다. 간단한 방법으로는 원래의 컬러 영상이 차지는 하는 RGB공간에서의 영역을 감싸는 최소의 박스를 찾고, 이것을 원하는 최종적으로 원하는 컬러수만큼의 박스로 분할을 하는 것이다. 그러나 박스를 어떨게 분할을 해야만 제대로 컬러를 나누고, 또한 효율적으로 할 수 있는가를 고려해야 한다. 분할된 박스의 크기가 너무 크면 제대로 된 대푯값을 부여하기가 힘들어지고, 너무 작게 만들면 원하는 수에 맞추기가 어렵다.

 

Median-Cut 양자화(quantization)에서 사용하는 방법은 박스의 가장 긴축을 기준으로 그 축으로  projection 된 컬러 히스토그램의 메디안 값을 기준으로 분할을 하여서 근사적으로 픽셀들을 절반 정도 되게 분리를 한다 (한축에 대한 메디안이므로 정확히 반으로 분리되지 않는다). 두 개의 박스가 이렇게 해서 새로 생기는데, 다시 가장 많은 픽셀을 포함하는 박스를 다시 위의 과정을 통해서 분할을 하게 된다. 이렇게 원하는 수의 박스를 얻을 때까지 위의 과정을 반복적으로 시행을 하게 된다.

 

여기서 원래의 컬러 값을 모두 이용하게 되면 계산에 필요한 히스토그램을 만들기 위해서 너무 많은 메모리를 사용하게 되고 이것이 또한 연산의 오버헤드로 작용하게 되므로 RGB 컬러 비트에서 적당히 하위 비트를 버리면, 초기의 RGB공간에서의 histogram의 크기를 줄일 수 있게 된다.(보통  하위 3-비트를 버려서, 각각 5-비트만 이용하여, 전체 컬러의 개수를 32x32x32= 32768로 줄인다)

 

이렇게 RGB공간에서의 컬러 분포가 몇 개의 대표적인 컬러(예:박스의 중앙값)로 줄어들면(양자화 과정:: 공간에 smooth 하게 분포한 것이 몇 개의 점으로 대체됨), 원본 영상의 각각의 픽셀에서의 대체 컬러 값은 원래의 컬러와 가장 유사한, 죽 RGB 공간에서 유클리디안 거리가 가장 작은 박스의 컬러로 대체하면 된다. 

 

그러나 너무 적은 수의 컬러로 줄이게 되면 인접 픽셀 간의 컬러 값의 차이가 눈에 띄게 나타나는 현상이 생기게 된다. 이러한 것을 줄이기 위해서는 디더링(dithering) 과정과 같은 후처리가 필요하게 된다.

 

**네이버 블로그에서 이전;

 

source code: dev.w3.org/Amaya/libjpeg/jquant2.c

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

Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Binary Image에서 Convex Hull  (0) 2021.01.06
Posted by helloktk

댓글을 달아 주세요

728x90

Union-Find 알고리즘을 이용하여 영역 분할 과정이다. 각각의 픽셀에서 4방향으로 연결된 픽셀이 속한 영역 merge 할 수 있는지를 시도한다. merge 조건은 현재 픽셀의 그레이 값과 인접 영역의 평균 그레이 값의 차이가 주어진 임계값보다 작은 경우다.

$$\text{merge 조건: }\quad | 그레이 값- \text{인접 영역 평균값}| < threshold$$

컬러 영상의 경우에는 RGB 채널에 조건을 부여하거나 gray level만 이용해서 판단할 수 있다. root node 수만큼의 분할 영역이 만들어지고, 임계값이 클수록 분할된 각 영역의 사이즈가 커진다.  

gray = 0.2989 * r + 0.5870 * g + 0.1140 * b ;

같은 평균 픽셀값을 가지고 있더라도 4-방향으로 서로 연결된 영역이 아니면 합병되지 못하고 서로 다른 영역으로 남는다. 이를 하나의 영역으로 만들려면 분할된 영역을 다시 검사를 하는 과정이 필요하다. 자세한 과정은 Statistical Region Merge(SRM) 알고리즘(kipl.tistory.com/98)에 잘 구현이 되어 있다.

int FindCompress ( int *p, int node ) { // find root node;
    if ( p[node] < 0 ) return node;     // root-node;
    else return p[node] = FindCompress ( p, p[node] ); // 찾음과 동시에 depth을 줄인다;      
}
// root가 r1인 집합(단일픽셀)을 root가 r2인 집합으로 합병을 시도한다;
// 합병이 되면, 전체 픽셀값의 합과 픽셀수를 update한다;
BOOL MergePixel(int r1, int r2, int thresh, int *sums, int *sizes, int *p) {
    double diff = sums[r2] / (double) sizes[r2] - sums[r1] / (double) sizes[r1];
    if ( fabs(diff) < thresh ) { // merge two sets;
        sums[r2]  += sums[r1];
        sizes[r2] += sizes[r1];
        p[r1]     = r2;
        return TRUE;
    }
    return FALSE;
}
// root 노드는 분할영역의 픽셀 갯수와 픽셀 값의 합을 저장한다.
// 처음 각 픽셀이 root 이고, 픽셀 갯수는 1, 픽셀값은 자신의 픽셀값을 저장;
BOOL UnionFindAverage ( BYTE *image, int width, int height, int thresh ){
    int npixs = width * height;
    int* p      = new int [npixs];    // parent table;
    int* sums   = new int [npixs];  
    int* sizes  = new int [npixs];
    // initial setting;
    for ( int k = 0; k < npixs; k++ ) {
        p[k]     = -1;                // 각 픽셀이 root임;
        sums[k]  = image[k];
        sizes[k] = 1;
    }
    // 4-connectivity:
    // 첫행; LEFT 픽셀이 속한 영역의 평균값과 차이가 thresh 이내이면 left로 합병함;
    for ( int x = 1; x < width; x++ ) {
        int left = FindCompress ( p, x - 1 );
        int curr = FindCompress ( p, x );
        MergePixel(curr, left, thresh, sums, sizes, p);
    }
    for ( int y = 1, pos = y * width; y < height; y++ ) {
        // x = 0; up 픽셀이 속학 영역과 합병 체크;
        int up = FindCompress ( p, pos - width );
        int curr = FindCompress ( p, pos );
        MergePixel(curr, up, thresh, sums, sizes, p);
        pos++;
        // UP-LEFT 픽셀 영역과 합병 체크;
        for ( int x = 1; x < width; x++, pos++ ) {
            int left = FindCompress ( p, pos  - 1 );
            int curr = FindCompress ( p, pos );
            // left와 합병이 되면 left가 up과 합병이 되는지 확인; 안되면 up과 합병시도;
            if (MergePixel(curr, left, thresh, sums, sizes, p)) curr = left;
            int up = FindCompress ( p, pos - width );
            if (curr != up) MergePixel(curr, up, thresh, sums, sizes, p);
        }
    }
    // 평균 이미지를 만듬(input is overwritten);
    for ( int k = 0; k < npixs; k++ ) {
        int pos = k;
        // find root node: p[root] = -1;
        while ( p[pos] >= 0 ) pos = p[pos];
        int a = ( int ) ( sums[pos] / ( double ) sizes[pos] + 0.5 );
        image[k] = a > 255 ? 255 : a;
    };
    delete [] p;
    delete [] sums;
    delete [] sizes;
    return TRUE;
};

네이버 블로그 이전
statistical region merge 알고리즘을 적용한 결과

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

Octree Quantization  (0) 2021.01.12
Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Binary Image에서 Convex Hull  (0) 2021.01.06
Kuwahara Filter  (2) 2020.12.28
Posted by helloktk

댓글을 달아 주세요

728x90

Otsu 알고리즘은 이미지를 이진화시키는데 기준이 되는 값을 통계적인 방법을 이용해서 결정한다. 같은 클래스(전경/배경)에 속한 픽셀의 그레이 값은 유사한 값을 가져야 하므로 클래스 내에서 픽셀 값의 분산은 되도록이면 작게 나오도록 threshold 값이 결정되어야 한다. 또 잘 분리가 되었다는 것은 클래스 간의 거리가 되도록이면 멀리 떨어져 있다는 의미이므로 클래스 사이의 분산 값은 커야 함을 의미한다. 이 두 가지 요구조건은 동일한 결과를 줌을 수학적으로 보일 수 있다.

이미지의 이진화는 전경과 배경을 분리하는 작업이므로 클래스의 개수가 2개, 즉, threshold 값이 1개만 필요하다. 그러나 일반적으로 주어진 이미지의 픽셀 값을 임의의 개수의 클래스로 분리할 수도 있다. 아래의 코드는 주어진 이미지의 histogram을 Otsu의 아이디어를 이용해서 nclass개의 클래스로 분리하는 알고리즘을 재귀적으로 구현한 것이다. 영상에서 얻은 히스토그램을 사용하여 도수를 계산할 수 있는 0차 cumulative histogram(ch)과 평균을 계산할 수 있는 1차 culmuative histogram(cxh)을 입력으로 사용한다. 

$$ th= \text {argmax} \left( \sigma^2_B = \sum_{j=0}^{nclass-1} \omega_j m_j^2 \right)$$

 

* Otsu 알고리즘을 이용한 이미지 이진화 코드: kipl.tistory.com/17

* 좋은 결과를 얻으려면 히스토그램에 적당한 필터를 적용해서 smooth하게 만드는 과정이 필요하다.

// 0 <= start < n;
double histo_partition(int nclass, double cxh[], int ch[], int n, int start, int th[]) {
    if (nclass < 1) return 0;
    if (nclass == 1) {
        int ws; double ms;
        if (start == 0) {
            ws = ch[n - 1];
            ms = cxh[n - 1] / ws;
        } else {
            ws = ch[n - 1] - ch[start - 1];             // start부터 끝까지 돗수;
            ms = (cxh[n - 1] - cxh[start - 1]) / ws;    // start부터 끝까지 평균값;
        }
        th[0] = n - 1;
        return ws * ms * ms;                            // weighted mean;
    }

    double gain_max = -1;
    int *tt = new int [nclass - 1];
    for (int j = start; j < n; j++) {
        int wj; double mj;
        if (start == 0) {
            wj = ch[j]; 
            mj = cxh[j];                    //mj = cxh[j] / wj;
        }
        else {
            wj = ch[j] - ch[start - 1];     //start부터 j까지 돗수;
            mj = (cxh[j] - cxh[start - 1]); //mj = (cxh[j] - cxh[start - 1]) / wj;
        }
        if (wj == 0) continue;
        mj /= wj;                           //start부터 j까지 평균;
        double gain = wj * mj * mj + histo_partition(nclass - 1, cxh, ch, n, j + 1, tt);
        if (gain > gain_max) {
            th[0] = j;
            for (int k = nclass - 1; k > 0; k--) th[k] = tt[k - 1];
            gain_max = gain;
        }
    }
    delete [] tt;
    return gain_max;
};

trimodal 분포의 분리;

class0: 0~th[0]

class1: (th[0]+1)~th[1],

class2: (th[1]+1)~th[2]=255

th[0]=103, th[1]=172  (th[2]=255)
th[0]=88, th[1]=176, th[2]=255

 

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

Median-Cut 컬러 양자화  (0) 2021.01.12
Union-Find 알고리즘을 이용한 영역분할  (0) 2021.01.11
Multilevel Otsu’s Thresholding  (0) 2021.01.09
Binary Image에서 Convex Hull  (0) 2021.01.06
Kuwahara Filter  (2) 2020.12.28
Moving Average을 이용한 Thresholding  (0) 2020.11.26
Posted by helloktk

댓글을 달아 주세요

728x90

사진 속에 포함된 물체를 판별하기 위해서는 보통 이미지를 이진화시켜 binary 이미지를 만들고, conneted component labeling 등을 이용하여 영상에서 겹치지 않는 다수의 물체들을 분리해 낸다. 분리된 물체의 외각을 감싸는 convex hull을 구하기 위해서는 우선 경계점들을 추적하여 폴리곤을 형성한 다음 이들의 convex hull 알고리즘을 이용함으로써 구할 수 있다. convex hull 알고리즘 중에서 chain-hull이나 Melkman 알고리즘은 점들이 일정한 규칙에 의해서 정렬되어 있다는 사실을 이용하여 빠르게 convex hull을 찾을 수 있게 한다.

 

이미지를 읽는 raster-scan 방식은 물체의 경계점들을 일정하게 정렬된 순서로 읽게 되므로(chain-hull에서도 동일한 순서가 들어오지만, 전체 점들을 다 얻은 다음에 다시 분류작업을 해야 하므로, 읽은 즉시 바로 convex  hull을 구성하는 알고리즘은 제공하지 못한다) on-site convex hull의 구성이 가능하게 된다. 아래의  코드는 binary raster image를 읽으면서 경계점들의 convex hull을 찾는다.

#define bg     0x00             /* foreground color (black) */ 
#define fg     0xFF             /* background color (white) */
//cross(AB, AC) > 0 if C is lying on the left side of edge(AB);
#define CCW(A,B,C) ((B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x))

/* find a convex hull of foreground object in a binary image; */ 
/* in some case L[0] = R[0], or L[ll] = R[rr] if first line or last line of object 
** is composed of a single point (수정본:: 배경 체크 부분);
** convex hull이 중복점을 갖는 경우를 허용하지 않을 때는 comment에 따라서 처리해야 한다.
*/ 
int RasterHull(BYTE *image, int w, int h, CPoint H[/* 2xh */]) { 
    CPoint *L = &H[0]; //stack of left-side hull; 
    CPoint *R=  &H[h]; //stack of right side hull; 
    int rr = -1, ll = -1, x; 
    CPoint q; 
    for (int y = 0 ; y < h; y++) { 
        for (x = 0; x < w; x++)
            if (image[y * w + x] != bg) break; 
        if (x == w) continue; 
        q = CPoint(x, y);
        while (ll > 0) {
            // q가 이전 edge의 오른쪽이면 convex hull의 꼭지점이 됨;
            if (CCW(L[ll - 1], L[ll], q) < 0) break; 
            else ll--; 
        } 
        L[++ll] = q; 
        for (x = w - 1; x >= 0; x--) //x=-1 never occurs; 
            if (image[y * w + x] != bg) break; 
        // 행에 싱글 픽셀 있는 경우 왼쪽과 오른쪽이 점이 중복이 되는 경우가
        // 생긴다. 이 중복을 없애고 싶으면
        if (x <= q.x) continue ;
        q = CPoint(x, y); 
        while (rr > 0) { 
            // q가 이전 edge의 왼편에 있으면 convex hull의 꼭지점이 됨;
            if (CCW(R[rr - 1], R[rr], q) > 0) break; 
            else rr--; 
        } 
        R[++rr] = q; 
    } 
    // 왼쪽 처음에 점과 오른쪽 처음 점유 중복될 가능성이 있다(싱글 픽셀)
    // 왼쪽 마지막 점과 오른쪽 마지막 점이 중복이 되는 경우도 있다(싱글 픽셀) 
    // convex hull에서 중복점을 완전히 제외하고 싶으면 이를 체크하는 부분이 필요;
    for (int i = 0; i <= ll; i++) H[i] = L[i]; //left hull;
    // remove bottom degeneracy; 
    if (H[ll] == R[rr]) rr--;
    for (int j = rr; j >= 0; j--) H[i++] = R[j];//right hull; 
    //remove top degeneracy;
    if (i > 1 && H[0] == H[i - 1]) i--;
    return (i); 
};

위 함수는 픽셀 값이 bg가 아닌 것은 모두 전경으로 처리하도록 써져 있는데, 만약 connected component를 구한 상태에서 특정 component의 convex hull을 구하고자 하는 경우에는 if-문안의 조건을 image[y * w + x] == index로 바꾸어야 한다 (여기서 index는 connected component label이다. 그리고 labeling 이미지는 일반적으로 정수 이미지이므로 그에 맞게 수정을 해야 한다.) 

 

data matrix 코드를 포함하는 이미지를 이진화시킨 후 connected component labeling을 하여 작은 blob을 제거하면 거의 코드만 남길 수 있다. 그 후 위의 알고리즘으로  convex hull(붉은색 라인)을 구성한 것을 보인 것이다.

C네이버 블로그에서 이전

Convex hull 알고리즘은 아래에서 찾을 수 있다.

Posted by helloktk

댓글을 달아 주세요

Kuwahara Filter

Image Recognition 2020. 12. 28. 18:57
728x90

주어진 픽셀을 한 코너로 하는 4개 block 중에서 분산이 가장 작은 블록의 평균으로 중앙 픽셀 값을 대체하는 비선형 filter다. 이 필터를 적용하면  균일한 영역에서  spike noise 픽셀은 주변의 영역의 픽셀 값으로 대체되므로 제거될 수 있고, edge 라인을 기준으로 수직방향으로는 거의 균일한 영역이므로 필터링을 하더라도 edge가 뭉개지지 않음을 알 수 있다. 즉, Kuwahara filters는 이미지의 노이즈를 제거하고 edge를(위치) 보존하는 특성을 가진다.  윈도의 크기가 5x5인 경우 4개의 3x3영역에서 평균과 분산을 계산해서 중앙 픽셀의 값을 결정한다. integral 이미지를 사용하면 연산 부담을 줄일 수 있다.

원본 이미지
7x7 필터 적용 결과

이 Kuwahara filter보다도 보다 향상된 filter는 symmetric nearest neighbour (SNN) filter가 있다. 이 filter는 기준 픽셀에서 대칭 지점에 있는 두 지점의 픽셀 값 중에서 중앙 픽셀 값에 가까운 것을 취해서 평균값을 취한다.

 

아래 코드를 컬러 이미지에 적용하기 위해서는 RGB 채널별로 따로 적용하면 된다. 경계 영역은 추가적인 처리가 필요하다.

#define GET_SUM(x,y)  sum[ (y) * width + (x)]
#define GET_SUM2(x,y) sum2[(y) * width + (x)]
#define BLK_SUM(x,y) ((GET_SUM(x,y) + GET_SUM(x-hs1,y-hs1) \
                        - GET_SUM(x-hs1,y) - GET_SUM(x,y-hs1)))
#define BLK_SUM2(x,y) ((GET_SUM2(x,y) + GET_SUM2(x-hs1,y-hs1) \
                        - GET_SUM2(x-hs1,y) - GET_SUM2(x,y-hs1)))
void integral_image(BYTE *image, int width, int height, int *sum);
void integral_image_sq(BYTE *image, int width, int height, int *sum2);
void kuwahara_filter ( BYTE *image, int width, int height, int size, BYTE *out ) {
    int *sum = new int [width * height] ;
    int *sum2 = new int [width * height];
    // integral image; for window mean;
    integral_image(image, width, height, sum);
    // integral image of square; for window variance;
    integral_image_sq(image, width, height, sum2);
    
    // do filtering;
    size = ((size >> 1) << 1) + 1; // make window size be an odd number;
    int offset = ( size - 1 ) >> 1;
    int hs1 = ( size + 1 ) >> 1;
    int nn = hs1 * hs1;    // # of pixels within the  window;
    int var[4], mean[4];
    for (int y = hs1; y < height - hs1; ++y) {
        for (int x = hs1; x < width - hs1; ++x) {
            mean[0] = BLK_SUM(x, y) / nn;
            mean[1] = BLK_SUM(x + offset,          y) / nn;
            mean[2] = BLK_SUM(x,          y + offset) / nn;
            mean[3] = BLK_SUM(x + offset, y + offset) / nn;
            var[0] = BLK_SUM2(x,                   y) / nn - mean[0] * mean[0];
            var[1] = BLK_SUM2(x + offset,          y) / nn - mean[1] * mean[1];
            var[2] = BLK_SUM2(x,          y + offset) / nn - mean[2] * mean[2];
            var[3] = BLK_SUM2(x + offset, y + offset) / nn - mean[3] * mean[3];
            int vmin = var[0], id = 0;
            for (int i = 1; i < 4; ++i)
                if (var[i] < vmin) {
                    vmin = var[i]; id = i;
                }
            int a = mean[id]; 
            out[y * width + x] = a > 255 ? 255: a < 0 ? 0: a;
        }
    }
    delete[] sum;
    delete[] sum2;
};
Posted by helloktk

댓글을 달아 주세요

  1. hgmhc 2020.12.28 23:50 신고  댓글주소  수정/삭제  댓글쓰기

    이를 이용하면, 용량은 낮추되 원본과 비슷한 형태로 보존이 가능한 것인가요?

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

댓글을 달아 주세요