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

평면 위에 점집합이 주어지고 이들을 잘 기술하는 직선의 방정식을 구해야 할 경우가 많이 발생한다. 이미지의 에지 정보를 이용해 선분을 찾는 경우에 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

히스토그램처럼 균일한 간격으로 주어진 데이터에서 피크 값의 위치를 찾기 위해서는 먼저 노이즈 제거를 위해 몇 번의 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

Adaptive threshold를 적용하는 데 있어서 윈도 계산의 로드를 줄이는 방법은 integral image을 이용하면 된다. 물론 메모리의 소요가 부가적으로 발생하지만, 요 근래의 스마트 기기에서는 메모리는 별로 문제가 안된다.

아래의 코드는 integral 이미지를 이용해서 moving 윈도 내의 픽셀 평균 (= local average)을 기준으로 영상을 이진화시키는 함수다 (정확히는 "평균값 - 3"이다. 여기서 3은 바코드 인식 open library인 zbar에서 쓰는 기준을 잡았다. zbar library에서는 moving average를 구해 임계값으로 사용하는데, 윈도가 움직이면서 나가는 픽셀과 들어오는 픽셀을 업데이트하는 과정이 정확히 구현이 되어 있지는 않다. 그렇지만 근사적으로는 맞게 구현되어 있으므로 코드는 대부분의 경우 원하는 데로 잘 동작을 한다. integral image를 이용하면 윈도가 이동에 따른 픽셀 정보를 업데이트하는 복잡한 과정이 필요 없이 integral image의 단순 합/차만 수행하면 된다)

"윈도 평균-3" 대신 윈도의 표준편차를 이용할 수 있다. 그러나 이 경우에는 합의 제곱에 대한 적분 영상이 하나 더 필요하고, 얼마의 편차를 허용할 것인지를 정해야 한다. 이 기준에 맞게 구현된 코드는 http://kipl.tistory.com/30에서 찾을 수 있다.

2차원 바코드가 아닌 일차원 바코드 영상을 이진화할 때는 이만큼 복잡한(?) 알고리즘을 쓸 필요가 없다. 일차원 바코드는 보통 한 scanline의 정보만으로도 인식이 가능하므로 라인 단위의 이진화를 시키면 충분히다. 이 경우도 moving average를 사용하면 매우 간단하게 adaptive 한 임계값을 구할 수 있다. scanline 기준이므로 integral image는 따로 필요하지 않다.

void makeIntegralImage(BYTE *image, int width, int height, int* intImage);
더보기
void makeIntegralImage(BYTE *image, int width, int height, int* intImage) {    
    intImage[0] = image[0]; 
    for (int x = 1; x < width; ++x)
        intImage[x] = intImage[x - 1] + image[x];
    //next line;
    image += width;
    for (int y = 1, offset = y * width; y < height; ++y, offset += width) {
        int linesum = 0;
        for(int x = 0; x < width; ++x) {
            linesum += image[x];
            intImage[offset + x] = intImage[offset - width + x] + linesum ;
        }
        //next line;
        image += width ;
    }
}
/*
** moving window의 중심에 해당픽셀을 놓을 필요는 없다; 
*/
void thresholdByIntegralImage(BYTE *image, int width, int height, int wsz, BYTE *matrix) { 
    std::vector<int> intImage(width * height);
    makeIntegralImage(image, width, height, &intImage[0]);
    const int winArea = wsz * wsz ;
    /* const int wsz = 10;*/
    for (int y = 0, offset = 0; y < height; y++, offset += width) {
        int top = y - (wsz >> 1) ;
        if (top < 0 ) top = 0;
        else if (top > height - wsz) top = height - wsz;
        int bottom = top + wsz - 1;
        // y-range = [top, bottom];
        for (int x = 0; x < width; x++) {
            int left = x - (wsz>>1);
            if (left < 0) left = 0;
            else if (left > width - wsz) left = width - wsz;
            int right = left + wsz - 1;
            // xrange = [left, right];
            //
            int sum1 = (left > 0  && top > 0) ? intImage[(top - 1) * width + left - 1] : 0;
            int sum2 = (left > 0) ? intImage[bottom * width + left - 1] : 0;
            int sum3 = (top > 0) ? intImage[(top - 1) * width + right] : 0;
            //
            int graySum = intImage[bottom * width + right] - sum3 - sum2 + sum1;
            // overflow ? 
            // Threshold T = (window_mean - 3); why 3?
            if ((image[offset + x] + 3) * winArea <= graySum)
                matrix[offset + x] = 0xFF; //inverted!
            else
                matrix[offset + x] = 0x00;
        }
    }
}

 

QR 코드가 인쇄된 지면에 그라데이션이 있어서 전역 이진화로는 코드의 분리가 쉽지 않다.

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

Least Square Estimation of Perspective Transformation  (4) 2012.02.15
Perspective Transformation  (2) 2012.02.14
Integral Image을 이용한 Adaptive Threshold  (0) 2012.02.04
Peak Finder  (1) 2012.02.02
QR-code: decoder  (0) 2012.01.26
QR-code: detector  (0) 2012.01.12
Posted by helloktk

댓글을 달아 주세요

728x90

제대로 segmented 된 그레이 영상은 원래의 영상이 나타내고자 하는 전경이 잘 표현이 된 것이다. 이 경우의 원래 영상과 segmented 된 영상은 높은 상관관계를 갖는다. 따라서, 세그먼트를 위한 임계값의 설정 기준으로 이 상관계수를 최대로 하는 임계값을 찾는 것도 좋은 방법 중의 하나가 될 수 있다.

여기서 사용할 상관계수는 원래의 영상(A)과 전경과 배경을 그들의 픽셀 평균값으로 대체한 segmented 된 영상(B) 간의 상관계수를 사용한다. 임계값이 $T$인 경우 세그먼트된 영상 B 

$$B(i,j) = \left\{\begin{array}{ll} m_0, & \text{if}~A(i,j) \le T\\ m_1, &\text{otherwise}\end{array}\right. $$

로 나타난다. 여기서 $m_0$는 배경 픽셀의 평균값이고, $m_1$은 전경 픽셀의 평균값이다. 이 값은 임계값 $T$에 따라 달라진다. 임계값이 높으면 $m_0$는 커지고, 반대로 $m_1$은 작아진다

 

임계값이 $T$일 때 배경 픽셀 비를 $p$, 전경 픽셀 비를 $q(=1- p)$라 하면 segmented된 영상 B는 각 영역에서의 픽셀 값을 평균으로 대체했으므로 원본 영상의 평균과 같다. 또한, 원본 영상의 분산은 임계값에 무관하게 일정한 값을 유지한다. 이를 정리하면,

$$E(A)=E(B)=m=\text{pixel mean}=p m_0 + q m_1$$

$$V(A)=\text{variance} =T\text{-independent} = \text{const}$$

$$V(B)=pm_0^2 + q m_1^2 - m^2 = pq (m_0 - m_1)^2$$

$$E(A,B)= p m_0^2 + q m_1^2 $$

$$E(A,B) - E(A) E(B) = V(B)$$ 이므로, 

\begin{align}\text{Correlation}(A,B) &=\frac{ {E(A,B)-E(A)E(B)} }{\sqrt{V(A)V(B)} } \\ &=\frac{\sqrt{pq(m_0 - m_1)^2 } }{\sqrt{V(A)} }\\ &\propto \sqrt{pq(m_0 -m_1)^2 }\\ &=\sqrt{\text{interclass variance}}\end{align}

, 원래의 그레이 영상 A와 전경과 배경 픽셀을 각각의 평균값으로 대체한 영상간의 상관계수는 전경과 배경 두 클래스 간의 분산이 최대일 때 가장 크게 나타난다. 이 기준은 Otsu 알고리즘에서 사용한 기준과 같다.

 

참고: Otsu Algorithm 구현 예.

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

Is Pow of 2  (0) 2012.02.13
Fixed-point RGB2Gray  (0) 2012.01.25
Otsu-알고리즘의 새로운 해석  (0) 2010.01.28
Object Orientation  (1) 2010.01.17
Bicubic Interpolation  (1) 2010.01.14
Bezier Curve을 이용한 Histogram Smoothing  (0) 2010.01.10
Posted by helloktk

댓글을 달아 주세요

728x90

영상처리에서 한 픽셀의 수정을 위해서 주변 픽셀의 정보를 요구하는 윈도 기반 필터들은 일반적 연산 비용이 큰 편이다. 한 변의 길이가 W인 윈도를 사용할 때  W^2 횟수의 픽셀을 참조해야 하므로 윈도가 클수록 그 비용은 제곱으로 증가한다. 선형 필터 중에는 x-방향과 y-방향 연산으로 각각 분리가 가능한 경우는 연산 비용이 필터의 크기에만 비례하도록 만들 수 있다. median 필터는 이런 분리가 안 되는 비선형 필터 중 하나다. 근사적으로 x-방향으로 median filtering을 하고, 그 결과를 y-방향으로 다시 median filtering을 하는 방법으로 연산을 줄이는 방법을 사용하기도 한다.   

윈도가 움직이면서 윈도 내 모든 픽셀이 다 바뀌는 것이 아니라 움직이는 방향에 수직인 가장자리 픽셀만 바뀐다는 사실을 이용하면 각 픽셀의 윈도에서 median을 찾는 작업을 할 필요가 없다. 예를 들면 scanline 방향(x-방향)으로 윈도를 움직이면서 median filter를 작용할 때, 윈도가 오른쪽으로 1 픽셀 움직이면 윈도의 왼쪽 가장자리의 수직 픽셀들은 새 윈도에서 사라지고, 기존 윈도의 오른쪽 수직 가장자리 앞의 픽셀들이 새로 들어온다. 따라서 각 스캔라인에서 처음 한 번만 윈도의 median을 찾으면 이후에는 윈도가 이동할 때 윈도를 나가는 픽셀과 새로 들어오는 픽셀 (2*W) 개에 대해서 이전 median과 비교만 하면 된다. 이 방법은 비교 횟수가 윈도 크기에 1차적으로 비례하므로 연산 부담을 많이 줄일 수 있다. 이 방법은 사각형 모양의 윈도를 가지는 다른 필터(mean filter, max-filter, min-filter,...)에 대해서도 쉽게 적용할 수 있다.

// boundary region도 처리할 수 있게 수정함; 2021-04-18;
// window size = wx * wy;
// median = argmin(hist[i] >= halfArea)
int RunningMedianFilter(BYTE* image, int w, int h, int wx, int wy, BYTE* out) {
    int hist[256], x;
    int wxhalf = wx >> 1;
    int wyhalf = wy >> 1;
    wx = (wxhalf << 1) + 1;  // size of window = odd number;
    wy = (wyhalf << 1) + 1;
    for (int y = 0, yw = 0; y < h; ++y, yw += w) {
        // calc available area;
        int wystart = max(0, y - wyhalf);
        int wystop  = min(h, y + wyhalf);
        int wysize  = wystop - wystart + 1;
        int wxstart  = 0;
        int wxstop   = wxhalf;
        int halfArea = (wxstop * wysize + 1) >> 1;
        // to avoid *w multiplication in y-step;
        wystart *= w; 
        wystop  *= w;
        // initial histogram of topleft window;
        memset(hist, 0, 256 * sizeof(int));
        for (int iy = wystart; iy <= wystop; iy += w) 
            for (int ix = wxstart; ix <= wxstop; ++ix) hist[image[iy + ix]]++;
        // find initial median;            
        int ltmed = hist[0];       // less_than_median;
        int med = 0;
        while (ltmed < halfArea) ltmed += hist[++med];  
        out[yw + 0] = med;
        // left edge rgn;
        for (x = 1; wxstop < wx; ++x) {
            ++wxstop;
            halfArea = (wxstop * wysize + 1) >> 1;
            for (int iy = wystart; iy <= wystop; iy += w) {
                int v = image[iy + wxstop];     // (x=wxstop) strip enters;
                ++hist[v];
                if (v <= med) ++ltmed;
            }
            while (ltmed >= halfArea) ltmed -= hist[med--];  
            while (ltmed < halfArea)  ltmed += hist[++med];  
            out[yw + x] = med;
        }
        // central rgn;
        for ( ; wxstop < w; ++x) {
            ++wxstop;
            for (int iy = wystart; iy <= wystop; iy += w) {
                int v = image[iy + wxstart];    // (x=wxstart) strip leaves;
                --hist[v];
                if (v <= med) --ltmed;
                v = image[iy + wxstop];         // (x=wxstop) strip enters;
                ++hist[v];
                if (v <= med) ++ltmed;
            }
            ++wxstart;
            while (ltmed >= halfArea) ltmed -= hist[med--];
            while (ltmed < halfArea)  ltmed += hist[++med];
            out[yw + x] = med;
        }
        // right edge rgn;
        for ( ; x <= w; ++x) {
            for (int iy = wystart; iy <= wystop; iy += w) {
                int v = image[iy + wxstart];  // (x=wxstart) strip leaves;
                --hist[v];
                if (v <= med) --ltmed;
            }
            ++wxstart;
            halfArea = ((wxstop - wxstart + 1) * wysize + 1) >> 1;
            while (ltmed >= halfArea) ltmed -= hist[med--];
            while (ltmed < halfArea)  ltmed += hist[++med];
            out[yw + x] = med;
        }
    }
    return 1;
};

wx=wy=7;

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

Bicubic Interpolation  (1) 2010.01.14
Bezier Curve을 이용한 Histogram Smoothing  (0) 2010.01.10
Running Median Filter  (0) 2010.01.07
Fant's Resampling  (0) 2008.12.17
Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (2) 2008.07.14
Posted by helloktk

댓글을 달아 주세요

728x90

이미지를 이진화시키기 위해서 여러 알고리즘이 사용된다. 그중 이미지 전체에 대해 하나의 임계값으로 이진화시키는 전역 이진화 알고리즘은 간단하고 빠르기 때문에 많이 이용이 된다. 그러나 이미지를 형성할 때 조명 조건이 균일하지 않은 경우에는 전역 이진화는 원하는 결과를 얻기가 힘들다. 이런 경우에는 각각의 픽셀 주위의 그레이 값을 참조하여 임계치를 결정하는 국소적 이진화 방법을 사용한다. 국소적 이진화에서 임계값을 추출하는 간단한 방법은 윈도 내의 평균값을 이용하면 된다. 좀 더 개선된 알고리즘은 평균값($m(x, y)$)을 참조하되, 편차($\sigma(x, y)$)를 한번 더 고려해 주는 것이다. 이렇게 하여 잡은 국소적 임계값은 다음과 같이 표현된다: 

$$T_{(x, y)} = m_{(x, y)} [1+ \text{factor}(\sigma_{(x, y)}-128)]$$

여기서 $128$은 그레이 값이 가질 수 있는 최대 편차를 의미한다. 편차가 $128$이면 단순 평균값으로 취한다는 의미가 된다. 그 외의 경우는 표준편차와 128의 차이(항상 음수다)에 비례하는 값으로 윈도 평균값을 offset 한 값을 임계치로 잡는다. $\text{factor}$는 일반적으로 정해지지 않고, 실험적으로 $[0.2, 0.5]$ 사이의 값이 취해진다. (문서처럼 배경이 흰색인 경우는 $\text{factor} > 0$이지만, 검정 배경에 흰색 글씨를 처리하는 경우는 음수의 값을 취하는 것이 맞다)
 
국소적인 이진화 알고리즘은 매 픽셀마다 윈도를 잡아서 계산해야 하므로 연산 비용이 많이 든다. 충분한 메모리를 갖춘 시스템의 경우에는 적분 이미지(integral image)를 이용하면 윈도 연산에 소요되는 비용을 대폭 줄일 수 있다..

국소적 이진화 알고리즘에서 윈도 크기와 $\text{factor}$를 결정하는 기준은 무엇일까? 이것은 해결하고자 하는 문제의 특성, 예를 들면 스캔된 문서를 이진화시키는 경우에는 윈도에 충분한 글자가 들어 있어야 한다... 등에 많이 의존한다.

void make_int_img12(BYTE *gray, int width, int height, *int intimage, int *intsqimg);

더보기
void make_int_img12(BYTE *gray, int width, int height, *int intimage, int *intsqimg) {
    // first row accumulation;
    intimage[0] = gray[0];
    for (int x = 1; x < width; ++x) {
        int a = gray[x] ;
        intimage[x] = intimage[x - 1] + a;
        intsqimg[x] = intsqimg[x - 1] + a * a;
    }
    for (int y = 1, pos = y * width; y < height; ++y) {
        int linesum = 0, linesqsum = 0 ;
        for (int x = 0; x < width; ++x, ++pos) {
            int a = gray[pos];
            linesum   += a;
            linesqsum += a * a;
            intimage[pos] = intimage[pos - width] + linesum ;
            intsqimg[pos] = intsqimg[pos - width] + linesqsum;
        }
    }
};
#define integral_image(x, y) (intimage[(y) * width + (x)])
#define integral_sqimg(x, y) (intsqimg[(y) * width + (x)])
//
void adap_binariztion(BYTE *gray, int width, int height, 
                      int w       /*window size = 15*/,
                      double k    /*factor           = 0.2*/,
                      BYTE *bimage) {
    int whalf = w >> 1; //half of adaptive window;
    int diff, sqdiff;
    // make integral image && square integral image; 
    // if image is sufficiently large, use int64 or floating point number;
    std::vector<int> intimage(width * height) ;
    std::vector<int> intsqimg(width * height) ;

    //make integral image and its square integral image;
    make_int_img12(gray, width, height, &intimage[0], &intsqimg[0]);  
    //algorithm main;
    for (int j = 0, pos = 0; j < height; j++) {
        for (int i = 0; i < width; i++, pos++) {
            // clip windows 
            int xmin = max(0, i - whalf);
            int ymin = max(0, j - whalf);
            int xmax = min(width - 1, i + whalf);
            int ymax = min(height - 1, j + whalf);
            int area = (xmax - xmin + 1) * (ymax - ymin + 1);
            // calculate window mean and std deviation;
            if (!xmin && !ymin) {     // origin
                diff   = integral_image(xmax, ymax);
                sqdiff = integral_sqimg(xmax, ymax);
            } else if (!xmin && ymin) { // first column
                diff   = integral_image(xmax, ymax) - integral_image(xmax, ymin - 1);
                sqdiff = integral_sqimg(xmax, ymax) - integral_sqimg(xmax, ymin - 1);
            } else if (xmin && !ymin){ // first row
                diff   = integral_image(xmax, ymax) - integral_image(xmin - 1, ymax);
                sqdiff = integral_sqimg(xmax, ymax) - integral_sqimg(xmin - 1, ymax);
            } else{ // rest of the image
                int diagsum    = integral_image(xmax, ymax) + integral_image(xmin - 1, ymin - 1);
                int idiagsum   = integral_image(xmax, ymin - 1) + integral_image(xmin - 1, ymax);
                diff           = diagsum - idiagsum;
                int sqdiagsum  = integral_sqimg(xmax, ymax) + integral_sqimg(xmin - 1, ymin - 1);
                int sqidiagsum = integral_sqimg(xmax, ymin - 1) + integral_sqimg(xmin - 1, ymax);
                sqdiff         = sqdiagsum - sqidiagsum;
            }
            // threshold = window_mean *( 1 + factor * (std_dev/128.-1));
            // 128 = max_allowed_std_deviation in the gray image;
            double mean = double(diff) / area;
            double std  = sqrt((sqdiff - double(diff) * diff / area) / (area - 1));
            double threshold = mean * (1.0 + k * ((std / 128.0) - 1.));
            if (gray[pos] < threshold) bimage[pos] = 0;
            else                       bimage[pos] = 255;
        }
    }   
};

사용자 삽입 이미지

 

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

Fant's Resampling  (0) 2008.12.17
Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (2) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Posted by helloktk

댓글을 달아 주세요

  1. 박영우 2011.08.26 13:01  댓글주소  수정/삭제  댓글쓰기

    글 잘 읽었습니다.
    좋은정보 주셔서 감사합니다.
    그런데 적응이진화를 적용하니 좌측상단에 정사각형 모양으로 점이 하니 찍히네요
    이것을 제가 할수는 없는건가요??

  2. 진만시 2020.04.13 13:04 신고  댓글주소  수정/삭제  댓글쓰기

    시간을 단축시킬만한 방법이 있을까요 ?

728x90

이미지에서 어떤 유용한 정보를 추출하기 위해서는 이미지가 담고 있는 객체들을 분리하는 작업을 해야 한다. 가장 단순한 것 방법 중의 하나가 이진화(binarization)이다. 이진화는 이미지를 픽셀 값에 따라 0과 1(또는 255)로 값만 가지는 이진 이미지로 변환하는 과정이다. 이진화 작업을 거치면 이미지가 담고 있는 객체를 배경에서 분리할 수 있다. 이때, 어떤 기준으로 픽셀을 분리하는가에 대한 문제가 생긴다. 기준이 되는 임계값(threshold value)의 설정에 대한 다양한 알고리즘이 알려져 있는데, 그중에서 통계적인 방법을 이용한 Otsu 알고리즘이 자연스러운 임계값을 준다.

Otsu 알고리즘은 classification 기법을 이용하고 있다. 임계값을 설정하는 데 있어 비용함수를 설정하고 그 비용함수의 최솟값을 주는 값으로 임계값을 취하는 방식이다. 그럼 어떻게 비용함수를 설정할 것인가? 이미지에서 나타나는 픽셀 값을 2개의 클래스로 분리할 때, 좋은 분리는 각각의 클래스에 속한 픽셀 값의 분포가 유사해야 한다. 즉, 같은 클래스에 들어 있는 픽셀 값의 분산이 작아야 한다는 의미다. 따라서 비용함수는 픽셀 수의 비율로 가중치를 준 각 클래스의 분산을 합산한 것이 될 것이고, 임계값은 이 비용함수를 최소화하는 픽셀 값이다.

     비용함수 = (가중치1 * 분산1) + (가중치2 * 분산2) <= 2개 클래스로 분리 시
                 =   q1 * V1 + q2 * V2 ;      

              q1 =  전체 이미지에서 클래스1에 해당하는 픽셀이 나타날 확률
              q2 =  전체 이미지에서 클래스2에 해당하는 픽셀이 나타날 확률
              V1 = 클래스1에서 픽셀 값의 분산.
              V2 = 클래스2에서 픽셀 값의 분산.

     임계값  -->  MIN ( 비용함수 )

이미지의 픽셀 값 분포는 히스토그램으로 표현되므로, 임계값은 히스토그램으로 분리하는 레벨 값이고, 클래스 1은 그 값보다도 작은 부분, 클래스 2는 큰 부분을 의미한다. 비용함수의 의미를 좀 더 살펴보기 위해서 식을 바꾸어서 적으면

   비용함수 = 전체 분산 - (가중치1*(전체평균 - 평균1)^2 + 가중치2*(전체평균 - 평균2)^2);
                 = V - (q1 * (m1 - m)^2  + q2 * (m2 - m)^2) ;
                         
              V = 전체 분산;
              m = 전체 평균,
              평균1 (m1) = 클래스1의 평균,
              평균2 (m2) = 클래스2의 평균,

여기서 q1*(m-m1)^2 + q2*(m-m2)^2는 클래스들의 분산이다. 전체 분산은 어떤 식으로 클래스를 분리하더라도 항상 일정한 값을 가지므로, 비용함수를 최소화하는 것은 클래스들의 분산을 최대화하는 것과 같다. 새로운 비용함수(엄밀한 의미로 이득함수다)를 이것으로 잡으면
 
             이득함수 = q1 * (m1 - m)^2 + q2 * (m2 - m)^2;
             임계값 --> MAX (이득함수)
 
새로운 이득함수는 약간의 계산을 하면 그 의미가 더 명확한 표현인
             
             이득함수 = q1 * q2 * (m1 - m2)^2 ;

로 쓸 수 있다. 즉, 클래스 분리하는 값은 두 클래스의 평균의 차이(가중치를 갖는)를 최대화시키는 값으로 주어진다.

이 알고리즘의 구현은 히스토그램의 각 레벨에 대해서 좌우를 각각 클래스 1과 2로 설정한 후 이득함수를 계산하여 최댓값을 업데이트하는 방식을 취한다. 0-번째 cumulative histogram과 1-번째 cumulative histogram을 사용하면 각 클래스의 가중치와 평균을 쉽게 계산할 수 있다. 
    cumulative_histogram_0th[k] = 0...k 까지 값이 나타날 확률.
    cumulative_histogram_1th[k]/
cumulative_histogram_0th[k] = 0...k까지 값의 평균.


Otsu 알고리즘은 2개의 클래스뿐만 아니라 여러 클래스로 히스토그램을 분리하도록 확장할 수 있다. 재귀 알고리즘을 이용하면 쉽게 구현할 수 있다. (see: kipl.tistory.com/258)

/* Otsu임계화 예: 설명을 위해 최적화는 하지 않은 것이다. 반드시 cumulative 히스토그램을 이용할 
** 필요는 없다:*/
/* 이득함수의 최대값을 주는 레벨이 연속적으로 여러개 나타나면 평균값을 취하도록 한다(2016.04.26)
*/ 
int OtsuThreshold(BYTE *src, int width, int height, BYTE *dst) {
    double hist[256] = {0.}, chist[256] = {0.}, cxhist[256] = {0.};
    int ntot = width * height;

    // make histogram ;
    for (int i = 0; i < ntot; i++) hist[src[i]] += 1.;
    // normalize;
    for (int i = 0; i < 256; i++) hist[i] /= ntot;

    // make 0-th and 1-st cumulative histogram;
    chist[0] = hist[0];
    cxhist[0] = 0;
    for (int i = 1; i < 256; i++) {
        chist[i] = chist[i - 1] + hist[i] ;               //0-th cumulative histogram ;
        cxhist[i] = cxhist[i - 1] + double(i) * hist[i] ; //1-st cumulative histogram ;
    };
    
    double gain_max = 0.;
    int thresh = 0;   
    double m = cxhist[255];                     //total mean = q1 * m1 + q2 * m2;
    int mul_count = 1;                          //number of degenerate maxima;
    for (int i = 0; i < 256; i++) {
        if (chist[i] == 0.) continue ;
        double q1 = chist[i] ;                  //weight1;
        double q2 = 1 - q1;                     //weight2;
        if (q2 == 0.) break;
        double m1 = cxhist[i] / q1;             //mean1 ;
        double m2 = (m - cxhist[i]) / q2;       //mean2 ;
        double gain = q1 * q2 * (m1 - m2) * (m1 - m2) ;
        if (gain_max < gain) {
            gain_max = gain; 
            thresh   = i;
            mul_count = 1;                      //reset mul_count=1;
        } else if (gain_max == gain)            //degenerate case;
            mul_count ++;
    }
    if (mul_count > 1) thresh = thresh + (mul_count - 1) / 2;    //2016.04.26;

    // threshold image;
    for (int i = 0; i < ntot; i++) dst[i] = (src[i] >= thresh) ? 0xFF : 0x00 ;

    return thresh;
}

 

사용자 삽입 이미지

 

히스토그램 (계산된 임계치는 100이다)

사용자 삽입 이미지

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

Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (2) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Otsu Algorithm  (6) 2008.05.30
Hough Transform  (2) 2008.05.22
Posted by helloktk

댓글을 달아 주세요

  1. 감사^^ 2008.09.10 01:21  댓글주소  수정/삭제  댓글쓰기

    세미나 준비하는 중에 막막했던 부분이었는데

    너무 잘 봤습니다.

    근데 한가지 궁금한점이 있어서요...

    두 클래스의 분산의 합에 대한 수식이 두 클래스의 평균의 차이로

    변하는 과정을 알고 싶은데요...혹시 알수 있을까요?

  2. helloktk 2008.09.11 06:22  댓글주소  수정/삭제  댓글쓰기

    전체평균=m=q1*m1 + q2*m2(정의로부터); q1 + q2=1;을 이용하면,
    q1*(m1-m)^2 + q2*(m2-m)^2;
    =q1*(m1-q1*m1-q2*m2) + q2*(m2-q1*m1-q2*m2)^2;
    =q1*((1-q1)*m1-q2*m2)^2 + q2*((1-q2)*m2-q1*m1)^2 ;
    =q1*(q2*m1-q2*m2)^2 + q2*(q1*m2-q1*m1)^2;
    =q1*q2^2*(m1-m2)^2 + q2*q1^2(m1-m2)^2;
    =q1*q2*(q1+q2)*(m1-m2)^2;
    =q1*q2*(m1-m2)^2;
    을 얻을 수 있습니다.

  3. 공대생 2010.05.29 22:58  댓글주소  수정/삭제  댓글쓰기

    혹시 보게 되시면요...

    저도 bmp 이미지를 받아

    히스토그램은..

    void DibHistogram(CDib& dib, float histo[256])
    {
    register int i,j;

    int w = dib.GetWidth();
    int h = dib.GetHeight();

    BYTE** ptr = dib.GetPtr();
    //Histogram 계산
    int temp[256];
    memset(temp, 0, sizeof(int)*256);
    for(j=0; j<h; j++)
    for(i=0; i<w; i++)
    {
    temp[ptr[j][i]]++;
    }
    //Histogram 정규화
    float area = (float)w*h;
    for(i=0; i<256; i++)
    histo[i] = temp[i]/area;
    }


    이진화는 현재,,,

    void DibBinarization(CDib& dib)
    {
    register int i, j;

    int w = dib.GetWidth();
    int h = dib.GetHeight();

    BYTE** ptr = dib.GetPtr();

    for(j=0; j<h; j++)
    for(i=0; i<w; i++)
    {
    ptr[j][i]=(ptr[j][i]>50) ? 255:0;
    }

    }


    이런식의 이진화 중이에요.

    근데 위에거 보고 따라하기 식으로 해도 안되더라구요...

    아직 초보자라서요

    소스 좀 가르쳐 주시면 감사하겠습니다(제꺼에 적용할수 있겠끔)

    부탁드립니다.

    E-mail : rockentop@cyworld.com

  4. helloktk 2010.06.02 11:48 신고  댓글주소  수정/삭제  댓글쓰기

    binirization함수시작부분에 아래 부분을 추가하고, 50 대신에 thresh를 기준으로 하면 됩니다.

    double chist[256], cxhist[256];
    memset(chist, 0, sizeof(chist));
    memset(cxhist, 0, sizeof(cxhist));

    chist[0] = histo[0];
    cxhist[0] = 0;
    for(i=1; i<256; i++){
    chist[i] = chist[i-1] + histo[i] ;
    cxhist[i] = cxhist[i-1] + double(i) * histo[i] ;
    };

    double cost_max = 0;
    int thresh = 0;
    double m = cxhist[255]; //total mean ;
    for(i=0; i<256; i++){
    if(chist[i]==0.) continue ;
    double q1 = chist[i] ; //weight1;
    double q2 = 1 - q1; //weight2;
    if(q2 == 0.) break;
    double m1 = cxhist[i] / q1; //mean1 ;
    double m2 = (m - cxhist[i]) / q2; //mean2 ;
    double cost = q1*q2*(m1-m2)*(m1-m2) ;
    if(cost_max < cost) {
    cost_max = cost;
    thresh = i;
    };
    }

  5. 학생 2010.11.08 13:32  댓글주소  수정/삭제  댓글쓰기

    안녕하세요. 님의 귀중한 정보 덕분에 어느정도의 해결을 했습니다.
    저도 위에 공대생님처럼 bmp파일을 읽어서 이진화를 하려고하는데
    저는 bmp파일의 색상값이 아니라 다른 값으로 이진화를 하려고했지만
    히스토그램을 그려본 결과 봉우리가 3개가 나타났습니다. 원하는 이진화값은
    1,2번 봉우리 사이의 최소값이 나와야하는데 2,3번 봉우리 사이의 최소값이
    앞의 최소값보다 더 작아서 원하는 이진화가 안되고 있습니다.
    대략적으로 보면 앞의 최소값은 약 80정도이고 뒤의 최소값은 187인데
    뒤의 최소값을 없애던지 아니면 두개의 임계값을 찾을 방법은 없는것인지요..

    • helloktk 2010.11.11 09:53 신고  댓글주소  수정/삭제

      정확히 뭘 이야기하는지 잘 이해가 안되지만, 아마 님의 히스토그램에서 봉우리가 3개가 나타나고, 원하는 값은 1-2봉사이의 최소값 지점인 것 같은데. 이것이 2-3봉사이의 최소값보다도 커서 원하는 값이 안 되는(??) 상황같은데..OTSU방법은 이 최소값이 중요한 것이 아니라 전체적인 봉우리들의 모양이 중요합니다. 그리고 봉우리가 여럿있을 떄 반드시 봉우리 사이의 계곡에서 최소인지점에서 나타난다고 할 수 는 없습니다. 따라서 오수 알고리즘을 쓰면, 이 알고리즘이 주는 값에서 정해야 하고, 결과가 원하는 것이 아니면, 이 알고리즘을 쓰면 안되고 다른 알고리즘으로 구현하는 것이 더 적합할 것입니다. 이 알고리즘이 반드시 가장 좋은 결과를 주지는 않습니다.