Retinex 알고리즘은 영상의 contrast를 향상하거나, sharpness를 증진시킬 때 많이 사용한다. 또 픽셀 값의 dynamic range가 큰 경우에 이것을 압축시켜 영상 데이터 전송에 따른 병목 현상의 해소에 이용할 수 있다. Retinex 알고리즘의 기본 원리는 입력 영상에 들어있는 배경 성분을 제거하는 것이다.

1. 배경 영상은 입력 영상의 평균적인 영상로 생각할 수 있는데, 이것은 적당한 스케일(필터 사이즈)의 가우시안 필터를 적용하여서 얻을 수 있다. 이 필터를 적용하면 입력 영상에서 필터 사이즈보다 작은 스케일은 무시하는 효과를 준다.
2. 입력 영상의 반사 성분(배경 조명에 무관한)은 입력 영상을 앞서 구한 배경 영상으로 나누면 된다
3. Retinex 출력은 이 반사 성분에 로그 값을 취한 것이다. 로그를 취함으로써 반사 성분이 분포 범위(=dynamic range)를 압축하는 효과를 얻는다.

- 입력 영상: $I(x, y)$ ;
- 가우시안 필터: $G_{\sigma }(x, y)=A\exp [-(x^2+y^2)/2\sigma^2]$
- 배경 영상: $(G_{\sigma}*I)(x, y) = \text{Gaussian convolution}$;

- Retinex 출력:
\begin{align}R(x,y;\sigma)&=\log[\text{반사 성분}] \\ &=\log[\text{입력 성분}/\text{배경 성분}]\\ &=\log [I(x, y)] - \log[ (G_{\sigma} * I)(x,y)]\end{align}

이처럼 하나의 스케일에 대해서 적용하는 경우를 SSR(single-scale retinex) 알고리즘이라고 한다. 컬러 영상의 경우에는 각각의 RGB 채널에 대해 알고리즘을 적용하면 된다. 


Retinex 영상을 구할 때 하나의 스케일이 아니라 다중 스케일에 대해서 적용한 retinex 영상을 적절한 가중치($\omega_\sigma$)를 주어서 더한 결과를 출력 영상으로 사용할 수 있다. 이 경우가 MSR(multi-scale retinex) 알고리즘이다. $$R_{MSR}(x,y) = \sum_{\sigma\in scales} \omega_{\sigma} R_{SSR}(x,y;\sigma) $$

이렇게 얻은 Retinex 출력 영상은 적당한 offset과 stretching을 하여서 픽셀 값이 [0,255] 구간에 있게 조절한다. 이 과정은 Retinex 출력 영상 픽셀의 평균값과 편차를 구하여 해결한다.

컬러 영상의 경우에 Retinex 출력 영상은 전체적으로 그레이화 되는 경향이 있어서 이것을 보완하기 위해서 아래의 추가적인 처리 과정을 더 거친다.  $$ {\tilde R}_{MSR}^{red}(x,y)=\log \left(\frac{C I_{red}}{ I_{red}+I_{green}+I_{blue}} \right) R_{MSR}^{red}(x,y)$$ $$ {\tilde R}_{MSR}^{grren}(x,y)=\log \left(\frac{CI_{green}}{ I_{red}+I_{green}+I_{blue}} \right) R_{MSR}^{green}(x,y)$$ $$ {\tilde R}_{MSR}^{blue}(x,y)=\log \left(\frac{C I_{blue}}{I_{red} +I_{green}+I_{blue}} \right) R_{MSR}^{blue}(x,y)$$

여기서, $C$는 상수이고, $I_{red}, I_{green}, I_{blue}$는 입력 영상의 RGB-채널이다.


*그레이 이미지 처리 결과:


*구현 코드는 다음을 참고: 그레이: http://kipl.tistory.com/33 
                                   컬러: http://blog.naver.com/helloktk/80039132534
* 기술적으로  log를 취하므로 입력 영상에 +1을 더해서 log(0)이 나오는 것을 방지해야 한다.
* convolution 된 영상도 1보다 작은 경우에 1로 만들어야 한다.
* 스케일이 큰 경우에 필터 사이즈가 매우 크므로 효과적인 convolution 알고리즘을 생각해야 한다. 보통 recursive 필터링을 하여서 빠르게 convolution결과를 얻는다.

 

728x90

'Image Recognition' 카테고리의 다른 글

Blur Detection  (0) 2010.05.25
Savitzky-Golay Smoothing Filter  (2) 2010.03.24
Gaussian Mixture Model & KMeans  (4) 2010.01.30
Image Morphing  (0) 2010.01.24
Fant's Algorithm  (0) 2010.01.22
Posted by helloktk
,
728x90

'Image Recognition' 카테고리의 다른 글

Similarity Transformation  (1) 2009.12.14
Eigenface  (0) 2009.12.12
Spline Based Snake  (0) 2008.08.15
Anisotropic Diffusion Filter  (0) 2008.08.11
Rolling Ball Transformation  (1) 2008.08.09
Posted by helloktk
,

 

vliet9504.pdf
다운로드

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을 적용하면 된다.


Q: How to find a optimal scale distribution for a given image?

#define SQR(x) ((x)*(x))
// 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 = 0;   
    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;
    
    double q2 = q * q;
    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];
}

// 멀티스케일 설정;
void retinex_scale_distribution(const int nscales/*=3*/, const int s/*=240*/,
                                double scales []) {

더보기
    //  ASSERT(nscales>=3);
    double size_step = double(s) / nscales;
    for (int i = 0; i < nscales; ++i)
        scales[i] = 2. + i * size_step;
}

// 가우시안  convolution(recursive);

void gausss_mooth (double *in, int size, int rowstride, double *out, double b [5]) {

더보기
    /* forward pass */
    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[0], size, &mean, &sig);
  
    double max_val = mean + 1.2 * sig;
    double min_val = mean - 1.2 * sig;    
    double range = max_val - min_val;
    if (!range) range = 1.0;
    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;
    double sigma[nfilter];
    double c[5];
    int default_scale = 240;
    retinex_scale_distribution(nfilter, default_scale, sigma);
    int csize = width * height;
    std::vector<double> in(csize), out(csize)
    memset(dst, 0, csize * sizeof(double));
    // scale-space gauss_smooth;
    double weight = 1.0 / double(nfilter); //equal-weight for each scale;
    for (int i = 0; i < nfilter; i++) {
        compute_coefs3(c, sigma[i]);
        // copy src to temporay buffer(=in);
        for (int pos = csize; 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++) 
            gausss_mooth(&in[x], height, width, &out[x], c);
        // 각 스케일에서 반사 성분을 누적(equal weight);
        for (int pos = csize; pos-- > 0;) 
            dst[pos] += weight * (log(src[pos] + 1.) - log(out[pos]));
    }
    double alpha = 128;
    double gain = 1;
    double offset = 0;
    for (int pos = csize; pos-- > 0;) {
        double logI = log(double(src[pos]) + 1.0);
    	dst[pos] = gain * ((log(alpha* (src[pos] + 1.)) - logI) * dst[pos]) + offset;
    }
    // scale to [0,255];
    rescale_range(&dst[0], csize); 
};

/** 2020.11.24.일 수정됨;**/

원본 이미지;

사용자 삽입 이미지

 

가우시안 Convolution 결과(sigma[0]=2);

 

가우시안 Convolution 결과(sigma[1]=82);

 

가우시안 Convolution 결과(sigma[2]=162);

 

결과 이미지(scale=240)

 

 
 
728x90

'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
Posted by helloktk
,