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

RGB 컬러 이미지의 gray level의 cdf를 이용해서 histogram equalization을 수행한다. 컬러 이미지의 gray level은 

$$gray = \frac{r + g+ b}{3}$$

으로 정의한다.

std::vector<BYTE> color_equalize_new(std::vector<BYTE> &rgb) {
    int hist[256] = {0};
    int chist[256], lut[256], partition[256 + 1];
    std::vector<BYTE> out = raster; // clone; 
    // pdf of gray level: g = (r + g + b) / 3
    for (int k = 0; k < rgb.size(); k += 3) {
        int a = rgb[k+0]; a += rgb[k+1]; a += rgb[k+2];
        ++hist[a/3];
    };
    // cdf;
    for (int i = 0, s = 0; i < 256; i++) {
        s += hist[i]; chist[i] = s;
    }
    
    /* assign equal number of pixels in each gray levels;*/
    int num_pixels = rgb.size()/3;
    double pixels_per_level = double(num_pixels) / 256;

    /* first and last in partition */
    partition[0]   = 0;
    partition[256] = 256;
    /* intermediate; */
    for (int j = 0, i = 1; i < 256; i++) {
        double desired = i * pixels_per_level;			
        while (chist[j + 1] <= desired) j++;
        /* chist[j] <= desired < chist[j+1];*/
        /* nearest interpolation */
        if ((desired - chist[j]) < (chist[j+1] - desired)) partition[i] = j;
        else partition[i] = j + 1;
    } 
    /* Find i s.t. partion[i] <= j < partition[i+1];*/
    for (int j = 0; j < 256; j++) {
        int i = 0;
        while (partition[i+1] <= j) i++;
        lut[j] = i;
    } 
    // needs hue preserving processes;
    for (k = rgb.size(); k-->0;)
        res[k] = lut[src[k]];
    return out;
}

std::vector<BYTE> color_equalize_HSV(std::vector<BYTE> &rgb) {
    int hist[256] = {0};
    int chist[256] = {0}, lut[256] = {0}, partition[256 + 1];
    std::vector<BYTE> out = raster; // cloning;
    std::vector<double> fH(rgb.size()/3);
    std::vector<double> fV(rgb.size()/3);
    std::vector<double> fS(rgb.size()/3);

    const int n = rgb.size()/3;
    double r, g, b;
    for (int k = 0; k < rgb.size(); k += 3) {
        b = rgb[k+0], g = rgb[k+1], r = rgb[k+2];
        RGBtoHSV(r / 255, g / 255, b / 255, fH[k], fS[k], fV[k]);
    };	
    // make histogram of V-color;
    for (int k = n; k-->0;)
        ++hist[int(fV[k] * 255)];

    // cdf;
    for (int i = 0, s = 0; i < 256; i++) {
        s += hist[i]; chist[i] = s;
    }
    /* assign equal number of pixels in each V-color;*/
    int num_pixels = rgb.size()/3;
    double pixels_per_level = double(num_pixels) / 256;

    /* first and last in partition */
    partition[0]   = 0;
    partition[256] = 256;
    /* intermediate; */
    for (int j = 0, i = 1; i < 256; i++) {
        double desired = i * pixels_per_level;			
        while (chist[j + 1] <= desired) j++;
        /* nearest interpolation */
        if ((desired - chist[j]) < (chist[j + 1] - desired)) partition[i] = j;
        else partition[i] = j + 1;
    } 
    /* fill equalization LUT */
    for (int j = 0; j < 256; j++) {
        int i = 0;
        while (partition[i + 1] <= j) i++;
        lut[j] = i;
    } 
    for (int k = n; k-->0;)
        fV[k]= lut[int(fV[k] * 255)] / 255.;

    for (int k = 0; k < rgb.size(); k += 3) {
        HSVtoRGB(fH[k], fS[k], fV[k], r, g, b);
        out[k+0] = int(b * 255); 
        out[k+1] = int(g * 255);
        out[k+2] = int(r * 255);
    }
    return out;
}

/* fR Red component, range: [0, 1]
** fG Green component, range: [0, 1]
** fB Blue component, range: [0, 1]
** fH Hue component, range: [0, 360]
** fS Hue component, range: [0, 1]
** fV Hue component, range: [0, 1] */
void RGBtoHSV(double fR, double fG, double fB, double& fH, double& fS, double& fV);
void HSVtoRGB(double fH, double fS, double fV, double& fR, double& fG, double& fB);

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

서로 다른 조건에서 생성이 된 두 개의 영상을 비교하기 위해서는 먼저 영상을 정규화시키어야 한다. 쉽게 생각할 수 있는 정규화의 방법은 두 개의 영상이 동일한 히스토그램을 같도록 만들면 된다. 어떻게 동일한 히스토그램을 만들 수 있는가? 간단히 생각할 수 있는 방법은 주어진 영상의 히스토그램을 전 픽셀에 걸쳐서 균일한 빈도를 보이는 히스토그램이 되도록 변환을 하는 것이다. 즉, 히스토그램을 평탄화시키는 것이다. 이 변환을 $F(x)$라고 하면, 입력 영상의 픽셀 값 $x$는 새로운 픽셀 값 $y = F(x)$로 바뀐다. 영상의 히스토그램은 픽셀 값의 빈도를 나타내므로 확률 밀도 함수로 해석할 수 있고 (픽셀 값을 연속 변수로 볼 때) 주어진 영상의 확률 밀도 함수를 $P(x)$라고 하자. 그러면 $P(x) dx$는 픽셀 값이 $[x, x+dx]$인 픽셀이 영상에서 나타날 확률을 의미한다. 변환 후에서는 임의의 픽셀 값이 균일해야 하므로 확률 밀도 함수는 상수 $a$이고, 픽셀 값의 구간이 $[0,255]$ 이므로 $a = 1/ 255$로 주어진다. 따라서, 변환이 확률을 보존하기 위해서는 $$P(x) dx = a dy$$을 만족해야 하므로 변환 후 픽셀 값 y와 변환 전 픽셀 값 x사이의 관계를 알 수 있다. $$\frac {dy}{ dx} = \frac {1}{a}  P(x) $$ 이 식을 적분하면, $$y = F(x) = \frac {1}{a} \int_0^x P(x') dx' =  255\times  \text {cumulative sum of } P(x)$$ 즉, 히스토그램을 평탄화시키는 변환은 원래 영상 히스토그램의 누적합을 전체 픽셀 수로 나눈 값,

$$k \rightarrow F(k)=255 \times \frac {\sum_{i=0}^{k} H [i] - H[0]}{ \sum_{i=0}^{255} H [i] - H[0]}$$으로 결정된다($0\to 0$이도록 $H[0]$을 뺀다). 영상의 히스토그램은 이산적이기 때문에 실제로 모든 픽셀 값에서 균일하게 나타나지 않지만, 구간 변환은 거의 일정하게 나타난다. 이 변환을 거친 영상의 누적 히스토그램(cumulative histogram)을 $y$-축 픽셀 값을 $x$-축으로 해서 그리면 거의 일직선으로 증가하는 모양을 나타낸다.

이 히스토그램의 평탄화는 픽셀 값이 좁은 영역에 몰려있는 영상을 전 영역에 분포하도록 한다. 인간은 영상의 밝기에 의해서 보다는 밝고 어두움의 변화의 크기에 의해서 인지도가 증가하게 되는데, 평탄화된 영상은 보다 쉽게 인식이 될 수 있는 구조를 가진다.

std::vector<BYTE> histogram_eq(const std::vector<BYTE> &src) {
    int hist[256] = {0}, cum[256], map[256];
    for (int k = src.size(); k-->0;) ++hist[src[k]];
    for (int k = 0, s = 0; k < 256; k++) {
    	s += hist[k]; cum[k] = s;
    }
    double factor = double(255) / (cum[255] - cum[0]) ;
    for (int k = 0; k < 256; k++) {
        int x = int(factor * (cum[k] - cum[0]));
        map[k] = x > 255 ? 255: x;
    }
    std::vector<BYTE> dst(src.size());
    for (int k = src.size(); k-->0;) dst[k] = map[src[k]];
    return dst;
}

동일한 물체나 풍경이라도 촬영에 사용된 센서의 종류나 조명 조건 때문에 사진이 더 어둡거나 반대로 더 밝게 표현될 수 있다. 이런 경우 사진을 잘 찍힌 사진과 비교하여 보정하려면 어떤 기준을 가지고 작업을 해야 할까? 동일한 내용의 영상이라면 단순히 히스토그램을 늘리거나 줄이는 과정을 (histogram stretching) 사용할 수 있다. 그러나 같은 물체나 풍경을 담지 않은 기준 영상의 히스토그램과 같은 내용의 히스토그램을 되도록 사진을 조작하고 싶으면 histogram stretching으로는 안된다. 

입력 영상의 주어진 픽셀 값이 기준 영상의 어떤 픽셀 값에 해당하게 변환을 해야 하는가를 결정해야 한다. 변환 $f$는 두 영상의 픽셀 값의 확률분포를 보존하도록 만들어진다. 입력 영상의 히스토그램 구간 $(i-1, i]$이 기준 영상의 히스토그램 구간 $(f(i-1), f(i)]$으로 변환될 때 각각 구간에서 픽셀이 나타날 확률이 동일해야 한다. 영상에서 픽셀 값은 이산적이므로 확률보다는 확률의 누적합을 이용하는 것이 더 편리하다. 입력 영상의 히스토그램 구간 $[0, i]$에서 픽셀이 나타날 누적 확률과 같은 누적 확률을 주는 기준 영상의 히스토그램 구간이 $[0, k=f(i)]$일 때 

$$i(영상)\rightarrow k = f(i) (기준영상)$$

으로 변환하면 된다. 이 과정은 히스토그램 stretching에도 그대로 적용이 된다. 누적 히스토그램(cumulative histogram)은 주어진 픽셀 값까지 확률 합을 쉽게 구할 수 있게 하는 일종의 lookup table이다. cumulative histogram의 $i$-번째 값은 픽셀 값이 $0$부터 $i$까지 히스토그램의 누적합이기 때문이다. 따라서 히스토그램 매칭(histogram matching, histogram specification)은 같은 cumulative histogram의 값을 갖는 픽셀 값 사이의 변환을 의미한다. 두  영상의 픽셀 수가 같지 않은 경우에는 cumulative histogram를 쓸 수 없다. 이 경우에는 전체 픽셀 수로 정규화시킨 cumulative distribution function(cdf)을 사용하면 된다.

히스토그램 평탄화(histogram equalization)는 각 그레이 값이 나타날 확률이 균일한 기준 영상을 (즉, 1차 cumulative histogram이 직선인 영상) 사용한 히스토그램 매칭에 해당한다. 따라서 기준 영상이 필요하지 않는다. 

//map을 lookup table로 이용해서 영상을 변환할 수 있다.
void matchHistogram(int srcHist[256], int refHist[256], int map[256]) {
    double cumSrc[256], cumRef[256];
    cumSrc[0] = srcHist[0];
    cumRef[0] = refHist[0];
    for (int i = 1; i < 256; i++) {
        cumSrc[i] = cumSrc[i - 1] + srcHist[i];
        cumRef[i] = cumRef[i - 1] + refHist[i];
    }
    // normalize ==> make CDF; 두 비교 영상이 같은 크기가 아닐 수 있으므로 필요하다;
    for (int i = 0; i < 256; i++) {
        cumSrc[i] /= cumSrc[255];
        cumRef[i] /= cumRef[255];
    }
    // 두 비교영상의 cumulant histogram이 가장 가까운(같은) 픽셀끼리 변환한다;
    for (int i = 0; i < 256; i++) {
        int k = 255;  
        while (cumRef[k] > cumSrc[i]) k--;
        // nearest interpolation; 2022.02.11
        if ((cumSrc[i] - cumRef[k]) < (cumRef[k + 1] - cumSrc[i])) map[i] = k;
        else map[i] = k + 1;
    }
};

실제 사용은 grey 이미지인 경우:

int histogramSpec(BYTE *src, int ws, int hs, BYTE *ref, int wr, int hr, BYTE *out) {
    int srcHist[256] = {0}, refHist[256] = {0}, map[256];
    makeHistogram(src, ws, hs, srcHist);
    makeHistogram(ref, wr, hr, refHist);
    matchHistogram(srcHist, refHist, map);
    for (int k = ws * hs; k-- > 0;)
        out[k] = map[src[k]];
    return 1;
}