방 안에 뜨거운 물체가 놓이면 그 열은 시간이 지나면서 사방으로 퍼져나간다. 열원에서 주위로 열이 시간에 따라 전달이 되는 프로세스를 기술하는 방정식이 열(전도)방정식이다. 이것을 이미지 필터링에서도 적용이 가능하다. 이미지에 나타난 각각의 노이즈를 열원으로 보면 주변으로 열의 전달은 열원의 온도를 낮추게 되는데 이를 이미지에서 노이즈를 제거하는 filtering 과정으로 생각할 수 있다. 그러나 열방정식은 방향성이 없는 등방성 확산(isotropic diffusion process)을 기술한다. 이미지에 열전도 방정식을 바로 적용하면 이미지의 블러링과 유사한 효과를 만든다. filtering를 하더라도 에지를 보존해야 이미지의 정보 손실이 적게 된다. 에지 영역에서는 gradient의 공분산 행렬($\Sigma$)의 고윳값의 차이가 크고, 노이즈나 corner 근방에서는 두 고윳값의 크기가 비슷하다. 따라서, 에지를 보존하기 위해서는 공간으로 확산하는 역할을 하는 이미지의 Laplace 항 ($\nabla^2 I = I_{xx} + I_{yy}$)을 gradient의 공분산 행렬의 고유 벡터와 고유치, 그리고 이미지의 hessian 행렬 등을 조합하여서 불변인 형태를 만들어야 한다. 여러 가능한 후보군 중에 하나의 예를 들면
여기서, $\bf U, V$는 그라디엔트 공분산 행렬의 두 고유 벡터이고, $s=\lambda_1 + \lambda_2$는 고유값의 합이다. $f_1$과 $f_2$의 함수 형태를 고유치가 큰 방향에서는 억제되고, 고유치가 작은 방향에서는 덜 억제가 되도록 적절히 선택하는 경우에 에지를 보존하는 필터링의 효과를 얻을 수 있다.
아래의 예는 이 식을 써서 300번 iteration을 한 결과이다. 원본 이미지는 첨부한 논문에서 추출하였다.
이미지를 획득할 때 조명 조건이 균일하지 않으면 이미지 처리 과정에서 복잡성이 증가하게 된다. 이런 경우에 배경의 불균일성을 제거하는 전처리 과정을 수행할 필요가 생기게 된다. 일반적으로 조명의 불균일성은 공간적으로 넓은 범위에서 천천히 일어나게 되는데, 이를 이용하면 이미지에서 배경을 좀 더 수월하게 제거할 수 있다. 보통 큰 윈도를 (전경의 물체의 크기에 비해서) 가지는 median filter를 연속적으로 적용하여 배경을 추출할 수 있지만, median filter는 일반적으로 계산 비용이 상당히 비싸다.
좀 더 저렴한 연산비용이 드는 방법으로는 이미지가 각 픽셀 좌표에서 하나의 값을 할당하므로 3차원 공간에서 곡면으로 생각할 수 있다는 점을 이용하면 된다. 일정하게 변하는 조명에 의한 배경은 3차원에서 평면을 형성하는 것으로 근사할 수 있다. 주어진 이미지 곡면에서 배경 평면을 추출하는 쉬운 방법은 최소자승법을 사용하여 이미지 곡면을 가장 잘 근사하는 평면을 찾는 것이다. 그러나 이 방법은 전경의 물체가 fitting에 있어 일종의 outlier 역할을 하게 된다. 전경에 해당하는 물체가 많은 이미지에서도 잘 동작하려면 보다 robust 한 fitting 알고리즘을 고안해야 한다.
3차원 공간에서 이미지 곡면을 고려할 때, 전경(밝은 색)이 이미지 영역에서 차지하는 비중이 작은 경우에는 그림처럼(A) 충분히 큰 공을 곡면 아래쪽에서 접하게 굴리면 공의 표면이 형성하는 곡면은 이미지의 작은 요철(전경 물체)에 영향을 거의 받지 않는다. 따라서 구르는 공에 의해 형성된 접곡면은 이미지의 배경으로 생각할 수 있다 (B는 전경이 어두운 경우).
이미지 곡면이 요철이 심하더라도 반지름이 충분히 큰 공을 사용하게 되면 부드러운 배경 표면을 얻을 수 있다. 공 표면과 이미지 표면과의 접촉 관계는 공의 표면에 놓인 patch의 좌표점에서 (미리 계산) 높이를 이미지의 픽셀 값과 비교하여 얻는다. 공의 반경이 큰 경우에는 연산을 줄이기 위해서 원본 이미지를 축소해서 작업을 하고, 다시 간단한 보간을 써서 원본 크기로 복원하면 충분하다.
다음 예는 밝은 배경에 어두운 전경인 이미지에 대해서 이 알고리즘을 적용한 경우다. 먼저, 원본 이미지의 negative 이미지를 만들어서 (어두운 배경에 밝은 전경) rolling ball transformation을 적용해서 배경을 얻었다.
#define SQR(x) ((x)*(x))
struct RollingBall {
std::vector<BYTE> data;
int patchwidth;
int shrinkfactor;
RollingBall(int radius) {
int artrimper;
if (radius <= 10) {
shrinkfactor = 1;
artrimper = 12; // trim 24% in x and y
} else if (radius <= 30) {
shrinkfactor = 2;
artrimper = 12; // trim 24% in x and y
} else if (radius <= 100) {
shrinkfactor = 4;
artrimper = 16; // trim 32% in x and y
} else {
shrinkfactor = 8;
artrimper = 20; // trim 40% in x and y
}
Build(radius, artrimper);
}
void Build(int ballradius, int artrimper) {
int sballrad = ballradius/shrinkfactor;
if (sballrad < 1) sballrad = 1;
int rsquare = SQR(sballrad);
int diam = sballrad * 2;
int xtrim = (artrimper * diam) / 100; // only use a patch of the rolling ball
patchwidth = diam - xtrim - xtrim;
int halfpatchwidth = sballrad - xtrim;
int ballsize = SQR(patchwidth + 1);
data.resize(ballsize);
for (int i=0; i<ballsize; i++) {
int xval = i % (patchwidth+1) - halfpatchwidth; //relative to patch center-x.
int yval = i / (patchwidth+1) - halfpatchwidth; //relative to patch center-y;
int temp = rsquare - SQR(xval) - SQR(yval);
if (temp >= 0)
data[i] = int(sqrt(temp));
else
data[i] = 0;
}
}
};
// processing을 빠르게 하기 위해서, ball의 반경이 큰 경우에는 원본이미지를
// 축소하여서(2,4,8,..) 처리한 후에, 처리 안된 점들은 interpolation을 하여서
// bacground image을 생성한다.
void RollBall(RollingBall& ball, //rolling ball object;
BYTE *smallImage, int swidth, int sheight,//shrinked source image;
BYTE *background) //shrinked background image
{
int wsz = ball.patchwidth ; //patch width(odd);
int hwsz = wsz >> 1;
BYTE *patch = &ball.data[0] ; //pre-calculated patch z-value;
int zctr = 0; // start z-center in the xy-plane
for (int y = 0; y < sheight; y++) {
int ystart = y - hwsz;
int yend = y + hwsz;
for (int x = 0; x < swidth; x++) {
// 현재의 패치내에서(중심(x-patchwidth/2, y-patchwidth/2)) 픽셀의 z값(그레이레벨)과
// zctr만큼 shift한 공의 윗표면에 놓인 패치의 z-값간의 차이를 구함. 그 값의 차이가
// 가장 작을 때 값만큼 공의 중심을 높이거나 내리면 공과 픽셀 surface가 밑에서
// 접하게 된다.
int zmin = 255;
int xstart = x - hwsz;
int xend = x + hwsz;
for (int yy = ystart, ballpt = 0; yy <= yend; yy++) {
for (int xx = xstart, imgpt = yy * swidth + xx; xx <= xend; xx++) {
if (xx >= 0 && xx < swidth && yy >= 0 && yy < sheight) {
int zdif = smallImage[imgpt] - (zctr + (patch[ballpt]));
if (zdif < zmin) zmin = zdif;
}
ballpt++; imgpt++;
}
}
if (zmin != 0) zctr += zmin;
// zmin<0 인 경우는 x의 오른쪽 반만 처리하면 된다(not yet);
for (int yy = ystart, ballpt = 0; yy <= yend; yy++) {
for (int xx = xstart, imgpt = yy * swidth + xx; xx <= xend; xx++) {
if (xx >= 0 && xx < swidth && yy >= 0 && yy < sheight) {
int zadd = zctr + (patch[ballpt]);
if (zadd > (sbackground[imgpt]))
sbackground[imgpt] = (BYTE)zadd;
}
ballpt++; imgpt++;
}
}
}
}
}
Mean-shift는 확률 밀도 함수의 극대점(mode)을 연속적으로 찾는 방법이다. 주어진 데이터(특징점)가 샘플링된 확률 밀도 함수의 평균과 분산을 알 수 없을 때 적용할 수 있는 non-parametric 알고리즘이다. Clustering, tracking, 또는 filtering 등에 사용될 수 있다. Mean-shift filter는 각 픽셀을 중심으로 하는 윈도 내의 모든 픽셀에 대해서 특징점 공간(= 픽셀 좌표 공간 + 픽셀 컬러 공간)에서 밀도가 극대인 특징점(특징점 공간에서 질량중심)을 찾고, 이 특징점에 해당하는 픽셀의 컬러로 윈도 중심 픽셀의 컬러를 대체한다.
아래 예제는 컬러 공간을 RGB에서 YIQ로 바꾸어서 알고리즘을 적용한 후에 그 결과를 RGB로 다시 변환하였다. 윈도를 크게 하면 알고리즘의 수행 속도가 매우 느려진다. 공간 윈도를 원에서 정사각형으로 바꾸어도 된다.
Clustering 알고리즘으로 볼 때 Mean-Shift 알고리즘의 특징은
단순한 알고리즘 구조를 갖지만 연산 비용이 매우 크다. 고차원 특징에 대해서는 부적합
클러스터 수에 대한 제한이 없다.(k-means와 비교)
클러스터 초기 위치 설정이 필요없다.(k-means와 비교)
outliers의 영향이 적다.(k-means와 비교)
유일한 매개변수는 윈도 크기다. 원도 크기를 미리 설정해야 한다.
// 구현 code
void RGB2YIQ(BYTE *color, int w, int h, float *ycomp, float *icomp, float *qcomp);
주어진 영상이 어떤 특정한 물체를 담고 있는지를 알아보려면, 그 물체의 template을 만들어서 영상 위에 겹친 후 픽셀 단위로 순차적으로 이동시키면서 일치하는가를 체크하면 된다. 일치의 여부는 보통의 경우 템플레이트와 영상의 상관계수를 측정하면 된다. 대상 영상이 이진 영상인 경우에는 distance map을 이용하면 좀 더 쉽게 매칭 문제를 해결할 수 있다. 에지 영상을 사용해서 만든 distance map은 이미지 영역의 한 지점에서 가장 가까이에 있는 에지 위치까지의 거리를 알려준다. 매칭하고자 template의 이진 영상을 담고 있는 이미지를 distance map에서 픽셀 단위로 차례로 옮기면서 template 위치에서 distance map이 알려주는 에지까지의 거리 합을 기록하면 총거리의 합이 가장 작은 지점에서 template와 가장 잘 맞는 형상을 찾을 수 있을 것이다.
이 메칭 기법(Chamfer Match)은 모델과 이미지가 서로 회전되어 있지 않아야 하고, 스케일 변형도 없는 경우에 좋은 결과를 얻는다. 스케일 변화가 있는 경우에는 스케일 공간에서 순차적으로 찾는 방법을 쓰면 된다. 모델 영상에 변화가 있는 경우에는 민감하게 반응하므로, 각각의 변화된 모델 영상을 따로 준비하여서 매칭 과정을 수행해야 한다. distance map만 구해지면 매칭 과정이 매우 단순하기 때문에 매칭 비용도 매우 저렴하다. 그리고 거리 값을 적당한 범위에 truncate 하는 것이 보다 robust 한 결과를 준다.
* procedure:
입력영상의 에지영상을 구함.
입력영상의 에지영상의 distance map을 구함.
매칭과정을 실행.
int GetAddressTable(BYTE *subImg, int w, int h, int stride/*=원본이미지 폭*/, int addr[]);
//계산을 보다 빠르게 하기 위해서 subimg에서 template 위치를(distimg에 놓인 경우) 미리 계산함.
int GetAddressTable(BYTE *subImg, int w, int h, int stride/*=원본이미지 폭*/, int addr[]) {
for(int y=0, k=0; y<h; y++) {
for(int x=0; x<w; x++){
if(subImg[y*w+x])
addr[k++]= y*stride + x ; //template가 dist이미지에서움직일 때 픽셀위치(상대적);
}
}
return (k); //# of template pixels<= w*h;
}
Matching code;
void ChamferMatch(float* distImg, int w1, int h1, //distance-map;
BYTE* subImg, int w2, int h2, //template-map;
float *match/*[w1 * h1]*/) { //match_score-map(filled with FLT_MAX)
int umax = w1 - w2;
int vmax = h1 - h2;
int uc = w2 / 2; //center_x of template-map;
int vc = h2 / 2; //center_y of template-map;
int *addr = new int [w2 * h2]; //max;
int n = GetAddressTable(subImg, w2, h2, w1, addr) ;
for (int v = 0; v <= vmax; v += 2){ //to speed-up;
for (int u = 0; u <= umax; u += 2){ //to speed-up;
double score=0;
#if (0)
for (int v1 = v, v2 = 0; v1 < h1 && v2 < h2; v1++, v2++){
int i1 = v1 * w1;
int i2 = v2 * w2;
for (int u1 = u, u2 = 0; u1 < w1 && u2 < w2; u1++, u2++) {
if (subImg[i2 + u2])
score += distImg[i1 + u1]; // truncate [0, th_dist];
}
}
#else
int offset = v * w1 + u; //start address;
int i = n;
while(i--) score += distImg[addr[i] + offset];
#endif
//at the center of subimg;
match[(uc + u) + (vc + v) * w1] = score;
}
}
delete[] addr ;
};
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]];
};
