히스토그램은 영상에서 각 그레이 값에 해당하는 픽셀의 수를 주는 일종의 이산적 함수로 생각할 수 있다. histogram에서 피크의 위치는 histogram을 연속적인 함수로 모델링하거나 또는 여러 개의 그룹으로 분리를 할 때 중요한 정보를 제공한다. 영상으로부터 얻은 histogram은 대부분 이웃하는 그레이 값 사이에서 smooth 하게 변하지 않기 때문에 피크를 찾는 작업을 하기 전에 미리 mean filter나 gaussian filter와 같은 smoothing 과정을 거친 후 사용한다. 여기서는 low-pass filter 대신에 histogram의 bin 인덱스와 bin 값을 컨트롤 포인트로 사용해서 만든 Bezier 곡선을 이용해서histogram을 smooth 한 곡선으로 근사하는 방법을 알아본다. 이 경우 Bezier 곡선은 255-차수의 곡선이 된다. 높은 차원의 Bezier 곡선 계산에 Berstein 함수를 사용하는 경우 truncation 등의 수치 에러 때문에 값이 불안정해지므로 De Casteljau's algorithm을 이용하여서 iterative 하게 값을 계산을 하면 된다.
// De Casteljau's algorithm (degree:=size(Q)-1)
double Bezier(int deg, double Q[], double t) {
for (int k = 0; k < deg; k++)
for (int j = 0; j < (deg - k); j++)
Q[j] = (1 - t) * Q[j] + t * Q[j + 1];
return Q[0];
}
void SmoothenHistogram (int hist[], int numLevels/* =256 */) {
std::vector<double> p(numLevels);
std::vector<int> hist2(numLevels);
// cloning;
for (int j = 0; j < numLevels; j++) hist2[j] = hist[j];
for (int j = 0; j < numLevels; j++) {
double t = double(j) / (numLevels - 1);
// control points {p}; calling of Bezier() modifies p's;
for (int i = 0; i < numLevels; i++)
p[i] = hist2[i];
hist[j] = int(Bezier(numLevels-1, &p[0], t) + 0.5);
}
};
영상처리에서 한 픽셀의 수정을 위해서 주변 픽셀의 정보를 요구하는 윈도 기반 필터들은 일반적 연산 비용이 큰 편이다. 한 변의 길이가 W인 윈도를 사용할 때 W^2 횟수의 픽셀을 참조해야 하므로 윈도가 클수록 그 비용은 제곱으로 증가한다. 선형 필터 중에는 x-방향과 y-방향 연산으로 각각 분리가 가능한 경우는 연산 비용이 필터의 크기에만 비례하도록 만들 수 있다. median 필터는 이런 분리가 안 되는 비선형 필터 중 하나다. 근사적으로 x-방향으로 median filtering을 하고, 그 결과를 y-방향으로 다시 median filtering을 하는 방법으로 연산을 줄이는 방법을 사용하기도 한다.
윈도가 움직이면서 윈도 내 모든 픽셀이 다 바뀌는 것이 아니라 움직이는 방향에 수직인 가장자리 픽셀만 바뀐다는 사실을 이용하면 각 픽셀의 윈도에서 median을 찾는 작업을 할 필요가 없다. 예를 들면 scanline 방향(x-방향)으로 윈도를 움직이면서 median filter를 작용할 때, 윈도가 오른쪽으로 1 픽셀 움직이면 윈도의 왼쪽 가장자리의 수직 픽셀들은 새 윈도에서 사라지고, 기존 윈도의 오른쪽 수직 가장자리 앞의 픽셀들이 새로 들어온다. 따라서 각 스캔라인에서 처음 한 번만 윈도의 median을 찾으면 이후에는 윈도가 이동할 때 윈도를 나가는 픽셀과 새로 들어오는 픽셀 (2*W) 개에 대해서 이전 median과 비교만 하면 된다. 이 방법은 비교 횟수가 윈도 크기에 1차적으로 비례하므로 연산 부담을 많이 줄일 수 있다. 이 방법은 사각형 모양의 윈도를 가지는 다른 필터(mean filter, max-filter, min-filter,...)에 대해서도 쉽게 적용할 수 있다.
// boundary region도 처리할 수 있게 수정함; 2021-04-18;
// window size = wx * wy;
// median = argmin(hist[i] >= halfArea)
int RunningMedianFilter(BYTE* image, int w, int h, int wx, int wy, BYTE* out) {
int hist[256], x;
int wxhalf = wx >> 1;
int wyhalf = wy >> 1;
wx = (wxhalf << 1) + 1; // size of window = odd number;
wy = (wyhalf << 1) + 1;
for (int y = 0, yw = 0; y < h; ++y, yw += w) {
// calc available area;
int wystart = max(0, y - wyhalf);
int wystop = min(h, y + wyhalf);
int wysize = wystop - wystart + 1;
int wxstart = 0;
int wxstop = wxhalf;
int halfArea = (wxstop * wysize + 1) >> 1;
// to avoid *w multiplication in y-step;
wystart *= w;
wystop *= w;
// initial histogram of topleft window;
memset(hist, 0, 256 * sizeof(int));
for (int iy = wystart; iy <= wystop; iy += w)
for (int ix = wxstart; ix <= wxstop; ++ix) hist[image[iy + ix]]++;
// find initial median;
int ltmed = hist[0]; // less_than_median;
int med = 0;
while (ltmed < halfArea) ltmed += hist[++med];
out[yw + 0] = med;
// left edge rgn;
for (x = 1; wxstop < wx; ++x) {
++wxstop;
halfArea = (wxstop * wysize + 1) >> 1;
for (int iy = wystart; iy <= wystop; iy += w) {
int v = image[iy + wxstop]; // (x=wxstop) strip enters;
++hist[v];
if (v <= med) ++ltmed;
}
while (ltmed >= halfArea) ltmed -= hist[med--];
while (ltmed < halfArea) ltmed += hist[++med];
out[yw + x] = med;
}
// central rgn;
for ( ; wxstop < w; ++x) {
++wxstop;
for (int iy = wystart; iy <= wystop; iy += w) {
int v = image[iy + wxstart]; // (x=wxstart) strip leaves;
--hist[v];
if (v <= med) --ltmed;
v = image[iy + wxstop]; // (x=wxstop) strip enters;
++hist[v];
if (v <= med) ++ltmed;
}
++wxstart;
while (ltmed >= halfArea) ltmed -= hist[med--];
while (ltmed < halfArea) ltmed += hist[++med];
out[yw + x] = med;
}
// right edge rgn;
for ( ; x <= w; ++x) {
for (int iy = wystart; iy <= wystop; iy += w) {
int v = image[iy + wxstart]; // (x=wxstart) strip leaves;
--hist[v];
if (v <= med) --ltmed;
}
++wxstart;
halfArea = ((wxstop - wxstart + 1) * wysize + 1) >> 1;
while (ltmed >= halfArea) ltmed -= hist[med--];
while (ltmed < halfArea) ltmed += hist[++med];
out[yw + x] = med;
}
}
return 1;
};
Histogram equalization(HE)은 주어진 이미지의 픽셀 분포가 모든 픽셀 값에서 같은 확률로 나타나도록 픽셀 값을 변환하여 이미지를 보다 잘 인식되게 만드는 영상 처리 과정이다. 이러한 픽셀 값의 균일한 분포는 엔트로피의 관점에서 보면 최대 엔트로피를 주는 분포이기도 하다. 그러나, HE는 원본 이미지의 밝기를 유지하지 않는다. 따라서, 원본 이미지의 밝기를 유지하면서 엔트로피를 최대화시키는 히스토그램 분포를 찾은 후 그것으로 변환(histogram specification)을 시도하는 방법을 고려하는 것이 보다 현실적일 수 있다. 그러한 목표 히스토그램을 $f(s)$ (연속적인 확률 밀도 함수(pdf)로 취급)라고 하면,
정규화된 조건 (픽셀 값을 $[0,255] \rightarrow [0,1]$ 로 rescale 함) $$f(s) \ge 0, \quad \int_0^1 f(s)ds = 1.$$
밝기 보존: $$\int_0^1 sf(s) ds = \mu=\text{pixel mean of input image}.$$
의 근이다. 오른쪽은 $λ$에 대해서 단조함수이므로 근이 하나만 존재한다. 원본 이미지에서 픽셀 평균을 계산하여 위 방정식에 대입하면 $λ$를 구할 수 있고, 목표 히스토그램을 얻을 수 있다. 목표 히스토그램의 cdf와 ($\text {chist}(x)=\int_0^x sf(s) ds$), 원본 이미지의 cdf를 구해서, 그 차이를 최소로 하는 픽셀 값의 대응관계를 찾으면 된다 (histogram specification):
$$x \rightarrow y = F(x) = \int_0^x f(s) ds$$
참고: Bright Preserving Histogram Equalization with Maximum Entropy: A Variational Perspective. C. Wang and Z.Ye, IEEE Trans. Consumer Electronics. V51. No4. (2005); // 알고리즘 적용 단계에서 픽셀 값의 범위를 [0,255] -> [0,1]로 변환해야 한다. // BPHEME의 cumulative histogram(continuous version)::integral of f(s) over s;
double cdf(double s, double mu, double lambda) {
if (mu == 0.5) return s ;
else return (exp(lambda * s) - 1) / (exp(lambda) - 1);
}
// histogram specification;1==>2 void hist_spec(double chist1[], //source cumulative histogram; double chist2 [], // target cumulative histogram; int n, //=256; int lut[]) // resultant mapping(1->2); {
for (int i = 0; i < n; i++) {
int j = n - 1;
while (chist2[j] > chist1[i]) j--;
lut[i] = j < 0 ? 0 : j
}
}
void BPHEME(BYTE* src, int width, int height, BYTE *dst) {
double hist[256] = {0}; // src(dst) histogram;
double chist[256], chist2[256]; // src(dst) cumulative histogram;
int lut[256]; // histogram specification mapping;
const int n = width * height;
make_hist(src, width, height, hist);
normalize_hist(hist, 256);
//cumulative histogram;
make_cumulative_hist0(hist, 256, chist);
// gray-mean; note, pixel range should be changed [0,255] -> [0,1]
double mu = hist_mean(hist, 256);
// determine lambda;
double lambda, th =1.e-15;
FindRoot(mu, th, lambda);
// entropy of src;
double entropy1 = hist_entropy(hist, 256);
// dst-cumulative
for (int i = 0; i < 256; i++) chist2[i] = cdf(double(i) / 255., mu, lambda);
// histogram-specification;
hist_spec(&chist[0], &chist2[0], 256, &lut[0]);
// make dst image;
for (int i = 0; i < n; i++) dst[i] = lut[src[i]];
};
이미지를 이진화시키기 위해서 여러 알고리즘이 사용된다. 그중 이미지 전체에 대해 하나의 임계값으로 이진화시키는 전역 이진화 알고리즘은 간단하고 빠르기 때문에 많이 이용이 된다. 그러나 이미지를 형성할 때 조명 조건이 균일하지 않은 경우에는 전역 이진화는 원하는 결과를 얻기가 힘들다. 이런 경우에는 각각의 픽셀 주위의 그레이 값을 참조하여 임계치를 결정하는 국소적 이진화 방법을 사용한다. 국소적 이진화에서 임계값을 추출하는 간단한 방법은 윈도 내의 평균값을 이용하면 된다. 좀 더 개선된 알고리즘은 평균값($m(x, y)$)을 참조하되, 편차($\sigma(x, y)$)를 한번 더 고려해 주는 것이다. 이렇게 하여 잡은 국소적 임계값은 다음과 같이 표현된다:
여기서 $128$은 그레이 값이 가질 수 있는 최대 편차를 의미한다. 편차가 $128$이면 단순 평균값으로 취한다는 의미가 된다. 그 외의 경우는 표준편차와 128의 차이(항상 음수다)에 비례하는 값으로 윈도 평균값을 offset 한 값을 임계치로 잡는다. $\text{factor}$는 일반적으로 정해지지 않고, 실험적으로 $[0.2, 0.5]$ 사이의 값이 취해진다. (문서처럼 배경이 흰색인 경우는 $\text{factor} > 0$이지만, 검정 배경에 흰색 글씨를 처리하는 경우는 음수의 값을 취하는 것이 맞다) 국소적인 이진화 알고리즘은 매 픽셀마다 윈도를 잡아서 계산해야 하므로 연산 비용이 많이 든다. 충분한 메모리를 갖춘 시스템의 경우에는 적분 이미지(integral image)를 이용하면 윈도 연산에 소요되는 비용을 대폭 줄일 수 있다.. 국소적 이진화 알고리즘에서 윈도 크기와 $\text{factor}$를 결정하는 기준은 무엇일까? 이것은 해결하고자 하는 문제의 특성, 예를 들면 스캔된 문서를 이진화시키는 경우에는 윈도에 충분한 글자가 들어 있어야 한다... 등에 많이 의존한다.
void make_int_img12(BYTE *gray, int width, int height, *int intimage, int *intsqimg);
