각 전경 픽셀에 대해서 모든 배경 픽셀까지 거리를 구한 다음 최솟값을 할당하면 된다. 픽셀 수가 n = w * h 개면 time complexity는 O(n^2)이다 (이미지 폭이나 높이의 4승이므로 연산량이 상당하다.) linear time Euclidean distance transform은 kipl.tistory.com/268.

ListPlot3D[ Reverse@ImageData[pic], PlotRange -> Full]

mathematica를 이용한 3D plot

int brute_force_edt(double *img, int w, int h) {
    int n = w * h;
    std::vector<int> list(n);
    int fg_cnt = 0;
    int bg_ptr = n - 1;               //배경 픽셀 위치는 끝에서 역으로 채움;
    //전경과 배경 픽셀 위치를 분리;
    for (int i = 0; i < n; ++i)
        if (img[i])                  // foreground;
            list[fg_cnt++] = i;
        else                         // background;
            list[bg_ptr--] = i;

    for (int i = 0; i < fg_cnt; ++i) {     // 전경 픽셀;
        int xi = list[i] % w, yi = list[i] / w;
        double d2min = DBL_MAX;
        for (int j = n; j-- > fg_cnt; ) { // 배경 픽셀까지 거리를 계산해서 최소값을 할당;
                                          // 배경이 list에 역순으로 저장됨으로 역방향 서치;
            int xj = list[j] % w, yj = list[j] / w;
            double dx  = xi - xj, dy = yi - yj;  //
            double dst = dx * dx + dy * dy;
            if (dst == 1) { 
            	d2min = 1; break;
            }
            if (d2min > dst) d2min = dst;
        }
        img[list[i]] = d2min;
    }
    return fg_cnt;
}
 
 
 
728x90

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

Interpolation Kernels  (0) 2021.05.05
Fowler Angle  (0) 2021.04.05
Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Posted by helloktk
,

segmentation result

1. 경계영역 처리를 보완해야 한다.

2. distance-map을 충분히 smooth하게 만들어야 한다: mean-filter 반복 적용.

3. 참고: imagej.net/Distance_Transform_Watershed의 구현보다 단순하지만 성능은 떨어진다.

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

 

 
 
728x90

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

Anisotropic Diffusion Filter  (0) 2024.02.23
Edge Preserving Smoothing  (0) 2024.02.14
Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (3) 2021.02.15
Fixed-Point Bicubic Interpolation  (1) 2021.01.19
Posted by helloktk
,

Distance Transform은 이미지 상의 각 픽셀에서 가장 가까이 놓인 물체(전경)의 경계까지 거리를 구하는 연산이다. 이진 영상의 경우 각 배경 픽셀에서 전경까지 거리를 구하는 과정이 된다. 이렇게 얻어진 거리 정보를 이용해서 영상 속 전경 물체의 skeleton을 구하거나 (OCR), 로봇의 팔 등을 자동 제어하는 데 있어서 움직이는 경로에 놓인 방해물의 위치를 예측할 수 있으며, 또한 물체까지의 경로 찾기(path finding) 등에 이용이 된다. distance transform은 사용되는 거리 측도(measure)에 따라 그 특성이 달라진다. 우선 가장 단순한 Euclidean 측도를 이용하도록 하자.

 

distance transform은 전경의 경계선을 모두 찾아서 배경의 각 픽셀에서 거리의 최솟값을 할당하면 된다. 그러나 경계가 많은 픽셀로 구성이 되는 이미지에서는 연산 비용이 매우 커진다. 여기서는 이미지의 크기에 선형으로 비례하는 알고리즘에 대해서 살펴보도록 하겠다. 픽셀 수에 선형이기 위해서는 알고리즘이 raster 스캔 방식을 따라야 한다.

 

Euclidean 측도를 사용하는 distance transform은 한 점 $(x, y)$에서 전경 픽셀 $(i, j)$까지 거리의 제곱이 $z = (x-i)^2 + (y-j)^2$로 표현된다. 각 전경 픽셀 꼭짓점으로 하는 paraboloid를 생각하면 distance transform은 주어진 배경 픽셀에서 볼 때 가장 밑에 위치한 paraboloid를 찾는 문제로 환원된다. 그러나 이차원 영역에서 이를 찾기는 쉽지 않은 작업이므로 이를 수평과 수직 방향으로 분리해서 적용할 수 알고리즘을 먼저 만들어야 한다. 우선 $x$-방향으로 스캔하면서 $x$-축으로 사영된 거리를 저장하면서 다시 $y$-축 방향으로 스캔하면서 저장된 $x$-방향의 거리의 제곱의 정보를 이용하면 된다. 이것은  $x$-방향으로 스캔할 때 거리의 제곱 값을 이미지의 그레이 값으로 저장하면 보다 쉽게 해결이 된다. 그러면 $y$-방향으로 스캔할 때는 꼭짓점이 이전에 $x$-방향 스캔 시 저장된 거리의 제곱만큼 offset 된 이차곡선들 중 맨 아래에 놓인 것을 찾으면 된다.  2차원 distance transform은 1차원의 distance transform을 행방향 먼저 수행하고, 그 결과를 다시 열방향으로 시행하면 얻을 수 있다.

 

1차원의 distance transform은 일직선에 놓인 각 전경 픽셀을 꼭짓점으로 하는 이차 곡선들 중에서 가장 아래에 놓인 곡선을 선택하여 그 값을 취하면 된다. 배경 픽셀은 고려대상이 아니므로 꼭지점이 배경 픽셀이면 이차곡선에는 매우 큰 (수학적으로는 무한대) offset을 준다:

$$\text {1-dim distance transform}(x)=\underset {p\in\text {pixels}}{\text {min}}  \left( (x-p)^2+ F(p) \right)$$

초기 offset $F(p)$ 값은 전경 픽셀은 0, 배경 픽셀은 +INFINITY로 설정한다. 전경의 경계에서 형성이 된 2차 곡선 가장 아래에 놓이므로 그것에 의해서 거리가 정해진다.

일반적인 F(q)의 값이 주어지는 경우 각 픽셀에서 전경까지 거리를 구하는 것은 그 픽셀에서 볼 때 가장 아래에 있는 이차곡선이 어느 것이고, 어는 영역까지 유효한지를 찾는 문제로 귀결이 된다. 픽셀 위치 $p$, $q$에 꼭짓점이 있는 두 이차곡선의 교점이 

$$s=\frac { (F(p) + p^2)-(F(q) + q^2) }{ 2(p-q)}$$

로 표현되므로  주어진 2차 곡선이 가장 낮게 되는 영역은 쉽게 찾을 수 있다. 아래 그림에서 $p<x \le s$ 에서는 꼭짓점이 $p$인 이차곡선이, $s <x <q$ 영역에서는 꼭짓점이 $q$인 이차곡선이 거리를 준다.

 

그리고 다음 그림은 일차원 이미지에 대한 distance transform의 결과이다: 1차원 이진영상에 대해 적용한 것으로 붉은 선은 이미지(0:전경, 255:배경)이고, 녹색 선은 거리 값을 나타낸다.

네이버 블로그 이전

이 일차원 distance transform을 각각의 행에 대해서 적용하여 거리를 구한 다음 그 결과를 영상의 픽셀 값으로 저장한다. 다시 각각의 열에 대해서 적용하는데, 이때 주어진 픽셀 $p$에서 영상 값이 $F(p)$의 역할을 하게 된다. 아래의 그림은 글씨에 대해서 skeleton을 구하기 위해서 distance transform을 적용하였다(글씨의 skeleton을 얻기 위해서 negate 영상에 대해서 적용했다: 흰색=배경, 검은색=전경)

** 참고: en.wikipedia.org/wiki/Distance_transform

** cs.brown.edu/people/pfelzens/dt/에서 관련 논문과 원 코드를 찾을 수 있다;

#define SQR(a)	((a) * (a))
/* distanceTransform of 1d function using the Euclidean distance */
struct para_envelop {
    int x;           // parabola의 꼭지점 좌표;
    double start;    // 현재 parabola가 가장 아래 놓인 구간의 시작점;
};
double para_intersection(int q, int p, double *data) {
    return ((data[q] + SQR(q)) - (data[p] + SQR(p))) / (q - p) / 2;
}
// dblimg 이미지의 초기화: 배경 = +INF, 전경 = 0 으로 설정한다.
void DistanceTransform_1D(double *data, int n, double* dist) {
    std::vector<para_envelop> env(n + 1);
    // 가장 낮게 형성이 되는 parabola의 꼭지점과 시작점(유효한 영역 시작점)을 env에 등록한다.
    env[0].x = 0; env[0].start = -FLT_MAX; //첫번째 parabola를 모든 영역에서 가장 낮다고 설정.
    env[1].start = +FLT_MAX;
    int k = 0;
    for (int x = 1; x < n; ++x) {
        // 현 parabola와 이전 오른쪽 가장 낮은 env와 교차점을 찾음;
        double s = para_intersection(x, env[k].x, data);
        // 교차점이 이전 env 꼭지점보다 왼쪽에 있으면 현재 parabola가 더 낮으므로 더 이전의 
        // env와 교차점 검사해서 꼭지점보다 오른쪽에 있을 때까지 계속 인덱스를 내린다.
        while (s <= env[k].start) s = para_intersection(x, env[--k].x, data);
        // env[k]까지는 가장 낮은 parabola꼭지점과 그 영역의 시작점이 기록되어 있음.
        // x > s 이후는 현 parabola가 가장 낮으므로 새로이 env등록을 함;
        env[++k].x = x ;    
        env[k].start = s ;
        env[k + 1].start = +FLT_MAX;
    }
    // 등록된 env 정보를 이용해서 거리정보를 업데이트 함;
    // {env[k]}는 0<=x<n사이에서 가장 낮은 parabola의 꼭지점과 그것이 성립하는 구간의 시작점
    // 정보를 담고 있음. 긱 픽셀에서 가장 낮은 env.start가 그 픽셀보다 왼쪽에 있는 env를 찾아
    // env의 꼭지점까지 거리를 구하면 된다.
    k = 0;
    for (int x = 0; x < n; ++x) {
        // x가 env.x 오른쪽, 다음 env.start의 왼쪽에 놓여야 한다;
        while (env[k + 1].start < x) k++;
        // x에서 가장 가까이 있는 전경 픽셀 위치가 env[k].x 이므로 
        // Euclidean 거리 정보를 업데이트; 
        dist[x] = SQR(x - env[k].x) + data[env[k].x];
    }
};
void DistanceTransform_2D(double *img, int width, int height) {
    int sz = max(width, height);
    std::vector<double> dblimg(sz), dist(sz);
    // transform along rows
    for (int y = 0; y < height; y++) {
        double *ptr = &img[y * width];
        for (int x = 0; x < width; x++)
            dblimg[x] = *ptr++ ? 0: INF;    //전경(white)=0, 배경(black)=INF;
        distanceTransform(&dblimg[0], width, &dist[0]);
        memcpy(&img[y * width], &dist[0], width * sizeof(double));
    };
    // transform along columns
    for (int x = 0; x < width; x++) {
        double *ptr = &img[x];  // copy column;
        for (int y = 0; y < height; y++, ptr += width)
            dblimg[y] = *ptr;
        distanceTransform(&dblimg[0], height, &dist[0]);
        ptr = &img[x];
        for (int y = 0; y < height; y++, ptr += width)
            *ptr = dist[y];
    }
};
 
728x90

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

Contrast Limited Adaptive Histogram Equalization (CLAHE)  (3) 2021.02.15
Fixed-Point Bicubic Interpolation  (1) 2021.01.19
FFT 알고리즘의 재귀적 구현  (0) 2021.01.14
Edge-Preserving Smoothing  (0) 2021.01.12
Octree Quantization  (0) 2021.01.12
Posted by helloktk
,

Chamfer Match

Image Recognition 2008. 8. 1. 13:11

주어진 영상이 어떤 특정한 물체를 담고 있는지를 알아보려면, 그 물체의 template을 만들어서 영상 위에 겹친 후 픽셀 단위로 순차적으로 이동시키면서 일치하는가를 체크하면 된다. 일치의 여부는 보통의 경우 템플레이트와 영상의 상관계수를 측정하면 된다. 대상 영상이 이진 영상인 경우에는 distance map을 이용하면 좀 더 쉽게 매칭 문제를 해결할 수 있다. 에지 영상을 사용해서 만든 distance map은 이미지 영역의 한 지점에서 가장 가까이에 있는 에지 위치까지의 거리를 알려준다. 매칭하고자 template의 이진 영상을 담고 있는 이미지를 distance map에서 픽셀 단위로 차례로 옮기면서 template 위치에서 distance map이 알려주는 에지까지의 거리 합을 기록하면 총거리의 합이 가장 작은 지점에서 template와 가장 잘 맞는 형상을 찾을 수 있을 것이다.

이 메칭 기법(Chamfer Match)은 모델과 이미지가 서로 회전되어 있지 않아야 하고, 스케일 변형도 없는 경우에 좋은 결과를 얻는다. 스케일 변화가 있는 경우에는 스케일 공간에서 순차적으로 찾는 방법을 쓰면 된다. 모델 영상에 변화가 있는 경우에는 민감하게 반응하므로, 각각의 변화된 모델 영상을 따로 준비하여서 매칭 과정을 수행해야 한다. distance map만 구해지면 매칭 과정이 매우 단순하기 때문에 매칭 비용도 매우 저렴하다. 그리고 거리 값을 적당한 범위에 truncate 하는 것이 보다 robust 한 결과를 준다.

* procedure:

  1. 입력영상의 에지영상을 구함.
  2. 입력영상의 에지영상의 distance map을 구함.
  3. 매칭과정을 실행.

사용자 삽입 이미지

int GetAddressTable(BYTE *subImg, int w, int h, int stride/*=원본이미지 폭*/, int addr[]);

더보기
//계산을 보다 빠르게 하기 위해서 subimg에서 template 위치를(distimg에 놓인 경우) 미리 계산함.
int GetAddressTable(BYTE *subImg, int w, int h, int stride/*=원본이미지 폭*/, int addr[]) {
    for(int y=0, k=0; y<h; y++) {
        for(int x=0; x<w; x++){
            if(subImg[y*w+x])
                addr[k++]= y*stride + x ; //template가 dist이미지에서움직일 때 픽셀위치(상대적);
        }
    }
    return (k);    //# of template pixels<= w*h;
}
Matching code;
void ChamferMatch(float* distImg, int w1, int h1,   //distance-map;
                  BYTE* subImg, int w2, int h2,     //template-map;
                  float *match/*[w1 * h1]*/) {      //match_score-map(filled with FLT_MAX)
    int umax = w1 - w2;
    int vmax = h1 - h2;
    int uc = w2 / 2;      //center_x of template-map;
    int vc = h2 / 2;      //center_y of template-map;
    int *addr = new int [w2 * h2]; //max;
    int n = GetAddressTable(subImg, w2, h2, w1, addr) ;
    for (int v = 0; v <= vmax; v += 2){           //to speed-up;
        for (int u = 0; u <= umax; u += 2){       //to speed-up;
            double score=0;
#if (0)
            for (int v1 = v, v2 = 0; v1 < h1 && v2 < h2; v1++, v2++){
                int i1 = v1 * w1;
                int i2 = v2 * w2;
                for (int u1 = u, u2 = 0; u1 < w1 && u2 < w2; u1++, u2++) {
                    if (subImg[i2 + u2])
                        score += distImg[i1 + u1]; // truncate [0, th_dist];
                }
            }           
#else 
            int offset = v * w1 + u; //start address;
            int i = n;
            while(i--) score += distImg[addr[i] + offset];
#endif
            //at the center of subimg;
            match[(uc + u) + (vc + v) * w1] = score; 
        }   
    }   
    delete[] addr ;
};


입력영상의 에지영상(모델이미지를 포함, 에지를 1픽셀로 조절하지 않았음)

 

사용자 삽입 이미지

모델영상:

사용자 삽입 이미지


 distance map에서 찾은 위치에 template를 오버랩함;

사용자 삽입 이미지

 

사용자 삽입 이미지
사용자 삽입 이미지
사용자 삽입 이미지

728x90

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

Rolling Ball Transformation  (1) 2008.08.09
Mean Shift Filter  (5) 2008.08.06
Retinex Algorithm  (2) 2008.07.26
RANSAC: Circle Fit  (0) 2008.07.21
KMeans Algorithm  (0) 2008.07.19
Posted by helloktk
,

(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 에서 옮긴 자료.
*/

728x90

'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
Posted by helloktk
,