Geometry & Recognition

Kuwahara filter는 주어진 픽셀을 한 코너로 하는 4개 block 중에서 분산이 가장 작은 블록의 평균으로 중앙 픽셀 값을 대체하는 비선형 filter다. 이 필터를 적용하면  균일한 영역에서  spike noise 픽셀은 주변의 영역의 픽셀 값으로 대체되므로 제거될 수 있고, edge 라인을 기준으로 수직방향으로는 거의 균일한 영역이므로 필터링을 하더라도 edge가 뭉개지지 않음을 알 수 있다. Kuwahara filter는 이미지의 노이즈를 제거하고 edge를(위치) 보존하는 특성을 가진다. 윈도의 크기가 5x5인 경우 4개의 3x3영역에서 평균과 분산을 계산해서 중앙 픽셀의 값을 결정한다. integral 이미지를 사용하면 block의 평균과 분산을 매번 새롭게 계산할 필요가 없으므로 연산 비용이 많이 감소한다.

Kuwahara filter보다 더 향상된 filter로는 symmetric nearest neighbour (SNN) filter가 있다. 이 filter는 기준 픽셀에서 대칭 지점에 있는 두 지점의 픽셀 값 중에서 중앙 픽셀 값에 가까운 것을 취해서 평균값을 취한다.

아래 코드는 gray 이미지에 대해 Kuwahara filter을 구현한 것으로 컬러 이미지의 경우는 RGB 채널별로 따로 적용하면 된다. 그리고 경계 영역에서는 처리를 하지 않도록 만들었다.

#define GET_SUM(x,y)  sum[ (y) * width + (x)]
#define GET_SUM2(x,y) sum2[(y) * width + (x)]
#define BLK_SUM(x,y) ((GET_SUM(x,y) + GET_SUM(x-hs1,y-hs1) \
                        - GET_SUM(x-hs1,y) - GET_SUM(x,y-hs1)))
#define BLK_SUM2(x,y) ((GET_SUM2(x,y) + GET_SUM2(x-hs1,y-hs1) \
                        - GET_SUM2(x-hs1,y) - GET_SUM2(x,y-hs1)))
void integral_image(BYTE *image, int width, int height, int *sum);
void integral_image_sq(BYTE *image, int width, int height, int *sum2);
void kuwahara_filter ( BYTE *image, int width, int height, int size, BYTE *out ) {
    std::vector<int> sum(width * height, 0);
    stt::vector<int> sum2(width * height, 0);
    // integral image; for window mean;
    integral_image(image, width, height, &sum[0]);
    // integral image of square; for window variance;
    integral_image_sq(image, width, height, &sum2[0]);
    
    // do filtering;
    size = ((size >> 1) << 1) + 1; // make window size be an odd number;
    int offset = ( size - 1 ) >> 1;
    int hs1 = ( size + 1 ) >> 1;
    int nn = hs1 * hs1;    // # of pixels within the  window;
    int var[4], mean[4];
    for ( int y = height - hs1; y-- > hs1;) {
        for ( int x = width - hs1; x-- > hs1; ) {
            mean[0] = BLK_SUM ( x, y ) / nn;
            mean[1] = BLK_SUM ( x + offset, y ) / nn;
            mean[2] = BLK_SUM ( x, y + offset ) / nn;
            mean[3] = BLK_SUM ( x + offset, y + offset ) / nn;
            var[0] = BLK_SUM2 ( x, y ) / nn                   - mean[0] * mean[0];
            var[1] = BLK_SUM2 ( x + offset, y ) / nn          - mean[1] * mean[1] ;
            var[2] = BLK_SUM2 ( x, y + offset ) / nn          - mean[2] * mean[2] ;
            var[3] = BLK_SUM2 ( x + offset, y + offset ) / nn - mean[3] * mean[3] ;
            int vmin = var[0], winner = 0;
            for ( int i = 1; i < 4; ++i )
                if ( var[i] < vmin )
                    vmin = var[i], winner = i;
            int a = mean[winner];
            out[y * width + x] = a > 255 ? 255 : a < 0 ? 0 : a;
        }
    }
};

이미지 처리에서 각도를 다루게 되면 삼각함수를 이용해야 한다. 일반적인 PC 상에서는 CPU의 성능이 크게 향상이 되어 있고, floating point 연산이 잘 지원이 되므로 큰 문제가 없이 삼각함수를 바로 이용할 수가 있다. 그러나 embeded 환경에서는 연산 속도가 그다지 빠르지 않고 floating point 연산은 에뮬레이션을 하는 경우가 대부분이다. 이런 조건에서 루프 내에서 연속적으로 삼각함수를 호출하는 것은 성능에 많은 영향을 주게 되므로 다른 방법을 고려해야 한다. 일반적으로 각도를 점진적으로 증가시키면 cosine이나 sine 값을 계산해야 할 필요가 있는 경우에는 아래와 같이 처리한다:

for (int i = 0; i < nsteps; ++i) {
    double angle = i * angle_increment;
    double cos_val = cos( angle ) ;
    double sin_val = sin( abgle ) ;
    // do something useful with cos_val, and sin_val;
    // .......................;
}

물론 미리 계산된 cosine/sine의 값을 가지고 처리할 수도 있다. 하지만 매우 작은 각도 단위로 계산하는 경우는 lookup table의 크기가 매우 커지게 되므로 직접 계산을 하는 것보다 이득이 없을 수 있다. 이런 경우에는 다음의  삼각함수의 관계식을 이용하면 매번 cosine이나 sine 함수의 호출 없이도 필요한 값을 얻을 수 있다;

$$\cos(\theta+\delta)-\cos(\theta)=- (2\sin^2(\delta/2) \cos(\theta) + \sin(\delta)\sin(\theta))= -(\alpha \cos(\theta) + \beta \sin (\theta));$$

$$\sin(\theta+\delta)-\sin(\theta)= -(2\sin^2(\delta/2) \sin(\theta) - \sin(\delta)\cos(\theta))= -(\alpha \sin(\theta) - \beta \cos (\theta));$$

여기서 $\alpha = 2 \sin^2(\delta/2)$, $\beta = \sin(\delta)$다. 처음에 한번 계산된 $\alpha$와 $\beta$ 값을 이용하면, 추가적인 계산이 없이 점진적으로 $\cos$과 $\sin$의 값을 계산할 수 있다.

    int i = 0 ;
    double cos_val = 1 ; // = cos(0);
    double sin_val = 0 ; // = sin(0);
    double beta = sin(angle_increment) ;
    double alpha = sin(angle_increment/2) ;
    alpha = 2 * alpha * alpha ; 
    do {
        // do something useful with cos_val and sin_val ;
        // ...........................;
        // calculate next values;

        double cos_inc = (alpha * cos_val + beta * sin_val) ;
        double sin_inc = (alpha * sin_val - beta * cos_val) ;
        cos_val -= cos_inc ;
        sin_val -= sin_inc ;
    } while (++i < nsteps);

이 방법은 이론상 정확한 방법이나 round-off error가 전파가 되어서 반복 회수가 커질수록 에러가 누적이 될 위험이 있다. 이를 방지하기 위해 중간에 카운터를 두어서 일정한 반복 이후에는 cosine/sine값 계산을 다시 하여 에러를 리셋시키면 된다.