Geometry & Recognition

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

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

$$B(i,j)=\{\begin{array}{ll}{m}_0,& \text{if}A(i,j)\le T\\ m_1,& \text{otherwise}\end{array}$$

로 나타난다. 여기서 $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)=p{m}_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{array}{rl}\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{array}$$

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

참고: Otsu Algorithm 구현 예.

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

윈도가 움직이면서 윈도 내 모든 픽셀이 다 바뀌는 것이 아니라 움직이는 방향에 수직인 가장자리 픽셀만 바뀐다는 사실을 이용하면 각 픽셀의 윈도에서 median을 찾는 작업을 할 필요가 없다. 예를 들면 scanline 방향($x$-방향)으로 윈도를 움직이면서 median filter를 작용할 때, 윈도가 오른쪽으로 1 픽셀 움직이면 윈도의 왼쪽 가장자리의 수직 픽셀들은 새 윈도에서 사라지고, 기존 윈도의 오른쪽 수직 가장자리 앞의 픽셀들이 새로 들어온다. 따라서 각 스캔라인에서 처음 한 번만 윈도의 median을 찾으면 이후에는 윈도가 이동할 때 윈도를 나가는 픽셀과 새로 들어오는 픽셀 ($2W$) 개에 대해서 이전 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;

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