타원도 기본적으로 원 그리기와 동일한 알고리즘을 이용하여 그릴 수 있다. 하지만 원과 달리 좀 더 고려할 부분이 있다. 전체 영역을 $45^\circ$씩 8개 영역으로 나눈 후 한 영역에서만 계산하고 나머지는 대칭성을 이용하는 방법을 타원에는 바로 적용이 어렵다. 원에서 $45^\circ$ 기준선은 접선의 기울기가 $\frac{dy}{dx} = -1$인 점과 원점을 잇는 선분인데 타원에서는 이 관계가 성립하지 않는다. 타원 방정식 $\frac{x^2 }{ a^2} +\frac{ y^2}{b^2} = 1$ 에서 접선의 기울기가 $-1$인 지점은

$$x = \frac{a^2 }{ \sqrt {a^2  + b^2 }} , \quad y = \frac{ b^2 }{ \sqrt {a^2  + b^2 }}$$

로 주어진다. 이 지점까지는 $x$를 독립변수로 잡아야 하고, 나머지 $y = 0$까지는(1 사분면에서) $y$를 독립변수로 잡아야 한다. 나머지 분면은 대칭성을 이용하면 바로 그릴 수 있다.

사용자 삽입 이미지

픽셀 $(x_k , y_k)$에  점이 찍히면 다음은 어느 픽셀에 점을 찍어야 할까? $0\ge \text{기울기}\ge-1$이므로 오른쪽 옆이나 또는 오른쪽 아래 픽셀이 될 것이다. 이를 판별하는 식을 만들기 위해서 다음 식을 고려하자(나눗셈을 피하기 위해서 타원식에 $a^2b^2$을 곱했다): 

$$F(x, y) = b^2 x^2 + a^2 y^2 - a^2 b^2$$

$(x_k, y_k )$가 타원 위의 점이면 $F=0$을 만족한다. 다음 점이 찍힐 후보 위치가 $(x_k+1, y_k)$ 또는 $(x_k+1, y_k-1)$이므로 중간 위치(Midpoint)에서 $F$의 부호를 기준으로 판별하도록 하자:

$$d_k = F(x_k+1, y_k-1/2) \\ d_k < 0  \quad \text {이면} \quad (x_k+1, y_k),\\   d_k > 0\quad \text {이면} \quad (x_k+1, y_k-1)$$

이후에 판별식의 업데이트는 

$$d_k <0~인 ~경우:\quad           d_{k+1} = F(x_k+2, y_k-1/2) = d_{k} + b^2 ( 2 x_{k+1} +1)$$

$$d_k \ge 0 ~인 ~경우: \quad  d_{k+1} = F(x_k + 2, y_k -3/2) = d_{k} + b^2 (2x_{k+1} + 1)- 2a^2  y_{k+1}$$

로 주어진다. 시작 픽셀에서 판별식의 값은

$$ x_0 = 0, y_0 = b \quad \Rightarrow\quad d_0 = F(1, b-1/2) = (4 b^2 + a^2 (1 - 4 b)) / 4$$

다음 문제는 $x$가 독립변수로써 얼마나 사용되고, 어느 지점에서 $y$가 독립변수가 되어야 하는가를 알아내는 것이다. 기울기 $dy /dx$ 가 0에서 -1 사이인 지점까지 $x$가 독립변수이므로, 우선 타원식을 미분하면

$$2b^2 x + 2a^2 y \frac{dy}{dx} = 0 \\ \Rightarrow b^2 x = -a^2 y \frac{dy}{dx} = a^2 y \left| \frac{dy}{dx}\right|  \le a^2 y, \quad -1 \le \frac{dy}{dx} \le 0 \\ \therefore  \quad b^2 x \le a^2 y$$

이 조건이 만족하는 동안 $x$를 독립변수로 사용하여 그리기를 한다.

 

$y$가 독립변수인 구간으로 넘어가자. 이 경우는 $x_0=a, ~y_0=0$에서 출발하면 위 영역에서 $x_k$와 $y_k$의 역할을 바꾼 방식이 된다. $(x_k, y_k)$에 점이 찍히면 다음 후보는 $(x_k, y_k+1)$ 또는 $(x_k-1, y_k+1)$이다. 중간점이 $(x_k-1/2, y_k+1)$이므로 판별식은

$$ d_k = F(x_k -1/2, y_k +1)$$

로 주어진다. 이 판별식의 업데이트는 

$$d_k <0~인 ~경우:\quad           d_{k+1} = F(x_k - 1/2, y_k+2) = d_{k} + a^2 ( 2 y_{k+1} +1)$$

$$d_k \ge 0 ~인 ~경우: \quad  d_{k+1} = F(x_k  -3 / 2, y_k +2) = d_{k} +a^2 (2x_{k+1} + 1)  - 2 b^2  x_{k+1}$$

로 주어진다. 그리고 시작 픽셀에서 판별식의 값은

$$ x_0 = a, y_0 = 0 \quad \Rightarrow\quad d_0 = F(a-1/2,0) = (4 a^2 + b^2 (1 - 4 a)) / 4$$

이다. 전체적인 코드는

// setPixel(CDC *pDC, int x, int y)에 (x, y), (-x, y), (x,-y), (-x,-y)를 
// 동시에 그리도록 구현이 되어야 한다. 그리고 중심이 원점이 아닌 경우에는 
// 평행이동 시키면 된다:
    int x = 0, y = b ;
    int aa = a * a ;
    int bb = b * b ;
    int d = (4 * bb + aa * (1 - 4 * b)) / 4 ;

    setPixel (pDC, x, y) ;     
    while ( bb * x < aa * y) {
        ++x ;
        if (d < 0) d +=  bb * (2 * x + 1);    // (x+1, y) 선택;
        else {                                // (x+1, y-1) 선택;
            --y ;
            d += bb * (2 * x + 1) - 2 * aa * y ;
        }
        setPixel(pDC, x, y);

    };
    //y-독립변수 구간;
    x = a, y = 0;
    d = (4 * aa + bb * (1 - 4 * a)) / 4;
    setPixel(pDC, x, y) ;
    while (bb * x > aa * y) {
        ++y ;
        if (d < 0) d += aa * (2 * y + 1);
        else {
            --x ;
            d += - 2 * bb * x + aa * (2 * y + 1);   
        }
        setPixel(pDC, x, y);
    }

사용자 삽입 이미지

물론 실수 연산을 하면, 아래의 코드로 타원을 그릴 수 있다. 접선의 기울기의 크기 $|dy/dx| = 1$ 인 점을 기준으로 하여 그 값이 1 보다 작은 지점에서는 $x$를 독립변수로, 큰 지점에서는 $y$를 독립변수로 잡으면,

          x-구간    ===>  b^2 * x < a^2 * y;
          y-구간    ===>  b^2 * x > a^2 * y;

이므로, 아래처럼 쓸 수가 있다:

    int aa = a * a, bb = b * b;
    x = 0; y = b;
    while (bb * x <= aa * y) {
        y = (int)((double) b / a * sqrt(aa - x * x) + 0.5);
        setPixel(pDC, x, y);
        ++x;
    }
    y = 0;
    while (bb * x > aa * y) {
        x = (int)((double) a / b * sqrt(bb - y * y) + 0.5);
        setPixel(pDC, x, y) ;
        ++y;
    }

last update: 2021.04.08;

728x90

'Computational Geometry' 카테고리의 다른 글

삼각형 채우기  (0) 2009.12.19
Triangulation을 이용한 Imge Effect  (0) 2009.12.18
Circle Drawing Algorithm  (1) 2008.06.03
Wu's Line Algorithm  (1) 2008.06.02
Brute Force Triangulation  (1) 2008.05.28
Posted by helloktk
,

원을 그리기 위해서는 원의 방정식을 이용하면 된다. 반지름이 $r$인 경우에 원점이 중심인 원은 $x^2 + y^2 = r^2$이므로 (중심이 원점이 아닌 경우는 단순 평행 이동만 하면 되기 때문에 따로 취급하지 않는다), $x$를 독립변수로 잡으면 $y= \sqrt {r^2 - y^2}$으로 주어져서 1/2 사분면을 그릴 수 있다. 그리고 부호를 반대로 하면 3/4 분면에 그려진다. 그러나 $x$가 $r$에 근접하는 경우 $y$의 변화율이 급격해져서 $y$값이 연속적으로 그려지지 않는다. 이 문제는 전체 4 사분면을 45도 단위로 8개 영역으로 쪼개서 $x$와 $y$의 역할을 적당히 바꾸면서 그리면 이러한 문제를 피할 수 있다. 그러나, 여전히 실수 연산을 수행해야 하는 문제가 남아 있다. 이를 피하는 방법을 고안하도록 하자. 힌트는 Bresenham의 라인 알고리즘에서 얻을 수 있다. 우선, 영역을 $x = 0$에서 $x = y$인 90-도에서 45-도 사이의 구간만 생각하자. 이 구간에서 $y$의 변화율의 크기는 항상 0에서 1 사이 이므로, 8-방향으로 연결성을 갖도록 그리려면 현재의 점이 $(x, y)$로 주어지는 경우

          (x, y) ---> (x+1, y)     변화율 = 0;
                      (x+1, y-1)   변화율 = 1      (절대값)

의 두 가지 경우만 있다.  그럼 이 두 점 중 어느 것을 선택을 해야 하는가? 이것을 판별하기 위한 식은 이 두 점의 중심인 $(x, y - 1/2)$이 원의 내부에 있는가 또는 외부에 있는가를 보고서 결정하면 된다.

          (x+1, y-1/2) 이  원의 내부점인 경우---->(x+1, y)  선택
                           원의 외부점인 경우---->(x+1, y-1) 선택;

하면 각각의 상황에서 보다 정확한 원에 가까운 점을 선택한 것이 된다. 다음 점에서도 동일한 과정을 거치면 된다.

사용자 삽입 이미지

그러면 $(x+1, y-1/2)$가 원의 내부점인지 외부점인지를 항상 계산해야 하는가? 매 스텝마다 이것을 계산을 할 필요는 없다. 왜냐면 처음 $(0, r-1/2)$에서 계산을 하고

          d = x^2+y^2-r^2 at (0, r-1/2) ;
             = (5-4*r)/4  ; -> 이런 식으로 표현을 하는 것은 정수 연산만으로 
                               하기 위해서 계산 우선 순위를 조절한 것임.

이것과, $F(x, y)=x^2-y^2-r^2$ 라고 할 때,

(1),  $(x+1, y)$를 선택할 때 새로운 $d$값은 $(x+2, y-1/2)$에서 계산하므로

         d_old = F(x+1, y-1/2) = (x+1)^2+(y-1/2)^2-r^2 ;
         d_new = F(x+2, y-1/2) = d_old + 2*(x+1) + 1;

(2),  $(x+1, y-1)$를 선택할 때, 새로운 $d$값은 $(x+2, y-1-1/2)$에서 계산하므로

        d_old = F(x+1, y-1/2) = (x+1)^2+(y-1/2)^2-r^2 ;
        d_new = F(x+2, y-3/2) = d_old + 2 * ((x+1) -(y-1)) + 1;

인 관계를 이용하면 쉽게 정수 연산만으로 $d$값을 업데이트할 수 있다. 따라서 전체적인 코드는

   int x = 0, y = r ;
   int d = (5 - 4 * r) / 4;
   setPixel (x, y) ; 
   while ( x < y ) {
        ++x ;
        if (d < 0)                            // (x+1, y) 선택;
            d +=  2*x + 1;                    // ++x와 연관해서 잘 살펴야 한다.
        else {                                // (x+1, y-1) 선택;
            --y ;
            d += 2 * (x - y) + 1;             // ++x,--y;
        }
        setPixel(x, y);
   };

8 방향 대칭성으로 인해 점 $(x, y)$를 그리면 아래의 점도 그려야 한다 (처음 점 $(0, r)$은 4방향만, 그리고 $x=y$인 경우도 마찬가지나 따로 분리하기보다는 그냥 그리는 것이 더 낫다).

   {(x,-y), (-x,y), (-x,-y), (y,x), (y,-x), (-y,x),(-y,-x)}

단, 원의 중심이 원점이 아닐 경우는 전체적으로 평행이동을 시키면 된다.

사용자 삽입 이미지

728x90

'Computational Geometry' 카테고리의 다른 글

Triangulation을 이용한 Imge Effect  (0) 2009.12.18
Ellipse Drawing Algorithm  (0) 2008.06.04
Wu's Line Algorithm  (1) 2008.06.02
Brute Force Triangulation  (1) 2008.05.28
Polygon Triangulation (II)  (0) 2008.05.26
Posted by helloktk
,

이미지 상의 한 점은 픽셀로 표현이 된다. 그러나 실제적으로 색을 표시하기 위해서 픽셀은 유한한 크기를 갖는다. 따라서 픽셀이 한 점을 표시한다는 말을 할 때는 실제상으로는 그 점은 픽셀의 중앙을 의미한다고 생각하면 된다. 예를 들어, 점(1, 3)에 색을 칠한다고 하면 엄밀한 의미로는 [0.5,1.5] x [2.5,3.5]로 주어지는 정사각형 영역을 채우는 것을 의미한다. 픽셀을 가지고 작업을 하는 경우는 항상 정수단위로 움직인다. 이런 사실 때문에 원치 않은 결과가 나오기도 한다. 기울기가 매우 작은 직선을 그리는 작업을 생각해 보자. 예를 들면 기울기가 2/10인 경우에 픽셀상에서의 작업은 0부터 4까지는 같은 y=0인 픽셀에, 5, 9까지는 y=1인 픽셀에 칠해야 한다. 실수 값을 다룰 수 없으므로 4와 5를 지나면서 갑자기 계단을 형성하게 된다. 

이러한 계단 현상은 눈에 보이는 그래픽의 품질을 저하시키는 결과를 초래하므로 이것을 피하는 선 그리기 방법이 있어야 한다. 간단히 생각할 수 있는 방법은 블러링을 주는 방법이다. 우선 기울기가 1보다 작은 경우만 생각하자. 나머지의 경우는 대칭을 생각하면 쉽게 확장이 가능하다. 이 상황에서는 x 좌표값이 독립변수의 역할을 하고, y가 종속변수다. 각각의 픽셀(x=정수)에서 y값을 계산하면 일반적으로 실수 값을 갖는다. 이 y값에 가까운 두 개의 정수 y값에 해당하는 픽셀에 색을 칠하는데, 이때 농도는 거리의 차이에 반하도록 잡는다.

     (x, y)  -> (x, int( y ) ),       명암 가중치 = 1 (y-int(y));

               -> (x, int( y ) + 1 ),  명암 가중치 = y int(y) ;

사용자 삽입 이미지

처음점과 끝점 좀 더 특별한 처리가 필요하다. 여기서도 역시 블러링의 효과를 주기 위해서 가중치를 주는데 y에 의한 가중치에 덧붙여 x에 의한 가중치도 준다. x가 정수 값을 취할 경우에 50%의 가중치를 주는 것을 기준으로 하면, x + 0.5를 기준으로 이것의 정수 값이 벗어나는 정도를 취하면 된다.  그리고  x의 경우는 가장 가까운 정수에서부터 시작한다.

     (x, y)     ->  (int(x + 0. 5), int(y)),  명암 가중치 = (1 (y - int(y))) * (1 - (x + 0.5 - int(x + 0.5)))

                  ->  (int(x + 0. 5), int(y) + 1), 명암 가중치 = (y - int(y)) * (1- ( x + 0.5 - int(x + 0.5)))

사용자 삽입 이미지

오른쪽 끝점은 x에서 오는 가중치는 왼쪽의 complement이므로 x+0.5-int(x+0.5)로 바꾸어야 한다.
아래의 그림은 이 알고리즘으로 그린 선분(위)과 Bresenham 알고리즘(아래)을 이용한 선분을 비교한 것이다( 2배로 확대시킨 그림이다)

사용자 삽입 이미지

이 알고리즘은 1991년도에 Wu에 의해서 제시되었다.

#define FPART(x) ((x) - int(x))     // fractional part of (x)
void DrawLine(double x1, double y1, double x2, double y2) {
    if (x2 < x1) {
    	swap(x1, x2); swap(y1, y2);
    }
    double dx = x2 - x1;
    double dy = y2 - y1;
    double grad = dy / dx ;
    if ( fabs(dx) > fabs(dy)) { // 수평에 가까운 선분;
        //(start point)
        int xs = int(x1 + 0.5);                 //x1 에 제일 가까운 정수 xs ;
        double ys = y1 + grad * (xs - x1);      //픽셀점(xs)에서 y값;
        int ixs = xs ;
        int iys = int(ys) ;
        double xgap = 1 - FPART(x1 + .5);   //x1에 의한 가중치
        setPixel(ixs, iys,     (1 - FPART(ys)) * xgap);
        setPixel(ixs, iys + 1, FPART(ys) * xgap);
        //(end point)
        int xe = int(x2 + 0.5);                 //x2에 제일 가까운 정수 xe;
        double ye = y2 + grad * (xe - x2);      //픽셀점(xe)에서의 y값
        xgap = FPART(x2 + .5);
        int ixe = xe; 
        int iye = int(y2);
        setPixel(ixe, iye,    (1 - FPART(ye)) * xgap);
        setPixel(ixe, iye + 1, FPART(ye) * xgap);
        
        //[ixs+1,..., ixe-1];
        double y = ys + grad ;            // 시작(xs) 다음점에서 y-값;
        for (int ix = ixs + 1; ix < ixe; ix++) {
            setPixel(ix, int(y),      1 - FPART(y));
            setPixel(ix, int(y) + 1,  FPART(y));
            y += grad ;
        }
    } else {
        // 이 경우는 x와 y의 역할을 바꾼다;생략;
    };
};

이 함수는 음수가 아닌 인자에 대해서만 성립한다. 그러나 일반적인 경우로도 쉽게 확장이 가능하다.
참고 :
1. http://en.wikipedia.org/wiki/Xiaolin_Wu's_line_algorithm
2. http://freespace.virgin.net/hugo.elias/graphics/x_wuline.htm
3. http://www.gamedev.net/reference/articles/article382.asp

728x90

'Computational Geometry' 카테고리의 다른 글

Ellipse Drawing Algorithm  (0) 2008.06.04
Circle Drawing Algorithm  (1) 2008.06.03
Brute Force Triangulation  (1) 2008.05.28
Polygon Triangulation (II)  (0) 2008.05.26
Polygon Triangulation  (4) 2008.05.25
Posted by helloktk
,

이미지에서 어떤 유용한 정보를 추출하기 위해서는 이미지가 담고 있는 객체들을 분리하는 작업을 해야 한다. 가장 단순한 것 방법 중의 하나가 이진화(binarization)이다. 이진화는 이미지를 픽셀 값에 따라 0과 1(또는 255)로 값만 가지는 이진 이미지로 변환하는 과정이다. 이진화 작업을 거치면 이미지가 담고 있는 객체를 배경에서 분리할 수 있다. 이때, 어떤 기준으로 픽셀을 분리하는가에 대한 문제가 생긴다. 기준이 되는 임계값(threshold value)의 설정에 대한 다양한 알고리즘이 알려져 있는데, 그중에서 통계적인 방법을 이용한 Otsu 알고리즘이 자연스러운 임계값을 준다.

Otsu 알고리즘은 classification 기법을 이용하고 있다. 임계값을 설정하는 데 있어 비용함수를 설정하고 그 비용함수의 최솟값을 주는 값으로 임계값을 취하는 방식이다. 그럼 어떻게 비용함수를 설정할 것인가? 이미지에서 나타나는 픽셀 값을 2개의 클래스로 분리할 때, 좋은 분리는 각각의 클래스에 속한 픽셀 값의 분포가 유사해야 한다. 즉, 같은 클래스에 들어 있는 픽셀 값의 분산이 작아야 한다는 의미다. 따라서 비용함수는 픽셀 수의 비율로 가중치를 준 각 클래스의 분산을 합산한 것이 될 것이고, 임계값은 이 비용함수를 최소화하는 픽셀 값이다.

     비용함수 = (가중치1 * 분산1) + (가중치2 * 분산2) <= 2개 클래스로 분리 시
                 =   q1 * V1 + q2 * V2 ;      

              q1 =  전체 이미지에서 클래스1에 해당하는 픽셀이 나타날 확률
              q2 =  전체 이미지에서 클래스2에 해당하는 픽셀이 나타날 확률
              V1 = 클래스1에서 픽셀 값의 분산.
              V2 = 클래스2에서 픽셀 값의 분산.

     임계값  -->  MIN ( 비용함수 )

이미지의 픽셀 값 분포는 히스토그램으로 표현되므로, 임계값은 히스토그램으로 분리하는 레벨 값이고, 클래스 1은 그 값보다도 작은 부분, 클래스 2는 큰 부분을 의미한다. 비용함수의 의미를 좀 더 살펴보기 위해서 식을 바꾸어서 적으면

   비용함수 = 전체 분산 - (가중치1*(전체평균 - 평균1)^2 + 가중치2*(전체평균 - 평균2)^2);
                 = V - (q1 * (m1 - m)^2  + q2 * (m2 - m)^2) ;
                         
              V = 전체 분산;
              m = 전체 평균,
              평균1 (m1) = 클래스1의 평균,
              평균2 (m2) = 클래스2의 평균,

여기서 q1*(m-m1)^2 + q2*(m-m2)^2는 클래스들의 분산이다. 전체 분산은 어떤 식으로 클래스를 분리하더라도 항상 일정한 값을 가지므로, 비용함수를 최소화하는 것은 클래스들의 분산을 최대화하는 것과 같다. 새로운 비용함수(엄밀한 의미로 이득함수다)를 이것으로 잡으면
 
             이득함수 = q1 * (m1 - m)^2 + q2 * (m2 - m)^2;
             임계값 --> MAX (이득함수)
 
새로운 이득함수는 약간의 계산을 하면 그 의미가 더 명확한 표현인
             
             이득함수 = q1 * q2 * (m1 - m2)^2 ;

로 쓸 수 있다. 즉, 클래스 분리하는 값은 두 클래스의 평균의 차이(가중치를 갖는)를 최대화시키는 값으로 주어진다.

이 알고리즘의 구현은 히스토그램의 각 레벨에 대해서 좌우를 각각 클래스 1과 2로 설정한 후 이득함수를 계산하여 최댓값을 업데이트하는 방식을 취한다. 0-번째 cumulative histogram과 1-번째 cumulative histogram을 사용하면 각 클래스의 가중치와 평균을 쉽게 계산할 수 있다. 
    cumulative_histogram_0th[k] = 0...k 까지 값이 나타날 확률.
    cumulative_histogram_1th[k]/
cumulative_histogram_0th[k] = 0...k까지 값의 평균.


Otsu 알고리즘은 2개의 클래스뿐만 아니라 여러 클래스로 히스토그램을 분리하도록 확장할 수 있다. 재귀 알고리즘을 이용하면 쉽게 구현할 수 있다. (see: kipl.tistory.com/258)

/* Otsu임계화 예: 설명을 위해 최적화는 하지 않은 것이다. 반드시 cumulative 히스토그램을 이용할 
** 필요는 없다:*/
/* 이득함수의 최대값을 주는 레벨이 연속적으로 여러개 나타나면 평균값을 취하도록 한다(2016.04.26)
*/ 
int OtsuThreshold(BYTE *src, int width, int height, BYTE *dst) {
    double hist[256] = {0.}, chist[256] = {0.}, cxhist[256] = {0.};
    int ntot = width * height;

    // make histogram ;
    for (int i = 0; i < ntot; i++) hist[src[i]] += 1.;
    // normalize;
    for (int i = 0; i < 256; i++) hist[i] /= ntot;

    // make 0-th and 1-st cumulative histogram;
    chist[0] = hist[0];
    cxhist[0] = 0;
    for (int i = 1; i < 256; i++) {
        chist[i] = chist[i - 1] + hist[i] ;               //0-th cumulative histogram ;
        cxhist[i] = cxhist[i - 1] + double(i) * hist[i] ; //1-st cumulative histogram ;
    };
    
    double gain_max = 0.;
    int thresh = 0;   
    double m = cxhist[255];                     //total mean = q1 * m1 + q2 * m2;
    int mul_count = 1;                          //number of degenerate maxima;
    for (int i = 0; i < 256; i++) {
        if (chist[i] == 0.) continue ;
        double q1 = chist[i] ;                  //weight1;
        double q2 = 1 - q1;                     //weight2;
        if (q2 == 0.) break;
        double m1 = cxhist[i] / q1;             //mean1 ;
        double m2 = (m - cxhist[i]) / q2;       //mean2 ;
        double gain = q1 * q2 * (m1 - m2) * (m1 - m2) ;
        if (gain_max < gain) {
            gain_max = gain; 
            thresh   = i;
            mul_count = 1;                      //reset mul_count=1;
        } else if (gain_max == gain)            //degenerate case;
            mul_count ++;
    }
    if (mul_count > 1) thresh = thresh + (mul_count - 1) / 2;    //2016.04.26;

    // threshold image;
    for (int i = 0; i < ntot; i++) dst[i] = (src[i] >= thresh) ? 0xFF : 0x00 ;

    return thresh;
}

 

사용자 삽입 이미지

 

히스토그램 (계산된 임계치는 100이다)

사용자 삽입 이미지

728x90

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

Bright Preserving Histogram Equalization with Maximum Entropy  (0) 2008.07.31
Adaptive Binarization  (2) 2008.07.14
Histogram Equalization  (0) 2008.06.22
FFT2D  (0) 2008.06.10
Hough Transform  (2) 2008.05.22
Posted by helloktk
,

주어진 점집합에서 세 점을 뽑아서 만든 삼각형의 외접원에 다른 점이 하나도 포함하지 않으면 triangulation의 기본 삼각형 cell이 된다. 주어진 점으로 만들 수 있는 삼각형의 총개수가 ${}_nC_3$ 이므로, 기본 삼각형을 찾기 위해서는 이들 각각의 삼각형에 대서 나머지 점을 가지고 incircle 테스트를 수행해야 한다. 따라서 이 알고리즘은 ${\cal O}(n^4)$의 스텝이 필요하게 된다.

/*brute force attack*/
foreach point p1
      foreach point p2 
             foreach point p3
                   foreach point p4 
                          if(incircle(p1,p2,p3,p4))
                                 iscell=false;
                                 break ;
                   endfor;
                   if(iscell) 
                           add(triangle(p1,p2,p3));
              endfor;
       endfor;
endfor;

세 점에 의해서 형성이 되는 외접원은 대수적으로 쉽게 구할 수 있다. 여기서는 좀 더 기하학적인 접근을 쓰면, 평면의 점은 

$$(x, y)\rightarrow (x, y, z=x^2 + y^2)$$

의 mapping에 의해서 3차원 paraboloid 곡면의 점으로 옮길 수 있다. paraboloid 위의 세 점이 형성하는 3차원에서 평면이 paraboloid를 절단하는 곡선을 $x-y$ 평면으로 정사영하면 원이 된다는 것을 쉽게 알 수 있다.(incircle 포스팅 참조). 따라서 주어진 점이 세 점의 외접원에 포함되는가를 테스트하는 작업은 이 점을 paraboloid로 올렸을 때의 점과 (paraboloid로 올려진) 외접원의 3점이 형성하는 3차에서의 평면과 관계를 조사하는 것으로 바꿀 수 있다.

주어진 점이 외접원에 포함되면 paraboloid로 변환된 점은 평면의 아래에 놓이고, 외접원 밖의 점이면 평면 위에 놓이게 된다. 물론 외접원 위의 점은 평면에 놓인다. 따라서 평면의 법선 벡터 구하고, 삼각형의 한 꼭짓점을 기준한 주어진 점의 변위 벡터와 내적을 구하면 내적의 값은 평면 위인지, 아래인지 또는 평면에 놓인 점인가에 따라서 부호가 달라진다. 평면의 수직 벡터를 고정하면(예제는 아래 방향: $n_z < 0$), 평면 위에 놓인 점과의 내적은 음수, 평면 아래에 놓인 점과의 내적은 양수가 되고, 평면의 점과의 내적은 0이다. 

주어진 세 점이 만드는 외접원 내부(and 경계)에 들어가는 점이 없으면 이 삼각형을 선택한다.

** 참고 : Computational Geometry in C(2nd Edition) by Joseph O'Rourke

// input  x[0], x[1],....,x[n-1],
// input  y[0], y[1],....,y[n-1];
// calculate z[0]=x[0]^2+y[0]^2, z[1]=x[1]^2+y[1]^2,...;
int dt4(int n, double x[], double y[]) {
    double *z = new double [n] ;
    for(int i = 0; i < n; i++) 
        z[i] = x[i] * x[i] + y[i] * y[i] ;

    int ntri = 0;
    /* For each triple (i,j,k) */
    for (int i = 0; i < n - 2; i++ )
        for (int j = i + 1; j < n; j++ )
            for (int k = i + 1; k < n; k++ )
                if ( j != k ) {
                    /* Compute normal to triangle (i,j,k)::  outter_product(j-i, k-i)*/
                    double nx = (y[j] - y[i]) * (z[k] - z[i]) - (y[k] - y[i]) * (z[j] - z[i]); 
                    double ny = (x[k] - x[i]) * (z[j] - z[i]) - (x[j] - x[i]) * (z[k] - z[i]);
                    double nz = (x[j] - x[i]) * (y[k] - y[i]) - (x[k] - x[i]) * (y[j] - y[i]);
                    
                    /* Only examine faces on bottom of paraboloid: nz < 0. */
                    int flag = (nz < 0);
                    if (flag) {
                        /* For each other point m */
                        for (int m = 0; m < n; m++) {
                            /* Check if m above (i,j,k)::i점을 기준으로 m 과 
                            ** normal 간의 내적으로 체크(내적이 양수이면 
                            ** m이 원의 내부에 완전히 들어 있는 경우가 된다. 
                            ** 0인 경우는 원주상에 놓인 경우이므로 배제한다
                            */
                            flag &= ((x[m]-x[i])*nx + (y[m]-y[i])*ny + (z[m]-z[i])*nz <= 0);
                        }
                    }
                    if (flag) {
                        ntri++;
                        // (i, j, k)의 외접원이 다른 점을 포함하지 않으므로 이 삼각형은 
                        // 삼각분할의 한 면을 형성하게 된다.
                        // addTriangle(tri(i, j, k));
                    }
                }
    delete[] z ;
    return ntri;
}

사용자 삽입 이미지

 

728x90

'Computational Geometry' 카테고리의 다른 글

Circle Drawing Algorithm  (1) 2008.06.03
Wu's Line Algorithm  (1) 2008.06.02
Polygon Triangulation (II)  (0) 2008.05.26
Polygon Triangulation  (4) 2008.05.25
Polygon Fill  (0) 2008.05.22
Posted by helloktk
,