적용 예는 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

댓글을 달아 주세요