Geometry & Recognition

제대로 segmented 된 그레이 영상은 원래의 영상이 나타내고자 하는 전경이 잘 표현이 된 것이다. 이 경우의 원래 영상과 segmented 된 영상은 높은 상관관계를 갖는다. 따라서, 세그먼트를 위한 임계값의 설정 기준으로 이 상관계수를 최대로 하는 임계값을 찾는 것도 좋은 방법 중의 하나가 될 수 있다.

여기서 사용할 상관계수는 원래의 영상(A)과 전경과 배경을 그들의 픽셀 평균값으로 대체한 segmented 된 영상(B) 간의 상관계수를 사용한다. 임계값이 $T$인 경우 세그먼트된 영상 B는 

$$B(i,j) = \left\{\begin{array}{ll} m_0, & \text{if}~A(i,j) \le T\\ m_1, &\text{otherwise}\end{array}\right. $$

로 나타난다. 여기서 $m_0$는 배경 픽셀의 평균값이고, $m_1$은 전경 픽셀의 평균값이다. 이 값은 임계값 $T$에 따라 달라진다. 임계값이 높으면 $m_0$는 커지고, 반대로 $m_1$은 작아진다. 

임계값이 $T$일 때 배경 픽셀 비를 $p$, 전경 픽셀 비를 $q(=1- p)$라 하면 segmented된 영상 B는 각 영역에서의 픽셀 값을 평균으로 대체했으므로 원본 영상의 평균과 같다. 또한, 원본 영상의 분산은 임계값에 무관하게 일정한 값을 유지한다. 이를 정리하면,

$$E(A)=E(B)=m=\text{pixel mean}=p m_0 + q m_1$$

$$V(A)=\text{variance} =T\text{-independent} = \text{const}$$

$$V(B)=pm_0^2 + q m_1^2 - m^2 = pq (m_0 - m_1)^2$$

$$E(A,B)= p m_0^2 + q m_1^2 $$

$$E(A,B) - E(A) E(B) = V(B)$$ 이므로, 

\begin{align}\text{Correlation}(A,B) &=\frac{ {E(A,B)-E(A)E(B)} }{\sqrt{V(A)V(B)} } \\ &=\frac{\sqrt{pq(m_0 - m_1)^2 } }{\sqrt{V(A)} }\\ &\propto \sqrt{pq(m_0 -m_1)^2 }\\ &=\sqrt{\text{interclass variance}}\end{align}

즉, 원래의 그레이 영상 A와 전경과 배경 픽셀을 각각의 평균값으로 대체한 영상간의 상관계수는 전경과 배경 두 클래스 간의 분산이 최대일 때 가장 크게 나타난다. 이 기준은 Otsu 알고리즘에서 사용한 기준과 같다.

참고: Otsu Algorithm 구현 예.

Beier & Neely Field Morphing Algorithm:

The Beier-Neely algorithm is a feature-based image-morphing method proposed by Beier and Neely in 1992 [3]. The basic idea of this method is to specify the correspondence between source and target images interactively using a set of line-segment pairs. The mapping of pixels from the source image to the target image is defined using these feature line segments.

모핑 결과: 치타의 이미지 영역밖의 픽셀값은 0으로 처리하였다.

이미지 resampling은 센서에서 받은 신호를 sampling 해서 만든 영상을 확대하거나 줄이는 경우, 또는 기하학적인 변환을 걸쳐서 다시 영상을 생성하는 과정이다. 이 과정에서 출력 영상의 픽셀 좌표 (x, y)에 대응하는 입력 영상의 픽셀 좌표 (u, v)가 요구되는데 일반적으로 이 좌표값은 정수로 주어지지 않는다. 따라서 입력 영상에서 픽셀 값을 얻기 위해서는 (u, v)에 가장 가까운 정수에 해당하는 지점의 픽셀 값을 사용하거나(nearest neighbor interpolation), 아니면 (u, v) 주변의 픽셀 값들을 보간하는 함수를 찾아서 필요한 픽셀 값을 추정할 수 있다: bilinear interpolation, bicubic interpolation,... 

실제 영상 데이터는 하드웨어적으로 스캔라인 방식으로 접근하므로 resampling 알고리즘도 스캔라인 방식에 맞추어서 수평방향으로 처리하고, 그 결과를 다시 수직방향으로 처리하는 separable 알고리즘이 요구된다. Fant 알고리즘은 스캔라인 방식이어서 하드웨어적으로 구현이 쉽도록 설계되었다. 

resampling 중에서도 특히 down sampling을 할 때는 입력 픽셀 정보의 손실로 인해서 나타나는 계단 현상이나 모아레 무늬와 같은 alias를 줄이는 알고리즘이 필요하다. Fant 알고리즘에서는 각각의 입력 픽셀들이 출력 픽셀을 만드는데 얼마나 기여하는지를 고려하여 그 기여만큼 가중치를 주어 합산한다. 이렇게 하면 입력 픽셀의 정보의 손실을 줄어 alias 현상을 억제할 수 있다.

Fant 알고리즘은 일반적으로 출력 영상에 해당하는 입력 픽셀의 좌표가 주어지는 역변환(inverse mapping)에 대해서 사용이 된다. (참고: 이곳에서는 입력 영상에 대응하는 출력 영상의 좌표가 주어지는 정변환(forward mapping)에 대한 Fant 알고리즘의 구현이 있다)

참고:

// dst--->src;
// index ++ ;
// pointer += stride ;
/* stride = horizontal_scan=3; vertical_scan=bytes_per_line;*/
void fant_resample24_inverse(const BYTE src[], const int src_len, const int src_stride,
                              BYTE dest[], const int dest_len, const int dest_stride, 
                              const float f[/*dest_len+1*/]) /*dest_pixel=>src_pixel map*/ 
{
    float inseg, outseg ;
    //dest scanline pixel projects into [0, src_len];
    if (f[dest_len] < 0 || f[0] > src_len) return ;
    //advance to;
    for (int x = 0; f[x] < 0; x++) ;
    int xl = x > 0 ? x - 1 :  x ;
    //
    for (x = dest_len; f[x] > src_len; x--) ;
    int xr = (x == dest_len) ? dest_len - 1 : x ;
    //
    float isf = f[xl + 1] - f[xl];  //inverse_size_factor;
    if (f[x] < 0) {
        inseg = 1;
        outseg = isf + f[xl] ;
    } else {
        inseg = 1 - (f[0] - int(f[0]));
        outseg = isf ;
    }
    //src_initial_pos
    int u = f[x] < 0 ? 0 : int(f[xl]) ;
    src += u * src_stride ;
    int bval = src[0];
    int gval = src[1];
    int rval = src[2];
    src += src_stride ;
    u++ ;
    //src_next_pos ;
    int bnext = src[0];
    int gnext = src[1];
    int rnext = src[2];
    src += src_stride ;
    u++;
    //dest_inital_pos;
    dest += xl * dest_stride ;
    float bsum = 0, gsum = 0, rsum = 0;
    for (x = xl ; x <= xr; ) {
        float binten = inseg * bval + (1 - inseg) * bnext;
        float ginten = inseg * gval + (1 - inseg) * gnext;
        float rinten = inseg * rval + (1 - inseg) * rnext;
        if (inseg < outseg) {
            bsum += binten * inseg ;
            gsum += ginten * inseg ;
            rsum += rinten * inseg ;
            //
            outseg -= inseg ;
            inseg = 1;
            //copy pixel values for next use;
            bval = bnext ;
            gval = gnext ;
            rval = rnext ;
            if (u < src_len) {   //dest를 채우기에 부족한 부분은 마지막에서 계속 반복 사용;
                bnext = src[0];  //여기서 끝낼수도 있음;
                gnext = src[1];
                rnext = src[2];
            }
            src += src_stride ;
            u++ ;
        } else {
            bsum += binten * outseg ; bsum /= isf ;
            gsum += ginten * outseg ; gsum /= isf ;
            rsum += rinten * outseg ; rsum /= isf ;
            dest[0] = (BYTE)min(max(bsum, 0), 255);
            dest[1] = (BYTE)min(max(gsum, 0), 255);
            dest[2] = (BYTE)min(max(rsum, 0), 255);
            dest += dest_stride ;
            //
            bsum = gsum = rsum = 0;
            inseg -= outseg ;
            x++ ;
            if (x == dest_len) break ;
            isf = f[x + 1] - f[x] ;
            outseg = isf ;
        }
    }
}

kipl.tistory.com/41

물체의 형상은 폴리곤이나 폴리곤의 집합으로 근사적으로 표현할 수 있다. 예를 들면 snake나 active shape model (ASM) 등에서 손 모양이나 얼굴의 윤곽, 또는 의료 영상 등에서 장기의 모양 등을 표현할 때 사용이 된다. 이러한 응용에서 주어진 형상을 기준으로 주어진 형상에 정렬을 시켜야 필요가 생긴다. 일반적으로 카메라를 써서 얻은 각 영상에서 추출한 정보들 사이에는 서로 사영 변환의 관계로 연결된다. 그러나 많은 경우에는 in-plane 변형만 고려해도 충분할 때가 많다. 이 경우에 가장 일반적인 형상의 변형은 affine 변환으로 표현된다. 회전(rotation), 평행 이동(translation), 크기 변환(scale transformation) 그리고 층 밀림(shear)을 허용하는 변환이다. 물론, 간단한 경우로는 shear를 제외할 수도 있고 (similarity transformation), 더 간단하게는 크기 변환을 제외할 수도 있다 (isometric transformation).

$N$개의 꼭짓점을 갖는 두 개의 형상 $S=\{(x_1, y_1), (x_2, y_2),..., (x_N, y_N) \}$, $S'=\{(x'_1, y'_1), (x'_2, y'_2),..., (x'_N, y'_N) \}$이 affine 변환에 의해서 연결이 되는 경우에 각 꼭짓점 사이의 관계는

\begin{align} x'_i &= a x_i  + b y_i + t_x \\ y'_i &= c x_i + d y_i + t_y, \quad (i=1,2,..., N);\end{align}

의 6개의 매개변수$(a, b, c, d, t_x, t_y)$에 의해서 기술이 된다(평행 이동: $x/y$축 방향 2개, 회전: 1개, shear: 1개, 스케일: $x/y$축 방향 2개). Affine 변환에 의해서 평행인 두 직선은 변환 후에도 평행인 관계를 유지한다.

꼭짓점 위치는 실제로 다양한 영상처리 과정에 의해서 얻어지므로 필연적으로 노이즈를 포함하게 되어서 일종의 랜덤 변수로 생각해야 한다. 주어진 랜덤 변수에서 최적으로 매개변수를 추출하기 위해 최소자승법을 이용한다. Affine 변환된 좌표와 실제 측정된 좌표 사이의 거리 차이를 최소화하는 매개변수를 찾도록 하자:

$$L=\sum_i \big| x'_i - a x_i - b y_i - t_x \Big|^2 + \big| y'_i - c x_i -d y_i - t_y\big|^2 $$

Affine변환을 규정하는 매개변수를 구하기 위해서는 L을 각 매개변수에 대해서 미분해서 극값을 가질 조건을 구하면 된다:

        ∂L/∂a = -2 * ∑ (x'_i - a * x_i - b * y_i - t_x) * x_i ;
        ∂L/∂b = -2 * ∑ (x'_i - a * x_i - b * y_i - t_x) * y_i ;
        ∂L/∂c = -2 * ∑ (y'_i - c * x_i - d * y_i - t_y) * x_i ;
        ∂L/∂d = -2 * ∑ (y'_i - c * x_i - d * y_i - t_y) * y_i ; 
        ∂L/∂t_x = -2 * ∑ (x'_i - a * x_i - b * y_i - t_x) ;
        ∂L/∂t_y = -2 * ∑ (y'_i - c * x_i - d * y_i - t_y); 

각 식을 0으로 놓아서 얻어지는 연립방정식을 행렬식으로 다시 정리하면,

$$\left[\begin{array}{ccc} S_{xx} & S_{xy} & S_x \\ S_{xy} & S_{yy} & S_y \\ S_x & S_y & N \end{array}\right]\left[ \begin{array}{ll} a & c \\ b & d\\ t_x & t_y \end{array} \right] = \left[\begin{array}{cc} S_{xx'} & S_{x y'} \\ S_{y x'} & S_{yy'} \\ S_{x'} & S_{y'}\end{array} \right]$$

여기서,
\begin{align} & S_{xx}= ∑ x^2, ~S_{yy} = ∑ y^2, ~S_{xy} = ∑ xy, \\ &S_x = ∑ x, ~S_y = ∑ y, ~S_{x'} = ∑ x', ~S_{y'} = ∑ y' \\ & S_{xx'} = ∑ xx', ~S_{xy'} = ∑ xy', ~S_{yx'} =∑ yx' \end{align} 이다.

// dst = (A,T)src;
//  [u]  = [ A0 A1 ][x] + A4
//  [v]  = [ A2 A3 ][y] + A5
//
BOOL GetAffineParameter(POINT *srcPts, POINT *dstPts, int n, double AT[6]) {
    double Sx, Sy, Sxx, Sxy, Syy;
    double Su, Sv, Sxu, Sxv, Syu, Syv ;
    double A[9], invA[9] ;
    double det ;
    Sx = Sy = Sxx = Sxy = Syy = 0;
    Su = Sv = Sxu = Sxv = Syu = Syv = 0;
    for (int i = 0; i < n; i++) {
        double x = srcPts[i].x, y = srcPts[i].y ;
        double u = dstPts[i].x, v = dstPts[i].y ;
        Sx += x;        Sy += y ;
        Sxx += (x * x); Sxy += (x * y); Syy += (y * y);
        Su += u;        Sv += v ;
        Sxu += (x * u); Sxv += (x * v); Syu += (y * u); Syv += (y * v);
    }
    A[0] = Sxx; A[1] = Sxy; A[2] = Sx;
    A[3] = Sxy; A[4] = Syy; A[5] = Sy;
    A[6] = Sx ; A[7] = Sy ; A[8] = n ;
    det = (A[0]*(A[4]*A[8]-A[5]*A[7])-A[1]*(A[3]*A[8]-A[5]*A[6])+A[2]*(A[3]*A[7]-A[4]*A[6])) ;
    if (det != 0.) {
        det = 1. / det; 
        invA[0] = (A[4]*A[8] - A[5]*A[7]) * det;
        invA[1] = (A[2]*A[7] - A[1]*A[8]) * det;
        invA[2] = (A[1]*A[5] - A[2]*A[4]) * det;
        invA[3] = (A[5]*A[6] - A[3]*A[8]) * det;
        invA[4] = (A[0]*A[8] - A[2]*A[6]) * det;
        invA[5] = (A[2]*A[3] - A[0]*A[5]) * det;
        invA[6] = (A[3]*A[7] - A[4]*A[6]) * det;
        invA[7] = (A[1]*A[6] - A[0]*A[7]) * det;
        invA[8] = (A[0]*A[4] - A[1]*A[3]) * det;
    }
    else return FALSE;

    AT[0] = invA[0] * Sxu + invA[1] * Syu + invA[2] * Su;
    AT[1] = invA[3] * Sxu + invA[4] * Syu + invA[5] * Su;
    AT[4] = invA[6] * Sxu + invA[7] * Syu + invA[8] * Su;
    AT[2] = invA[0] * Sxv + invA[1] * Syv + invA[2] * Sv;
    AT[3] = invA[3] * Sxv + invA[4] * Syv + invA[5] * Sv;
    AT[5] = invA[6] * Sxv + invA[7] * Syv + invA[8] * Sv;
    return TRUE ;
};

아래의 그림은 지문에서 얻은 특징점을 가지고 변환을 한 것이다. 밑에 그림이 기준 template (붉은 점)이고 윗 그림은 이 기준  template와 입력된 지문의 특징점(노란 점+ 녹색점) 사이에 서로 메칭이 되는 특징점(노란색)을 찾고, 그것을 기준으로 두 지문 영상 간의 affine 파라미터를 찾아서 기준 template을 변환시킨 것이다. 이렇게 하면 새로 찾은 특징점 중에서 기준 template에 없는 특징점(녹색점)을 발견할 수 있고, 이 특징점을 기준 template에 추가하여서 좀 더 넓은 범위를 커버할 수 있는 template을 만들 수 있다. 물론 추가된 녹색점이 신뢰할 수 있는 것인가에 대한 판단을 하기 위해서는 추가적인 정보가 더 요구된다.

컬러 이미지가 있을 때, 몇 가지 색상이 쓰였는가를 찾아내는 것은 간단한 문제가 아니다. 왜냐면 24비트 RGB 컬러 영상에서 최대로 가능한 컬러의 수는 256*256*256=16777216나 되어서 간단히 히스토그램을 이용해서 셀 문제가 아니다. 정수의 히스토그램을 이용한다고 하더라도 대략 64MB의 메모리가 필요하다.

그런데 한 영상에는 픽셀 수 이상의 컬러를 담을 수 없으므로, 실제로 필요한 메모리는 정수 배열로 하면, 4*(이미지 폭)*(이미지 높이) 정도만 필요로 한다. 정수 배열을 준비하여 컬러 값을 옮기고 나서 정렬 알고리즘을 이용해서 정수 배열을 크기 순서대로 정렬하면 컬러의 수와 도수를 구할 수 있다(적은 사이즈의 이미지/메모리가 작을 때 유리)

물론 STL의 map을 이용해서 간단히 컬러 히스토그램을 구현할 수 있다 (이것은 아무래도 모기를 잡는데 도끼를 휘두르는 것 같은 느낌이다). 이들 방법을 동원해서 프로그램을 하면 큰 이미지에 대해서는 그다지 빠르다는 느낌이 없다.
         std::map <DWORD, int> table ; 
         .... table [color]++;

만약에 사용된 컬러의 도수가 필요하지 않고 사용된 색상의 총 수에만 관심이 있다면 좀 더 빠른 방법을 찾을 수 있다. 색상이 사용되었는지 안 되었는지 표시하는 데는 1비트의 정보만 기록하면 되므로 총 16777216개의 색상이 사용되었는지 안되었는지를 기록하기 위해서는 1 바이트가 8개의 비트를 저장할 수 있으므로 16777216/8 = 2MB의 메모리만 필요로 한다. 이 값은 이미지의 크기에 상관없이 일정하다. 

아래 소스는 사용된 색상의 수를 총 수를 얻는 법을 구현한 것이다. 주어진 컬러 값을 8로 나누면 비트 정보를 기록할 바이트 배열의 위치가 나오고, 컬러 값을 8로 나눈 나머지가 그 바이트에서 비트의 위치이므로 해당 비트가 체크가 안되어 있으면 새로운 컬러를 찾은 것이므로 비트 마스크를 이용해서 비트를 기록하면 된다.

static BYTE bitmasks[] = {0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80};//1바이트에서 비트 위치를 표시하기 위한 마스크
int CountColors(BYTE *rgbimage, int w, int h) {
    std::vector<BYTE> (1024 * 1024 * 2, 0);
    int count = 0; 
    int k = w * h ;
    while (k--) {
        DWORD color = (*(DWORD*)rgbimage) & 0xFFFFFF;  //Little endian machine ;
        DWORD idx = color >> 3;
        DWORD bit_pos = color - (idx << 3);
        BYTE *a = &table[idx];
        if (!((*a) & bitmasks[bit_pos])) {//if bit is not set;
            *a |= bitmasks[bit_pos] ; //set bit;
            count++ ;
        };
        rgbimage += 3; //rgb
    }
    return count ;
} ;