주어진 히스토그램 데이터를 두 개의 가우시안 분포의 혼합으로 모델링하여 데이터를 분리하기 위한 decision boundary 값 (threshold value)을 expectation maximization 알고리즘을 적용하여 구한다. 

E-step: compute responsibility (of class 2; for class 1, 1-gamma_i)

M-step: compute the weighted means, variances and mixing probability

log-likelihood: 

void estimGaussParams(double data[], int start, int end, double *mean, double *var) {

코드보기

void initGuess(double data[], int n, double mean[], double var[], double *mixprob) {

코드보기

#define PI (4.0 * atan(1.))

double gaussDist(double x, double mean, double var) {   

코드보기

double responsibility2(double x, double mean[], double var[], double mixprob) {   

코드보기

double weightedMeanVar(double data[], double gamma[], int n, double mean[], double var[]) {

코드보기

#define EPSILON  1e-6

// Expectation Maximization algorithm applied to Two component Gaussian Mixture Model;

double emTwoCompGMM(double data[], int n) {

    double mean[2], var[2], mixprob;

    double *gamma = new double [n];     // responsibilities for class 2;

    initGuess(data, n, mean, var, &mixprob);

    

    // begin algorithm;

    while (1) {

        // E-step;

        for (int i = 0; i < n; i++) 

            gamma[i] = responsibility2(i, mean, var, mixprob);


        double old_mixprob mixprob;


        // M-step;

        mixprob = weightedMeanVar(data, gamma, n, mean, var);

        TRACE("mixing probability= %f\n", mixprob);


        // check convergence(usually loglikelihood is tested);

        if (fabs(mixprob - old_mixprob) < EPSILON)

            break;

    }

    // estimate decision boundary;

    for (int k = n - 1; k >= 0; k--) 

        if (gamma[k] < 0.5) break;


    delete [] gamma;

    return (2 * k + 1) / 2.; // = average of {k, k+1};

}



저작자 표시 비영리 변경 금지
신고
Posted by helloktk

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

영상의 주어진 그레이 값이 기준영상의 어떤 그레이 값으로 바뀌어야 하는가를 결정해야 하는데, 이때 기준은 주어진 영상의 i번째 그레이 값까지 픽셀이 나타날 확률합과 같은 확률합을 주는 기준영상의 그레이 값이 k일 때 

i(영상) -> k(기준영상)

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



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


//map을 lookup table로 이용해서 영상을 변환할 수 있다.

void matchHistogram(int srcHist[256], int refHist[256], int map[256]) {

    double cumSrc[256], cumRef[256];

    cumSrc[0] = srcHist[0];

    cumRef[0] = refHist[0];

    for (int i = 1; i < 256; i++) {

        cumSrc[i] = cumSrc[i - 1] + srcHist[i];

        cumRef[i] = cumRef[i - 1] + refHist[i];

    }

    // normalize ==> make CDF; 두 비교영상이 같은 크기가 아닐 수 있으므로 필요하다;

    for (int i = 0; i < 256; i++) {

        cumSrc[i] /= cumSrc[255];

        cumRef[i] /= cumRef[255];

    }

    // 두 비교영상의 cumulant histogram이 가장 가까운(같은) 픽셀끼리 변환한다;

    for (int i = 0; i < 256; i++) {

        int k = 255;  

        while (cumRef[k] > cumSrc[i]) k--;

        map[i] = k < 0 ? 0 : k;

        // 좀 더 smooth한 interpolation이 필요함;

    }

};

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

int histogramSpec(BYTE *src, int ws, int hw, BYTE *ref, int wr, int hr, BYTE *out) {

코드보기


저작자 표시 비영리 변경 금지
신고
Posted by helloktk

Tarjan의 Union-Find algorithm을 이용한 이미지상의 Connected Component Labeling 코드이다. 이미지의 배경은 검정(0)으로 하고 전경은 0이 아닌 값이면 된다. 
이 코드 자체는 좀 더 최적화를 시킬 수 가 있다. 예를 들면, 4-방향연결에서 LEFT와 TOP이 전경이면, 이 코드에서는 UNION을 두 번 호출한다. 그러나 생각해보면, 현위치의 라벨링을 LEFT와 같게 만들고 TOP과 UNION을 한 번만 호출하여도  TOP, LEFT,  현위치가 같은 parent를 갖음을 알 수 있다.

 

#define CCL_LEFT         (-1)
#define CCL_TOP          (-w)
#define CCL_TOPLEFT   (-1 - w)
#define CCL_TOPRIGHT  (1 - w)

#define CCL_BG            (-1)

 

#define UNION(x, y, label) {                    \
    int r0 = (x), r1 = (y), p;                        \
    while((p = label[r0]) != r0) { r0 = p; }    \
    r0 = p ;                                              \
                                                            \
    while((p = label[r1]) != r1) { r1 = p; }     \
    r1 = p ;                                              \
    if (r0 > r1)  label[r0] = r1;                     \
    else          label[r1] = r0;                     \
}                  

// #include <vector>  
int ConnectedComponentLabel(BYTE *image, int w, int h, int *out)
{
    /* 
    ** Computes "Connected Components" using Tarjan's Union-Find algorithm;
    ** the result is returned in the same buffer as gray_level image with values
    ** equal to the label of the component.
    ** Returns: the number of connected components 
    */

    /* if image !=0, it's FG */

    std::vector<int> label(w * h) ; 
    BYTE *q = &image[0] ;
    for (int i = 0, curpos = 0; i < h; i++) {
        for (int j = 0; j < w; j++, q++, curpos++) {
            if (*q) {
                label[curpos] = curpos;             
                if ((i > 0) && q[CCL_TOP]) {
                    UNION(curpos, curpos + CCL_TOP, label);         
                }
                if ((j > 0) && q[CCL_LEFT]) {
                    UNION(curpos, curpos + CCL_LEFT, label);
                }
#if defined (_EIGHTCONN_)
                if ((j+1 < w) && (i > 0) && q[CCL_TOPRIGHT]) {
                    UNION(curpos, curpos + CCL_TOPRIGHT, label);
                }
                if ((j   > 0) && (i > 0) && q[CCL_TOPLEFT]) {
                    UNION(curpos, curpos + CCL_TOPLEFT, label);
                }
#endif // _EIGHTCONN_
            }      
            else {
                label[curpos] = CCL_BG;     //BACKGROUND;
            }
        }
    }
    int curlab = 1; //출력용.1부터 시작한다; 0 은 bg;

    /* 라벨을 순차적으로 정리한다; 출력을 정수형으로 하는 경우에는 실제로 out이 필요없고, 

    ** label를 그대로 내보내도 된다. 단, CCL_BG 값을 갖는 픽셀을 0 으로 처리해야 한다;*/
    for (int pos = 0 ; pos < w * h; pos++) {
        int r = label[pos];
        if (r == CCL_BG ) {         // background ;
            out[pos] = 0;
        }
        else if (r == pos) {
            out[pos]  = curlab;
            label[pos] = curlab++;
        } 
        else { 
            label[pos] = label[r];
            out[pos]  = label[r];
        }
    }  
    return (curlab-1);
} ;

//http://blog.naver.com/helloktk/80026521024 에서 옮김.

                          <CCL을 이용하여 DataMatrix-code를 검출하는 과정>

 

 


20000 성분의 체스보드에 대한 4-방향 연결테스트; 출력은 정수형으로 바꾸어야 한다; 물론 8-방향 연결로 하면 1개 성분만 나온다.

 

 

저작자 표시 비영리 변경 금지
신고
Posted by helloktk

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


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



이제 원의 방정식은 xc , yc, d에 대한 1차 식으로 정리되었다. 삼각형의 외접원의 방정식은 세 꼭지점을 통과하므로 세 꼭지짐 A(x1, y1), B(x2, y2), C(x3, y3)을 위의 식에 넣으면 3개의 방정식을 얻는다. 이 방정식을 행렬로 표현하면


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


determinant를 정리하면,


을 얻을 수 있다. 이 식들의 분모는 세 꼭지점에 동일직선 위에 있지 않으면 0이 아니다. 한 꼭지점과 다른 꼭지점간의 좌표의 차이 4개 (x2-x1, y2-y1, x3-x1, y3-y1)와 좌표의 합 4개(x1+x2, y1+y2, x1+x3, y1+y3)를 계산하여서 중심점의 좌표를 얻을 수 있다. 

//

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

int circumcenter2(CPoint P, CPoint Q, CPoint R, double *xc, double *yc) {

double A = Q.x - P.x;

double B = Q.y - P.y;

double C = R.x - P.x;

double D = R.y - P.y;

double E = A * (P.x + Q.x) + B * (P.y + Q.y);

double F = C * (P.x + R.x) + D * (P.y + R.y);

double G = 2. * ( A * D - B * C);

if (G != 0.) { // 세 점이 일직선에 놓이지 않는 경우; 이 경우만 외접원이 정의된다;

*xc = (D * E - B * F) / G;

*yc = (A * F - C * E) / G;

return 1;

}

else return 0;

}

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

저작자 표시 비영리 변경 금지
신고
Posted by helloktk

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


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


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


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

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


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


따라서, 외접원의 중심은  



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


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


따라서, 

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


이므로, 이 값이 0이 아니기 위해서는 세 꼭지점이 일직선상에 놓여있지 않으면 된다.


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


또는, 성분으로 표현하면



이 설명은 공간상에 놓인 삼각형에 대해서도, 삼각형이 놓인 평면을 xy평면으로 생각하면 적용된다. 단지 조심해야 할 것은 m, n벡터를 구하는 것이다.


int circumcenter(CPoint A, CPoint B, CPoint C, double *xc, double *yc){

double ax = A.x - C.x ;

double ay = A.y - C.y ;

double bx = B.x - C.x ;

double by = B.y - C.y ;

double asq = ax * ax + ay * ay;

double bsq = bx * bx + by * by;

double ccw = ax * by - ay * bx;

if (ccw != 0.) { // 세 점임 일직선 위에 있지 않는 경우; 이 경우만 외접원이 정의됨;

*xc = C.x + (by * asq - ay * bsq) / (2. * ccw) ;

*yc = C.y + (-bx * asq + ax * bsq) / (2. * ccw) ;

return 1;

}

else return 0;

}



저작자 표시 비영리 변경 금지
신고
Posted by helloktk

세 개의 점 A, B, C가 만드는 삼각형의 외접원의 반지름을 구해보자. 


우선 세점이 일직선상에 놓이지 않기 위해서는 


이어야 한다. 

그림의 삼각형에서 각(BOC)가 각(BAC)의 두배이므로,


이 각도를 이용하면 삼각형의 면적은 세 변의 길이(a, b, c)와 외접원의 반지름(R)로 표현된다:


또한, 삼각형의 면적이 세 점이 직선상에 있는가를 테스트 하는 값으로 아래 식처럼 쓰여지므로,


외접원의 반지름은 다음과 같이 구할 수 있다:


#define SQR(x) ((x)*(x))

double circumradius(CPoint A, CPoint B, CPoint C) {

double ax = C.x - B.x; 

double ay = C.y - B.y;

double bx = A.x - C.x; 

double by = A.y - C.y;

double crossab = ax * by - ay * bx;

if (crossab != 0.) { 

double a = sqrt(SQR(ax) + SQR(ay));

double b = sqrt(SQR(bx) + SQR(by)); 

double cx = B.x - A.x; 

double cy = B.y - A.y;       

double c = sqrt(SQR(cx) + SQR(cy));

return (0.5 * a * b * c/fabs(crossab));

else

return ( -1.0 ); 

}


저작자 표시 비영리 변경 금지
신고
Posted by helloktk

conic section은 직교좌표계에서 (x, y)에 대한 2차식으로 쓰여진다;



이 conic section의 파라미터 {a,b,c,d,e,f}를 이용해서 타원의 중심과 장축 및 단축의 길이, 그리고 회전각에 대한 공식을 구해보자. 2차 항을 행렬을 써서 표현하면



따라서, 적절한 회전변환을 하면 두 좌표의 곱으로 주어지는 xy-항을 없앨 수 있다.  회전변환은 행렬의 determinant를 보존하고, 그 값은 행렬의 두 eigenvalue의 곱으로 쓰여진다.



따라서, 회전후의 방정식이 타원의 방정식(원점이 이동된)을 기술하기 위해서는 det > 0 이어야 한다 (회전시킨 후 식에서 x^2항과 y^2의 항의 계수는 각각 eigenvalue로 주어지므로 같은 부호를 가져야 한다.)



conic section F가 타원을 기술한다면, 다음과 같이 평행이동 및 회전(rigid transformation)변환을 해서 타원의 중심과 회전각도를 얻을 수 있다



이 식을 대입해서  xy-항이 없어지는 조건과 1차 항이 없어지는 조건을 사용하면


그리고, x와 y의 제곱의 계수는


이 eigenvalue는 위에서 구한 회전각도를 이용해서 표현이 가능하다:



나머지 상수항은



로 주어짐을 알 수있다. 따라서 타원의 각축의 길이는



로 주어진다.

// conic_params(a, b, c, d, e, f);

//ellipse_params(major_axis, minor_axis, center_x, center_y, tilt_angle);

bool conic_to_ellipse(double conic_param[6], double ellipse_param[5]) {

    const double a = conic_param[0];

    const double b = conic_param[1];

    const double c = conic_param[2];

    const double d = conic_param[3];

    const double e = conic_param[4];

    const double f = conic_param[5];

    //get ellipse orientation

    const double theta = 0.5 * atan2(b, a - c);


    //get scaled major/minor axes

    const double ct = cos(theta);

    const double st = sin(theta);

    const double ap = a * ct * ct + b * ct * st + c * st * st;

    const double cp = a * st * st - b * ct * st + c * ct * ct;


    //get translations

    const double cx = (2 * c * d - b * e) / (b * b - 4 * a * c);

    const double cy = (2 * a * e - b * d) / (b * b - 4 * a * c);


    //get scale factor

    const double val = a * cx * cx + b * cx * cy + c * cy * cy;

    const double scale_inv = val - f;


    if (scale_inv / ap <= 0 || scale_inv / cp <= 0) {

        //TRACE("Error! ellipse parameters are imagina\n");

        return 0;

    }


    ellipse_param[0] = sqrt(scale_inv / ap);

    ellipse_param[1] = sqrt(scale_inv / cp);

    ellipse_param[2] = cx;

    ellipse_param[3] = cy;

    ellipse_param[4] = theta;

    return 1;

};

저작자 표시 비영리 변경 금지
신고
Posted by helloktk

다섯개의 점들 (x1,y1), (x2,y2), (x3,y3), (x4,y4), (x5,y5)에 대해서 conic eq의 왼쪽값(대수적인 거리^2)이 최소가 되도록 하는 conic parameter를 결정하는 식은 다음과 같다: 

이 식은 해는 다음의 방정식을 만족시키는 해를 구하는 것과 같다:

또는 행렬식으로 표현하면,

의 형태로 쓰여진다. 6개의 conic-parameters를 결정해야 하므로, A는 마지막 행을 0으로 채워서 6x6 행렬로 만들어서 사용한다.  따라서 A(->A6x6)는 singular 행렬이 된다. A를 다음과 같이 singular value decomposition(SVD)를 할 때

가장 작은 eigenvalue(w)에 해당하는 V의 열벡터가 에러를 최소로 하는 conic parameter를 준다 (행렬 V의 각 열들이 6차원 벡터공간의 기저를 형성한다. U는 orthogonal 행렬이어서 오른쪽에서 곱해지는 벡터의 크기를 변화시키지 않고 단지 회전만 시킨다. 따라서 가장 작은 고유값에 해당하는 V의 열벡터가 |A.x|를  최소화시킨다)



// 주어진 점에서 타원까지의 대수적 거리: conic-eq의 좌변의 값;

static double fitting_error(double conic_params[6], double x, double y) {

더보기

//

static int ransac_sample_number(int ndata, int ngood, double prob, int ntry) {

더보기

int ransac_ellipse_fitter(std::vector<CPoint>& points, int width, int height,

                          double ellipse_params[5],/*(a, b, cx, cy, theta)*/

                          std::vector<int> &best_index)

{

    const int npts = points.size(); 

    if (npts < 5) {

        return 0;

    }

    // 안정적인 conic fitting을 위해서 데이터를 (-1,1)구간으로 normalize함;

    std::vector<CfPt> nor_pts(npts);

    double ellipse_scale; //원점에서  타원점들까지의 평균거리;

    CfPt ellipse_center;  //타원의 중심 예측;

    normalize_points(points, nor_pts, &ellipse_scale, &ellipse_center);

    //

    std::vector<int> inliers_index(npts);

    best_index.resize(npts);

    // inverse_scale을 기준으로 만들어졌음;

    const double dist_threshold = sqrt(3.84) / ellipse_scale;

    double best_params[5];

#define MAX_ITERATION (1500)

    int sample_num = 1000; //number of sample

    int ransac_iter = 0;

    int max_inliers = 0;

    while (sample_num > ransac_iter) {

        int rand_idx[5];

        //5 점을 랜덤하게 선택;

        random_quintet((npts - 1), rand_idx);

  //5 점을 이용해서 conic parameter estimate;

double conic_params[6];

get_conic_params(nor_pts, rand_idx, conic_params);


        //conic parameter를 이용해서 inlier후보들을 모음;

        int ninliers = 0;

        for (int i = 0; i < npts; i++) {

            double error = fitting_error(conic_params, nor_pts[i].x, nor_pts[i].y);

            if (fabs(error) < dist_threshold) {

                inliers_index[ninliers] = i;

                ninliers++;

            }

        }


        if (ninliers > max_inliers) {

            if (conic_to_ellipse(conic_params, ellipse_params)) {

                get_real_params(ellipse_params, ellipse_scale, ellipse_center, ellipse_params);

                // 주어진 영역안의 정상적인 타원만 선택함;

                if (regular_ellipse(ellipse_params, width, height, 0.5, 2.0)) {

                    max_inliers = ninliers;

                    //이 단계에서 inliers를 이용해서 다시 타원을 예측한다;

                    //if (ninliers >= 6) {

                    // std::vector<double> pts(2 * ninliers);

                    // for (int i = 0; i < ninliers; i++) {

                    // pts[2 * i + 0] = points[inliers_index[i]].x ;

                    // pts[2 * i + 1] = points[inliers_index[i]].y ;

                    // }

                    // EllipseFit(&pts[0], ninliers, conic_params);

                    // for (int i = 0; i < 5; i++) best_params[i] = conic_params[i];

                    //}

                    // if (ninliers == 5) 

                    for (int i = 0; i < 5; i++) best_params[i] = ellipse_params[i];

                    for (int i = 0; i < ninliers; i++) best_index[i] = inliers_index[i];

                    sample_num = ransac_sample_number(npts, ninliers, 0.99, 5);

                }

            }

        }

        ransac_iter++;

        if (ransac_iter > MAX_ITERATION) {

            TRACE("Error! ransac_iter exceeds!\n");

            max_inliers = 0;

            break ;

        }

    }


    if (best_params[0] > 0 && best_params[1] > 0)

        for (int i = 0; i < 5; i++) ellipse_params[i] = best_params[i];

    else

        max_inliers = 0;


    best_index.resize(max_inliers);

    return best_index.size();

}

저작자 표시 비영리 변경 금지
신고
Posted by helloktk

1. Derivative Based Algorithms:

* f = ∑ | I(x+1, y) - I(x, y)|,    only for   | I(x+1, y) - I(x,y)| > threshold;

* f = ∑ ( I(x+1, y) - I(x, y))^2,  only for  ( I(x+1, y) - I(x, y))^2 > threshold;

* f = ∑ ( I(x+2, y) - I(x, y))^2,  only for  ( I(x+1, y) - I(x, y))^2 > threshold;

* f = ∑ ( SobelX(x, y)^2 + SobelY(x, y)^2);

* f = ∑ ((C * I)(x, y))^2,          where, C = [ -1 - 4 -1; -4 20 -4; -1 -4 -1];

* f = ∑(|LaplaceX(x, y)| + |LaplaceY(x, y)|)

* f = (1/MN)∑ (GaussianDerivativeX(x,y)^2 + GaussianDerivativeY(x,y)^2),

                    where scale = d/2sqrt(3), d = dimension of the smallest feature;


2. Statistical Algorithms:

* Variance Measure ;

* Normalized Variance Measure ;

* Auto-correlation Measure:  f = ∑ I(x,y).I(x+1,y) - ∑ I(x,y).I(x+2,y)

* Standard Deviation-based Correlation: ∑ I(x,y).I(x+1,y) - H.W.*Average(I)


3. Histogram-based Algorithms:

* Range Algorithm: f = max{i|histo(i) > 0} - min{i| histo(i) > 0};

* Entropy Algorithm: f = - ∑ p(i) log2 p(i),    where p(i) = histo(i) / H.W;


4. Other Algorithms:

* Threshold Contents: f = ∑ I(x, y),       only for I(x, y) >= threshold;

* Thresholded Pixel Count: f = ∑ Step(threshold - I(x,y));

* Image Power: f = ∑ I(x, y)^2,       only for I(x, y) >= threshold;



Ref: Dynamic evaluation of autofocusing for automated microscopic analysis of blood smear and pap smear, J. Microscopy,227, 15(2007).                           


저작자 표시 비영리 변경 금지
신고

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

Union Find Connected Component Labeling  (0) 2012.11.01
RANSAC Ellipse Fitting  (0) 2012.10.07
Autofocus Algorithm  (0) 2012.06.03
Statistical Region Merging  (2) 2012.03.25
Savitzky-Golay Filter의 주파수 특성  (0) 2012.03.03
2차원 Savitzky-Golay Filters 응용  (0) 2012.02.28
Posted by helloktk
http://www.lix.polytechnique.fr/~nielsen/tpami04-nn.pdf
Abstract—This paper explores a statistical basis for a process often described in computer vision: image segmentation by region merging following a particular order in the choice of regions. We exhibit a particular blend of algorithmics and statistics whose segmentation error is, as we show, limited from both the qualitative and quantitative standpoints. This approach can be efficiently approximated in linear time/space, leading to a fast segmentation algorithm tailored to processing images described using most common numerical pixel attribute spaces. The conceptual simplicity of the approach makes it simple to modify and cope with hard noise corruption, handle occlusion, authorize the control of the segmentation scale, and process unconventional data such as spherical images. Experiments on gray-level and color images, obtained with a short readily available C-code, display the quality of the segmentations obtained. 

영상를 보면 인간은 비슷한 픽셀값을 모아서 몇 개의 영역으로 쉽게 나누어서 인식을 한다.  역으로, 우리가 보는  실제영상이 여러 개의 균일한 영역으로 나누어진 기본영상에 추가된 랜덤노이즈에 의해서 만들어졌다고 생각해 보자. 기본영상의 균일한 영역에서 랜덤변수는 일정구간(Q)에서 값을 취하여, 영역의 픽셀값과 더해져서 실제영상의 픽셀값 (0, ... g - 1 = 255)을 만든다. 이렇게 만들어진 실제영상에서의 두 영역의 픽셀평균값의 차이와 기본영상에서 랜덤변수에 의한 통계적 기대값의 차이는 주어진 0 < δ < 1 에 다음과 같은 하한을 갖음을 보일 수 있다:

따라서, 이 하한보다 작은 경우에서 두 영역은 하나로 인식될 수 있다.


#region count = 576; 많은 영역이 1픽셀로 구성이 되었있다;
#region count > 1 = 302;
#region count > 2 = 222;
#region count > 3 = 179;
#region count > 4 = 140;



저작자 표시 비영리 변경 금지
신고

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

RANSAC Ellipse Fitting  (0) 2012.10.07
Autofocus Algorithm  (0) 2012.06.03
Statistical Region Merging  (2) 2012.03.25
Savitzky-Golay Filter의 주파수 특성  (0) 2012.03.03
2차원 Savitzky-Golay Filters 응용  (0) 2012.02.28
webcam용 QR code detector  (0) 2012.02.19
Posted by helloktk
관련 논문들
1: http://www.hpl.hp.com/techreports/2010/HPL-2010-109.pdf 

2: http://inst.eecs.berkeley.edu/~ee123/fa11/docs/SGFilter.pdf
저작자 표시 비영리 변경 금지
신고
Posted by helloktk
Savitzky-Golay 필터는 일차원의 데이터에 대해서 일종의 이동평균을 취하는 경우와 동일하게 동작하는 필터이지만, 추정하는 지점의 주변의 모든 점에 동일한 가중치를 주는 방식(이동평균)을 택하지 않고, 그들을 보간하는 다항식을 최소자승법으로 찾아서 해당 지점의 값을 추정하는 방식을 택한다(frequency domain에서 분석을 하면 Savitzky-Golay 필터의 특성, 예를 들면, 피크의 위치나 등이 잘 유지되는 점등, 좀 더 다양하게 볼 수 있다). 이 필터를 쓰기 위해서는 몇 차의 다항식과 얼마의 윈도우 크기를 사용해야 하는지 설정을 해야 한다. (다항식의 찻수가 정해지면 최소의 윈도우 크기가 정해진다).
동일한 방식으로 이차원에 대해서도 Savitzky-Golay를 적용할 수 있다. 이 경우에는 다항식이 (x,y)의 2변수 함수 (2차원 평면에서의 정의되는 곡면)로 주어질 것이다. 이차원의 경우도 국소적인 필터로 사용하여서 영상의 smoothing 필터로 사용할 수 있지만, 필터의 윈도우를 영상전체로 잡아서, 즉 영상을 구성하는 픽셀값을 전 영역에서 보간하는 곡면을 찾을 수도 있다. 이렇게 찾은 곡면은 만들어진 영상의 배경조명이 균일하지 않는 경우에 이 추정된 곡면을 이용하면, 조명에 의한 효과를 예측할 수 있고, 배경조명이 보정된 영상을 만들어서 영상의 인식에 도움을 받을 수 있다. (문서인식에서 문서를 스캔할 때 생기는 균일하지 않은 배경이나, 2차원 코드 인식에서 배경의 추정등 다양한 부분에서 사용할 수 있다. 물론 간단한 경우에는 배경의 변화를 균일하게 기울어진 평면으로 근사를 하여서 추정할 수 있다)  

간단한 경우가 3차 다항식으로 영상을 보간하는 경우:

I(x, y) = a00
         + a
10*x + a01*y
         + a
20*x2 + a11*x*y + a02*y2
         + a
30*x3 + a21*x2*y + a12*x*y2 + a03*y3 
                                                      (x, y) ∈ image 


다항식은 x = [a00, a10,..., a03]T 의 10개의 필터계수를 추정하면 얻어진다. 추가적으로 Savitzky-Golay을 이용하면 영상의 미분값을 쉽게 구할 수 있다. 로컬버전의 필터인 경우에 필터적용값은 윈도우의 중심인 (x, y) = (0, 0)에서 다항식 값인 a00이다. 이 지점에서 x-방향의 편미분값은 a10, y방향의 편미분 값은 a01, 식으로 미분값을 구할 수 있다.

필터의 계수 x는 최소자승법 적용하면 얻어질 수 있다.  위의 다항식에  N (= width * height) 개의 픽셀로 구성된 영상의 각각의 픽셀에서의 좌표와 픽셀값을 대입하면,  N개의 식을 얻는다. 이것을 행렬식으로 쓰면

A.x = b

A는 N x 10 의 행렬로 각 행은 픽셀의 좌표로 구해진다:

A = [1, x0,  y0,  x02,  x0*y0,  y02,  x03,  x02*y0,  x0*y02,  y03]
     [1, x1,  y1,  x12,  x1*y1,  y12,  x13,  x12*y1,  x1*y12,  y13] 
     [1, x2,  y2,  x22,  x2*y2,  y22,  x23,  x22*y2,  x2*y22,  y23] 
     [1, .....................................................................] 
     [1, .....................................................................] 
                       ......................................
     [1, .....................................................................] 

여기서, 영상을 읽을 때, i-번째의 픽셀 위치가 (xi, yi) 로 주어진 경우다.

b 는 N-(열)벡터로 각각의 픽셀 위치에서 픽셀 값을 나타내는 벡터이다:

b = [ I(x0, y0) ]
     [ I(x1, y1) ] 
     [ I(x2, y2) ] 
     [ ............] 
     [ ............]
      ..............
     [.............]  


최소자승법을 적용하면, 추정된 다항식의 계수벡터 x는 

x = (A
T .A)-1.AT.b

임을 알 수 있다.

이렇게 추정된 2차원 곡면은 영상에서 추정된 배경의 픽셀 값 분포를 의미한다. 문자인식의 예를 들면, 보통의 경우에 흰 배경에 검정색 활자를 인식한다. 스캔된 영상에 검정색 활자들 때문에 추정된 곡명은 일반적으로 주어진 픽셀이 만드는 곡면보다도 낮게 된다. 픽셀 값이 추정된 곡면보다 더 낮은 픽셀들은 보통 검정색 문자들을 의미하므로, 이 차이의 평균값을 구하면, 대략적으로 어떤 픽셀이 배경에 속하는지 (곡면과 차이가 평균보다 작고, 또한 픽셀 값이  곡면의 아래에 놓인 경우), 아니면 문자영역인지(곡면과 차이가 평균보다 크고, 픽셀 값이 곡면의 아래에 놓인 경우)를 구별할 있게 된다.  
이제 이 정보들을 이용해서 추정을 다시 하는데 이번에는 1차 추정에서 글자영역으로 분류된 픽셀을 제외하고 배경을 추정하면 좀 더 정확한 배경을 기술하는 곡면을 얻을 수 있다.

로컬버전 필터로 사용할 때는 1차원에서와 마찬가지로 필터계수를 lookup table로 만들어서 사용할 수 있으나, 전영역을 대상으로 할 때는 행렬의 크기가 매우 커져서 연산도 많아진다. 

""

""

""


저작자 표시 비영리 변경 금지
신고
Posted by helloktk
웹캠의 영상에 들어오는 QR코드의 FinderPattern (topLeft, topRight, bottomLeft ==> 붉은색 십자)와 Alignment pattern(==>초록색 색자)을 찾고, 이들을 이용해서 perspective 변환을 구해서 code의 bounding box을 찾는다. Alignment Pattern을 못찾는 경우에는 bottomRight의 코너(==>초록색 십자)를 찾고, 그마저도 실패하면, FinderPattern 3개을 이용하여서 Affine 변환으로 bounding box을 찾는다.

웹  카메라:  이미지 형식: RGB24만 지원 (마이크로소프트의 vx-1000으로 테스트함)
                  이미지 크기: 640x480만 지원
                  영상을 보여주는 Callback함수 내에서 알고리즘을 호출하여서
                                       빠른 컴퓨터에서는 마크가 보이지 않을 수 있음.

테스트용이므로 상업적사용을 허용하지 않음.
 



실행화일 (detector + decoder 포함):


저작자 표시 비영리 변경 금지
신고
Posted by helloktk
32비트 머신에서 float 형의 변수는 4바이트의 메모리 공간을 차지한다. 즉, int 와 같은 메모리가 할당이 된다. 그리고, 4바이트의 최상위 비트가 1로 세팅이 되면, 이 float형의 변수는 음수를 의미한다. 따라서, float 형의 변수의 절대값을 구하고 싶으면, 메모리에 접근해서 4바이트를 얻어 온 후에, 최상위비트만 0으로 만들고 나머지는 그대로 두어야 한다. 따라서 비트마스크를 이용해서 이 과정을 수행하고 싶으면, and 연산의 비트마스크를 

01111111 11111111 11111111 11111111 = 2^31 - 1 = 0x7F FF FF FF 

처럼 잡으면 된다.

float  fast_abs(float fx) {
    int  ix = *(int *)&fx ;   // 4바이트 얻어 오기;
    ix &= 0x7FFFFFFF;    // MSB 지우기
    return *(float *)&ix ;   // float 형으로 되돌려 줌;
}

 
저작자 표시 비영리 변경 금지
신고

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

삼각형의 외접원 : 외접원의 반지름  (0) 2012.10.13
Ellipse Parameters  (0) 2012.10.13
float 타입 변수의 절대값은?  (0) 2012.02.17
x 보다 크거나 같은 가장 작은 2^n ?  (0) 2012.02.13
Is Pow of 2  (0) 2012.02.13
Fixed-point RGB2Gray  (0) 2012.01.25
Posted by helloktk

두 영상 사이의 perspective 변환은 8개의 매개변수(a,b,c,d,e,f,g,h)에 의해서 다음 식처럼 기술이 된다. (see, http://kipl.tistory.com/86)

또는, 

따라서, 매개변수를 찾기 위해서는 두 영상에서 서로 대응하는 점이 4개 이상 주어져 야 한다. 만약에 N개의 대응점들이 주어진 경우에


위의 식을 각각의 대응점에 넣어서 정리하면 아래의 행렬식을 얻을 수 있다.(좌변행렬의 마지막 열은 전부 -부호가 들어가야 한다) 
 

또는 간단히 

로 쓸 수 있다. 그러나 실제 대응점을 찾을 때 들어오는 noise로 인해서 실제 데이터를 이용하는 경우에는 정확히 등호로 주어지지 않는다. 따라서, 좌변과 우변의 차이의 제곱을 최소로 만드는 x를 찾아야 할 것이다.


에 대해서 미분을 하여서,

를 만족시키는 극값 x*을 구하면 된다. 는 8x8의 대칭행렬이어서 대각화가 가능하므로 역행렬을 구할 수 있다 (주어진 점들 중 한 직선 위에 놓이지 않는 점이 4개 이상이 있어야 한다). 따라서, 최소제곱해는 다음과 같이 주어진다:

.

저작자 표시 비영리 변경 금지
신고

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

2차원 Savitzky-Golay Filters 응용  (0) 2012.02.28
webcam용 QR code detector  (0) 2012.02.19
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
Posted by helloktk

한 평면에서 다른 평면으로 변환 중에서 직선의 직선성을 유지하는 변환은 2차원의 perspective 변환이다. 이 변환의 부분인 어파인 변환은 평행한 두 직선의 평행성을 그대로 유지하면서 변환시킨다. 따라서 이 perspective 변환은 사각형을 사각형으로 변환시키는 특성도 있다. 물론 사각형을 다른 사각형으로 변환시키는 다른 이차원 변환인 bilinear 변환이 있으나 일반적으로 직선의 직선성은 보전하지 못한다. 이 직선성의 보존은 매우 중요한 특성이다. 카메라도 일종의 perspective 변환기로 영상을 맺을 때도 직선은 그대로 직선으로 나타난다.(FOV가 큰 카메라는 렌즈왜곡이 심해서 보존이 안된다) 

평면에서의 변환을 다룰 때는 2x2행렬보다는 3x3 행렬을 이용하는 것이 더 편리하다. 이렇게 하면 평면에서 점의 이동을 행렬의 요소로 넣어서 생각할 수 있다.

(ex) affine 변환:
x = a11 * u + a21 * v + tu 
y = a12 * u + a22 * v + tv;
==> 
[x]   [a11   a21  tu] [u]
[y] = [a12   a22  tv] [v]
[1]   [   0      0   1] [1]

그리고 이것은 perspective 변환이 선형변환임을 명시적으로 보여주어서, 직선성이 보존된다는 사실이 자명해진다. 3x3행렬로 표현할 때, 평면의 좌표는 (x, y ,1)^T 처럼 3번째 좌표의 값은 항상 1로 고정한다( homogeneous coordinate) 
 
카메라로 물체를 촬영할 때, 가까운 거리에서 촬영을 하던, 먼 거리에서 촬영을 하던 두 영상은 크기차이만 있는 동일한 모양의 물체상을 만들어 낸다. perspective 변환은 3차원에 놓인 평면에서 평면으로 변환으로 생각할 수 있는데, 크기의 차이만 있는 경우에 같은 것으로 본다. 3차원에서 행렬변환은 9개의 매개변수에 의해서 기술이 되는데, 전체적인 크기의 차이를 무시하므로 1개 매개변수가 줄어들어서 8개의 매개변수로 표현이 된다.  
perspective 변환을 아래처럼 쓰면 변환된 좌표의 3번쨰 성분은 일반적으로 1이 아니다. 3번째 좌표 w을 구한 후에 이 값으로 x, y를 나누어서 생각하면 된다.

[x]   [ a11  a21 a31 ] [u]
[y] = [ a12  a22 a32 ] [v] 
[w]   [ a13  a23 a33 ] [1]             (a33 = 1)
 
x  = a11 * u +  a21 * v + a31     ==>  x =  x / w;
y  = a12 * u +  a22 * v + a32     ==>  y =  y / w;
w  = a13 * u +  a23 * v + a33  

perspective 변환행렬 [aij]는 4개의 점에 대응하는 출력영상에서의 4개의 점이 주어지면 8개의 방정식을 만들 수 있고, 이것을 이용해서 계수를 구할 수 있다. 그러나, 8차 방정식의 근의 공식이 없는 관계로 일반적으로 수치해석적으로 해결해야 한다. 그리고, 주어진 4점이 (입력 또는 출력) 일직선상에 있으면 답을 구할 수 없고, 그 중에 3개만 일직선상에 있는 경우에는 이 변환은 평행성을 보존하는 affine 변환이 된다.(affine은 6개의 매개변수로 기술되고, 평행이동을 빼면, 4개의 매개변수가 남는데, 4차 방정식의 근의 공식이 있으므로 답을 적을 수 있다)

다행이 정사각형에서 사변형으로 변환은 수치해석에 의존하지 않고도 답을 적을 수 있다.
(0,0) -> (x0, y0);               
(1,0) -> (x1, y1);                         
(1,1) -> (x2, y2);
(0,1) -> (x3, y3);
==>
denom = (x1 - x2) * (y3 - y2) - (x3 - x2) * (y1 - y2);     
a11 = x1 - x0 + a13 * x1 ;
a21 = x3 - x0 + a23 * x3 ;
a31 = x0 ;
a12 = y1 - y0 + a13 * y1;
a22 = y3 - y0 + a23 * y3;
a32 = y0;
a13 = ((x0-x1+x2-x3)*(y3-y2) - (x3-x2)*(y0-y1+y2-y3)) / denom;
a23 = ((x1-x2)*(y0-y1+y2-y3) - (x0-x1+x2-x3)*(y1-y2)) / denom;
a33 = 1.0;


따라서 일반적인 사변형에서 사변형으로의 변환은 

사변형1 --> 정사각형 --> 사변형2


의 2단계 변환의 곱으로 주어진다. 사변형에 정사각형으로 변환은 정사각형에서 사변형으로 변환의 역변환이므로 역행렬을 구해야 하나, 이것 보다는 수치적으로 안정적인 adjoint행렬을 이용하는 것이 낳다(adjoint을 쓰면 determinant로 나누기를 할 필요가 없다). 이것은 perspective변환에서 항상 좌표를 3번째 좌표로 나누어서 사용하기 때문에 가능하다.


 


저작자 표시 비영리 변경 금지
신고

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

webcam용 QR code detector  (0) 2012.02.19
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
Posted by helloktk
양의 정수  x가 주어질 때, 이보다 크거나 같은 가장 작은 2^n의 숫자를 찾아보자. 왜? FFT 때문이다. 물론,  n = (int)ceil(log(x) / log(2.)) 로 계산하거나, 
                    int a = 1;
                    while (a < x) a <<= 1; 
로 찾을 수 있다;

2^(n-1) < x <= 2^n 사이에 있으면, 원하는 답은 2^n이다. x=2^n 인 경우를 제외하고 모두 n번째 비트가 1로 세팅이 된다. 따라서, 통일시키기 위해서 1을 빼면, x-1 은 n 번째 비트가 항상 1로 주어진다.  이제, 남은 일은 n-1 번째에서 0번째까지 모든 비트를 1로 채우면 된다. 그 결과에 1을 더하면 원하는 2^n을 얻는다.  n-번째 비트가 1이고 이를 이용해서 하위비트를 채우기 위해서는 차례로 >> 연산을 한 후에 or 연산을 하면 된다.

a = x - 1
a                   = 1xxxxxxxxxxxxxxxxxxxxxxxx
a >> 1             = 01xxxxxxxxxxxxxxxxxxxxxxx
a >> 2             = 001xxxxxxxxxxxxxxxxxxxxxx 
.......

31번의 >> 연산을 한 후에 or 연산을 하면 n이하의 모든 비트가 1로 채워진다(31은 최대로의 경우다. 사전에 답을 모르므로 31번 >>까지 해야 한다). 이것은 너무 낭비다.  잘 살펴보면,

a | (a>>1) = 11xxxxxxxxxxxxxxxxxxxxxxxxx

형태이어서, 상위 두 자리 비트가 이미 세팅이 되어 있으므로, 이 결과를 다시  >> 2 하여서 이전 결과와 or연산을 하면 상위 4자리의 비트가가 채워지고, 또다시 이 결과를 >>4 하고 난 후에 직전 결과와 or 연산을 하면 상위 8자리의 비트가 1로 채워진다. 이런식으로 하면 >>16까지 하면 32자리를 모두 커버할 수 있다.

따라서, 전체적인 알고리즘은 (x > 0)
 

   x-- ; 
   x |= (x >> 1);  //상위 2 자리 채움
   x |= (x >> 2);  //상위 4자리
 채움 
   x |= (x >> 4);  //상위 8자리
 채움 
   x |= (x >> 8);  //상위 16자리
 채움 
   x |= (x >> 16);//상위 32자리
 채움 

   x++;
   return x;


단, 32비트 머신에서만이다. 64비트 프로그래밍에서는 >>32 도 추가로 해주어야 한다.
저작자 표시 비영리 변경 금지
신고

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

Ellipse Parameters  (0) 2012.10.13
float 타입 변수의 절대값은?  (0) 2012.02.17
x 보다 크거나 같은 가장 작은 2^n ?  (0) 2012.02.13
Is Pow of 2  (0) 2012.02.13
Fixed-point RGB2Gray  (0) 2012.01.25
A Fast 2D Otsu Thresholding Algorithm Based on Improved Histogram  (0) 2010.11.17
Posted by helloktk
FFT를 수행하려면 이미지의 폭이나 높이가 2의 지수승인 경우가 가장 단순하게 주어진다. 따라서  주어진 정수 x가 2의 지수승인가 판별할 수 있어야 한다. x가 양의 정수이고, 2^n으로 표현이 된다면, 2진수로 나타낼 때, (n+1) 번째 비트만 1이고 (0부터 센다), 나머지 비트는 모두 0이다. 그리고 x-1은  n 번째에서 0번째까지 모든 비트가 1이 된
다.  
    x                    x(2진수)            x-1
    1                      1                        0
    2                    10                         1
    4                   100                       11
    8                  1000                     111
     ................................................................
이 표를 보면, x와 x-1사이에는 겹치는 비트가 없다. 따라서 두 수를  and 연산을 하면 0 이 되는 경우에는 2의 지수승이고, 그 이외의 경우에는 0이 아님을 알 수 있다. x = 2의 지수승  판별은    
 

                           return  x & (x - 1) == 0         /* x != 0 인 정수*/
 

인가를 보면 된다. 

그런데 이 판별식은 x = 0 인 경우에는 성립이 안된다. x-1 = -1 이므로 32비트 자리 전부가 1로 채워지므로 x & (x-1) = 0 이어서 2의 지수승으로 판별한다. 따라서 0을 제외하는 방법을 찾아야 한다. (물론 함수 인자에서 양수로 제한을 하면 되지만 폼이 안난다). 음수는 최상위비트가 1로 채워지므로 이 사실을 이용하자. 최상위 비트를 1로 만들려면,

~0U                    =111111111111111..11111111(32개) 
~0U>>1              = 011111111111111..11111111
~(~0U>>1)          =100000000000000..00000000 
~(~0U>>1)|x       = x의 최상위 비트를 항상 1로 채워준다(음수 일때는 자동으로 만족)
                             나머지 비트는 그대로 둔다.

따라서 이 값과 x-1을 and 연산을 하면 0 이하인 수가 들어오면 연산이 결과를 항상 0이 아니게 된다.


                          return  !((~(~0U>>1)|x) & (x - 1))          /*x = 정수 */    
 
0 하나를 예외처리하기 위해서 너무 많은 과정을 거치는 것이 아닌가? 
저작자 표시 비영리 변경 금지
신고
Posted by helloktk
Adaptive threshold 방법을 적용하는 데 있어서 윈도우 계산의 로드를 줄이는 방법은 integral image을 이용하면 된다(물론 메모리의 소요가 부가적으로 발생하지만, 현재의 기기에서는 메모리는 별로 문제가 안된다).

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

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

2차원 바코드가 아닌 일차원 바코들 영상을 이진화할 때는 이 만큼 복잡한(?) 알고리즘을 쓸 필요가 없다. 일차원 바코드는 한 스캔라인의 정보만으로도 보통 인식이 가능하므로 라인단위의 이진화면 충분히다. 이 경우에도 이동평균을 사용하면 매우 간단하게, 그리고 adaptive한 임계값을 구할수 있는데, 라인기준이므로 적분영상이 따로 필요하지 않다.

void makeIntegralImage(BYTE *image, int width, int height, int* intImage) {

더보기

/*
** 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;
            }
        }
    }
}



저작자 표시 비영리 변경 금지
신고

'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

Peak Finder

Image Recognition 2012.02.02 21:41
주어진 데이터를 받아서 극대점과 극소점을 찾는다. 극점은 이전 극점과 일정한 차이(delta)가 날 때만 취한다. 맨 처음과 나중에 들어오는 데이터가 극점일 때, 별도의 처리가 필요하다.

BOOL peakFinder(std::vector<double>& data, double delta, /*IN*/
                          std::vector<int>& maxPosition, /*OUT*/
                          std::vector<int>& minPosition) /*OUT*/
{
    int maxPos = -1;
    int minPos = -1;
    bool lookForMax = true;
    double maxVal = -FLT_MAX ;       
    double minVal = FLT_MAX ;
    for (int i = 0; i < data.size(); i++) {
        double currVal = data[i] ;
        if (currVal > maxVal) {
            maxVal = currVal ;
            maxPos = i;
        } else if (currVal < minVal) {
            minVal = currVal ;
            minPos = i;
        }
        //
        if (lookForMax) {
            // 극대값을 찾는 상태이고, 현재값이 maxVal보다 충분이 아래면,
            // 이 maxVal이 극대값임 & 현재위치가 극소값의 후보가 됨;
            if (currVal < (maxVal - delta)) {
                maxPosition.push_back(maxPos);
                minVal = currVal ;
                minPos = i;
                lookForMax = false; //다음번에는 최소값을 찾는다;
            }
        } else {
            // 극소값을 찾는 상태이고, 현재값이 minVal보다 충분히 위면,
            // 이 minVal이 극소점임 & 현재위치가 극대값 후보가 됨;
            if (currVal > (minVal + delta)) {
                minPosition.push_back(minPos);
                maxVal = currVal;
                maxPos = i;
                lookForMax = true; //다음번에는 최대값을 찾는다;
            }
        }
 
    }
    return TRUE;
} 

데이터 구간이 (0,255) 히스트그램:

delta = 5;

""


delta  = 10;

""

delta = 15 ;

""

  

 
저작자 표시 비영리 변경 금지
신고

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

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
Adaboost  (0) 2010.12.28
Posted by helloktk

이미지에서 코드의 영역을 검출한 후에 는 코드영역의 비트값을 읽어들여서 디코딩하는 과정을 거쳐야 최종적으로 QR-code에 embedding 정보를 읽을 수 있다. 검출된 코드 영역은 역 perspective 변환을 거치면 정사각형의 코드영역으로 변환된다. 그러나 이 과정은 실제적인 설계에서는 불필요한 과정이다. 이미 주어진 perspective 변환을 이용하면, 이진화된 원래의 이미지에서 바로 코드를 재구성할 수 있다. 아래의 그림은 이 과정을 빠르게 수행하기 위해서 코드영역을 블럭단위로 구분하는 마스크 이미지를 만든 것을 보여준다.

마스크 이미지: 컬러는 각각의 블럭의 라벨링을 표현하기 위해서 사용했다.

이 마스크를 이용하면 각 블럭의 비트값을 majority test에 의해서 쉽게 결정할 수 있다.
재구성된 코드블럭:


이제 나머지는 QR-code의 스펙을 참고하면 (Reed-Solomon error-correction code는 영상처리 관점에서 약간 벗어나므로 취급하지 않음) 인코딩된 정보를 얻을 수 있다:

 

예제 코드의 정보는

저작자 표시 비영리 변경 금지
신고

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

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
Adaboost  (0) 2010.12.28
Blur Detection  (0) 2010.05.25
Posted by helloktk
TAG QR code

matlab에서 rgb 컬러이미지를 gray 이미지로 바꿀 때, 호출하는 함수가 rgb2gray이다. 이 함수는 주어진 RGB값에 다음과 같은 weight를 주어서 명암값을 얻는다:

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

이 계산을 fixed-point 버전으로 바꾸어 보자. 정밀도를 유지하기 위해서는 소수점 이하 4자리까지 유지해야 하는데, 이것은 weight에 각각 10000을 곱한 값을 사용한 후에 다시 10000으로 나누면 된다 (10000을 곱하더라도 r, g, b가 0-255사이의 값이므로 최대로 255 * 10000 보다 작은 값이 나와서 32-비트 정수 범위 내에 있게 된다)

gray = (2989 * r + 5870 * g + 1140 * b) / 10000;

여기서 좀 더 개선을 할 수 있다. 10000으로 나누는 과정을 shift 연산으로 바꾸면 된다. 정밀도를 보존하기 위해서 10000보다 큰 2의 power의 수를 찾으면 2^14 = 16384 이 주어진다. 따라서 10000 대신에 2^14을 곱한 weight를 쓰고, 계산 결과를 오른쪽으로 14만큼 shift 연산을 하면 된다.

0.2989 * 2^14 = 4897.1776;
0.5870 * 2^14 = 9617.408;
0.1140 * 2^14 = 1867.776;

gray = (4897 * r + 9617 * g + 1868 * b) >> 14;

weight의 총 합 = 4897 + 9617 + 1868 < 2^14 이어서 gray 값은 255를 넘지 않고, 중간 계산은 32-비트 정수 범위 내에서만 이루어 진다.

정말 빠를까?
저작자 표시 비영리 변경 금지
신고
Posted by helloktk
이제는 내 땀이 들어갔던 핫코드가 잊혀지고 QR코드가 대세인 모양이다...
우연히, ZXing library를 보고, c#(cpp-코드도 있음) ->cpp 코드로 바꾸어 본 것이다.

영상에서 코드를 검출하기 위한 전단계로 영상을 이진화 시켜야 한다. QR-code는 대체로 흰 배경에 검정색으로 인쇄가 되고 코드가 있는 영역은 검정과 흰색이 일정한(?) 비율로 섞여있다는 특징을 이용하면, 좀더 용이하게 이진화를 시킬 수 있게 된다. 전체 영상을 8x8 크기의 블럭으로 나누고, 각각의 블럭을 이진화시키는데, 이진화의 임계값은 주변 25개의 블럭의 평균값을 이용한다(코드가 없는 영역은 픽셀값이 균일할 수 있으므로 이런 경우에 예외적으로 처리하는 방법이 필요하다). 약각, empirical 한 수치가 들어간다.


// 한 8x8블럭의 검정(코드)과 하얀(배경)으로 구분짓는 값을 결정한다;
// 밝은 픽셀과 어두운 픽셀값간의 값의 차이가 24이상이면, 블럭평균값을 반환하고,
// 아니면, 이 블럭은 하얀색(배경)으로 인식하고, 이 경우에 임계값은 배경의 가장 어두운
// 값의 절반을 이용한다.
void calculateBlackPoints(BYTE *luminances, int stride,
                          int *blackPoints, int subWidth, int subHeight)
{
    int lineStart = 0;
    for (int y = 0; y < subHeight; y++)
    {
        for (int x = 0; x < subWidth; x++)
        {
            int sum = 0;
            int min = 255;
            int max = 0;
            for (int yy = 0; yy < 8; yy++)
            {
                int offset = ((y << 3) + yy) * stride + (x << 3);
                for (int xx = 0; xx < 8; xx++)
                {
                    int pixel = luminances[offset + xx] & 0xff;
                    sum += pixel;
                    if (pixel < min)
                    {
                        min = pixel;
                    }
                    if (pixel > max)
                    {
                        max = pixel;
                    }
                }
            }
           
            int average = (max - min > 24)?(sum >> 6):(min >> 1);
            blackPoints[lineStart + x] = average;
        }
        // move to next line;
        lineStart += subWidth;
    }
}  
// 주어진 임계값을 이용해서 한 개의 8x8 블럭내의 픽셀들을 이진화시킨다.
void threshold8x8Block(BYTE * luminances, int stride,
                       int xoffset, int yoffset,
                       int threshold,
                       BYTE *matrix)
{
    int offset = yoffset * stride + xoffset;
    for (int y = 0; y < 8; y++)
    {
        matrix[offset + 0] = luminances[offset + 0] < threshold ? 0x00 : 0xff ;
        matrix[offset + 1] = luminances[offset + 1] < threshold ? 0x00 : 0xff ;
        matrix[offset + 2] = luminances[offset + 2] < threshold ? 0x00 : 0xff ;
        matrix[offset + 3] = luminances[offset + 3] < threshold ? 0x00 : 0xff ;
        matrix[offset + 4] = luminances[offset + 4] < threshold ? 0x00 : 0xff ;
        matrix[offset + 5] = luminances[offset + 5] < threshold ? 0x00 : 0xff ;
        matrix[offset + 6] = luminances[offset + 6] < threshold ? 0x00 : 0xff ;
        matrix[offset + 7] = luminances[offset + 7] < threshold ? 0x00 : 0xff ;
        //moveto next line;
        offset += stride;
    }
}
// 각각의 8x8 블럭은 주변의 25개 블럭의 그레이값의 평균을 이용해서 이진화를 시도한다;
// 가장자리 특히 오른쪽과 아래쪽 0-7픽셀은 처리가 안될 수 도 있다.
void  calculateThresholdForBlock(BYTE *luminances, int stride,
                                 int *blackPoints, int subWidth, int subHeight,
                                 BYTE *matrix)
{
    for (int y = 0; y < subHeight; y++)
    {
        for (int x = 0; x < subWidth; x++)
        {
            int left = (x > 1)?x:2;
            left = (left < subWidth - 2)?left:subWidth - 3; // center-x;
            int top = (y > 1)?y:2;
            top = (top < subHeight - 2)?top:subHeight - 3;  // center-y;
            int sum = 0;
            // real top-line pointer;
            int *blackRow = &blackPoints[(top - 2) * subWidth];
            for (int yy = - 2; yy <= 2; yy++)
            {
                // int *blackRow = &blackPoints[(top + yy) * subWidth];
                sum += blackRow[left - 2];
                sum += blackRow[left - 1];
                sum += blackRow[left];
                sum += blackRow[left + 1];
                sum += blackRow[left + 2];
                // moveto nextline;
                blackRow += subWidth;
            }
            int average = sum / 25;
            threshold8x8Block(luminances, stride, x << 3, y << 3, average, matrix);
        }
    }
}
// test module;
void binarizeEntireImageHybrid(CRaster& raster, CRaster& out)
{
    CSize sz = raster.GetSize();
    int width = sz.cx ;
    int height = sz.cy ;
    int subWidth = width >> 3 ;
    int subHeight = height >> 3 ;
    int *blackPoints = new int [subWidth * subHeight];
    calculateBlackPoints((BYTE *)raster.GetDataPtr(), (int)raster.GetBytesPerLine(),
                            blackPoints, subWidth, subHeight);
    calculateThresholdForBlock((BYTE *)raster.GetDataPtr(), (int)raster.GetBytesPerLine(),
                                blackPoints, subWidth, subHeight, (BYTE *)out.GetDataPtr());
    delete [] blackPoints;
}


-이 단계후에는 QR코드의 finderpattern(영상의 검정/하얀 네모박스로 둘러싸인 패턴)를 찾으면 된다. 이 패던은 영상의 라인을 스캔하면 검점-하얀-검정-하얀-검점이 폭의 비가 1 : 1 : 3 : 1 : 1으로 나타난다는 QR코드의 특성을 이용하면 쉽게 찾을 수 있다.



아래 그림(반전을 함)은 이것을 구현하여서 finder pattern의 중심(빨간 십자)를 찾는 결과를 표시한 것이다. 이 부분의 구현은 코드가 좀 길므로(그러나 쉬움) ZXing라이브러를 직접 참고하면 된다. 이 작업은 ZXing 라이브러리에 몇가지 문제점과 비효율적인 부분이 있어서 코드를 새롭게 작성했음을 밝힌다.

 


-finder pattern의 세 중심을 찾은 후에는 추가로 bottom-right 영역의 alignment pattern (W-B-W=1:1:1)을 찾으면 QR코드가 perspective view에 의해서 코드가 왜곡이 된 정보를 구할 수 있다. 그러나, 이 alignment pattern은 finder pattern 보다도 찾기가 쉽지는 않다. 이 경우에는 finder pattern의 3점을 이용해서 왜곡정보를 얻어야 하는데, 이 경우에 가능한 기하학적인 변환은 affine변환으로 근사를 시켜야 한다 (직선의 평행성을 보존하는 2d-변환으로 일반적인 직사각형을 평형사변형으로 바꾸고, 6개의 매개변수에 의해서 결정이 된다).
alignment pattern의 중심이 얻어지면 4개의 중심점으로 perspective변환의 매개변수를 얻을 수 있다.

-alignment pattern의 중심(녹색십자)을 찾는 결과이다: finder pattern 3개 (bottomLeft, topLeft, topRight)와 alignment pattern이 화면상에서 시계방향으로 정렬이 된 사변형의 꼭지점을 이룰 때, 정상적인 코드 이미지이고, 반대로 반시계방향으로 주어지면 뒷면에서 찍은 코드 이미지이므로, 디코딩시 이 점을 고려해야 한다.


- finder pattern과 alignenment pattern(이것을 찾을 수 없는 경우에는 affine변환을 이용한다)의 정보를 이용해서 코드영역을 inverse perspective 변환을 통해서 rectify한 결과을 얻을 수 있다:


저작자 표시 비영리 변경 금지
신고

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

Peak Finder  (1) 2012.02.02
QR-code : decoder  (0) 2012.01.26
QR-code : detector  (0) 2012.01.12
Adaboost  (0) 2010.12.28
Blur Detection  (0) 2010.05.25
Savitzky-Golay smoothing filter  (2) 2010.03.24
Posted by helloktk

Adaboost

Image Recognition 2010.12.28 15:17


회색으로 칠해진 바닥에 검정색 점과 흰색 점들이 분포해 있다. 이제 이 점들을 몇 단계의 간단한 직선가르기로 분류를 하려고 한다. 먼저 바닥에 임의로 선을 그려서 바닥의 점이 선의 오른쪽에 (직선의 방향벡터를 생각하면 오른쪽의 의미가 명확해진다) 위치하면 흰색으로 분류하는 규칙을 마련하자.

그런데 이 직선가르기에 의한 분류는 "직선의 오른쪽에 검정색의 점이 있거나, 직선의 왼쪽에 흰색의 점이 있을 때 점들을 잘못 분류하는 에러를 일으킨다. 따라서 에러를 최소화시키기 위해서 여러 가지 직선을 그려보아서 분류에러가 최소인것을 분류기로 선택한다. 이렇게 만든 직선분류기는 대부분의 경우에 단순한 동전던지기에 의한 랜덤분류방법(50% 에러)보다는 좀 더 낳은 결과를 낸다.

다음 단계의 직선 가르기를 할 때는 이전 단계의 분류규칙에 의한 잘못을 참조하여서 만들어야 한다. 간단한 방법은 이전에 제대로 분류된 흰색 점은 약간 덜 희게 하고, 제대로 분류된 검정점도 덜 검게 하면, 잘못 분류된 흰색과 검정색의 점만 두드려지게 보이게 되어 다음번 직선 가르기를 할 때 이전 단계에서 잘못 분류한 점들이 더 눈에 잘 보이게 되므로 더 많이 참고를 하게 된다(기술적으로는 잘못 분류된 점들이 에러에 기여하는 가중치를 상대적으로 더 크게 만드는 것이다). 이러한 과정을 반복적으로 하면 횟수만큼의 직선분류기를 가지게 되고, 각각의 분류기에 의한 분류에러도 가지게 된다.

이들 분류기를 선형결합을 하면 개별적인 분류기보다도 훨씬 분류에러가 적은 강력한 분류기를 만들 수 있는데, 이때, 각각의 직선분류기에 주는 가중치는 그 분류기가 일으키는 에러를 참조해서 만든다. Adaboost 알고리즘에서는 직선분류기가 ε
의 에러를 주면 최종분류에 기여하는 정도는 α = log (1-ε)/ε 로 잡는다. 

저작자 표시 비영리 변경 금지
신고

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

QR-code : decoder  (0) 2012.01.26
QR-code : detector  (0) 2012.01.12
Adaboost  (0) 2010.12.28
Blur Detection  (0) 2010.05.25
Savitzky-Golay smoothing filter  (2) 2010.03.24
Watershed Algorithm Code  (0) 2010.03.19
Posted by helloktk

A Fast 2D Otsu Thresholding Algorithm Based on Improved Histogram
Ningbo Zhu, Gang Wang, Gaobo Yang, Weiming Dai
Pattern Recognition Letters 16 (1995) 653-666

Based on the gray-level histogram of an image, Otsu’s thresholding method takes the variance between clusters as the criterion to select the optimal threshold. This method has been cited as an effective thresholding technique and widely used in real thresholding tasks. However, there are always uncertain
disturbing factors in practical applications. These disturbing factors probably make one-dimensional (1D) histogram fail to possess obvious vally between two peaks. As a result, Otsu method can not obtain optimal segmentation results. To solve this problem, Liu et. al.[5] extended 1D histogram into
two-dimensional(2D) histogram, which utilizes not only gray value distribution of an image, but also the dependency of pixel in its neighborhood. This method greatly improved the segmentation results, specially in low SNR condition.

그레이 이미지의 2차원 히스토그램은 원래 이미지의 픽셀값 분포를 x-축으로, 주변 픽셀값의 평균을 이용해서 만든 이미지에서의 픽셀값의 분포를 y-축으로 하여서 만든 일종의 2차원 영상이다(256x256). 
원본이미지에서 전경 픽셀(검정)은 주변 픽셀을 평균한 값도 전경에 비슷한 값을 가져야 한다. 그러나 그 픽셀이 배경(흰색)에 나타난 검정색의 노이즈라면, 주변 픽셀의 평균값은 흰색에 가까운 값을 가지게 된다. 반대로, 배경픽셀이면, 주변 픽셀의 평균값도 배경에 비슷한 값을 가져야 하고, 만약에 그것이 전경에 놓인 흰색잡음이었다면, 주변픽셀의 평균값은 검정으로 나타날 것이다(전경과 배경의 경계 근처에 놓인 픽셀들도 마찬가지다). 따라서 이 2차원 히스토그램을 보면 잡음이 많은 경우에 어떤 픽셀값들을 노이즈로 처리할 수 있는가를 보다 명확하게 알 수 있다.
위 2차원 히스토그램의 x-축은 원본, y-축은 주변픽셀을 평균이미지의 픽셀값 분포(0~255=L)를 표현한다. A 영역에 해당하는 픽셀은 전경(검정)픽셀로 생각할 수 있고, C-영역에 해당하는 픽셀은 배경(흰색)으로 생각할 수 있다. B영역에 해당하는 픽셀은 배경에 놓인 검정색 노이즈이고, D-영역은 전경에 놓인 흰색 노이즈 또는 전경과 배경의 경계지점이라고 생각할 수 있다(전경이 배경보다도 작은 영역을 차지하는 경우에)

이미지크기 : 512x512;

픽셀평균을 구하는 윈도우 크기를 5x5로 할 때, 2차원 히스토그램(돗수를 log-스케일로 바꾸어서 표현한 것임. 붉은색일수록 돗수가 많음을 의미한다):
저작자 표시 비영리 변경 금지
신고
Posted by helloktk
"Optimizing Connected Component Labeling Algorithms"
Kesheng Wu, Ekow Otoo and Arie Shoshani
Lawrence Berkeley National Laboratory, Berkeley, CA, USA

ABSTRACT
This paper presents two new strategies that can be used to greatly improve the speed of connected component labeling algorithms. To assign a label to a new object, most labeling algorithms use a scanning step that examines some of its neighbors. The first strategy exploits the dependencies among the neighbors to reduce the number of neighbors examined. When considering 8-connected components in a 2D image, this can reduce the number of neighbors examined from four to one in many cases. The second strategy uses an array to store the equivalence information among the labels. This replaces
the pointer based rooted trees used to store the same equivalence information. It reduces the memory required and also produces consecutive final labels. Using an array based instead of the pointer based rooted trees speeds up the connected component labeling algorithms by a factor of 5 100 in our tests on random binary images.

Keywords: Connected component labeling, Union-Find, optimization


저작자 표시 비영리 변경 금지
신고
Posted by helloktk

Image Partial Blur Detection and Classication,  http://www.cse.cuhk.edu.hk/~leojia/all_final_papers/blur_detect_cvpr08.pdf

Blur Detection for Digital Images Using Wavelet Transform, http://bit.kuas.edu.tw/~jihmsp/2010/vol1/JIH-MSP-2010-01-003.pdf
저작자 표시 비영리 변경 금지
신고

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

QR-code : detector  (0) 2012.01.12
Adaboost  (0) 2010.12.28
Blur Detection  (0) 2010.05.25
Savitzky-Golay smoothing filter  (2) 2010.03.24
Watershed Algorithm Code  (0) 2010.03.19
Retinex 알고리즘  (10) 2010.02.03
Posted by helloktk

The Savitzky–Golay method essentially performs a local polynomial regression (of degree k) on a series of values (of at least k+1 points which are treated as being equally spaced in the series) to determine the smoothed value for each point.
""

Savitzky–Golay 필터는 일정한 간격으로 주어진 데이터들이 있을 때(이들 데이터는 원래의 정보와 노이즈를 같이 포함한다), 각각의 점에서 주변의 점들을 가장 잘 피팅하는 k-차의 다항식을 최소자승법을 이용해서 찾아서 그 지점에서의 데이터 값을 결정하는 필터이다. 이 필터는 주어진 데이터에서의 극대나 극소, 또는 봉우리의 폭을 상대적으로 잘 보존한다.(주변 점들에 동등한 가중치를 주는 Moving Average Filter와 비교해 볼 수 있다).

간단한 예로, 2차의 다항식과 5개의 데이터 점

{(-2, d0), (-1, d1), (0, d2), (1, d3), (2, d4)}

을 이용해서 중앙에서의 값을 결정하는 방법을 살펴보자. 사용하려고 하는 다항식은

p[x] = a0 + a1*x + a2*x2

이다. 다항식의 계수들은 다항식의 값과 실제 데이터의 값과의 차이를 최소화시키는 것들로 잡아야 한다.. 즉, 최소자승의 원리를 적용하여서 구하면 된다. 계산된 다항식의 값과 실제 데이터 값 사이의 차의 제곱을 구하면:

L = |a0-2*a1+4*a2-d0|2+|a0 -a1+ a2-d1|2+|a0-d2|2+|a0+a1+a2-d3|2
                          
+|a0+2*a1+4*a2-d4|2

이 식의 a0, a1, a2에 대한 극값을 가질 조건은

5*a0 + 10*a2 = d0 + d1 + d2 + d3 + d4

10*a1 = -2d0– d1 + d3+ 2d4

10*a0 + 34*a2 = 4d0 + d1+ d3 + 4d4        

이 식을 만족시키는 a
0를 구하면, 필터의 출력값 (원도우 중앙에서 값)이 결정된다.

필터출력 = a0 = (-3*d0 + 12*d1 + 17*d2 + 12*d- 3*d4) / 35;

위에서 계수 a0, a1, a2를 결정하는 방정식은 행렬로 정리하면 아래의 식과 같이 표현할 수 있다.

""

좌변의 5행3열 행렬을 A, a = [a0, a1, a2]T, d = [d0, d1, d2, d3, d4]T로 놓으면, 이 행렬방정식은 A.a = d 형태로 쓸 수 있고, 양변에 AT를 곱해서 왼쪽을 대칭행렬로 바뀌면

(AT.A).a = AT.d

따라서 해는,

a = (AT.A)-1.(AT.d)

이 식은 임의의 k-차 다항식을 이용한 경우에도 사용할 수 있다. 이 경우에 행렬 (ATA) 는 (k+1)x(k+1)의 대칭행렬이 된다. 행렬 A는 다항식의 찻수와 피팅에 사용이 될 데이터의 구간의 크기가 주어지면 정해지므로,  윗 식에서 (AT.A)-1.AT의 첫행 (a0을 d로 표현하는 식의 계수들)을 구해서 프로그램상에서는 lookup table를 만들어서 사용할 수 있다. 아래표는 mathematica 를 이용해서 윈도우 크기가 7 (7개 점)인 경우 2차 다항식을 사용할 때 계수를 구하는 과정이다.

 

""

2차 다항식일 때, 같은 방식으로 다양한 윈도우 크기에 따른 계수를 구할 수 있다.
*크기(n)에 따른 필터값 결정계수 (중앙에 대해 좌우대칭이다);

n=5;  W[] = {-3, 12, 17, 12, -3};
n=7;  W[] = {-2, 3,  6, 7, 6, 3. -2};
n=9;  W[] = {-21, 14. 39, 54, 59, 54, 39, 14, -21};
n=11; W[] = {-36, 9, 44, 69, 84, 89, 84, 69, 44, 9, -33};

필터출력 =   ∑ W[i]d[i] / ∑ W[i]

저작자 표시 비영리 변경 금지
신고

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

Adaboost  (0) 2010.12.28
Blur Detection  (0) 2010.05.25
Savitzky-Golay smoothing filter  (2) 2010.03.24
Watershed Algorithm Code  (0) 2010.03.19
Retinex 알고리즘  (10) 2010.02.03
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Posted by helloktk

적용 예는 http://kipl.tistory.com/1 
#include <vector>
#include <queue>
#define WMASK   -2
#define WINIT   -1
#define WSHED    0
#define FICTION ((UINT)(INT)-1)
/**
*/
void watershed_sorting_step(const std::vector<BYTE> & img,
                            int &hmin, int &hmax, UINT fr[256],
                            std::vector<UINT> & sort)
{
    UINT freq[ 256];
    /* Obtain frequency distribution */
    memset( freq, 0, 256*sizeof(UINT));
    for (int i = 0; i < img.size(); i++)
        freq[img[i]]++;
    memcpy( fr, freq, sizeof(freq));
   
    /* Calculate hmin & hmax */
    hmin = 0;  
    while ( hmin < 256 && freq[hmin] == 0) hmin++;
    hmax = 255;
    while ( hmax > 0 && freq[hmax] == 0)   hmax--;
   
    /* Obtain cumulative frequency distribution */
    for ( i = 1; i < 256; i++)
        freq[i] += freq[i-1];
   
    /* allocate and fill an array of sorted positions */
    for (int pos = 0; pos < img.size(); pos++) {
        int val = img[pos];
        // freq[100]은 cumulative histogram이므로 100값까지 같은 픽셀의 총수다.
        // freq[100]-freq[99]이 정확히 100값을 갖는 픽셀의 총수이다;
        // [ 전체 이미지                                          ]
        // 0----------------->1-------->2-------->........>255----]
        //이런식으로 파티션이 되어있으므로, 특정한 값이 한 번 나타나는 경우에
        int num = --freq[val] ;
        sort[num] = pos;
    };
    //sort-[000000000111112222222233333....................255255255]
    //각각의 위치는 그 값이 이미지에서 나타난 순서대로, 그 픽셀의 위치를 나타낸다.
}
#define NN(p) ((p)-width)   
#define SS(p) ((p)+width)
#define EE(p) ((p)+1)
#define WW(p) ((p)-1)
#define NW(p) ((p)-width-1)
#define NE(p) ((p)-width+1)
#define SW(p) ((p)+width-1)
#define SE(p) ((p)+width+1)
//
void Watershed8UC(CRaster & raster, CRaster & res)
{
    int neighborhood = 8;
    int N, width, height, maxx, maxy;
    int hmin, hmax;
 
    UINT maxdist = 1;
    UINT freq[256];
    UINT pbis[8];

    /* Initialization */
    CSize sz = raster.GetSize();
    width = sz.cx;
    height = sz.cy;
    maxx = width - 1;
    maxy = height - 1;
    N = width * height;

    /* Convert our image to non-aligned data structure */
    std::vector<BYTE> inp(N) ;
    for (int y = 0; y < height; y++) {
        BYTE *p = (BYTE*)raster.GetLinePtr(y) ;
        memcpy(&inp[y * width], p, width);
    }

    /* Obtain sorted array of pixel positions */
    std::vector<UINT> sorted(N);
    watershed_sorting_step(inp, hmin, hmax, freq, sorted);

    /* Check the correctness of sorting */
    if (1) {
        BYTE last = 0;
        for (int i = 0; i < N; i++)
            if (inp[sorted[i]] < last)
                TRACE("%s: incorrectly sorted", "check");
            else
                last = inp[sorted[i]];
    }

    /* create and initialize output data array */
    std::vector<int> out(N, WINIT);

    /* create and initialize distance work image */
    std::vector<UINT> distance(N, 0);
   
    /* create and initialize the queue */
    std::queue<UINT> Q;

    /* other initialization */
    int label = 0;
    int i = 0, p;

    /* main loop */
    for (int h = hmin; h <= hmax; h++) {
        /*=== geodesic SKIZ of level h-1 inside level h */
        /*=== for every pixel p such that image(p) == h */
        int n = freq[h];
        if (neighborhood == 4) {
            while (n--) {
                UINT p = sorted[i++];
                if (inp[p] != h)
                    TRACE("sort assertion failed: %d(%d) != %d", inp[p], p, h);
                int x = p % width;
                int y = p / width;
                out[p] = WMASK;
                if ((x < maxx && out[EE(p)] >= WSHED) ||
                    (x > 0    && out[WW(p)] >= WSHED) ||
                    (y > 0    && out[NN(p)] >= WSHED) ||
                    (y < maxy && out[SS(p)] >= WSHED))
                {
                    distance[p] = 1;
                    Q.push(p);
                }
            }
        } else {                  /* neighborhood == 8 */
            while (n--) {
                UINT p = sorted[i++];
                int x = p % width;
                int y = p / width;
                out[p] = WMASK;
                if ((x < maxx && (out[EE(p)] >= WSHED || (y > 0 && out[NE(p)] >= WSHED) || (y < maxy && out[SE(p)] >= WSHED))) ||
                    (x > 0    && (out[WW(p)] >= WSHED || (y > 0 && out[NW(p)] >= WSHED) || (y < maxy && out[SW(p)] >= WSHED))) ||
                    (y > 0    &&  out[NN(p)] >= WSHED) || (y < maxy && out[SS(p)] >= WSHED))
                {
                    distance[p] = 1;
                    Q.push(p);//FIFO_ADD(p);
                }
            }
        }
        int dist = 1;
        Q.push(FICTION);
        while (1) {
            UINT p=Q.front(); Q.pop();
            if (p == FICTION) {
                if (Q.empty())
                    break;
                Q.push(FICTION);
                dist++;
                if (dist > maxdist)
                    maxdist = dist;
                p=Q.front();Q.pop();
            }
            if (Q.empty())
                TRACE("assertion failed: !FIFO_EMPTY");
            int x = p % width;
            int y = p / width;
            int npbis = 0 ;
            if (neighborhood == 4) {
                if (x > 0)
                    pbis[npbis++] = WW(p);
                if (x < maxx)
                    pbis[npbis++] = EE(p);
                if (y > 0)
                    pbis[npbis++] = NN(p);
                if (y < maxy)
                    pbis[npbis++] = SS(p);
            } else {            /* neighborhood == 8 */
                if (x > 0) {
                    pbis[npbis++] = WW(p);
                    if (y > 0)
                        pbis[npbis++] = NW(p);
                    if (y < maxy)
                        pbis[npbis++] = SW(p);
                }
                if (x < maxx) {
                    pbis[npbis++] = EE(p);
                    if (y > 0)
                        pbis[npbis++] = NE(p);
                    if (y < maxy)
                        pbis[npbis++] = SE(p);
                }
                if (y > 0)
                    pbis[npbis++] = NN(p);
                if (y < maxy)
                    pbis[npbis++] = SS(p);
            }
            /*=== for every pixel p' belonging to Ng(p) */
            while (--npbis >= 0) {
                UINT pbs = pbis[npbis];
                if (distance[pbs] < dist && out[pbs] >= WSHED) {
                    if (out[pbs] > 0) {
                        if (out[p] == WMASK || out[p] == WSHED)
                            out[p] = out[pbs];
                        else if (out[p] != out[pbs])
                            out[p] = WSHED;
                    } else if (out[p] == WMASK)
                        out[p] = WSHED;
                } else if (out[pbs] == WMASK && distance[pbs] == 0) {
                    distance[pbs] = dist + 1;
                    Q.push(pbs);
                }
            }
        }                       /* while(1) */
        /*=== checks if new minima have been discovered */
        n = freq[h];
        i -= n;
        while (n--) {
            UINT p = sorted[i++];
            distance[p] = 0;
            if (out[p] == WMASK) {
                label++;
                Q.push(p);
                out[p] = label;
                while (!Q.empty()){
                    UINT pp=Q.front();Q.pop();
                    int x = pp % width;
                    int y = pp / width;
                    int npbis = 0;
                    if (neighborhood == 4) {
                        if (x > 0)
                            pbis[npbis++] = WW(pp);
                        if (x < maxx)
                            pbis[npbis++] = EE(pp);
                        if (y > 0)
                            pbis[npbis++] = NN(pp);
                        if (y < maxy)
                            pbis[npbis++] = SS(pp);
                    } else {    /* neighborhood == 8 */
                        if (x > 0) {
                            pbis[npbis++] = WW(pp);
                            if (y > 0)
                                pbis[npbis++] = NW(pp);
                            if (y < maxy)
                                pbis[npbis++] = SW(pp);
                        }
                        if (x < maxx) {
                            pbis[npbis++] = EE(pp);
                            if (y > 0)
                                pbis[npbis++] = NE(pp);
                            if (y < maxy)
                                pbis[npbis++] = SE(pp);
                        }
                        if (y > 0)
                            pbis[npbis++] = NN(pp);
                        if (y < maxy)
                            pbis[npbis++] = SS(pp);
                    }
                    while (--npbis >= 0) {
                        UINT pbs = pbis[npbis];
                        if (out[pbs] == WMASK) {
                            Q.push(pbs);
                            out[pbs] = label;
                        }
                    }
                }
            }
        }
    }

    if (!Q.empty())//FIFO_EMPTY)
        TRACE("$s: queue is not empty - can't be", "");

    /* Convert the result to a suitable form */
    res.SetDimensions(width, height, 8, 1) ;
    for (y = 0, p = 0; y < height; y++) {
        for (int x = 0; x < width; x++, p++) {
            // int p = x + width * y;
            if (out[p] == WMASK) TRACE("%s: %d,%d has mask", "", x, y);
            if (out[p] == WINIT) TRACE("%s: %d,%d has init", "", x, y);
            if (out[p] == WSHED) continue;
            if (neighborhood == 4) {
                if ((x < maxx && out[EE(p)] > WSHED && out[EE(p)] < out[p]) ||
                    (x > 0    && out[WW(p)] > WSHED && out[WW(p)] < out[p]) ||
                    (y > 0    && out[NN(p)] > WSHED && out[NN(p)] < out[p]) ||
                    (y < maxy && out[SS(p)] > WSHED && out[SS(p)] < out[p]))
                    out[p] = WSHED;
            } else {
                if ((x < maxx && ((out[EE(p)] > WSHED && out[EE(p)] < out[p]) ||
                    (y > 0    && out[NE(p)] > WSHED   && out[NE(p)] < out[p])||
                    (y < maxy && out[SE(p)] > WSHED   && out[SE(p)] < out[p]))) ||
                    (x > 0    && ((out[WW(p)] > WSHED && out[WW(p)] < out[p])||
                    (y > 0    && out[NW(p)] > WSHED   && out[NW(p)] < out[p])||
                    (y < maxy && out[SW(p)] > WSHED   && out[SW(p)] < out[p]))) ||
                    (y > 0    && out[NN(p)] > WSHED   && out[NN(p)] < out[p]) ||
                    (y < maxy && out[SS(p)] > WSHED   && out[SS(p)] < out[p]))
                    out[p] = WSHED;
            }
        }
    }
    for (y = 0; y < height; y++)
    {
        BYTE *p = (BYTE*) res.GetLinePtr(y) ;
        for (int x = 0; x < width; x++) {
            UINT a=out[x + width * y] ;
            p[x] = a;//((a == WSHED) ? 255 : 0);
        }
    }
    return ;
}

저작자 표시 비영리 변경 금지
신고

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

Blur Detection  (0) 2010.05.25
Savitzky-Golay smoothing filter  (2) 2010.03.24
Watershed Algorithm Code  (0) 2010.03.19
Retinex 알고리즘  (10) 2010.02.03
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Image Morphing  (0) 2010.01.24
Posted by helloktk
TAG Watershed


Retinex 알고리즘은 영상의 contrast를 향상시키거나, sharpness를 증진시킬 때 많이 사용한다. 또한 픽셀값의 dynamic range가 큰 경우에 이것을 압축시켜 영상 데이터 전송에 따른 병목현상의 해소에 이용할 수 있다.

Retinex 알고리즘의 기본원리는 입력영상에 들어있는 배경성분을 제거하는 것이다.
1. 배경영상는  입력영상의 평균적인 영상로 생각할 수 있는데, 이것은 적당한 스케일(필터사이즈)의 가우시안 필터를 적용하여서 얻을 수 있다. 이 필터를 적용하면 입력영상에서 필터사이즈 보다도 작은 스케일은 무시하는 효과를 준다.
2. 입력영상의 반사성분(배경조명에 무관한)은 입력영상을 앞서 구한 배경영상으로 나누면 된다
3. Retinex 출력은 이 반사성분에 로그값을 취한 것이다. 이 로그값을 취하므로써 반사성분이 분포범위(=dynamic range)를 압축하는 효과를 얻는다( 원래 로그함수는 큰 수의 연산을 작은 수의 연산으로 바꾼다)
- 입력영상 = I(x, y) ;
- 가우시안 필터 = G(x, y) = A exp(-(x2 + y2 )/2σ2)
- 배경영상 = G(x, y)* I(x, y) = convolution;
- Retinex 출력 
   R(x, y) = log(반사성분)
             = log(입력영상/배경영상)
             = log(I(x, y)) - log(G(x, y)* I(x, y));


이처럼 하나의 스케일에 대해서 적용하는 경우를 SSR(single-scale retinex)알고리즘이라고 한다. 컬러영상의 경우에는 각각의 RGB채널에 대해서 알고리즘을 적용하면 된다.

Retinex 영상을 구할 때 하나의 스케일이 아니라 여러가지 스케일에 대해서 적용한 retinex 영상을 적절한 가중치를 주어서 합한 결과를 출력영상으로 사용할 수 있다. 이 경우가 MSR(multi-scale retinex)알고리즘이다.

  Rmsr(x,y) = i wiRssr(x,y) = i (가중치) * SSR.

이렇게 얻은  Retinex 영상은 적당한 offset과 stretching을 하여서 픽셀값이 [0,255]구간에 있게 조절한다. 이 과정은 Retinex 출력영상 픽셀의 평균값과 편차를 구하여서 분석할 수 있다.

컬러영상의 경우에 Retinex 출력영상은 전체적으로 그레이화되는 경향이 있어서 이것을 보완하기 위해서 아래의 추가적인 처리과정을 더 거친다.

   R'msr(x,y)_red = log(C* Ir/(Ir + Ig + Ib)) * Rmsr(x,y)_red ;
   R'msr(x,y)_gre = log(C* Ig/(Ir + Ig + Ib)) * Rmsr(x,y)_gre ;
   R'msr(x,y)_ble = log(C* Ib/(Ir + Ig + Ib)) * Rmsr(x,y)_blu ;


여기서, C는 상수값이고, Ir, Ig, Ib는 입력영상의 RGB-채널이다.

*그레이 이미지 처리 결과:


*구현 코드는 다음을 참고: 그레이: http://kipl.tistory.com/33 
                               컬러: http://blog.naver.com/helloktk/80039132534
* 기술적으로  log를 취하므로 입력영상에 +1을 더해서 log(0)이 나오는 것을 방지해야 한다.
* convolution된 영상도 1보다 작은 경우에 1로 만들어야한다.
* 스케일이 큰 경우에 필터사이즈가 매우 크므로 효과적인 convolution 알고리즘을 생각해야 한다. 보통 recursive 필터링을 하여서 빠르게 convolution결과를 얻는다.

저작자 표시 비영리 변경 금지
신고

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

Savitzky-Golay smoothing filter  (2) 2010.03.24
Watershed Algorithm Code  (0) 2010.03.19
Retinex 알고리즘  (10) 2010.02.03
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Image Morphing  (0) 2010.01.24
Fant's Algorithm  (0) 2010.01.22
Posted by helloktk