Geometry & Recognition

픽셀값의 분포가 좁은 구간에 한정된 영상은 낮은 명암 대비(contrast)를 가지게 된다. 이 경우 픽셀값을 가능한 넓은 구간에 재분포하도록 변환시키면 전체적으로 훨씬 높은 contrast를 가지는 영상을 얻을 수 있는데 이러한 기법을 histogram equalization이라고 한다. Histogram equalization은 픽셀 분포의 cdf(cumulative distribution function)가 가능한 직선이 되도록 변형을 시킨다. 그러나 어떤 경우에는 과도한 contrast로 인해서 영상의 질이 더 안 좋아지는 경우가 발생하거나 영상의 일부 영역에서는 오히려 contrast가 감소하는 경향이 발생할 수 있다. 이를 개선하기 위해는 cdf의 과도한 변형을 막고, equaliztion을 국소적으로 다르게 적용할 수 있는 adaptive 알고리즘을 고려해야 한다. 픽셀값의 분포를 연속적 변수로 생각하면 히스토그램 $h(x)$은 픽셀값의 pdf(probability density function)를 형성하고, 히스토그램의 최댓값과 최솟값을 각각 $m=\text{min}(h)$, $M=\text{max} (h)$라면 

$$ m \le \frac{1}{255}\int_0^{255} h(x) dx \le M$$

Histogram equalization은 $m \approx M$이 되도록 $h(x)$을 변환을 시키는 과정에 해당한다.  Histogram 함수는 cdf함수의 미분값이므로 cdf 곡선의 접선의 기울기에 해당하고, histogram equation은 cdf 함수의 기울기가 모든 픽셀값에서 일정하게 조정하는 과정이 된다. 과도한 contrast를 갖는 영상을 만들지 않기 위해서 일정한 cdf의 기울기를 요구하지 않고 일정범위를 넘어서는 경우만 제한을 두도록 하자. 제한값은 cdf의 평균변화율

$$\overline{s}=\frac{1}{255}\int_0^{255} h(x) dx$$

을 기준으로 선택하면 된다.  제한 기준으로 $slope \times \overline{s}$를 선택하면 변형된 histogram은 

$$ h_\text{modified} (x) = \left\{  \begin{matrix}  slope \times \overline{s} & \text{if} ~~h(x) \ge slope \times \overline{s} \\h(x) & \text{otherwise} \end{matrix}   \right.$$

확률을 보존하기 위해 기준을 넘는 픽셀 분포는 초과분을 histogram의 모든 bin에 고르게 재분배를 하는 과정이 필요하다.

$$h_\text{modified}(x) \underset{\text{redistribution of excess}}{\Longrightarrow}\tilde{h}_\text{modified}(x)$$

이 새로이 변환된 histogram에 대해 equalization을 적용하여 윈도 중심 픽셀의 값을 설정한다. 이 알고리즘을 contrast limited histogram equalization(CLAHE)라고 부른다. 보통 CLAHE 알고리즘에서는 각각의 픽셀에 대해 윈도를 선택하지 않고, 원본 영상을 겹치는 타일영역으로 분할한 후 각각의 타일에서 제한된 histogram equalization을 변환을 구한다. 그리고 겹치는 인접 타일에서 변환의 선형보간을 이용하면 타일과 타일 사이에서 계단현상을 막을 수 있다. 여기서는 각각의 픽셀에 대해서 슬라이딩 윈도를 설정하는 방법을 사용한다. 윈도 크기는 영상에서 관심 대상을 충분히 커버할 정도의 크기여야 한다. 또한 $slope$은 보통 $3\sim 4$정도를 선택한다. $slope\to 1$인 경우는 원본 영상과 거의 변화가 없고($h(x) \approx \overline{s}$), $slope\gg 1$이면 재분배되는 픽셀이 없으므로 마지막 단계에서 equalization이 global equalization과 동일하게 된다.

static const int halfwindow	= 63;
static const double slope	= 3;
void clahe2(BYTE **input, int width, int height, BYTE **output) {
    const int bins = 255;
    int hist[bins+1], clippedHist[bins+1];
    for (int y = 0; y < height; ++y) {
        int top = max(0, y - halfwindow);
        int bottom = min(height-1, y + halfwindow);
        int h = bottom - top + 1;

        /* 픽셀 access를 줄이기 위해 sliding 윈도우 기법을 사용함;*/
        memset(hist, 0, sizeof(hist));
        int right0 = min(width-1, halfwindow);	// x=0일 때 오른쪽;
        for (int yi = top; yi <= bottom; ++yi)
            for (int xi = 0; xi < right0; ++xi) // x=right는 다음 step에서 더해짐;
                ++hist[input[yi][xi]];

        for (int x = 0; x < width; ++x) {
            int left = max(0, x - halfwindow );
            int right = x + halfwindow;
            int w = min(width-1, right) - left + 1;
            int npixs = h * w;	// =number of pixels inside the sliding window;

            int limit = int(slope * npixs / bins + 0.5);  
            // slope >>1 -->hist_eq;
            /* 윈도우가 1픽셀 이동하면 왼쪽 열은 제거; */
            if (left > 0)
                for (int yi = top; yi <= bottom; ++yi)
                    --hist[input[yi][left-1]];						

            /* 윈도우가 움직이면 오른쪽 열을 추가함; */
            if (right < width)
                for (int yi = top; yi <= bottom; ++yi)
                    ++hist[input[yi][right]];						

            /* clip histogram and redistribute excess; */
            memcpy(clippedHist, hist, sizeof(hist));
            int excess = 0, excessBefore;
            do {
                excessBefore = excess;
                excess = 0;
                for (int i = 0; i <= bins; ++i) {
                    int over = clippedHist[i] - limit;
                    if (over > 0) {
                        excess += over;
                        clippedHist[i] = limit;
                    }
                }

                int avgExcess = excess / (bins + 1);
                int residual  = excess % (bins + 1);// 각 bin에 분배하고 남은 나머지;
                for (int i = 0; i <= bins; ++i)	// 먼저 전구간에 avgExcess를 분배;
                    clippedHist[i] += avgExcess;

                if (residual != 0) {
                    int step = bins / residual;	// 나머지는 일정한 간격마다 1씩 분배;
                    for (int i = 0; i <= bins; i += step)
                        ++clippedHist[i];
                }
            } while (excess != excessBefore);

            /* clipped histogram의 cdf 구성;*/
            int hMin = bins;
            for (int i = 0; i < hMin; ++i)
                if (clippedHist[i] != 0) hMin = i;
            int cdfMin = clippedHist[hMin];
            
            int v = input[y][x]; //현 위치의 픽셀 값;
            int cdf = 0;
            for (int i = hMin; i <= v; ++i) // cdf(현 픽셀);
                cdf += clippedHist[i];

            int cdfMax = cdf;
            for (int i = v + 1; i <= bins; ++i) // 총 픽셀;
                cdfMax += clippedHist[i];

            // 현픽셀 윈도의 clipped histogram을 구하고, 
            // 이 clipped histogram을 써서 equaliztion을 함;
            output[y][x] = int((cdf - cdfMin)/double(cdfMax - cdfMin) * 255);
        }
    }
}

https://kipl.tistory.com/291

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;
    }
}

https://kipl.tistory.com/608