Geometry & Recognition

영상에는 보통 우리가 관심이 있는 물체와 배경을 동시에 담고 있다. 그런데 어떤 영상의 경우 전체 픽셀 값의 다이내믹 레인지(dynamic range)가 매우 좁아서, 달리 이야기하면 영상의 히스토그램을 살펴볼 때 픽셀 값의 분포가 매우 좁은 영역에 한정이 되어있으면 물체와 배경의 구별이 잘 안되는 영상이 되게 된다. 이러한 영상의 경우 몇 가지 영상 처리 방법을 이용하면 좋은 contrast를 갖는 영상을 만들 수 있다.

Histogram Equalization(HE)은 이러한 기법 중의 하나로 좁은 영역에 분포하는 픽셀 값을 가능한 모든 범위에 골고루 분포하도록 재배치하여서 영상의 contrast를 증대시키는 방법이다. Histogram Equalization에서 사용하는 픽셀 값의 mapping은 

(1) 출력 영상의 픽셀 값 범위는 전 그레이 레벨이어야 하고 

(2) 각각의 그레이 레벨에 해당하는 픽셀의 수가 거의 일정하여야 한다: 픽셀 값이 이산적이어서 완벽하게 상수가 될 수는 없다, 좀 더 정확히는 출력 영상의 히스토그램 누적합이 그레이 레벨에 선형 비례해야 한다는 조건으로 결정된다.

입력 영상에서 값이 k 인 픽셀은 출력 영상에서 어떤 그레이 값을 가져야 할까? 조건에 따르면 입력 영상의 k 그레이 레벨이 출력 영상의 j 그레이 레벨로 매핑이 된다면, 출력 영상에서는 0부터 j까지 그레이 레벨에 해당하는 픽셀 수는 모든 그레이 레벨에서 픽셀 분포가 균일하다는 조건 때문에 

로 주어진다. 이 값과 입력 영상에서 0부터 k까지 그레이 레벨을 갖는 픽셀의 수가 같으면 된다.

입력 영상에서 그레이 값이 0부터 k까지인 픽셀 수는 

총 픽셀 수는 cumulative histogram[255]로 주어지므로 원하는 매핑 j = T [k]는

이 매핑에 의하면 입력 영상의 픽셀이 가지는 가장 큰 그레이 값은 출력 영상에는 255에 해당하게 된다. 출력 영상은 정의에 의해서 cumulative histogram이 직선적으로 증가하는 형태를 가지게 된다. 큰 영상의 경우 이 매핑을 Lookup 테이블로 저장하면 연산의 횟수를 줄일 수 있다.

note)
1. 0차 cumulative histogram을 이용하면 0부터 k까지의 histogram의 합은 한 번의 호출로 얻어진다.
2. 0차 cumulative histogram의 마지막 bin의 값은 전체 픽셀 수의 합이 된다.

// histogram의 합은 cumulative histogram을 구성하면 된다;
void HistogramEqualize(BYTE *src, int width, int height, BYTE* dst) {
    int histogram[256] = {0}, cumhist[256];
    for (int k = width * height; k-->0;) ++histogram[src[k]];
    
    // make 0-st order cumulative histogram;
    cumhist[0] = histogram[0];
    for (int k = 1; k < 256; k++) cumhist[k] = cumhist[k-1] + histogram[k] ;
    
    int map[256];
    double v = 255.0 / cumhist[255]; 	// cumhist[255] = number of pixels;
    for (int k = 0; k < 256; k++) {
    	int v = int(v * cumhist[k]);
        map[k] = v > 255 ? 255: v;
    }
    
    for (int k = width * height; k-->0;) dst[k] = map[src[k]];
};

평탄화 전의 히스토그램(붉은색) 및 누적 히스토그램(연두색)

평탄화 후의 히스토그램(붉은색) 및 누적 히스토그램(연두색)

• Forward transformation:
  $$X_{k} =  \sum_{n = 0} ^{N-1} x_n e^{ -2i  \pi k n /N }, ~\quad k=0,1,..., N-1;$$
• Backward transformation:
  $$x_n = \frac{1}{N}\sum_{k = 0} ^{N-1} X_k  e^{ + 2i  \pi k n /N }, ~\quad n=0,1,..., N-1;$$
• Butterworth Filter:
  $$H(x, y) = \frac{1}{1+\left( \frac{x^2 +y^2}{H_c^2 }\right)^N} ,\quad N=\text{order},~H_c =\text{cut-off};$$
• 주어진 수가 $2^n$인가?(kipl.tistory.com/84). 입력 데이터 개수가 2의 거듭제곱이 아니면 부족한 부분은 0으로 채워서 $2^n \ge \text{데이터 개수}$ 크기가 되도록 만든다. 

Cooley–Tukey FFT algorithm:  

/* 1차원 FFT:: dir=1 for forward, -1 for backward;
** x = real-component, and y = imaginary componet.
** nn = length of data: 2의 거듭제곱이어야 한다.
*/
int fft1d(int dir, int nn, double x[], double y[]) {
    /* Calculate the number of points */
    int m = 0;
    while (nn > 1) {
        nn >>= 1;  m++;
    }
    nn = 1 << m;
    
    /* Do the bit reversal */
    int i2 = nn >> 1;
    int j = 0;
    for (int i = 0; i < nn - 1; i++) {
        if (i < j) { // swap(x[i], x[j]), swap(y[i], y[j]);
            double tx = x[i], ty = y[i];
            x[i] = x[j];  y[i] = y[j];
            x[j] = tx;    y[j] = ty;
        }
        int k = i2;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
    
    /* Compute the FFT */
    double c1 = -1.0;
    double c2 = 0.0;
    int l2 = 1;
    for (l = 0; l < m; l++) {
        int l1 = l2;
        l2 <<= 1;
        double u1 = 1.0;
        double u2 = 0.0;
        for (int j = 0; j < l1; j++) {
            for (int i = j; i < nn; i += l2) {
                int i1 = i + l1;
                double t1 = u1 * x[i1] - u2 * y[i1];
                double t2 = u1 * y[i1] + u2 * x[i1];
                x[i1] = x[i] - t1;
                y[i1] = y[i] - t2;
                x[i] += t1;
                y[i] += t2;
            }
            double z = u1 * c1 - u2 * c2;
            u2 = u1 * c2 + u2 * c1;
            u1 = z;
        }
        c2 = sqrt((1.0 - c1) / 2.0);
        if (dir == 1) c2 = -c2;
        c1 = sqrt((1.0 + c1) / 2.0);
    }
    
    /* Scaling for forward transform */
    if (dir == 1) {
        for (i = 0; i < nn; i++) {
            x[i] /= double(nn);
            y[i] /= double(nn);
        }
    }
    return 1;
};

/* 2차원 FFT. data = (2 * nn) x mm 크기의 버퍼;
** 각 행은 nn개의 실수 성분 다음에 nn개의 허수 성분이 온다.
** nn, mm은 각각 행과 열은 나타내며, 2의 거듭제곱이어야 한다.
** isign = 1 for forward, -1 for backward. 
** copy = buffer of length 2*mm ;
*/
void fft2d(double* data, int nn, int mm, int isign, double* copy) { 
    int stride = nn << 1; //row_stride;
    /* Transform by ROWS for forward transform */
    if (isign == 1) {
        int index1 = 0;
        for (int i = 0; i < mm; i++) {
            // real = &data[index1], imag = &data[index1 + nn];
            fft1d(isign, nn, &data[index1], &data[index1 + nn]);
            index1 += stride;
        }
    }
    
    /* Transform by COLUMNS */
    for (int j = 0; j < nn; j++) {
        /* Copy pixels into temp array */
        int index1 = j ;
        int index2 = 0;
        for (i = 0; i < mm; i++) {
            copy[index2] = data[index1];
            copy[index2 + mm] = data[index1 + nn];
            index2++;
            index1 += stride;
        }
        
        /* Perform transform */
        fft1d(isign, mm, &copy[0], &copy[mm]);
        
        /* Copy pixels back into data array */
        index1 = j;
        index2 = 0;
        for (i = 0; i < mm; i++) {
            data[index1] = copy[index2];
            data[index1 + mm] = copy[index2 + mm];
            index2++;
            index1 += stride;
        }
    }
    
    /* Transform by ROWS for inverse transform */
    if (isign == -1) {
        int index1 = 0;
        for (i = 0; i < mm; i++) {
            fft1d(isign, nn, &data[index1], &data[index1 + nn]);
            index1 += stride;
        }
    }
}

void butterworth( double * data, int Xdim, int Ydim, 
                 int Homomorph, int LowPass,
                 double Power, double CutOff, double Boost);

// H * F;
void butterworth(double * data, int Xdim, int Ydim, 
                 int Homomorph, int LowPass,
                 double Power, double CutOff, double Boost) 
{
    double CutOff2 = CutOff * CutOff;
    /* Prepare to filter */
    int stride = Xdim << 1;  // width of a sigle row(real + imag);    
    int halfx = Xdim >> 1;
    int halfy = Ydim >> 1;
    double *rdata = &data[0];   //real part;
    double *idata = &data[Xdim];//imag part;
    /* Loop over Y axis */
    for (int y = 0; y < Ydim; y++) {
        int y1 = (y < halfy) ? y: y - Ydim;
        /* Loop over X axis */
        for (int x = 0; x < Xdim; x++) {
            int x1 = (x < halfx) ? x: x - Xdim;
            /* Calculate value of Butterworth filter */
            double filter;
            if (LowPass)
                Filter = (1 / (1 + pow((x1 * x1 + y1 * y1) / CutOff2, Power)));
            else if ((x1 != 0) || (y1 != 0))
                Filter = (1 / (1 + pow( CutOff2 / (x1 * x1 + y1 * y1), Power)));
            else
                Filter = 0.0;
            if (Homomorph)
                Filter = Boost + (1 - Boost) * Filter;
            /* Do pointwise multiplication */
            rdata[x] *= Filter;
            idata[x] *= Filter;
        };
        rdata += stride; //go to next-line;
        idata += stride; 
    }
};

Butterworth Low Pass filter 적용의 예(order=2, cutoff=20)

PowerSpectrum 변화

이미지에서 어떤 유용한 정보를 추출하기 위해서는 이미지가 담고 있는 객체들을 분리하는 작업을 해야 한다. 가장 단순한 것 방법 중의 하나가 이진화(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이다)

주어진 영상에서 선분을 찾고자 하면 어떻게 해야 할까. 우선 에지 영상을 만들어 선분을 강조하고, 선분이 하나만 있는 경우에는 에지 데이터를 가지고 line-fitting 알고리즘을 적용해서 해당 선분을 기술하는 파라미터를 추출할 수 있다. 그러나 영상이 선분을 여러 개 포함하는 경우에는 이 방법을 이용하기 어렵다(불가능한 것은 아니지만 노이즈 등에 의한 영향을 고려해야 하므로 복잡한 과정을 거쳐야 한다). 평면에서 선분은 원점까지 거리와 그것의 수직(수평)인 방향만 주어지면 결정이 된다. 직선에서 원점까지 거리를 $r$, 법선이 $x$-축과 이루는 각도가 $θ$면

$$r = \cos(θ)  x + \sin(θ)  y;$$

의 형태로 주어진다. 즉, $(r,θ)$ 한쌍이 주어지면 직선 한 개가 정의된다 (주어진 $(r,θ)$ 값이 만드는 직선을 $a  x + b  y = c$  꼴로 쓰면 계수는 $a = \cos θ$, $b = \sin θ$, $c = r$로 표현된다). 

이 직선 방정식은 $rθ$-평면의 stripe=$[0,∞) \times  [0,2π]$을 $xy$-평면으로 보내는 변환으로도 생각할 수 있다. 이 변환은 $rθ$-평면의 한 점을 $xy$-평면의 한 직선으로 보낸다. 역으로는 $xy$-평면의 한 점은 $rθ$-평면의 한 곡선으로 변환이 된다. 직선 상의 점들은 같은 원점 거리=$r$, 같은 기울기=$\theta$를 가지므로, 직선 위의 각 점들이 $rθ$-평면에 만드는 곡선들은 공통의 교점을 가지게 된다. 이 교점의 위치가 $xy$-평면에서 직선을 기술하는 파라미터를 준다.

Hough Transform은 이 변환의 특성을 이용한 것이다. 이미지에서 에지에 해당하는 점들을 위의 변환에 의해서 $rθ$-평면의 곡선으로 보내는 것이다 (프로그램적으로는, 에지점 좌표 $(x, y)$가 주어지면 $θ=0 \rightarrow π$까지 일정 간격으로 증가시키면서 곡선 $r=x \cos(θ) + y\sin(θ)$의 값을 구해서 메모리 상의  $[r,θ]$ 지점의 도수를 1씩 증가시킨다). 만약 에지에 해당되는 점들이 직선의 관계를 가지게 되면, 각 점들에 해당하는 $rθ$-평면에서 곡선들은 한 교점을 형성하게 될 것이다. 긴 선분은 $rθ$-평면의 한 점에서 더 많은 곡선들의 교점을 형성하게 된다. 따라서, $rθ$-평면에서 이러한 교점의 히스토그램을 구하여, 교점 수가 일정 이상 누적이 된 경우를 취하면, $xy$-이미지 평면에서 일정한 길이 이상을 갖는 선분을 골라낼 수 있다.

그러면 왜 $rθ$-평면으로 변환을 생각하는 것일까? 직선이 $y= a x+b$의 형태로 표시되므로 기울기와 $y$축과의 교점을 기술하는 $ab$-공간을 이용할 수도 있지만, 이 경우에 $y$-축에 평행인 직선의 경우에 기울기 $a$ 값이 무한히 커지는 경우가 발생하여 저장공간 할당에 문제가 생긴다. 실제적인 문제에서 되도록이면 compact 한 공간으로 변환을 고려해야 하는 것이다. $rθ$-공간으로의 변환은 유한한 이미지의 경우 항상 유한한 $rθ$-공간의 영역으로 변환된다.

아래의 그림은 $x-y=-5$인 직선상의 세 점 $(5,10)$, $(10,15)$, $(15,20)$에 해당하는 $rθ$-평면상의 세 곡선을 보인 것이다: $r = 5 \cos(θ) + 10\sin(θ)$, $r = 10 \cos(θ) + 15  \sin(θ)$, $r = 15 \cos(θ) + 20 \sin(θ)$. 이 세 곡선은 $(r,θ) = (5/\sqrt{2}, 3 \pi /4)$ 지점에서 만남을 알 수 있다. 이 만나는 지점을 구하면, 거꾸로 세 점이  $x-y=-5$ 인 직선 위의 점들이었다는 것을 추정할 수 있다.

std::vector<Line> HoughTransformLine(BYTE* image, int width, int height, /*background=0*/
                                     double rho/*=2*/, 
                                     double theta/*=Pi/180.*/,
                                     int threshold/*=20*/) 
{
    /* use cosine and sine look-up tables*/
    const double irho = 1 / rho;
    const int numangle = (int) (Pi / theta);
    const int numrho = (int) (((width + height) * 2 + 1) / rho);
    const int rwidth = (numrho + 2);
    std::vector<int> accum((numangle + 2) * (numrho + 2), 0);
    double ang; int n;
    for (int  y = 0; y < height; y++ ){
        for (int x = 0; x < width; x++ ){
            if ( image[y * width + x] != 0 ){ // only for foreground pixels;
                for (ang = 0, n = 0; n < numangle; ang += theta, n++ ) {
                    double rho = (x * cos(ang) + y * sin(ang)) * irho ;
                    int r = rho >= 0 ? int(rho + .5) : (-int(-rho + .5));//round to nearest int;
                    // accum을 이차원배열로 생각하고, 1 픽셀의 border를 두는 형태로 한다.
                    // 이는 아래의 local-maxima를 찾을때 경계를 고려할 필요가 없어서 유리하다.
                    r += (numrho - 1) / 2;
                    accum[(n + 1) * rwidth + r + 1]++;
                }
            }
        }
    }
    //find local maxima;
    std::vector<Line> vecLine;
    for (int  r = 0; r < numrho; r++ ) {
        for (int  n = 0; n < numangle; n++ ) {
            int base = (n + 1) * rwidth + r + 1;
            //test whether it is local-maximum(4-way);
            if ( accum[base] > threshold &&
                 accum[base] > accum[base - 1] && accum[base] > accum[base + 1] &&
                 accum[base] > accum[base - rwidth] && accum[base] > accum[base + rwidth] )
            {               
                Line line;
                line.rho = (r - (numrho - 1) * .5) * rho;
                line.angle = n * theta;
                vecLine.push_back( line );
            }
        }
    }
    return vecLine;
}

아래 그림에서 첫 번째는 소스 이미지이고, 이 이미지에서 $rθ$-평면으로 Hough Transform 한 것이 두 번째 것이다(rho=2, theta=Pi/360). 원본 이미지에는 8개의 선분이 있고, 4 개씩 같은 방향을 갖고 있다. Hough Transform 된 이미지에서 이러한 특징을 볼 수 있다(가로축이 $θ$이고 세로축이 $r$이다. $r$ 축은 가운데가 $r=0$이다). 누적이 많이 된 부분이 두 군데의 동일한 θ에서 나타나고 있다. 그러나 결과에서 8개의 피크를 분리하기는 쉽지 않다. 위의 코드는 각각의 점에서 4방향으로 체크하여 극대 값을 찾고 있으나, 항상 잘 동작하는 것은 아니다.

local-maxima를 좀 더 잘 잡고 싶으면, 위처럼 주변의 4점만 체크할 것이 아니라, 윈도를 좀 더 크게 잡고, 그 윈도 내에서 최댓값을 찾아보는 것도 한 방법이 된다.

/**
** http://blog.naver.com/helloktk/80051779331
*/