Geometry & Recognition

1. gradient 대신 negate시킨 distance map을 사용한다.
2. 경계영역 처리를 보완해야 한다.
3. distance-map을 충분히 smooth하게 만들어야 한다: mean-filter 반복 적용.
4. 참고: https://kipl.tistory.com/66의 구현보다 단순하지만 성능은 떨어진다.

#define insideROI(x, y) ((x) >= 0 && (x) < w && (y) >= 0 && (y) < h)
#define BACKGND 0xFF
int watershed(float *distMap, int w, int h, int *rgnMap) {
    const int xend = w - 1, yend = h - 1;
    const int nx[] =  {-1, 0, 1, -1, 0, 1, -1, 0, 1};
    const int ny[] =  {-1, -1, -1, 0, 0, 0, 1, 1, 1};
    // processing image;
    std::vector<float>   wbuffer(w * h);
    std::vector<float* > procImg(h); 
    for (int y = 0; y < h; y++)
        procImg[y] = &wbuffer[y * w];
        
    // presmooth; 10 iterations;
    mean_filter(distMap, w, h, 10, procImg); 

    // 입력 영상에서 경계는 global maximum으로 설정;
    for (int y = 0; y < h; y++ ) procImg[y][0] = procImg[y][xend] = BACKGND;
    for (int x = 0; x < w; x++ ) procImg[0][x] = procImg[yend][x] = BACKGND;
    // find local minima for each pixel;
    int minCnt = 1;
    for (int y = 0, pos = 0; y < h; y++) {
        for (int x = 0; x < w; x++, pos++) {
            float minVal = procImg[y][x];
            if (minVal == BACKGND) {
                rgnMap[pos] = 1;	// background(white);
                continue;
            }
            int minIdx = 4; //current position;
            for (int k = 0; k < 9; k++) {
                int xx = x + nx[k], yy = y + ny[k];
                if (insideROI(xx, yy) && procImg[yy][xx] <= minVal) {
                    minVal = procImg[yy][xx];
                    minIdx = k;
                }
            }
            if (minIdx == 4) rgnMap[pos] = ++minCnt; // center(x, y) is a new local minimum;
            else             rgnMap[pos] = -minIdx;  // direction of local-min =  "-(dir)"
        }
    }
    // follow gradient downhill for each pixel;
    // reset rgnMap to have the id's of local minimum connected with current pixel;
    for (int y = 1, pos = y * w; y < yend; y++ ) {
        pos++;
        for (int x = 1; x < xend; x++, pos++ ) {
            int minpos = pos;
            while ( rgnMap[minpos] <= 0 ) { //it is not a local minimum.
                switch ( rgnMap[minpos] ) {
                case  0: minpos += -w - 1; break; // top-lef = downhill direction;
                case -1: minpos += -w;	   break; // top
                case -2: minpos += -w + 1; break; // top-right;
                case -3: minpos--;         break; // left;
                case -5: minpos++;         break; // right;
                case -6: minpos += w - 1;  break; // bot-left;
                case -7: minpos += w;      break; // bot;
                case -8: minpos += w + 1;  break; // bot-right;
                }
            }
            // speed-up: use a stack to record downhill path.
            //assign the id of a local min connected with current pixel;
            rgnMap[pos] = rgnMap[minpos]; 
        }
        pos++;
    }
    // rgnMap는 각 픽셀과 연결되는 국소점의 id을 알려주는 lookup table이다;
    return ( minCnt );
}

적용 예는 http://kipl.tistory.com/1         

std::vector<int> sort_pixels(std::vector<BYTE> &img,   std::vector<int> &sorted) ;

// img의 픽셀 위치를 gray값이 커지는 순서로 정렬하고, cumulative histogram chist를 반함함.

// chist은 sorted img에서 특정한 gray값을 갖는 픽셀이 포함된 위치 구간을 알 수 있게 한다.

// 8-way connectivity case only;
// last update: 2024-7-22;
#define MASK -2
#define WSHED 0
#define INIT -1
#define RANGE(i,j,h,w) (((i)>=0) && ((j)>=0) && ((i)<h) && ((j)<w))
// main procedure;
int watershed8(std::vector<BYTE> &img, int width, int height, 
               std::vector<int*> &lab)
{
    int cur_label = 0 ;
    std::queue<int> q;
    std::vector<int> dmap(img.size(), 0);  
    std::vector<int> sorted(img.size());
    std::vector<int> chist = sort_pixels(img, sorted); 
    int hmin = -1, hmax = 255;
    while (!chist[++hmin]);
    while (chist[--hmax]==chist[255]); hmax++; 
    // prepare label image;
    for (int y = 0; y < height; y++)
        for (int x = 0; x < width; x++) lab[y][x] = INIT;
    //
    int *plab = lab[0]; //1-dim representative of label image
    for (int h = hmin; h <= hmax; h++) {
        int *bin_left = &sorted[0] + (h==0 ? 0: chist[h-1]);
        int *bin_right = &sorted[0] + chist[h];
        for (int *curr = bin_left; curr < bin_right; curr++) {
            int pos = *curr;
            plab[pos]=MASK;      
            int x=pos % width;
            int y=pos / width;      
            for (int k=-1;k<=1;k++)
                for (int m=-1;m<=1;m++) 
                    if (RANGE(y+k, x+m, height, width)) 
                        if ((k!=0) || (m!=0))
                            if ((lab[y+k][x+m]>0) || (lab[y+k][x+m]==WSHED)) {
                                q.push(pos);
                                dmap[pos]=1;
                                k=1; break;
                            }
        }      

        int cur_dist=1;
        q.push(-1);
        while (1) {
            int pos = q.front(); q.pop();
            if (pos==-1) {
                if (q.empty()) break;
                else {
                    q.push(-1);
                    cur_dist++;
                    pos = q.front(); q.pop();
                }
            }
            int x=pos % width;
            int y=pos / width;
            for (int k=-1;k<=1;k++)
                for (int m=-1;m<=1;m++) {
                    if (RANGE(y+k, x+m, height, width)) 
                        if ((k!=0) || (m!=0)) {
                            int nn = (y+k)*width + (x+m);// nearest neightbor;
                            if ((dmap[nn]<cur_dist) && ((plab[nn]>0) || (plab[nn]==WSHED))) { 
                                if (plab[nn]>0) {
                                    if ((plab[pos]==MASK) || (plab[pos]==WSHED))
                                        plab[pos]=plab[nn];
                                    else if(plab[pos]!=plab[nn])
                                        plab[pos]=WSHED;
                                }
                                else if(plab[pos]==MASK) plab[pos]=WSHED;
                            } else if((plab[nn]==MASK) && (dmap[nn]==0)) {
                                dmap[nn]=cur_dist+1;
                                q.push(nn);
                            }
                        }
                }
        }  // end while(1)	    	  
        for (int *curr = bin_left; curr < bin_right; curr++) {
            int pos = *curr;
            dmap[pos]=0;	
            if (plab[pos]==MASK) {
                cur_label++;
                q.push(pos);
                plab[pos] = cur_label;
                while (!q.empty()) {
                    int pos1=q.front(); q.pop();
                    int x=pos1 % width;
                    int y=pos1 / width;
                    for (int k=-1;k<=1;k++)
                        for (int m=-1;m<=1;m++)
                            if (RANGE(y+k, x+m, height, width)) 
                                if ((k!=0) || (m!=0)) {
                                    int nn = (y+k)*width + (x+m);
                                    if (plab[nn]==MASK) {
                                        q.push(nn);
                                        plab[nn]=cur_label;
                                    }
                                }
                }
            }
        }
    }  //endfor h    	    	
    
    // convert the result(labeled rgns + boundary) to a suitable form;
    mark_boundary8(lab, width, height);
    return cur_label;
}

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

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

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

(4) final result; watershed 결과를  mask로 이용해서 이미지를 분할한 것이다. 겹친 cell이 분리되었다.

다른 예:

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

* 구현: https://kipl.tistory.com/66

https://kipl.tistory.com/299