Geometry & Recognition

참고:  ndevilla.free.fr/median/median/src/optmed_bench.c

/*
 * Optimized median search on 9 values
 */
#define PIX_SORT(a,b) { if ((a) > (b)) PIX_SWAP((a), (b)); }
#define PIX_SWAP(a,b) { BYTE temp = (a); (a) = (b); (b) = temp; }
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]) ;
}

BYTE opt_med25(BYTE p[]) { 
    PIX_SORT(p[0], p[1]) ;   PIX_SORT(p[3], p[4]) ;   PIX_SORT(p[2], p[4]) ;
    PIX_SORT(p[2], p[3]) ;   PIX_SORT(p[6], p[7]) ;   PIX_SORT(p[5], p[7]) ;
    PIX_SORT(p[5], p[6]) ;   PIX_SORT(p[9], p[10]) ;  PIX_SORT(p[8], p[10]) ;
    PIX_SORT(p[8], p[9]) ;   PIX_SORT(p[12], p[13]) ; PIX_SORT(p[11], p[13]) ;
    PIX_SORT(p[11], p[12]) ; PIX_SORT(p[15], p[16]) ; PIX_SORT(p[14], p[16]) ;
    PIX_SORT(p[14], p[15]) ; PIX_SORT(p[18], p[19]) ; PIX_SORT(p[17], p[19]) ;
    PIX_SORT(p[17], p[18]) ; PIX_SORT(p[21], p[22]) ; PIX_SORT(p[20], p[22]) ;
    PIX_SORT(p[20], p[21]) ; PIX_SORT(p[23], p[24]) ; PIX_SORT(p[2], p[5]) ;
    PIX_SORT(p[3], p[6]) ;   PIX_SORT(p[0], p[6]) ;   PIX_SORT(p[0], p[3]) ;
    PIX_SORT(p[4], p[7]) ;   PIX_SORT(p[1], p[7]) ;   PIX_SORT(p[1], p[4]) ;
    PIX_SORT(p[11], p[14]) ; PIX_SORT(p[8], p[14]) ;  PIX_SORT(p[8], p[11]) ;
    PIX_SORT(p[12], p[15]) ; PIX_SORT(p[9], p[15]) ;  PIX_SORT(p[9], p[12]) ;
    PIX_SORT(p[13], p[16]) ; PIX_SORT(p[10], p[16]) ; PIX_SORT(p[10], p[13]) ;
    PIX_SORT(p[20], p[23]) ; PIX_SORT(p[17], p[23]) ; PIX_SORT(p[17], p[20]) ;
    PIX_SORT(p[21], p[24]) ; PIX_SORT(p[18], p[24]) ; PIX_SORT(p[18], p[21]) ;
    PIX_SORT(p[19], p[22]) ; PIX_SORT(p[8], p[17]) ;  PIX_SORT(p[9], p[18]) ;
    PIX_SORT(p[0], p[18]) ;  PIX_SORT(p[0], p[9]) ;   PIX_SORT(p[10], p[19]) ;
    PIX_SORT(p[1], p[19]) ;  PIX_SORT(p[1], p[10]) ;  PIX_SORT(p[11], p[20]) ;
    PIX_SORT(p[2], p[20]) ;  PIX_SORT(p[2], p[11]) ;  PIX_SORT(p[12], p[21]) ;
    PIX_SORT(p[3], p[21]) ;  PIX_SORT(p[3], p[12]) ;  PIX_SORT(p[13], p[22]) ;
    PIX_SORT(p[4], p[22]) ;  PIX_SORT(p[4], p[13]) ;  PIX_SORT(p[14], p[23]) ;
    PIX_SORT(p[5], p[23]) ;  PIX_SORT(p[5], p[14]) ;  PIX_SORT(p[15], p[24]) ;
    PIX_SORT(p[6], p[24]) ;  PIX_SORT(p[6], p[15]) ;  PIX_SORT(p[7], p[16]) ;
    PIX_SORT(p[7], p[19]) ;  PIX_SORT(p[13], p[21]) ; PIX_SORT(p[15], p[23]) ;
    PIX_SORT(p[7], p[13]) ;  PIX_SORT(p[7], p[15]) ;  PIX_SORT(p[1], p[9]) ;
    PIX_SORT(p[3], p[11]) ;  PIX_SORT(p[5], p[17]) ;  PIX_SORT(p[11], p[17]) ;
    PIX_SORT(p[9], p[17]) ;  PIX_SORT(p[4], p[10]) ;  PIX_SORT(p[6], p[12]) ;
    PIX_SORT(p[7], p[14]) ;  PIX_SORT(p[4], p[6]) ;   PIX_SORT(p[4], p[7]) ;
    PIX_SORT(p[12], p[14]) ; PIX_SORT(p[10], p[14]) ; PIX_SORT(p[6], p[7]) ;
    PIX_SORT(p[10], p[12]) ; PIX_SORT(p[6], p[10]) ;  PIX_SORT(p[6], p[17]) ;
    PIX_SORT(p[12], p[17]) ; PIX_SORT(p[7], p[17]) ;  PIX_SORT(p[7], p[10]) ;
    PIX_SORT(p[12], p[18]) ; PIX_SORT(p[7], p[12]) ;  PIX_SORT(p[10], p[18]) ;
    PIX_SORT(p[12], p[20]) ; PIX_SORT(p[10], p[20]) ; PIX_SORT(p[10], p[12]) ;
    return (p[12]);
}

/*---------------------------------------------------------------------------
   Function :   kth_smallest()
   In       :   array of elements, # of elements in the array, rank k
   Out      :   one element
   Job      :   find the kth smallest element in the array
   Notice   :   use the median() macro defined below to get the median. 
                Reference:
                  Author: Wirth, Niklaus 
                   Title: Algorithms + data structures = programs 
               Publisher: Englewood Cliffs: Prentice-Hall, 1976 
    Physical description: 366 p. 
                  Series: Prentice-Hall Series in Automatic Computation 
 ---------------------------------------------------------------------------*/
#define SWAP(a, b) {BYTE t = (a); (a) = (b); (b) = t; }
BYTE kth_smallest(BYTE a[], int n, int k) {
    int l = 0, m = n - 1 ;
    while (l < m) {
        BYTE x = a[k] ;
        int i = l, j = m ;
        do {
            while (a[i] < x) i++ ;
            while (x < a[j]) j-- ;
            if (i <= j) {
                SWAP(a[i], a[j]) ;
                i++ ; j-- ;
            }
        } while (i <= j) ;
        if (j < k) l = i ;
        if (k < i) m = j ;
    }
    return a[k] ;
}
#define median(a, n) kth_smallest(a, n, (((n) & 1) ? ((n) / 2) : (((n) / 2) - 1)))

지문에서 ridge의 방향(orientation)은 gradient에 수직한 방향이다(그런데 벡터인 gradient와는 달리 ridge의 방향은 모호함이 있다. 시계방향 또는 반시계방향으로 90도 회전이 모두 동일한 ridge의 방향이다).  gradient의 방향각이 $\alpha$일 때 수직인 ridge의 방향각은 $\theta =\frac{\pi}{2} + \alpha$로 정의하자. 방향각을 주어진 gradient 성분 $\nabla I = (g_x, g_y)$로 표현하기 위해서

$$ \sin 2\alpha = 2 \sin \alpha \cos \alpha = 2 \frac{g_x}{\sqrt{g_x^2 + g_y^2}}\frac{g_y}{\sqrt{g_x^2 + g_y^2}}$$

$$\cos 2\alpha = \cos^2 \alpha - \sin^2\alpha = \frac{g_x^2}{g_x^2 + g_y^2} -\frac{g_y^2}{g_x^2 +g_y^2}$$ 

임을 이용하자. $$\tan 2\alpha = \frac{ 2g_x g_y}{ g_x^2 - g_y^2}$$로 쓸 수 있으므로 ridge의 방향각은

$$ \theta =\frac{\pi}{2} + \frac{1}{2} \arctan \frac{2g_x g_y}{g_x^2 - g_y^2}$$

임을 알 수 있다. $(g_x, g_y) \rightarrow (-g_x, -g_y)$이더라도 ridge의 방향은 동일하므로 $\theta$의 범위는 $0\le \theta < \pi$로 제한된다.

이미지의 한 지점에서 ridge의 방향은 그 점을 중심으로 하는 적당한 크기의 윈도 평균값으로 추정한다. gradient 성분은 Sobel 연산자나 Prewitt 연산자를 convolution해서 얻을 수 있다. 따라서 지문 이미지 상의 픽셀 좌표 $(i,j)$에서 ridge의 방향은

$$\theta_{ij} = \frac{\pi}{2} +\frac{1}{2}\arctan \frac{2 G_{xy}}{G_{xx} - G_{yy}}$$

$$ G_{xy} = \sum_{W} \text{gradX}_{ij} \times \text{gradY}_{ij}$$

$$ G_{xx} = \sum_{W} \text{gradX}_{ij} \times \text{gradX}_{ij}$$

$$ G_{yy} = \sum_{W} \text{gradY}_{ij} \times \text{gradY}_{ij}$$

아래 결과는 5x5 크기의 Gaussian filter를 적용한 지문 이미지에 17x17 크기의 윈도를 겹치지 않게 잡아서 평균을 계산한 결과다. $G_{xx}+G_{yy}$가 최대값의 $20\%$보다 작은 영역은 방향을 표시하지 않았다 (red:0~45, yellow:45~90, magenta: 90~135, cyan: 135~180).

int RidgeOrientaion(BYTE *img, int w, int h, int bsz/*17*/, double *omap) {
    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};
    const int xmax = w - 1, ymax = h - 1;
    std::vector<int> gradx(w * h, 0);
    std::vector<int> grady(w * h, 0);
    for (int y = 1, pos = w; y < ymax; y++) {
        pos++; //x = 0;
        for (int x = 1; x < xmax; x++, pos++) {
            int sx = 0, sy = 0; 
            for (int k = 0; k < 9; k++) {
                int v = img[pos + nn[k]];
                sx += sobelX[k] * v;
                sy += sobelY[k] * v;
            }
            gradx[pos] = sx;
            grady[pos] = sy;
        }
        pos++; // x = xmax;
    }
    const int ny = h / bsz;
    const int nx = w / bsz;
    for (int iy = 0; iy < ny; iy++) {
        for (int ix = 0; ix < nx; ix++) {
            int sxy = 0, sxx = 0, syy = 0; 
            int pos = (iy * w + ix) * bsz; //starting position of block(ix, iy)
            int *dx = &gradx[pos];
            int *dy = &grady[pos];
            for (int jj = 0; jj < bsz; jj++) {
                for (int ii = 0; ii < bsz; ii++) {
                    int gx = dx[jj * w + ii];
                    int gy = dy[jj * w + ii];
                    sxy += gx * gy;
                    sxx += gx * gx;
                    syy += gy * gy;
                }
            }
            omap[iy * nx + ix] = 0.5 * (PI + atan2(double(2.0 * sxy), double(sxx - syy)));
        }
    }
    //...identify fingerprint regions and draw orientation map;
    return 1;
}

더보기

Neighborhood Map;

[0 1 2]

[7 8 3]

[6 5 4]

int Thinning_2pass(BYTE *image, int w, int h) {
    const int xmax = w - 1, ymax = h - 1;
    const int nn[9] = {-w - 1,- w, -w + 1, 1, w + 1, w, w - 1, -1, 0};//clockwise;
    const BYTE FG = 255, BG = 0;
    bool *flag = new bool [w * h];
    int pass = 0, ok = 0;
    int nb[9];
    while (!ok) {
        ok = 1; pass = (pass + 1) % 2;
        for (int i = w * h; i-->0; ) flag[i] = false;
        for (int y = 1, pos = w; y < ymax; y++) {
            pos++;//x=0;
            for (int x = 1; x < xmax; x++, pos++) {
                if (image[pos] == FG) { //fg;
                    // condition 1; 
                    int count = 0;
                    for (int k = 0; k < 8; k++) 
                        if (image[pos + nn[k]] == FG) count++;
                    if (count >= 2 && count <= 6) {
                        for (int k = 0; k < 8; k++) nb[k] = image[pos + nn[k]];
                        nb[8] = nb[0]; //cyclic;
                        // condition 2;
                        int trans = 0;
                        for (int k = 0; k < 8; k++)
                            if (nb[k] == BG && nb[k + 1] == FG) trans++;
                        if (trans == 1) {
                            // condition3: top&&left=bg || bot=bg || right=bg
                            if (pass == 0 && (nb[3] == BG || nb[5] == BG ||	
                                (nb[1] == BG && nb[7] == BG))) {
                                    flag[pos] = true; ok = 0;
                            } else { // condition4: bot&&right=bg || top=bg || left=bg 
                                if (pass == 1 && (nb[1] == BG || nb[7] == BG ||	
                                    (nb[3] == BG && nb[5] == BG))) {
                                        flag[pos] = true; ok = 0;
                                }
                            }
                        }
                    }//(2<=count<=6);
                }
            }//for_x;
            pos++;//x = w - 1 skip;
        } //for_y;
        // remove flaged pixels;
        for (int y = 1, pos = w; y < ymax; y++) {
            pos++;//x = 0;
            for (int x = 1; x < xmax; x++, pos++) 
                if (flag[pos]) image[pos] = BG;
            pos++; //x=w-1;
        }
    }
    delete [] flag;
    return 1;
}

Contrast Limited Adaptive Histogram Equalization (CLAHE). CLAHE는 영상의 평탄화 (histogram equalization) 과정에서 contrast가 과도하게 증폭이 되는 것을 방지하도록 설계된 adaptive algorithm의 일종이다. CLAHE algorithm은 전체 영상 영역을 바로 적용하지 않고 tile이라 불리는 작은 영역으로 영상을 분할하여 tile 별로 적용한다. 그리고 인접하는 tile 간에 bilinear 보간을 사용함으로써 tile 경계에서 급격한 변화가 생기는 것을 방지한다. gray 영상에 CLAHE 알고리즘을 적용하면 영상의 contrast을 개선해 준다. 그리고 컬러 영상에도 적용을 할 수 있는데, 보통 luminance 채널에만 적용한 결과(HSV로 변환 후 V에 대해서만 적용)가 RGB 각 채널별로 적용하는 것보다 좋은 결과를 준다.

아래는 RGB 영상의 각 채널에 CLAHE 알고리즘을 적용한 결과이다.

rgn: 히스토그램을 구하는 이미지 상의 영역;

tile: adaptive HE가 적용되는 영역

첫 행/열 tile의 폭/높이 = rgn 폭/높이의 절반

마지막 tile 크기 = 이미지의 나머지 폭/높이 

cliplimit = 각 영역의 히스토그램 슬롯당 들어가는 최대 픽셀 수

/* "Graphics Gems IV" 의 code를 수정함;
** 2023-01-13 updated: 이미지 사이즈가 rgn 사이즈의 정수배가 되지 않을 때 문제 해결;*/
int CLAHE ( BYTE* image, int width, int height,
            BYTE gMin/*0*/, BYTE gMax/*255*/, int rgn_nx/*16*/, int rgn_ny/*16*/,
            int bins/*256*/, double fcliplimit )
{
    std::vector<int> map_array(rgn_nx * rgn_ny * bins);
    int rgn_xsz		= width / rgn_nx;
    int rgn_ysz		= height / rgn_ny; 
    int rgn_area	= rgn_xsz * rgn_ysz;
    int clipLimit	= int( fcliplimit * ( rgn_xsz * rgn_ysz ) / bins );
    clipLimit = ( clipLimit < 1 ) ? 1 : clipLimit;
    //  
    /* calculate greylevel mappings for each contextual region */
    for (int iy = 0; iy < rgn_ny; iy++) {
        for (int ix = 0; ix < rgn_nx; ix++) {
            int *hist = &map_array[bins * ( iy * rgn_nx + ix )];
            int start = (iy * rgn_ysz) * width + ix * rgn_xsz;
            RgnHistogram (&image[start], width, rgn_xsz, rgn_ysz, hist, bins);
            ClipHistogram ( hist, bins, clipLimit );
            //convert hist to cuimulative histogram normalized to [gMin=0, gMax=255]
            MapHistogram ( hist, gMin, gMax, bins, rgn_area );
        }
    }

    /* bilinearInterp greylevel mappings to get CLAHE image */
    /* 첫열/행에 해당하는 타일은 rgn의 절반 크기만 처리하므로 ix(iy)는 0~rgn_nx(rgn_ny)까지임*/
    int szx, szy, ixl, ixr, iyt, iyb;
    BYTE *first_tile = &image[0]; // 같은 iy를 갖는 첫 타일의 주소;
    for (int iy = 0; iy <= rgn_ny; iy++ ) {
        if ( iy == 0 ) {				
            // 첫 타일은 절반만;
            szy = rgn_ysz / 2; iyt = iyb = 0;
        } else if ( iy == rgn_ny ) {
            // 마지막 타일은 나머지 전부;
            szy = height - ((rgn_ny-1)*rgn_ysz + rgn_ysz/2);
            iyt = rgn_ny - 1; iyb = iyt;
        } else {
            szy = rgn_ysz; iyt = iy - 1; iyb = iyt + 1;
        }
        BYTE *ptile = &first_tile[0]; //각 타일의 시작 주소;
        for (int ix = 0; ix <= rgn_nx; ix++ ) {
            if ( ix == 0 ) {	
                //첫 타일은 절반만;
                szx = rgn_xsz / 2; ixl = ixr = 0;
            } else if ( ix == rgn_nx ) {		
                //마지막 타일은 나머지 전부;
                szx = width - ((rgn_nx-1)*rgn_xsz + rgn_xsz/2);
                ixl = rgn_nx - 1; ixr = ixl;
            } else {
                szx = rgn_xsz; ixl = ix - 1; ixr = ixl + 1;
            }
            // cumulative histogram data;
            int *LU = &map_array[bins * ( iyt * rgn_nx + ixl )];
            int *RU = &map_array[bins * ( iyt * rgn_nx + ixr )];
            int *LB = &map_array[bins * ( iyb * rgn_nx + ixl )];
            int *RB = &map_array[bins * ( iyb * rgn_nx + ixr )];
            BilinearInterp (ptile, width, LU, RU, LB, RB, szx, szy );
            ptile += szx;			  
        }
        first_tile += szy * width;
    }
    return 0;						 
}

void ClipHistogram ( int* hist, int gLevels, int clipLimit ) {
    int excess = 0;
    for (int i = 0; i < gLevels; i++ ) { /* calculate total number of excess pixels */
        int diff =  hist[i] - clipLimit;
        if ( diff > 0 ) excess += diff;	 /* excess in current bin */
    };

    /* clip histogram and redistribute excess pixels in each bin */
    int incr = excess / gLevels;		/* average binincrement */
    int upper =  clipLimit - incr;		/* bins larger than upper set to cliplimit */
    for (int i = 0; i < gLevels; i++ ) {
        if ( hist[i] > clipLimit ) hist[i] = clipLimit; /* clip bin */
        else {
            if ( hist[i] > upper ) {			/* high bin count */
                excess -= hist[i] - upper;
                hist[i] = clipLimit;
            } else {							/* low bin count */
                excess -= incr;
                hist[i] += incr;
            }
        }
    }

    while ( excess ) { /* redistribute remaining excess  */
        int start = 0;
        while ( excess && start < gLevels) {
            int step = gLevels / excess;
            if ( step < 1 ) step = 1;		 /* step size at least 1 */
            for (int i = start; i < gLevels && excess; i += step) 
                if (hist[i] < clipLimit) { 
                    hist[i]++; excess--;
                }
            start++;
        }
    }
}
void BilinearInterp ( BYTE * image, int stride, 
                    int* LU, int* RU, int* LB,  int* RB,
                    int rgn_xsz, int rgn_ysz ) {
    int area = rgn_xsz * rgn_ysz; /* Normalization factor */
    for (int y = 0, yrev = rgn_ysz; y < rgn_ysz; y++, yrev--) {
        BYTE *line = &image[y * stride];
        for (int x = 0, xrev = rgn_xsz; x < rgn_xsz; x++, xrev-- ) {
            int v = line[x]; /* get histogram bin value */
            v = (yrev * ( xrev * LU[v] + x * RU[v] ) + y * ( xrev * LB[v] + x * RB[v] ) ) / area;
            line[x] = v > 255 ? 0xFF: v < 0 ? 0x00: v;
        }
        //image += stride - rgn_xsz;  // goto the beginning of rgn's next-line;
    }
}
void RgnHistogram ( BYTE* image, int stride, int rgn_xsz, int rgn_ysz, 
                        int* hist, int gLevels) {
    // image is pointing the beginnign of a given region;
    for (int i = 0; i < gLevels; i++ ) hist[i] = 0;
    for (int y = 0; y < rgn_ysz; y++ ) {
        BYTE *line = &image[y * stride];
        for (int x = 0; x < rgn_xsz; x++)
            hist[line[x]]++;
    }
}
// transform hist[] to a normalized cumulative histogram;
void MapHistogram ( int* hist, BYTE gMin, BYTE gMax, int gLevels, int rgn_area ){
    double scale = double( gMax - gMin ) / rgn_area;
    int iMin = int(gMin);
    for (int i = 0, sum = 0; i < gLevels; i++ ) {
        sum += hist[i];
        int a = int(iMin + sum * scale );
        hist[i] = a > gMax ? gMax: a < gMin ? gMin: a;
    }
}

n이 양수일 때 n 과 n-1은 짝홀/홀짝의 짝이므로 2의 파워가 아니면 이진수로 표현할 때 한 곳에서 비트 차이가 난다. 따라서 AND연산을 하면 0이 아니다.  그런데 n이 2의 파워이면 최상위 비트만 1이고, n-1은 n의 최상위 비트를 제외한 하위비트 모두가 1로 구성되므로 AND 연산을 하면 0이다.

예:   n=6 = 110,  6-1=101 -->  6 & 5 = 100;

       n=8 = 1000, 8-1=0111 --> 8 & 7=0000

// n이 양수이고, n-1과 겹치는 비트가 없으면 된다;
bool IsPowerOf2(int n) {
    return n > 0 && (n & (n - 1)) == 0;
}

n>1 일 떄 2로 나누기를 계속할 때 1을 제외한 홀수가 나오면 안됨; 

bool IsPowerOf2(int n) {
    if (n == 1) return true;
    if (n == 0 || n & 1) return false;
    return IsPowerOf2(n >> 1);
} // n == 0;

x보다 크기 않은 2^i을 구한 후 x와 같은가 체크;

bool IsPowerof2(int x){
    int i = 0;
    while ((1 << i) < x) i++;
    if (x == (1 << i)) return true;
    return false;
}

2의 파워이면 log_2(n) = 자연수 이므로 ceil과 floor가 같다.

bool IsPowerOf2(int n) {
   if (n == 0) return false;
   return (ceil(log2(n)) == floor(log2(n)));
}

주어진 자연수보다 작지 않은 최소의 2^n;

int NextPowerOf2(int n) { //32-bit;
    n--;
    n |= n >> 1;
    n |= n >> 2;
    n |= n >> 4;
    n |= n >> 8;
    n |= n >> 16;
    n++;
    return n;    
} // NextPowerOf2(5) -> 8; NextPowerOf2(8) -> 8;