Geometry & Recognition

영상의 히스토그램($h[z]$)이 bimodal로 주어지는 경우 적절한 threshold 값을 선택해서 전경과 배경을 분리할 수 있다. 전경을 대표하는 픽셀 값을 $z_f$, 배경을 대표하는 픽셀 값을 $z_b$라면 이진화 후 정규화된 히스토그램은 

$$ \tilde{h}[z] =p_b \delta_{z, z_{b}} + p_f \delta_{z, z_{f}}$$

로 표현된다. $p_b$은 배경에 해당하는 픽셀 비율이고, $p_f$은 전경에 해당하는 픽셀 비율이다.

threshold 값을 어떻게 선택하면 이진화된 영상의 히스토그램이 원 영상의 히스토그램의 특성을 최대한 담게 할 수 있을까? 이에 대한 기준으로 두 히스토그램의 $n$차 moment가 같은 값을 갖도록 하자. 주어진 미지수가 $p_b$, $p_f$, $z_b$, $z_f$이 있으므로 최소한 4개의 moment가 같도록 만들어야 한다. 가장 낮은 찾수의 moment로부터 시작해서 3차까지 4개의 moments가 같다는 조건에서 아래의 식들을 얻을 수 있다.

$$ \text{0-차 moment: }  m_0 \equiv \sum_{z=0}^{255} h[z] = p_b  + p_f = 1$$

$$ \text{1-차 moment: }  m_1 \equiv \sum_{z=0}^{255}z h[z] = p_b z_b + p_f z_f$$

$$ \text{2-차 moment: }  m_2 \equiv \sum_{z=0}^{255}z^2 h[z] = p_b z_b^2 + p_f z_f^2$$

$$ \text{3-차 moment: }  m_3 \equiv \sum_{z=0}^{255}z^3 h[z] = p_b z_b^3+ p_f z_f^3$$

원 영상의 moment $m_0$, $m_1$, $m_2$, $m_3$을 계산해서 풀면

\begin{gather} c_0 = \frac{m_3 m_1 - m_2^2 }{ m_0 m_2 - m_1^2 } ,\quad  c_1 = \frac{m_1 m_2 - m_0 m_3}{ m_0 m_2 - m_1^2 } \\ z_b = \frac{1}{2} \left(-c_1 - \sqrt{ c_1^2 - 4 c_0} \right) \\  z_f = \frac{1}{2}\left( -c_1 + \sqrt{ c_1^2 -4 c_0} \right) \\ p_b = \frac{z_f - m_1}{z_f - z_b} \\ p_f = 1 - p_b\end{gather}

따라서 threshold 값

$$ \sum_{z=0}^{T-1} h[z] = p_b  $$

을 만족하는 $T$을 선택하면 된다.

Ref: W. Tsai, "Moment-preserving thresholding: a new approach," Computer Vision, Graphics, and Image Processing, vol. 29, pp. 377-393, 1985.

int MomentsPreseving_threshold(int histogram[256]) {
    int tot = 0;
    for (int i = 0; i < 256; i++)
        tot += histogram[i];
    //normalised histogram
    double hist[256];
    for (int i = 0; i < 256; i++)
        hist[i] = double(histogram[i]) / tot;
    /* moments calculation: zero moment is 1 by defintion*/
    double m0 = 1, m1 = 0, m2 = 0, m3 = 0;
    for (int i = 0; i < 256; i++ ) {
        double h = hist[i];
        m1 += i * h;
        m2 += i * i * h;
        m3 += i * i * i * h;
    }
    double det = m0 * m2 - m1 * m1;
    double c0 = (m1 * m3 - m2 * m2) / det;
    double c1 = (m2 * m1 - m3 * m0) / det;
    double zb = 0.5 * (-c1 - sqrt (c1 * c1 - 4.0 * c0));
    double zf = 0.5 * (-c1 + sqrt (c1 * c1 - 4.0 * c0));
    double pb = (zf - m1) / (zf - zb);  
    double s = 0;
    for (int i = 0; i < 256; i++) {
        s += hist[i];
        if (s > pb)
            return i; // threshold
    }
    return 0;
}

$$threshold = \underset{0\le t < 256}{\text{argmin}} \left[-p_f(t) \log(\mu_f(t) )- p_b (t) \log( \mu_b (t)) \right] $$

where

$$ p_b (t) =\int_0^t h(z)dz, \quad  p_f(t) = \int_t^{255} h(z)dz$$

$$ \mu_b(t) = \frac{1}{p_b(t)} \int_0^{t} z h(z) dz,\quad \mu_f(t) = \frac{1}{p_f(t)} \int_t^{255} z h(z)dz $$

Ref: Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum Cross Entropy Thresholding"Pattern Recognition Letters, 18(8): 771-776

double MCE_threshold(int hist[256]) {
    int chist[256], cxhist[256];
    chist[0] = hist[0]; cxhist[0] = 0;
    for (int i = 1; i < 256; i++) { 
        chist[i] = hist[i] + chist[i - 1];
        cxhist[i] = i * hist[i] + cxhist[i - 1];
    }
    int num = chist[255];
    double mean = double(cxhist[255]) / num;
    /* Initial estimate */
    double threshold = mean;
    while (1) {
        double old_thresh = threshold;
        int t = int(old_thresh + .5);
        /* background */
        int bgnum = chist[t];
        int bgsum = cxhist[t];
        double bgmean = bgnum == 0 ? 0: double(bgsum) / bgnum;
        /* foreground */
        int fgnum = num - bgnum;
        int fgsum = cxhist[255] - bgsum;
        double fgmean = fgnum == 0 ? 0: double(fgsum) / fgnum;
        threshold = (bgmean - fgmean) / (log(bgmean) - log(fgmean));
        // new thresh is a simple round of theta;
        ASSERT(threshold >= 0);
        if (fabs(threshold - old_thresh) < 0.5)
           break;
    }
    return threshold;
}

이미지를 균일한 컬러 영역으로 분할하는 것은 압축이나 segmentation 등에 매우 유용하다. 4개의 child node 가운데 parent node와 차이가 많이 나지 않는 것은  부분적으로 parent node에 merge 시키는 방법도 있다.

#include "QuadTree.h"
#include <queue>
#include <math.h>
QuadNode::QuadNode() {
    for (int i = 0; i < 4; i++) children[i] = NULL;
    rect = CRect(0, 0, 0, 0);
    hasChildren = false;
    mean = -1;
    stdev = 0;
}
QuadNode::QuadNode(CRect& rc) {
    for (int i = 0; i < 4; i++) children[i] = NULL;
    rect = rc;
    hasChildren = false;
    mean = -1;
    stdev = 0;
}
void QuadNode::Subdivide() {
    hasChildren = true;
    children[0] = new QuadNode(CRect(rect.TopLeft(), rect.CenterPoint()));
    children[1] = new QuadNode(CRect(rect.left, rect.CenterPoint().y, rect.CenterPoint().x, rect.bottom));
    children[2] = new QuadNode(CRect(rect.CenterPoint(), rect.BottomRight()));
    children[3] = new QuadNode(CRect(rect.CenterPoint().x, rect.top, rect.right, rect.CenterPoint().y));
}
// quadtree에 할당된 메모리를 정리함; root node는 별도로 정리.
void QuadNode::Clear() {
    if (this->hasChildren)
        for (int i = 0; i < 4; i++)
            if (children[i]) {
                children[i]->Clear();
                delete children[i];
            }
}
QuadTree::QuadTree(CRect rc, CRaster& image) {
    rect = rc & image.GetRect(); // rect는 이미지 안에.
    raster = image;
}
void QuadTree::Build() {
    root = new QuadNode();
    root->rect = this->rect;
    std::queue<QuadNode *> TreeQueue;
    TreeQueue.push(root);
    while (!TreeQueue.empty()) {
        QuadNode* node = TreeQueue.front();
        TreeQueue.pop();
        CRect rc = node->rect;
        if ((rc.Width() > 2 && rc.Height() > 2) && !IsHomogenous(node)) {
            node->Subdivide();
            for (int i = 0; i < 4; i++) 
                TreeQueue.push(node->children[i]);
        }
    }
}
// subdivision 조건: 편차, 최대-최소 차이등의 criteria를 사용할 수 있다.
#define THRESH 15
bool QuadTree::IsHomogenous(QuadNode* node) {
    CRect rc = node->rect;
    double s = 0, ss = 0;
    int minval = 255, maxval = 0;
    for (int y = rc.top; y < rc.bottom; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y) + rc.left;
        for (int x = rc.left; x < rc.right; x++) {
            int a = *p++; 
            s += a; ss += a * a;
            maxval = max(maxval, a); minval = min(minval, a);
        }
    }
    int n = rc.Width() * rc.Height();
    node->mean = s / n;
    node->stdev = sqrt((ss - s * s / n) / n);
    return (node->stdev >= THRESH) ? false: true;
    //return (maxval - minval) > THRESH ? false: true;
};
// leaf node를 찾아 list에 넣음;
void QuadTree::FindLeaf(QuadNode* node) {
    if (node == NULL) return;
    if (!node->hasChildren) {
        LeafList.push_back(node);
        return;
    }
    for (int i = 0; i < 4; i++) 
        FindLeaf(node->children[i]);
}

사용 예:

void test_main(CRaster& raster, CRaster& out) {
    ASSERT(raster.GetBPP() == 8);  // gray input;
    QuadTree qtree(raster.GetRect(), raster);
    qtree.Build();
    qtree.FindLeaf(qtree.root);
    out.SetDimensions(raster.GetSize(), 24); //RGB-color;
    // segment image;
    for (int i = 0; i < qtree.LeafList.size(); i++) {
        QuadNode *node = qtree.LeafList[i];
        int a = int(node->mean + .5);
        FillRegion(out, node->rect, RGB(a, a, a));  //gray average;
        // draw rectangle with red color;
        DrawRectangle(out, node->rect, RGB(0, 0, 0xFF));
    }
    qtree.root->Clear();
    delete qtree.root;
}

여러 가지 함수를 Fouries series로 재구성할 때 적은 수의 부분합으로도 함수의 원래 모양이 잘 재현되는 경우가 있는 반면에 아무리 많은 항을 더하더라도 재현되지 않는 경우도 있다. 왜 이런 차이가 생기는 것일까? 급격한 모양의 변화가 없는 연속함수에 비해 중간에 값이 점프하는 불연속점이 있는 함수는 Fourier series로 전개를 했을 때 함수의 원래 모양이 잘 재현되지 않는다. 불연속점 전후로 값이 큰 쪽에서는 overshoot가 발생하고, 값이 작은 쪽에서는 undershoot가 발생한다. 이 돌출부는 Fourier series의 부분합의 개수를 늘리더라도 없어지지 않는다. 이처럼 불연속점에서 Fouries series에서 overshoot 또는 undershoot가 발생하는 현상을 Gibbs phenomenon이라 한다. 이는 불연속점에서 직교 고유 함수 전개를 할 때도 나타나는 일반적인 현상이다. 

구체적으로 다음 식으로 표현이 되는 square wave

\begin{gather}f(x)=\left\{ \begin{array}{cc}1, &  0<x <\pi \\ 0 , & x=0, \pm \pi \\ -1, & -\pi <x <0\end{array}\right. \end{gather}

을 아래처럼 홀수의 harmonics을 순차적으로 더해서 근사하는 경우를 살펴보자. 

\begin{gather}f_{k}(x) = \frac{4}{\pi} \sum_{j=1}^{k} \frac{ \sin (2j-1)x}{2j-1}=\frac{4}{\pi}\left(\sin (x) + \frac{\sin 3x}{3} + \cdots + \frac{\sin(2k-1)x}{2k-1}\right) .  \end{gather}

$k\to\infty$이면 Fourier 전개에 해당한다. 이 합은 $f(x)$의 불연속점인 $x=0, \pm\pi$을 제외한 영역에서 square wave $f(x)$로 균일하게 수렴한다. 이는 $k$를 증가시키면 불연속점 근방을 제외한 영역에서 부분합 $f_{k}(x)$의 그래프가 square wave의 그래프에 한없이 가깝게 접근함을 의미한다. 그러나 harmonics의 개수를 증가시키더라도 불연속 근방에서 돌출부의 overshoot나 undershoot의 세기는 변하지 않는다. 다만 $k$를 증가시키면 그 폭은 점점 불연속점을 향해 좁아진다. 불연속점$(x=0)$을 기준으로 좌우에서 첫째로 나타나는 극대/극소값의 크기는 $|f(0)+ \frac{f(0^+)-f(0^-)}{\pi} \int_0^\pi \frac{\sin(x)}{x} dx|$로 계산된다. 이 값과 $f(0^\pm)$과의 차이의 절반인 

$$\frac{1}{\pi} \int_{0}^{\pi} \frac {\sin (x)}{x} dx -\frac {1}{2}= 0.0894898\dots$$가 불연속점 $x=0$에서 overshoot나 undershoot의 세기에 해당한다. 이 값은 Wilbraham-Gibbs 상수로 불린다.