Retinex 알고리즘은 영상의 밝기와 시각적으로 인지된 감각 사이에는 로그의 관계를 가진다는 사실과, 영상의 밝기는 실제의 밝기인 반사 성분과 조명에 의한 성분의 곱으로 주어진다는 실험적인 사실에 근거하여 영상에서 조명 성분을 줄이고, 반사 성분만을 나타냄으로써 영상의 콘트라스트를 증대시키고자 하는 시도이다.
$$B(x, y)=R(x, y)\times I(x, y)$$
where
$$ B(x, y)=\mbox {observed image}$$
$$ R(x, y)=\mbox {perceived reflectance(조명 성분)}$$
$$ I(x, y)=\mbox {perceived illumination(반사 성분)}$$
따라서,
$$\log R(x, y)= \log B(x, y) - \log I(x, y)$$
조명 성분을 구하기 위해서는 영상의 디테일에 의존하지 않는 큰 스케일에서의 영상 성분만 보면 된다 (low frequency 성분). 따라서 조명 성분은 원영상을 적당한 스케일($s$)로 blurring 한 영상에서 얻을 수 있다.
$$I\sim \text{Gaussfilter}(s)* B.$$
물론, 여러 스케일을 조합을 하여서 사용을 할 수 있다. 멀티 스케일을 적용하는 경우에 retinex 알고리즘에서 반사 성분은 (각 스케일의 기여에 같은 가중치를 준 경우)
$$\text{Dst}(x, y) = \sum_{s\in \text{scales}} \left(\log [\text{Src}(x, y)] - \log [ (\text{Gaussfilter}(s) * \text{Src})(x, y)] \right)$$
이다.
아래 코드는 출력 영상의 dynamic range를 $[m-1.2\sigma, m+1.2\sigma]$로 한정한 후, 픽셀 값을 $[0.255]$로 stretching 하였다. 픽셀 값에 $\log$를 취할 때는 $\log(0)$을 막기 위해서 $\mbox+1$을 하였다. 컬러 영상은 각각의 RGB 성분에 대해서 그레이 영상의 retinex process을 적용하면 된다.// recursive 가우시안 필터의 계수 계산;
void compute_coefs3(double c[5], double sigma) {
/* "Recursive Implementation of the gaussian filter.",
** Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
*/
double q;
if (sigma >= 2.5)
q = 0.98711 * sigma - 0.96330;
else if ((sigma >= 0.5) && (sigma < 2.5))
q = 3.97156 - 4.14554 * sqrt (1 - 0.26891 * sigma);
else
q = 0.1147705018520355224609375;
const double q2 = q * q;
const double q3 = q * q2;
c[0] = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));
c[1] = ( (2.44413*q)+(2.85619*q2)+(1.26661 *q3));
c[2] = ( -((1.4281*q2)+(1.26661 *q3)));
c[3] = ( (0.422205*q3));
c[4] = 1.0 - (c[1]+c[2]+c[3])/c[0];
};// Multiscale 설정;
void retinex_scale_distribution(const int nscales/*=3*/, const int s/*=240*/,
double scales [])
{
// ASSERT(nscales>=3);
const double size_step = double(s) / nscales;
for (int i = 0; i < nscales; ++i)
scales[i] = 2. + i * size_step;
}// recursive Gaussian convolution;
void gauss_smooth (double *in, int size, int rowstride, double *out, double b[5]) {
/* forward pass */
const int bufsize = size + 3;
std::vector<double> w1(bufsize, 0), w2(bufsize, 0);
w1[0] = in[0]; w1[1] = in[0]; w1[2] = in[0];
for (int i = 0, n = 3; i < size ; i++, n++)
w1[n] = b[4] * in[i*rowstride] + \
(b[1] * w1[n-1] + b[2] * w1[n-2] + b[3] * w1[n-3]) / b[0];
/* backward pass */
w2[size+0] = w1[size+2];
w2[size+1] = w1[size+2];
w2[size+2] = w1[size+2];
for (int i = size - 1, n = i; i >= 0; i--, n--)
w2[n] = out[i*rowstride] = b[4] * w1[n] + \
(b[1] * w2[n+1] + b[2] * w2[n+2] + b[3] * w2[n+3]) / b[0];
}// 이미지의 평균과 표준편차 계산;
void image_statistics(double *img, int size /*=width*height*/,
double &mean, double &std) {
double s = 0, ss = 0;
for (int i = size; i-- > 0;) {
double a = img[i];
s += a;
ss += SQR(a) ;
} ;
mean = s / size ;
std = sqrt((ss - s * s / size) / size);
}
//
void rescale_range(double *data, int size) {
double mean, sig;
image_statistics(data, size, mean, sig);
const double max_val = mean + 1.2 * sig;
const double min_val = mean - 1.2 * sig;
double range = max_val - min_val;
if (!range) range = 1;
range = 255 / range ;
// change the range;
for (int i = size; i-->0;)
data[i] = (data[i] - min_val) * range;
}// 메인;
void retinex_process(double *src, int width, int height, double *dst) {
const int nfilter = 3;
const double weight = 1.0 / double(nfilter); //equal-weight for each scale;
const int default_scale = 240;
double sigma[nfilter];
retinex_scale_distribution(nfilter, default_scale, sigma);
std::vector<double> in(width*height);
std::vector<double> out(width*height, 0);
// scale-space gauss_smooth;
double c[5];
for (int i = 0; i < nfilter; i++) {
compute_coefs3(c, sigma[i]);
// copy src to temporay buffer(=in);
for (int pos = width*height; pos-->0; )
in[pos] = double(src[pos]) + 1;
// (1) horizontal convolution(stride = 1 for grey);
for (int y = 0; y < height; y++)
gauss_smooth(&in[y * width], width, 1, &out[y * width], c);
// (2) vertical convolution(stride = width for grey);
in.swap(out); //copy out to in;
for (int x = 0; x < width; x++)
gauss_smooth(&in[x], height, width, &out[x], c);
// 각 스케일에서 반사 성분을 누적(equal weight);
for (int pos = width*height; pos-->0;)
dst[pos] += weight * (log(src[pos] + 1) - log(out[pos]));
}
const double alpha = 128;
const double gain = 1;
const double offset = 0;
for (int pos = width*height; pos-->0;) {
double logI = log(double(src[pos]) + 1);
dst[pos] = gain * ((log(alpha* (src[pos] + 1.)) - logI) * dst[pos]) + offset;
}
// scale to [0,255];
rescale_range(dst, width*height);
};/** 2020.11.24.일 수정됨;**/
원본 이미지;
가우시안 Convolution 결과(sigma[0]=2);
가우시안 Convolution 결과(sigma[1]=82);
가우시안 Convolution 결과(sigma[2]=162);
결과 이미지(scale=240)
'Image Recognition' 카테고리의 다른 글
| Mean Shift Filter (5) | 2008.08.06 |
|---|---|
| Chamfer Match (0) | 2008.08.01 |
| RANSAC: Circle Fit (0) | 2008.07.21 |
| KMeans Algorithm (0) | 2008.07.19 |
| Robust Line Fitting (0) | 2008.07.08 |
