'Watershed'에 해당되는 글 2건

  1. 2010.03.19 Watershed Algorithm Code
  2. 2008.05.21 Watershed Algorithm 적용의 예 (2)

적용 예는 http://kipl.tistory.com/1 
#include <vector>
#include <queue>
#define WMASK   -2
#define WINIT   -1
#define WSHED    0
#define FICTION ((UINT)(INT)-1)
/**
*/
void watershed_sorting_step(const std::vector<BYTE> & img,
                            int &hmin, int &hmax, UINT fr[256],
                            std::vector<UINT> & sort)
{
    UINT freq[ 256];
    /* Obtain frequency distribution */
    memset( freq, 0, 256*sizeof(UINT));
    for (int i = 0; i < img.size(); i++)
        freq[img[i]]++;
    memcpy( fr, freq, sizeof(freq));
   
    /* Calculate hmin & hmax */
    hmin = 0;  
    while ( hmin < 256 && freq[hmin] == 0) hmin++;
    hmax = 255;
    while ( hmax > 0 && freq[hmax] == 0)   hmax--;
   
    /* Obtain cumulative frequency distribution */
    for ( i = 1; i < 256; i++)
        freq[i] += freq[i-1];
   
    /* allocate and fill an array of sorted positions */
    for (int pos = 0; pos < img.size(); pos++) {
        int val = img[pos];
        // freq[100]은 cumulative histogram이므로 100값까지 같은 픽셀의 총수다.
        // freq[100]-freq[99]이 정확히 100값을 갖는 픽셀의 총수이다;
        // [ 전체 이미지                                          ]
        // 0----------------->1-------->2-------->........>255----]
        //이런식으로 파티션이 되어있으므로, 특정한 값이 한 번 나타나는 경우에
        int num = --freq[val] ;
        sort[num] = pos;
    };
    //sort-[000000000111112222222233333....................255255255]
    //각각의 위치는 그 값이 이미지에서 나타난 순서대로, 그 픽셀의 위치를 나타낸다.
}
#define NN(p) ((p)-width)   
#define SS(p) ((p)+width)
#define EE(p) ((p)+1)
#define WW(p) ((p)-1)
#define NW(p) ((p)-width-1)
#define NE(p) ((p)-width+1)
#define SW(p) ((p)+width-1)
#define SE(p) ((p)+width+1)
//
void Watershed8UC(CRaster & raster, CRaster & res)
{
    int neighborhood = 8;
    int N, width, height, maxx, maxy;
    int hmin, hmax;
 
    UINT maxdist = 1;
    UINT freq[256];
    UINT pbis[8];

    /* Initialization */
    CSize sz = raster.GetSize();
    width = sz.cx;
    height = sz.cy;
    maxx = width - 1;
    maxy = height - 1;
    N = width * height;

    /* Convert our image to non-aligned data structure */
    std::vector<BYTE> inp(N) ;
    for (int y = 0; y < height; y++) {
        BYTE *p = (BYTE*)raster.GetLinePtr(y) ;
        memcpy(&inp[y * width], p, width);
    }

    /* Obtain sorted array of pixel positions */
    std::vector<UINT> sorted(N);
    watershed_sorting_step(inp, hmin, hmax, freq, sorted);

    /* Check the correctness of sorting */
    if (1) {
        BYTE last = 0;
        for (int i = 0; i < N; i++)
            if (inp[sorted[i]] < last)
                TRACE("%s: incorrectly sorted", "check");
            else
                last = inp[sorted[i]];
    }

    /* create and initialize output data array */
    std::vector<int> out(N, WINIT);

    /* create and initialize distance work image */
    std::vector<UINT> distance(N, 0);
   
    /* create and initialize the queue */
    std::queue<UINT> Q;

    /* other initialization */
    int label = 0;
    int i = 0, p;

    /* main loop */
    for (int h = hmin; h <= hmax; h++) {
        /*=== geodesic SKIZ of level h-1 inside level h */
        /*=== for every pixel p such that image(p) == h */
        int n = freq[h];
        if (neighborhood == 4) {
            while (n--) {
                UINT p = sorted[i++];
                if (inp[p] != h)
                    TRACE("sort assertion failed: %d(%d) != %d", inp[p], p, h);
                int x = p % width;
                int y = p / width;
                out[p] = WMASK;
                if ((x < maxx && out[EE(p)] >= WSHED) ||
                    (x > 0    && out[WW(p)] >= WSHED) ||
                    (y > 0    && out[NN(p)] >= WSHED) ||
                    (y < maxy && out[SS(p)] >= WSHED))
                {
                    distance[p] = 1;
                    Q.push(p);
                }
            }
        } else {                  /* neighborhood == 8 */
            while (n--) {
                UINT p = sorted[i++];
                int x = p % width;
                int y = p / width;
                out[p] = WMASK;
                if ((x < maxx && (out[EE(p)] >= WSHED || (y > 0 && out[NE(p)] >= WSHED) || (y < maxy && out[SE(p)] >= WSHED))) ||
                    (x > 0    && (out[WW(p)] >= WSHED || (y > 0 && out[NW(p)] >= WSHED) || (y < maxy && out[SW(p)] >= WSHED))) ||
                    (y > 0    &&  out[NN(p)] >= WSHED) || (y < maxy && out[SS(p)] >= WSHED))
                {
                    distance[p] = 1;
                    Q.push(p);//FIFO_ADD(p);
                }
            }
        }
        int dist = 1;
        Q.push(FICTION);
        while (1) {
            UINT p=Q.front(); Q.pop();
            if (p == FICTION) {
                if (Q.empty())
                    break;
                Q.push(FICTION);
                dist++;
                if (dist > maxdist)
                    maxdist = dist;
                p=Q.front();Q.pop();
            }
            if (Q.empty())
                TRACE("assertion failed: !FIFO_EMPTY");
            int x = p % width;
            int y = p / width;
            int npbis = 0 ;
            if (neighborhood == 4) {
                if (x > 0)
                    pbis[npbis++] = WW(p);
                if (x < maxx)
                    pbis[npbis++] = EE(p);
                if (y > 0)
                    pbis[npbis++] = NN(p);
                if (y < maxy)
                    pbis[npbis++] = SS(p);
            } else {            /* neighborhood == 8 */
                if (x > 0) {
                    pbis[npbis++] = WW(p);
                    if (y > 0)
                        pbis[npbis++] = NW(p);
                    if (y < maxy)
                        pbis[npbis++] = SW(p);
                }
                if (x < maxx) {
                    pbis[npbis++] = EE(p);
                    if (y > 0)
                        pbis[npbis++] = NE(p);
                    if (y < maxy)
                        pbis[npbis++] = SE(p);
                }
                if (y > 0)
                    pbis[npbis++] = NN(p);
                if (y < maxy)
                    pbis[npbis++] = SS(p);
            }
            /*=== for every pixel p' belonging to Ng(p) */
            while (--npbis >= 0) {
                UINT pbs = pbis[npbis];
                if (distance[pbs] < dist && out[pbs] >= WSHED) {
                    if (out[pbs] > 0) {
                        if (out[p] == WMASK || out[p] == WSHED)
                            out[p] = out[pbs];
                        else if (out[p] != out[pbs])
                            out[p] = WSHED;
                    } else if (out[p] == WMASK)
                        out[p] = WSHED;
                } else if (out[pbs] == WMASK && distance[pbs] == 0) {
                    distance[pbs] = dist + 1;
                    Q.push(pbs);
                }
            }
        }                       /* while(1) */
        /*=== checks if new minima have been discovered */
        n = freq[h];
        i -= n;
        while (n--) {
            UINT p = sorted[i++];
            distance[p] = 0;
            if (out[p] == WMASK) {
                label++;
                Q.push(p);
                out[p] = label;
                while (!Q.empty()){
                    UINT pp=Q.front();Q.pop();
                    int x = pp % width;
                    int y = pp / width;
                    int npbis = 0;
                    if (neighborhood == 4) {
                        if (x > 0)
                            pbis[npbis++] = WW(pp);
                        if (x < maxx)
                            pbis[npbis++] = EE(pp);
                        if (y > 0)
                            pbis[npbis++] = NN(pp);
                        if (y < maxy)
                            pbis[npbis++] = SS(pp);
                    } else {    /* neighborhood == 8 */
                        if (x > 0) {
                            pbis[npbis++] = WW(pp);
                            if (y > 0)
                                pbis[npbis++] = NW(pp);
                            if (y < maxy)
                                pbis[npbis++] = SW(pp);
                        }
                        if (x < maxx) {
                            pbis[npbis++] = EE(pp);
                            if (y > 0)
                                pbis[npbis++] = NE(pp);
                            if (y < maxy)
                                pbis[npbis++] = SE(pp);
                        }
                        if (y > 0)
                            pbis[npbis++] = NN(pp);
                        if (y < maxy)
                            pbis[npbis++] = SS(pp);
                    }
                    while (--npbis >= 0) {
                        UINT pbs = pbis[npbis];
                        if (out[pbs] == WMASK) {
                            Q.push(pbs);
                            out[pbs] = label;
                        }
                    }
                }
            }
        }
    }

    if (!Q.empty())//FIFO_EMPTY)
        TRACE("$s: queue is not empty - can't be", "");

    /* Convert the result to a suitable form */
    res.SetDimensions(width, height, 8, 1) ;
    for (y = 0, p = 0; y < height; y++) {
        for (int x = 0; x < width; x++, p++) {
            // int p = x + width * y;
            if (out[p] == WMASK) TRACE("%s: %d,%d has mask", "", x, y);
            if (out[p] == WINIT) TRACE("%s: %d,%d has init", "", x, y);
            if (out[p] == WSHED) continue;
            if (neighborhood == 4) {
                if ((x < maxx && out[EE(p)] > WSHED && out[EE(p)] < out[p]) ||
                    (x > 0    && out[WW(p)] > WSHED && out[WW(p)] < out[p]) ||
                    (y > 0    && out[NN(p)] > WSHED && out[NN(p)] < out[p]) ||
                    (y < maxy && out[SS(p)] > WSHED && out[SS(p)] < out[p]))
                    out[p] = WSHED;
            } else {
                if ((x < maxx && ((out[EE(p)] > WSHED && out[EE(p)] < out[p]) ||
                    (y > 0    && out[NE(p)] > WSHED   && out[NE(p)] < out[p])||
                    (y < maxy && out[SE(p)] > WSHED   && out[SE(p)] < out[p]))) ||
                    (x > 0    && ((out[WW(p)] > WSHED && out[WW(p)] < out[p])||
                    (y > 0    && out[NW(p)] > WSHED   && out[NW(p)] < out[p])||
                    (y < maxy && out[SW(p)] > WSHED   && out[SW(p)] < out[p]))) ||
                    (y > 0    && out[NN(p)] > WSHED   && out[NN(p)] < out[p]) ||
                    (y < maxy && out[SS(p)] > WSHED   && out[SS(p)] < out[p]))
                    out[p] = WSHED;
            }
        }
    }
    for (y = 0; y < height; y++)
    {
        BYTE *p = (BYTE*) res.GetLinePtr(y) ;
        for (int x = 0; x < width; x++) {
            UINT a=out[x + width * y] ;
            p[x] = a;//((a == WSHED) ? 255 : 0);
        }
    }
    return ;
}

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

Blur Detection  (0) 2010.05.25
Savitzky-Golay smoothing filter  (2) 2010.03.24
Watershed Algorithm Code  (0) 2010.03.19
Retinex 알고리즘  (11) 2010.02.03
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Image Morphing  (0) 2010.01.24
Posted by helloktk
TAG Watershed

댓글을 달아 주세요

(1) cell image : 원본 영상을 이진화(Otsu 알고리즘)시킨  결과이다. 두 군데서 셀이 겹쳤다. 단순히 connected component labelling을 적용하여서는 이것을 분리 할 수 없다.

사용자 삽입 이미지

(2) distance transform : distance 변환 결과에 블러링을 추가한 결과이다. distance 변환은 셀 외부는 셀로부터의 거리를, 셀내부는 배경으로부터의 거리의 음의 값을 취하고, 전체적으로 다시 리스케일링한 것이다.  블러링은 워터쉐드 알고리즘이 보다 정확히 동작하는데 필요하다.

사용자 삽입 이미지

(3) watershed segmentation: 분할된 영역의 라벨이 나온다(경계라벨=0). 이 라벨을 가지고  false coloring을 한 결과이다. 이 알고리즘은 "The Watershed Transform: Definitions, Algorithms and Parallelization Strategies", Jos B.T.M. Roerdink and Arnold Meijster에 따라서 구현이 된것이다. 픽셀연결성은 8방향을 이용하였다.

사용자 삽입 이미지

(4) final cell segmentation; watershed 결과를  마스크로 이용하여서 cell이미지를 분할한 것이다. 겹친 cell들이 분리되었다.

사용자 삽입 이미지

다른 예:

사용자 삽입 이미지


사용자 삽입 이미지


사용자 삽입 이미지


/**
** http://blog.naver.com/helloktk/80051779331 에서 옮긴 자료.
*/


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

Gaussian Mixture Model  (2) 2008.06.07
Rasterizing Voronoi Diagram  (0) 2008.05.26
RANSAC Algorithm  (0) 2008.05.24
Contour Tracing  (0) 2008.05.22
Gausssian Scale Space  (0) 2008.05.22
Watershed Algorithm 적용의 예  (2) 2008.05.21
Posted by helloktk

댓글을 달아 주세요

  1. 2013.06.26 11:24  댓글주소  수정/삭제  댓글쓰기

    비밀댓글입니다

  2. helloktk 2013.06.27 23:28 신고  댓글주소  수정/삭제  댓글쓰기

    컴퓨터가 망가지는 통에 요 몇년간 작업한 내용은 다 사라진 상태입니다. 이 예제는 여러가지 알고리즘을 조합해서 쓰긴 했지만, 이 블로그와 네이버블로그를 보면 아마 필요한 부분은 다 찾을 수 있을 것입니다.