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

 

728x90

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

Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Harris Corner Detector  (0) 2022.04.07
Valley emphasis Otsu threshold  (0) 2022.02.23
Posted by helloktk
,

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

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

Moment-preserving Thresholding  (0) 2022.05.29
Minimum Cross Entropy Thresholding  (0) 2022.05.29
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
,

Harris response:

$$ R=\text{Det}(M)-k\text{Trace}(M)^{2}~\quad k \in [0.04,0.06]$$

참고: http://www.ipol.im/pub/art/2018/229/article_lr.pdf

number of corners=86
nonmaximal suppression를 하기 전 conrner의 response map

double PI = 4.0 * atan(1.0);
double Gaussian(double x, double y, double sigma) {
    double sigma2 = sigma * sigma;
    double t = (x * x + y * y) / (2 * sigma2);
    double n = 1.0 / (2.0 * PI * sigma2);
    return n * exp( -t );
};

struct Corner {
    double x, y;
    double strength;
    int    flag;
};
int HarrisCornerDetector(BYTE *img, int w, int h,
                    double sigma/*=1.2*/, double kvalue/*=0.04~0.06*/,
                    double minMeasure/*=1e-1*/,
                    std::vector<Corner>& corners)
{
    // Sobel gradient 3x3
    const int nn[] = { -w - 1, -w, -w + 1, -1, 0, 1, w - 1, w, w + 1}; // 3x3 window
    const int sobelX[] = { -1, 0, 1, -2, 0, 2, -1, 0, 1};
    const int sobelY[] = { -1, -2, -1, 0, 0, 0, 1, 2, 1};
    const int xmax = w - 1, ymax = h - 1;
    std::vector<double> Gx(w * h, 0);
    std::vector<double> Gy(w * h, 0);
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                int v = img[pos + nn[k]];
                sx += sobelX[k] * v;
                sy += sobelY[k] * v;
            }
            Gx[pos] = sx / (4 * 255);
            Gy[pos] = sy / (4 * 255);
        }
        pos++; // skip x=xmax;
    }
    // precompute the coefficients of the Gaussian filter
    const int radius = int(2 * sigma);
    const int window = 1 + 2 * radius;
    std::vector<double> gFilter(window * window);
    for (int j = -radius, pos = 0; j <= radius; j++)
        for (int i = -radius; i <= radius; i++, pos++)
            gFilter[pos] = Gaussian(i, j, sigma);

    // compute the Harris measure for each pixel of the image
    std::vector<double> harrismap(w * h, 0);
    double hmax = 0;
    for (int y = 0, pos = 0; y < h; y++) {
        for (int x = 0; x < w; x++, pos++) {
            // Convolve gradient with Gaussian filter:
            double filtersum = 0, sxx = 0, syy = 0, sxy = 0;
            for (int yk = y - radius, wpos = 0; yk <= y + radius; yk++) {
                for (int xk = x - radius; xk <= x + radius; xk++, wpos++) {
                    if (yk < 0 || yk >= h) continue;
                    if (xk < 0 || xk >= w) continue;
                    // Gaussian weight
                    double gauss = gFilter[wpos];
                    filtersum += gauss;
                    // convolution;
                    int loc = yk * w + xk; //image location;
                    sxx += gauss * Gx[loc] * Gx[loc];
                    syy += gauss * Gy[loc] * Gy[loc];
                    sxy += gauss * Gx[loc] * Gy[loc];
                }
            }
            sxx /= filtersum;
            syy /= filtersum;
            sxy /= filtersum;
            // Harris corner measure = Det(M) - kvalue.(Trace(M))^2
            double hmeasure = sxx * syy - sxy * sxy - kvalue * (sxx + syy) * (sxx + syy);
            if (hmeasure <= 0) continue;
            if (hmeasure > hmax) hmax = hmeasure;
            harrismap[pos] = hmeasure;
        }
    }
    // log scale
    double scale = 255 / log(1 + hmax);
    for (int pos = w * h; pos-- > 0;)
        harrismap[pos] = scale * log(1 + harrismap[pos]);
    // Nonmaximal suppression: keep only local maxima;
    minMeasure *= 255; //255 = the maximum value of Harris measure(in log scale);
    const int nn8[] = { -w - 1, -w, -w + 1, 1, w + 1, w, w - 1, -1}; // 8-neighbors;
    for (int y = 1, pos = w; y < ymax; y++) {
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            // thresholding: Harris measure > epsilon
            double hmeasure = harrismap[pos];
            if (hmeasure <= minMeasure) continue;
            bool isLocalMax = true;
            for (int k = 0; k < 8; k++) {
                double hk = harrismap[pos + nn8[k]];
                if (hk >= hmeasure) {
                    isLocalMax = false;
                    break;
                }
            }
            if (isLocalMax) {
                Corner c;
                c.x = x, c.y = y, c.strength = hmeasure, c.flag = 1;
                corners.push_back(c);
            }
        }
        pos++; //skip x=w-1;
    }

    // remove corners too close to each other: keep the highest measure
    double minDistance = 8;
    std::vector<Corner>::iterator iter = corners.begin();
    while (iter != corners.end()) {
        for (int i = 0; i < corners.size(); i++) {
            Corner& c1 = corners[i];
            if (c1.flag == 0) continue;
            if ((c1.x == iter->x && c1.y == iter->y) ) continue;
            if (sqrt((c1.x-iter->x)*(c1.x-iter->x) + (c1.y-iter->y)*(c1.y-iter->y)) > minDistance) 
                continue;
            if (c1.strength < iter->strength) continue;
            iter->flag = 0;    // iter = corners.erase(iter); berased = true;
            break;
        }
        iter++;
    }
    return 1;
};
728x90

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

Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Valley emphasis Otsu threshold  (0) 2022.02.23
Image rotation by FFT  (0) 2022.02.18
FFT를 이용한 영상의 미분  (0) 2022.02.12
Posted by helloktk
,

Otsu 이진화에서 threshold 값의 선택은 전경과 배경 사이의 분산을 최대로 하는 픽셀 값이다. 영상의 히스토그램이 bimodal인 경우는 잘 동작하지만, unimodal인 영상의 경우는 제대로 처리가 안된다. Valley emphasis 방법은 Otsu 방법을 사용하면서 histogram profile의 valley에 가중치를 더해주는 방식으로 threshold 값을 찾는다.

\begin{gather} \text{Ostu's method: } ~~~\text{threshold}= \underset{0\le t< 256}{\text{argmax}}~\left(\omega_1(t) m_1(t)^2 + \omega_2(t) m_2^2(t)\right)\\  \text{valley-emphasis: }~~\text{threshold}=\underset{0\le t < 256} {\text{argmax}}~  (1-p_t)\left(\omega_1(t) m_1(t)^2 + \omega_2(t) m_2^2(t)\right) \end{gather}

Bimodal 분포의 영상은 Otsu와 거의 같은 결과를 준다. 두 개 이상의 클래스로 분리도 쉽게 구현이 된다.

 

참고: http://msn.iecs.fcu.edu.tw/report/data/ori_paper/2006-12-21/2006-Automatic%20thresholding%20for%20defect%20detection.pdf

영상은 참고 논문에서 추출함. 중간=Otsu, 마지막=valley-emphasis

 

 

int valley_emphasis_Otsu(int hist[], int N) {
    int istart = 0, iend = N - 1;
    while (!hist[istart]) istart++;
    while (!hist[iend])   iend--;
    double st = 0, sxt = 0;           // total pdf, total cdf;
    for (int i = istart; i <= iend;i++) { 
        st += hist[i]; 
        sxt += double(i) * hist[i];
    }
    int winner = istart; 
    double maxgain = 0;
    double s0 = 0, sx0 = 0;
    for(int i = istart; i <= iend; i++) {
        s0  += hist[i];				// 1-pdf;
        sx0 += double(i) * hist[i];		// 1-cdf;
        double s1  = st - s0;			// 2-pdf
        double sx1 = sxt - sx0;			// 2-cdf
        double m0  = sx0 / s0;			// E(X1);
        double m1  = sx1 / s1;			// E(X2);
        double gain = (1 - hist[i] / st) * (s0 * m0 * m0 + s1 * m1 * m1);
        if (gain > maxgain) {
            maxgain = gain; 
            winner = i;
        }
    }
    return winner;
};

 

728x90

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

Quadtree Segmentation  (0) 2022.05.21
Harris Corner Detector  (0) 2022.04.07
Image rotation by FFT  (0) 2022.02.18
FFT를 이용한 영상의 미분  (0) 2022.02.12
SVD Fitting  (0) 2022.02.07
Posted by helloktk
,

평면을 한 점 $(x, y)$이 원점을 축으로 하는 회전에 의해서 $(x', y')$ 위치로 옮겼다면 두 점 사이의 관계는 

$$ \begin{pmatrix} x' \\ y'\end{pmatrix} = \begin{pmatrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta\end{pmatrix}\begin{pmatrix} x\\ y\end{pmatrix} = \mathbf{R}_\theta \begin{pmatrix} x\\ y\end{pmatrix}$$

로 표현된다. 여기서 $\theta$는 시계방향으로 회전한 각이다. 이 회전 변환은 3번의 shear 변환의 곱으로 다음처럼 쓸 수 있다.

$$ \mathbf{R}_\theta = \mathbf{S}_x. \mathbf{S}_y . \mathbf{S}_x$$

$$ \mathbf{S}_x =\begin{pmatrix} 1 & \tan (\theta/2) \\0 &1\end{pmatrix}, ~~\mathbf{S}_y =\begin{pmatrix}1 & 0 \\ -\sin \theta & 1\end{pmatrix}$$

 

이제 Fourier변환과 shear 변환 사이의 관계를 알아보자. 일반적인 2차원 Fourier 변환은

\begin{gather} F(u, v)  =\iint   f(x, y) e^{-2\pi i (ux + vy)}dxdy={\cal F}[ f(x, y)] \\= {\cal F_x}\left[{\cal F_y} [f(x, y)]\right]={\cal F_y}\left[  {\cal  F_x}  [f(x, y)]\right] = F_1 (u) F_2 (v)\end{gather}

여기서 $\cal F_{x}, F_y$는 각각 $x$, $y$-방향 Fourier 변환이고 순서에 상관이 없다.

 

이제 $f(x, y)$에 $x$-방향으로 $a$ 만큼 shear 변환을 한 결과가 $g(x, y)$일 때,

$$ g(x, y) = f(x + ay, y)$$

로 쓸 수 있고, 여기에 $x$ 방향으로 Fourier 변환을 하면

$$ G_1(u, y) = {\cal F_x}[g(x, y)]= {\cal F_x}[f(x+ay, y)] = F_1(u, y) e^{-2\pi i ua y} =e^{-2\pi i  ua y}  {\cal F_x}[ f(x, y)]$$ 

임을 알 수 있다. 즉, shear 변환된 함수의 Fourier 변환은 원함수의 Fourier변환에 shear 파라미터에 의존하는 위상 요소 $e^{-2\pi i  ua y}$만 곱해주면 쉽게 구할 수 있음을 알 수 있다. 또한 shear 변환에 의해서 스펙트럼이 위상만 변화고 크기에는 변화가 생기지 않음도 알 수 있다.

참고: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.190.4473&rep=rep1&type=pdf

 

아래의 코드는 이미지의 중앙을 기준으로 회전을 한 경우를 보여준다: $\theta = 60^\circ$. 중간 단계에서 artifact를 없애기 위해 원본 이미지의 가장자리를 $0$으로 채워 두 배 크기로 만든 후 수행한다. 임의의 지점에 대한 회전도 shear 변환 코드를 조금 손보면 쉽게 구현할 수 있다.

중간단계: 

$x$ 방향 shear 변환 결과
추가적인 $y$ 방향 shear 변환 결과

int fft1d(int dir, int nn, double x[], double y[]); /* https://kipl.tistory.com/22 */
#define HORIZONTAL 0
#define VERTICAL   1
/* perform shear transform along "axis" direction; */
void shearImage(int axis, double shear,
                int width, int height, double *input, double *output) 
{
    int size, halfsz, othersize, halfothersz;
    double TWOPI = 8.0 * atan(1.);
    /* size is the size of the image in the direction of the axis of shear. 
    ** othersize the size of the image in the orthogonal direction 
    */
    if (axis == HORIZONTAL) {
        halfsz = (size = width) >> 1;
        halfothersz = (othersize = height) >> 1;
    }
    else {
        halfsz = (size = height) >> 1;
        halfothersz = (othersize = width) >> 1;
    }

    std::vector<double> re_sig(size, 0);
    std::vector<double> im_sig(size, 0);
    for (int i = 0; i < othersize; i++){
        if (axis == HORIZONTAL)
            memcpy(&re_sig[0],  &input[i * size], size * sizeof(double));
        else
            for (int j = 0; j < size; j++) re_sig[j] = output[j * othersize + i];

        std::fill(im_sig.begin(), im_sig.end(), 0);  // im: padding 0;
        fft1d(1, size, &re_sig[0], &im_sig[0]);
        double shift = - (i - halfothersz - .5) * shear  * TWOPI;
        // No need to touch j = 0;
        for (int j = 1; j < halfsz; j++) {
            double cc = cos(j * shift / size);
            double ss = sin(j * shift / size);
            double reval = re_sig[j];
            double imval = im_sig[j];
            re_sig[j] = cc * reval - ss * imval;
            im_sig[j] = ss * reval + cc * imval;
            re_sig[size - j] =  re_sig[j];
            im_sig[size - j] = -im_sig[j];
        }
        re_sig[halfsz] = cos(shift/2) * re_sig[halfsz] - sin(shift/2) * im_sig[halfsz];
        im_sig[halfsz] = 0.;

        fft1d(-1, size, &re_sig[0], &im_sig[0]);
        if (axis == HORIZONTAL)
            memcpy(&output[i * size], &re_sig[0], size * sizeof(double));
        else
            for (int j = 0; j < size; j++) output[j * othersize + i] = re_sig[j];
    }
};
 
728x90

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

Harris Corner Detector  (0) 2022.04.07
Valley emphasis Otsu threshold  (0) 2022.02.23
FFT를 이용한 영상의 미분  (0) 2022.02.12
SVD Fitting  (0) 2022.02.07
Color Histogram Equalization  (4) 2022.02.07
Posted by helloktk
,