728x90

Statistical region merging은 이미지의 픽셀을 일정한 기준에 따라 더 큰 영역으로 합병하는 bottom-up 방식의 과정이다. 두 영역 $R_1$과 $R_2$가 하나의 영역으로 합병이 되기 위해서는 두 영역의 평균 픽셀값의 차이가 

$$ g \sqrt{ \frac{\ln(2/\delta)}{2Q}  \Big( \frac{1}{|R_1|}+ \frac{1}{|R_2|}\Big) }$$

를 넘지 않아야 한다. $g=256$으로  gray 레벨의 갯수를 의미하고, $|R_i|$는 $R_i$ 영역에 포함된 픽셀 수를 나타낸다. $\delta$는 작은 수로 이미지의 픽셀 수의 제곱에 반비례한다. 보통 $\delta = 1/(6\times \text{width}\times \text{height})^2$로 선택한다. $Q$는 이미지의 통계적인 복잡성을 정량화하는 양으로 이 알고리즘에서는 외부에서 설정이 되는 값이다. 낮은 $Q$값을 선택하면 분할된 영역의 수가 작아지고(undersegmentation), 반대로 높은 $Q$ 값을 입력하면 분할된 영상에 너무 많은 영역이 나타나게 된다(oversegmentation).

Ref:

https://en.wikipedia.org/wiki/Statistical_region_merging

http://www.lix.polytechnique.fr/~nielsen/Srmjava.java

 

Srm::Srm(int width, int height, BYTE *image) {
    this->width  = width;
    this->height = height;
    int n        = width * height;
    this->count = new int[n];
    this->Ravg  = new float [n];
    this->Gavg  = new float [n];
    this->Bavg  = new float [n];
    this->image = image;
    // disjoint sets with n pixels;
    this->UF = new Universe(n);
    // initialize to each pixel a leaf region;
    for (int i = 0, pos = 0; i < n; i++, pos += 3) {
        count[i] = 1;
        Bavg[i] = image[pos    ];
        Gavg[i] = image[pos + 1];
        Ravg[i] = image[pos + 2];
    }
    this->Q = 32;		// adjustable.					    
    this->g = 256.0;
    this->logDelta = 2. * log(6.0 * n);
}
bool Srm::Predicate(int reg1, int reg2) {
    double dR = (Ravg[reg1] - Ravg[reg2]); dR *= dR;
    double dG = (Gavg[reg1] - Gavg[reg2]); dG *= dG;
    double dB = (Bavg[reg1] - Bavg[reg2]); dB *= dB;
    double logreg1 = min(g, count[reg1]) * log(1.0 + count[reg1]);
    double logreg2 = min(g, count[reg2]) * log(1.0 + count[reg2]);
    double factor = g * g / (2.0 * Q);
    double dev1 = factor * (logreg1 + logDelta) / count[reg1] ;
    double dev2 = factor * (logreg2 + logDelta) / count[reg2] ;
    double dev = dev1 + dev2;
    return ( (dR < dev) && (dG < dev) && (dB < dev) );
}
void Srm::Merge(int root1, int root2) {
    if (root1 == root2) return;
    int w1 = count[root1], w2 = count[root2];
    int root = UF->Union(root1, root2);
    //update the merged region;
    count[root] = w1 + w2;
    double count_sum = w1 + w2;
    Ravg[root] = (w1 * Ravg[root1] + w2 * Ravg[root2]) / count_sum;
    Gavg[root] = (w1 * Gavg[root1] + w2 * Gavg[root2]) / count_sum;
    Bavg[root] = (w1 * Bavg[root1] + w2 * Bavg[root2]) / count_sum;
}
Edge* Srm::Pairs(int nedge) {
    // 4-connectivity;
    int ymax = height - 1, xmax = width - 1;
    Edge* edgeList = new Edge[nedge];
    int cnt = 0;
    for (int y = 0; y < ymax; y++) {
        for (int x = 0; x < xmax; x++) {
            int pos = y * width + x;
            int b1 = image[3 * pos + 0];
            int g1 = image[3 * pos + 1];
            int r1 = image[3 * pos + 2];
            //right: x--x
            edgeList[cnt].r1 = pos;     //current
            edgeList[cnt].r2 = pos + 1; //right
            int bdiff = abs(b1 - image[3 * (pos + 1) + 0]);
            int gdiff = abs(g1 - image[3 * (pos + 1) + 1]);
            int rdiff = abs(r1 - image[3 * (pos + 1) + 2]);
            edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff) ;
            //below: x
            //       |
            //       x
            edgeList[cnt].r1 = pos;
            edgeList[cnt].r2 = pos + width;
            bdiff = abs(b1 - image[3 * (pos + width) + 0]);
            gdiff = abs(g1 - image[3 * (pos + width) + 1]);
            rdiff = abs(r1 - image[3 * (pos + width) + 2]);
            edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
        }
    }
    //x=width-1;
    for (int y = 0; y < ymax; y++) {
        int pos = y * width + (width - 1); // (x,y) = (width-1, y)
        // x
        // |
        // x
        edgeList[cnt].r1 = pos;
        edgeList[cnt].r2 = pos + width;
        int bdiff = abs((int)image[3 * pos + 0] - image[3 * (pos + width) + 0]);
        int gdiff = abs((int)image[3 * pos + 1] - image[3 * (pos + width) + 1]);
        int rdiff = abs((int)image[3 * pos + 2] - image[3 * (pos + width) + 2]);
        edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
    }
    //y=height-1;
    for (int x = 0; x < xmax; x++) {
        int pos = (height - 1) * width + x;      //(x,y)=(x, height-1);
        //right; x--x
        edgeList[cnt].r1 = pos;
        edgeList[cnt].r2 = pos + 1;
        int bdiff = abs((int)image[3 * pos + 0] - image[3 * (pos + 1) + 0]);
        int gdiff = abs((int)image[3 * pos + 1] - image[3 * (pos + 1) + 1]);
        int rdiff = abs((int)image[3 * pos + 2] - image[3 * (pos + 1) + 2]);
        edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
    }
    return edgeList;
}
int Srm::Segment() {
    // 4-connectivity 
    int nedge = 2 * (width - 1) * (height - 1) + (height - 1) + (width - 1);
    Edge* edgeList = Pairs(nedge);
    BucketSort(edgeList, nedge);
    for (int i = 0; i < nedge; i++) {
        int root1 = UF->Find(edgeList[i].r1);
        int root2 = UF->Find(edgeList[i].r2);
        if ((root1 != root2) && (Predicate(root1, root2)))
            Merge(root1, root2);
    }
    delete [] edgeList;
    int rgn_count = 0;
    for (int node = width * height; node-- > 0;)
        if (UF->IsRoot(node)) rgn_count++;
    return rgn_count;
}
// sorting with buckets; returns an ordered edgeList;
void BucketSort(Edge* &edgeList, int n) {
    int hist[256] = {0}, chist[256];
    for (int i = 0; i < n; i++) hist[edgeList[i].diff]++;
    // cumulative histogram
    chist[0] = 0;  // Note, chist[0] ne hist[0];
    for (int i = 1; i < 256; i++)
        chist[i] = chist[i - 1] + hist[i - 1];

    Edge *ordered = new Edge [n];
    for (int i = 0; i < n; i++)
        ordered[chist[pair[i].diff]++] = pair[i];        
    delete[] edgeList;
    edgeList = ordered;
}

 

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

Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Image Inpainting  (0) 2022.04.19
Harris Corner Detector  (0) 2022.04.07
Posted by helloktk

댓글을 달아 주세요

728x90

영상의 히스토그램($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;
}

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

Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Image Inpainting  (0) 2022.04.19
Harris Corner Detector  (0) 2022.04.07
Posted by helloktk

댓글을 달아 주세요

728x90

$$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 $$

중간=Otsu, 마지막=MCE

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

 

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

Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Image Inpainting  (0) 2022.04.19
Harris Corner Detector  (0) 2022.04.07
Posted by helloktk

댓글을 달아 주세요

728x90

이미지를 균일한 컬러 영역으로 분할하는 것은 압축이나 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;
}

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

Moment-preserving Thresholding  (0) 2022.05.29
Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Image Inpainting  (0) 2022.04.19
Harris Corner Detector  (0) 2022.04.07
Valley emphasis Otsu threshold  (0) 2022.02.23
Posted by helloktk

댓글을 달아 주세요

728x90

Inpainting is a conservation process where damaged, deteriorated, or missing parts of an artwork are filled in to present a complete image.

ref: https://www.inf.ufrgs.br/~oliveira/pubs_files/inpainting.pdf

image source: https://github.com/Mugichoko445/Fast-Digital-Image-Inpainting

static const double _a = 0.073235;
static const double _b = 0.176765;
static const double Ker[]={ _a, _b, _a, _b, 0, _b, _a, _b, _a};
// mask는 color별로 지정할 수 있음;
//
void inPainting(CRaster& raster, CRaster& mask, CRaster& result) {
    CSize sz = raster.GetSize();
    double bsum = 0, gsum = 0, rsum = 0;
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            bsum += *p++;
            gsum += *p++;
            rsum += *p++;
        }
    }
    bsum /= (sz.cx * sz.cy);
    gsum /= (sz.cx * sz.cy);
    rsum /= (sz.cx * sz.cy);
    result = raster; // average color;
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)result.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            *p++ = int(bsum);
            *p++ = int(gsum);
            *p++ = int(rsum);
        }
    }
    std::vector<double> red(sz.cx * sz.cy);
    std::vector<double> blu(sz.cx * sz.cy);
    std::vector<double> grn(sz.cx * sz.cy);
    for (int y = 0, k = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *m = (BYTE *)mask.GetLinePtr(y);
        BYTE *a = (BYTE *)result.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++, k++) {
            blu[k] = (*m / 255.) * *p++ + (1 - *m++ / 255.) * *a++; 
            grn[k] = (*m / 255.) * *p++ + (1 - *m++ / 255.) * *a++;
            red[k] = (*m / 255.) * *p++ + (1 - *m++ / 255.) * *a++;
        }
    }

    int maxIters = 100;
    // color 별로 convolution;
    static int nn[] = {-sz.cx - 1, -sz.cx, -sz.cx + 1, -1, 0, 1, sz.cx - 1, sz.cx, sz.cx + 1};
    std::vector<double> bbuf(sz.cx * sz.cy);
    std::vector<double> gbuf(sz.cx * sz.cy);
    std::vector<double> rbuf(sz.cx * sz.cy);
    for (int iter = 0; iter < maxIters; ++iter) {
        // blu, grn, red를 임시 버퍼에 복사를 한 후; 
        // 복사본에 대해서 convolution 수행;
        std::copy(blu.begin(), blu.end(), bbuf.begin());
        std::copy(grn.begin(), grn.end(), gbuf.begin());
        std::copy(red.begin(), red.end(), rbuf.begin());
        for (int y = 1; y < sz.cy - 1; ++y) {	
            BYTE *m = (BYTE *)mask.GetLinePtr(y) + 3;
            double *pb = &bbuf[y * sz.cx + 1];
            double *pg = &gbuf[y * sz.cx + 1];
            double *pr = &rbuf[y * sz.cx + 1];
            for (int x = 1; x < sz.cx - 1; ++x){
                if (!*m++) { //blue
                    double s = 0;
                    for (int i = 0; i < 9; i++)
                        s += pb[nn[i]] * Ker[i]; 
                    blu[y * sz.cx + x] = s;
                }
                if (!*m++) { //green;
                    double s = 0;
                    for (int i = 0; i < 9; i++)
                        s += pg[nn[i]] * Ker[i]; 
                    grn[y * sz.cx + x] = s;
                }
                if (!*m++) { //red;
                    double s = 0;
                    for (int i = 0; i < 9; i++)
                        s += pr[nn[i]] * Ker[i]; 
                    red[y * sz.cx + x] = s;
                }
            }
        }
    }

    for (int y = 0, i = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)result.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++, i++) {
            *p++ = int(bbuf[i]);
            *p++ = int(gbuf[i]);
            *p++ = int(rbuf[i]);
        }
    }
}

 

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

Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Image Inpainting  (0) 2022.04.19
Harris Corner Detector  (0) 2022.04.07
Valley emphasis Otsu threshold  (0) 2022.02.23
Image rotation by FFT  (0) 2022.02.18
Posted by helloktk

댓글을 달아 주세요