Geometry & Recognition

RANSAC 알고리즘을 써서 주어진 2차원 점집합에서 원을 추정한다. 원을 만들기 위해서는 최소한 3점이 필요하고, 또 일직선에 있지 않아야 한다. 이렇게 만들어진 원은 세 점을 꼭짓점으로 하는 삼각형의 외접원(circumcircle)이다. 주어진 외접원에서 크게 벗어나지 않는 inliers를 찾는다(추가로 이 inliers에 대해 최소자승법으로 원의 중심과 반지름을 다시 구해서 보다 정밀하게 추정하는 과정을 넣을 수도 있다). 무작위로 선택된 세 점에 대해 위의 과정을 반복 시행해서 구한 원 중에서 가장 많은 inliers를 포함하는 원을 결과로 사용한다.

// 참고: http://en.wikipedia.org/wiki/RANSAC

// 2024.5.31 재작성;
double dist_deviate(const CfPt& pts, double cparam[3]) {
    double dx = pts.x - cparam[0];
    double dy = pts.y - cparam[1];
    return fabs(hypot(dx, dy) - cparam[2]);
}
int circumcircle(CfPt pts[3], double cparam[3]) {
    double x1 = pts[0].x, x2 = pts[1].x, x3 = pts[2].x;
    double y1 = pts[0].y, y2 = pts[1].y, y3 = pts[2].y;
    double bax = x2 - x1, bay = y2 - y1;
    double cax = x3 - x1, cay = y3 - y1;
    double E = bax * (x1 + x2) + bay * (y1 + y2);
    double F = cax * (x1 + x3) + cay * (y1 + y3);
    double G = 2. * (bax * (y3 - y2) - bay * (x3 - x2));
    if (G == 0.) return 0;    //error;
    //assert(fabs(G)>small_epsilon); //to prevent collinear or degenerate case;
    cparam[0] = (cay * E - bay * F) / G; //cx;
    cparam[1] = (bax * F - cax * E) / G; //cy;
    cparam[2] = hypot(cparam[0]-x1, cparam[1]-y1); //rad;
    return 1;
};
int num_sampling3(double prob_fail, double inlier_ratio) {
    return int(log(prob_fail)/log(1-pow(inlier_ratio, 3))); 
}
std::vector<int> Ransac_CircleFit(std::vector<CfPt>& points, double circle_param[3]) {
    if (points.size() < 3)
        return std::vector<int> (); //return null_vector;

    CfPt center; double inv_scale;
    // normalize input points for the sake of numerical stability;
    std::vector<CfPt> nor_pts = normalize(points, inv_scale, center);
    // distance threshold;
    double distance_thresh = sqrt(double(points.size())) * inv_scale;
    //ransac
    int sample_num = 1000;	//number of sample
    int ransac_count = 0;
    const double prob_fail = 0.01;
    double best_cparam[3] = {0};
    std::vector<int> best_inliers;
    while (sample_num > ransac_count) {
        // pick random 3 indices:[0,points.size()-1];
        int triple[3];
        random_triple(points.size()-1, triple);
        CfPt selected[3];
        for (int i = 0; i < 3; i++) 
            selected[i] = nor_pts[triple[i]];
        // circumcircle of 3 points;
        if (circumcircle(selected, circle_param)) {
            // find inliers;
            std::vector<int> inliers;
            inliers.reserve(points.size());
            for (int i = nor_pts.size(); i-->0;) {
                // error measure = algebric distance;
                double deviation = dist_deviate(nor_pts[i], circle_param);
                if (fabs(deviation) < distance_thresh)
                    inliers.push_back(i);
            }
            if (inliers.size() > best_inliers.size()) {			
                // update sampling_num;
                sample_num = num_sampling3(prob_fail, double(inliers.size())/points.size());
                // update best_inliers;
                best_inliers.swap(inliers);
                // update best circle param;
                for (int i = 0; i < 3; i++) 
                    best_cparam[i] = circle_param[i];
            }
        }
        if (++ransac_count > 1500) {
            TRACE("error! ransac_count exceed!\n");
            break;
        }
    }
    // recover original coordinate and scale;
    denormalize(best_cparam, best_cparam, inv_scale, center);
    if (best_cparam[0] > 0 && best_cparam[1] > 0) {
        for (int i = 0; i < 3; i++)
            circle_param[i] = best_cparam[i];
        TRACE("circle_found(%d, %d)\n", sample_num, ransac_count);
        // more accurate estimation needed at this stage;
    } else 
        best_inliers.clear();
    return best_inliers;
}

https://kipl.tistory.com/207

https://kipl.tistory.com/357

KMeans 알고리즘은 주어진 데이터들을 사전에 정해진 수만큼의 클러스터로 자동 분류하는 간단한 클러스터링 알고리즘의 일종이다. 이 알고리즘의 주된 아이디어는 각 클러스터의 중심을 어떻게 하면 잘 잡는가이다. 좋은 선택은 각각의 중심이 가능한 서로 멀리 떨어지게 잡아야 한다. 일단 클러스터 중심이 설정되면 각각의 데이터에 대해서 가장 가까이에 있는 클러스터에 할당시키면 된다. 다음으로 클러스터의 중심을 그 클러스터에 속한 데이터의 무게중심으로 재설정하고, 다시 각각의 데이터를 클러스터에 할당하는 작업을 반복한다. 이 과정은 클러스터의 중심이 더 이상 이동하지 않을 때까지 반복 시행하면 된다. 결과적으로 이 알고리즘은 클러스터의 중심에서 그 클러스터에 할당된 데이터와의 거리의 제곱을 그 클러스터의 모든 점과 전체 클러스터에 대해서 합한 값을 최소화시키는 클러스터 중심과 데이터의 분할을 찾는 것이 된다.

1. $k$개의 클러스터의 중심을 설정한다.
2. 설정된 중심을 가지고, 각 데이터를 분할한다.
3. 분할된 데이터를 이용하여서 클러스터의 중심을 재설정한다.
4. 중심이 변화가 없을 때까지 2,3 과정을 반복한다.

비용 함수의 관점에서 보면 KMeans 클러스터링은 다음 함수를

$$ \text{cost}=\sum_{k \in \rm classes} \sum_{i \in \rm class_k} \big| \text{data}_k(i) - \text{mean}_k \big|^2 $$

을 최소화시키는 분할을 하는 것이다. 그러나,  KMeans는 항상 최적의 비용 함숫값을 보장하지는 않는다. 그리고 알고리즘의 결과는 초기값의 설정에 의존하여서 다른 결과를 낼 수가 있다. 보다 큰 단점은 클러스터의 수를 사전에 정해주어야 한다는 사실인데, 잘못된 클러스터 개수의 설정은 성능에 크게 영향을 미친다. 자동의 클러스터의 수를 설정하면서 클러스터링을 할 수 있는 간단한 것으로는 isodata 알고리즘이 있다.

typedef double * vector_t ;
double kmeans_label(
             const vector_t *obs,             // observations;
             const int n_obs,                   /* in # of vectors */            
             const vector_t *mean,           // mean vectors;
             const int n_mean,                /* # of mean vectors */
             int* label,                           // cluster labeling;
             const int veclen)                 // dimension of feature;
{
    double sqerr=0;
    for (int i = 0; i < n_obs; i++) {
        vector_t c = obs[i];
        /* Get an estimate of best distance (min_dist) and (min_idx) */
        int min_idx = label[i];
        vector_t m = mean[min_idx];
        double min_sqdist = 0 ;
        for (int l = 0; l < veclen; l++) {
            double t = m[l] - c[l]; min_sqdist += t * t;
        }
        for (int j = 0; j < n_mean; j++) {
            vector_t m = mean[j];
            // save time by checking:: dist < min_sqdist=previous min_val;
            double sqdist = 0;
            for (l = 0; (l < veclen) && (sqdist < min_sqdist); l++) {
                double t = m[l] - c[l]; sqdist += t * t;
            }
            if (sqdist < min_sqdist) {
                min_sqdist = sqdist; min_idx = j;
            }
        }
        label[i] = min_idx;
        sqerr += min_sqdist;
    }
    return sqerr ;
};

int kmeans_update(
              const vector_t *obs, //observations;
              const int n_obs,       //# of observations;     
              vector_t *mean,       // mean vectors;
              const int n_mean,    // # of mean vectors;
              int *label,               // cluster labelling;
              const int veclen      // dimension of features;
              ) {
    std::vector<int> cnt(n_mean);
    for (int i = 0; i < n_mean; i++)
        for (int l = 0; l < veclen; l++) mean[i][l] = 0.0;
        
    for (i = 0; i < n_obs; i++) {
        ASSERT((0 <= label[i]) && (label[i] < n_mean));  
        vector_t m = mean[label[i]];
        cnt[label[i]]++;  
        vector_t c = obs[i];
        if (c == NULL)
            TRACE("No observations for %u, but expected up through %u", i, n_obs-1);
        for (int l = 0; l < veclen; l++) m[l] += c[l];
    }
    for (i = 0; i < n_mean; i++) {
        int j = cnt[i];
        if (j != 0) for (int l = 0; l < veclen; l++) mean[i][l] /= (double) j;
        else {
            TRACE("Empty cluster %u\n", i);
            return FALSE;
        }
    }
    return TRUE;
}

double kmeans(
       const vector_t *obs,     // observations;
       const int n_obs,          /* # of observations */
       vector_t *mean,           /* initial set of means */
       const int n_mean,        /* # of means (should be k_mean?) */
       int *label,                   /* cluster labeling*/
       const int veclen,         /* vector length of means and corpus */    
       const double min_conv_ratio,
       const int max_iter
    ) {
    double p_sqerr = DBL_MAX;
    int ret;
    double sqerr = kmeans_label( obs, n_obs,mean, n_mean, &label[0], veclen);
    double conv_ratio = (p_sqerr - sqerr) / p_sqerr;
    for (int i = 0; (i < max_iter) && (conv_ratio > min_conv_ratio); i++) {
        TRACE("kmeans iter [%u] %e ...\n", i, conv_ratio);
        ret = kmeans_update(obs, n_obs, mean, n_mean, &label[0], veclen );
        if (ret != TRUE)
            return (double)ret;
        p_sqerr = sqerr;
        sqerr = kmeans_label(  obs, n_obs,mean, n_mean, &label[0], veclen);
        conv_ratio = (p_sqerr - sqerr) / p_sqerr;
    }
    TRACE("kmeans n_iter %u sqerr %e conv_ratio %e\n", i, sqerr, conv_ratio);   
    return sqerr;
}

이미지를 이진화시키기 위해서 여러 알고리즘이 사용된다. 그중 이미지 전체에 대해 하나의 임계값으로 이진화시키는 전역 이진화 알고리즘은 간단하고 빠르기 때문에 많이 이용이 된다. 그러나 이미지를 형성할 때 조명 조건이 균일하지 않은 경우에는 전역 이진화는 원하는 결과를 얻기가 힘들다. 이런 경우에는 각각의 픽셀 주위의 그레이 값을 참조하여 임계치를 결정하는 국소적 이진화 방법을 사용한다. 국소적 이진화에서 임계값을 추출하는 간단한 방법은 윈도 내의 평균값을 이용하면 된다. 좀 더 개선된 알고리즘은 평균값($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;
        }
    }   
};