평균이 $\lambda$인 Poisson 분포를 가지는 random number를 만드는 레시피는 [0,1]에서 추출된 균일한 random number 수열이 $u_1, u_2, u_3,...$로 주어질 때, $u_1  u_2 u_3....u_{k+1}\le  e^{-\lambda}$을 처음 만족하는 $k$를 찾으면 됨이 Knuth에 의해서 알려졌다.

int PoissonValue(double lambda) { // lambda = mean value
    double L = exp(-lambda);
    int k = 0;
    double p = 1;
    do {
        k++;
        // generate uniform random number u in [0,1] and let p ← p × u.
        p *= double(rand()) / RAND_MAX;
    } while (p > L);
    return (k - 1);
}

이미지에서 각각의 픽셀 값은 영상신호를 받아들이는 ccd 센서에서 노이즈가 포함된 전기신호(광자 수)의 평균값이 구현된 것으로 볼 수 있다. 따라서 영상에 Poisson 노이즈를 추가하는 과정은 역으로 주어진 픽셀 값을 평균으로 가지는 Poisson 분포의 램덤 변수를 픽셀값으로 대체하면 된다.(물론 이미지에 무관하게 노이즈를 추가할 수도 있다. 예를 들면 모든 픽셀에 대해서 같은 평균값을 갖는 포아송 노이즈를 주는 경우다)

void AddPoissonNoise(CRaster& raster, CRaster& noised) {
    if (raster.IsEmpty() || raster.GetBPP() != 8) return;
    CSize sz = raster.GetSize();
    noised = raster;
    srand(unsigned(time(0)));
    for (int y = 0; y < sz.cy; y++) {
        BYTE *p = (BYTE *)raster.GetLinePtr(y);
        BYTE *q = (BYTE *)noised.GetLinePtr(y);
        for (int x = 0; x < sz.cx; x++) {
            int a = PoissonValue(*p++);
            *q++ = a > 0xFF ? 0xFF: a;
        }
    }
}

아래 예제는 cos(t)로 표현되는 신호에 Poisson noise를 주는 경우를 보여준다. Poisson noise 함수는 신호값을 평균으로 주는 event의 횟수를 찾아주므로 음수의 신호는 바로 다룰 수 없다. 따라서 원신호에 +1 만큼 더해 준 신호를 기준으로 한다. 그런데 cos(t) + 1의 범위가 [0,2]이므로 이산분포인 Poisson 분포를 적용하기 위해서는 적당한 스케일링이 필요하다. 예를 들면 cos(t) + 1이 신호의 절대값을 측정한 것이 아닌 상대적인 세기를 나타내고, 실제 신호의 최대값은 50으로 주어졌다면 Poisson noise을 만들 때 평균값이 주어진 신호 함수의 50배인 경우에 해당하는 값을 입력으로 넣은 후 출력값을 다시 50으로 나누어주면 된다.

void PoissonNoiseEx1D() {
    FILE *fp = fopen("poisson.txt", "w");
    for (double t = 0; t < 10; t += 0.01) {
        double y = cos(t) + 1.;
        fprintf(fp, "%f %f %f\n", t, y, (double)PoissonValue(y * 50) / 50.);
    }
    fclose(fp);
}

Poisson 잡음은 random noise처럼 단순히 더해지는 잡음과 달리 신호값이 클 때 잡음도 더 커짐을 알 수 있다. 그러나 SNR은 신호값이 클수록 더 좋아진다.

728x90

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

Fowler Angle  (0) 2021.04.05
Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Posted by helloktk
,

blue cell: obstacles, black cell = start( or end)

$$\text{distance measure}:~\text{Manhattan distance}=|x_2 - x_1 | + |y_2 - y_1|$$

Grassfire 알고리즘을 써서 최단경로를 찾는 문제는 DFS로도 구현할 수 있으나 비효율적이고(왜 그런지는 쉽게 알 수 있다) BFS로 구현하는 것이 더 적합하다. 

void grassfire(CPoint q, int **map, int w, int h, int **dist) {
    // left-top-right-bottom;
    const int dx[] = {-1,  0, 1, 0};
    const int dy[] = { 0, -1, 0, 1};
    for (int y = 0; y < h; y++) 
        for (int x = 0; x < w; x++) 
            dist[y][x] = INF;  //unvisited cells;

    std::queue<CPoint> Q;
    dist[q.y][q.x] = 0;     //start( or end) position: distance = 0;
    Q.push(q);
    while (!Q.empty()) {
        CPoint p = Q.front(); Q.pop();
        int distance = dist[p.y][p.x];
        // 4-way search;
        for (int i = 0; i < 4; i++) {
            CPoint q = CPoint(p.x + dx[i], p.y + dy[i]);
            if (q.x < 0|| q.y < 0|| q.x >= w|| q.y >= h) continue;
            if (map[q.y][q.x] == 0 && dist[q.y][q.x] == INF) {
                dist[q.y][q.x] = distance + 1;
                Q.push(q);
            }
        }
    }
};
// back tracking;
CPoint back_track(CPoint p, int **dist, int w, in h) {
    // left-top-right-bottom;
    const int dx[] = {-1,  0, 1, 0};
    const int dy[] = { 0, -1, 0, 1};
    int depth = dist[p.y][p.x]; 
    if (--depth < 0) return p;
    for (int i = 0; i < 4; i++) {
        CPoint q = CPoint(p.x + dx[i], p.y + dy[i]);
        if (q.x < 0 || q.y < 0 || q.x >= w || q.y >= h) continue; // out of ROI;
        else if (dist[q.y][q.x] == depth)
            return q;
    }
    return p; // never hit;
}

728x90

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

Brute-Force Euclidean Distance Transform  (0) 2021.03.14
Poisson Noise  (0) 2021.03.06
Image Sharpness  (0) 2021.02.25
Selection Sort  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Posted by helloktk
,

segmentation result

  1. gradient 대신 negate시킨 distance map을 사용한다.
  2. 경계영역 처리를 보완해야 한다.
  3. distance-map을 충분히 smooth하게 만들어야 한다: mean-filter 반복 적용.
  4. 참고: 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
,
#define PIX_SORT(a,b) { if ((a) > (b)) PIX_SWAP((a), (b));} 
#define PIX_SWAP(a,b) { BYTE tmp = (a); (a) = (b); (b) = tmp;} 
BYTE opt_med9(BYTE p[9]) { 
    PIX_SORT(p[1], p[2]); PIX_SORT(p[4], p[5]); PIX_SORT(p[7], p[8]); 
    PIX_SORT(p[0], p[1]); PIX_SORT(p[3], p[4]); PIX_SORT(p[6], p[7]); 
    PIX_SORT(p[1], p[2]); PIX_SORT(p[4], p[5]); PIX_SORT(p[7], p[8]); 
    PIX_SORT(p[0], p[3]); PIX_SORT(p[5], p[8]); PIX_SORT(p[4], p[7]); 
    PIX_SORT(p[3], p[6]); PIX_SORT(p[1], p[4]); PIX_SORT(p[2], p[5]); 
    PIX_SORT(p[4], p[7]); PIX_SORT(p[4], p[2]); PIX_SORT(p[6], p[4]); 
    PIX_SORT(p[4], p[2]); return(p[4]); 
}
double ImageSharpness(BYTE *img, int w, int h) {
    const int npixs = w * h;
    const int xend = w - 1, yend = h - 1;
    const int nn[] = {-w - 1, -w, -w + 1, -1, 0, 1, w - 1, w , w + 1};
    const int sobelX[] = {-1, 0, 1, -2, 0, 2, -1, 0, 1};
    const int sobelY[] = {-1, -2, -1, 0, 0, 0, 1, 2, 1};
    //median filter;
    BYTE arr[9];
    std::vector<BYTE> medimg(npixs, 0);
    for (int y = 1, pos = y * w; y < yend; y++) {
        pos++;
        for (int x = 1; x < xend; x++, pos++) {
            for (int k = 0; k < 9; k++) arr[k] = img[pos + nn[k]];
            medimg[pos] = opt_med9(arr);
        }
        pos++;
    }
    // Tenenbaum gradient;
    double sharpness = 0;
    for (int y = 1, pos = y * w; y < yend; y++) {
        pos++;
        for (int x = 1; x < xend; x++, pos++) {
            double gx = 0, gy = 0;
            for (int k = 0; k < 9; k++) {
                int v = medimg[pos + nn[k]];
                gx += sobelX[k] * v;
                gy += sobelY[k] * v;
            }
            sharpness += gx * gx + gy * gy;
        }
        pos++;
    }
    return sharpness;
}
728x90

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

Poisson Noise  (0) 2021.03.06
Grassfire Algorithm  (0) 2021.03.05
Selection Sort  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Insertion Sort  (0) 2021.02.24
Posted by helloktk
,

void selection_sort(BYTE data[], int n) {
    for (int i = 0; i < n - 1; i++) {
        int jmin = i;
        for (int j = i + 1; j < n; j++)
            if (data[j] < data[jmin]) jmin = j;
        if (jmin != i) SWAP(data[jmin], data[i]);
    }
}
728x90

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

Grassfire Algorithm  (0) 2021.03.05
Image Sharpness  (0) 2021.02.25
Bubble Sort  (0) 2021.02.24
Insertion Sort  (0) 2021.02.24
Optimized Median Search  (0) 2021.02.24
Posted by helloktk
,