Geometry & Recognition

Histogram equalization(HE)은 주어진 이미지의 픽셀 분포가 모든 픽셀 값에서 같은 확률로 나타나도록 픽셀 값을 변환하여 이미지를 보다 잘 인식되게 만드는 영상 처리 과정이다. 이러한 픽셀 값의 균일한 분포는 엔트로피의 관점에서 보면 최대 엔트로피를 주는 분포이기도 하다. 그러나, HE는 원본 이미지의 밝기를 유지하지 않는다. 따라서, 원본 이미지의 밝기를 유지하면서 엔트로피를 최대화시키는 히스토그램 분포를 찾은 후 그것으로 변환(histogram specification)을 시도하는 방법을 고려하는 것이 보다 현실적일 수 있다. 그러한 목표 히스토그램을 $f(s)$ (연속적인 확률 밀도 함수(pdf)로 취급)라고 하면,

1. 정규화된  조건 (픽셀 값을 $[0,255] \rightarrow [0,1]$ 로 rescale 함)
   $$f(s) \ge 0, \quad \int_0^1 f(s)ds = 1.$$
2. 밝기 보존: 
   $$\int_0^1 sf(s) ds = \mu=\text{pixel mean of input image}.$$

이러한 제약 조건에서 아래의 최대 엔트로피 조건을 만족시켜야 한다.

$$\underset{f}{\text{argmax}}\int_0^1 \Big[ -f(s) \ln f(s)\Big] ds.$$

위 조건를 만족시키는 $f(s)$는 Lagrange-multiplier를 이용하고, 변분 방법을 쓰면 쉽게 구해진다: 다음의 범함수 $J [f]$

$$J=-\int_0^1  f(s) \ln f(s) ds + \lambda_1  \left( \int_0^1 f(s) ds-1\right) + \lambda \left(\int_0^1 s f(s) ds -\mu\right)$$

의 극값을 구하면 된다.

$$\frac{\delta J}{\delta f}=0$$

에서

$$ -\log f - 1 + \lambda_1 f +\lambda  s = 0$$

이를 정리하면,

$$f(s)=e^{\lambda s } e^{\lambda_1-1}$$

을 얻고, 제약 조건을 적용하면 ($ e^{\lambda_1 -1}=\lambda/(e^{\lambda}-1))$

$$f(s) = \left\{\begin{array}{ll} 1,& \mu = 0.5 \\ \frac{\lambda e^{\lambda s}}{e^\lambda -1},& \mu \ne 0.5 \end{array} \right.$$

그리고 $λ$는 

$$ \mu = \frac{\lambda e^\lambda - e^\lambda +1}{\lambda (e^\lambda -1)},\quad(\mu \ne 0.5)$$

의 근이다. 오른쪽은 $λ$에 대해서 단조함수이므로 근이 하나만 존재한다.  원본 이미지에서 픽셀 평균을 계산하여 위 방정식에 대입하면 $λ$를 구할 수 있고, 목표 히스토그램을 얻을 수 있다. 목표 히스토그램의 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]];
};

// Root finding procedure ::
// define F(s, mu);
double F(double s, double mu) {

더보기

    if (s == 0)
        return (mu - 0.5);
    double ev = exp(s);
    return mu - (s * ev - ev + 1.) / (s * (ev - 1.)) ;
};

//derivative of F(s, mu) w.r.t. s;
double DF(double s, double mu) {   

더보기

    if (s == 0)
        return -1.0 / 12.0;
    double ev = exp(s);
    double ss = s * s;
    return -(ev * (ev - ss - 2.) + 1.) / (ss * (ev - 1.)*(ev - 1.));
}

// Find root of F(s,mu) based on Newton-Raphson method.
int FindRoot(double mu, double threshold, double& s) {     
     int iter = 0, max_iter = 100;
     s = 0. ; // initial guess; problem specific!
     for (iter = 0; iter < max_iter; ++iter) {
          double sold = s;
          // ASSERT(DF)
          s = s - F(s, mu) / DF(s, mu) ;
          if (fabs(s - sold) < threshold)
              break ;
     } ;
     TRACE("terminating iters=%d\n", iter);
     return (iter < max_iter) ;
};

histogram equaliztion 결과;

BPHEME  결과:

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)

RANSAC 알고리즘을 써서 주어진 2차원 점집합에서 원을 추정한다. 원을 만들기 위해서는 최소한 3점이 필요하고, 또 일직선에 있지 않아야 한다. 이렇게 만들어진 원은 세 점을 꼭짓점으로 하는 삼각형의 외접원이다. 주어진 외접원에서 크게 벗어나지 않는 inliers를 찾고, 이 inliers에 대해 최소자승법으로 원의 중심과 반지름을 다시 추정한다. 무작위로 선택된 세 점에 대해 위의 과정을 반복 시행해서 구한 원 중에서 inliers의 벗어남 편차가 최소인 것을 결과로 선택한다.
// 참고: http://en.wikipedia.org/wiki/RANSAC 
//

// [0, max_num) 사이에서 3개의 숫자를 무작위로 선택함;
void GetRandomTriplet(int max_num, int triplet [3]);

더보기

void getRandTriple(int max_num, int triplet[3]){
    // [0, max_num] select 3 numbers;
    int index = 0;
    while (index < 3) {
        BOOL newone = TRUE;
        int r = (int)((double)(rand()) / RAND_MAX * max_num);
        for (int i = 0; i < index; ++i)
            if (r == triplet[i]) { newone = FALSE; break; }
        if (newone) triplet[index++] = r;
    }   
}

// 세 점의 외접원을 구함;
int CircumCircle(double tx[3], double ty[3], double cparams[4]) ;

더보기

int circumCircle(double tx[3], double ty[3], double cparam[4]) {
    double x1 = tx[0], x2 = tx[1], x3 = tx[2];
    double y1 = ty[0], y2 = ty[1], y3 = ty[2];
    double bax = x2 - x1, bay = y2 - y1;
    double cax = x3 - x1, cay = y3 - y1;
    double E = bax * (x1 + x2) + bay * (y1 + y2);
    double F = cax * (x1 + x3) + cay * (y1 + y3);
    double G = 2. * (bax * (y3 - y2) - bay * (x3 - x2));
    if (G == 0.) return 0; //error;
    //assert(fabs(G)>small_epsilon); //to prevent collinear or degenerate case;
    cparams[0] = (cay * E - bay * F) / G; //cx;
    cparams[1] = (bax * F - cax * E) / G; //cy;
    cparams[2] = hypot(cparams[0]-x1, cparam[1]-y1); //radius;
    return 1;
}

// 3x3 선형 방정식을 푼다: A * x = b;
bool SolveLinearEQ3x3(double A [9], double bb [3], double x [3]) ;

더보기

bool SolveLinearEQ3x3(double A[9], double bb[3], double x[3]) {
    double invA[9];
    double det=(A[0]*(A[4]*A[8]-A[5]*A[7])-A[1]*(A[3]*A[8]-A[5]*A[6])+A[2]*(A[3]*A[7]-A[4]*A[6])) ;
    if (det != 0.) {
        det = 1./det;
        invA[0] = (A[4] * A[8] - A[5] * A[7]) * det;
        invA[1] = (A[2] * A[7] - A[1] * A[8]) * det;
        invA[2] = (A[1] * A[5] - A[2] * A[4]) * det;
        invA[3] = (A[5] * A[6] - A[3] * A[8]) * det;
        invA[4] = (A[0] * A[8] - A[2] * A[6]) * det;
        invA[5] = (A[2] * A[3] - A[0] * A[5]) * det;
        invA[6] = (A[3] * A[7] - A[4] * A[6]) * det;
        invA[7] = (A[1] * A[6] - A[0] * A[7]) * det;
        invA[8] = (A[0] * A[4] - A[1] * A[3]) * det;
        //
        x[0] = invA[0] * bb[0] + invA[1] * bb[1] + invA[2] * bb[2];
        x[1] = invA[3] * bb[0] + invA[4] * bb[1] + invA[5] * bb[2];
        x[2] = invA[6] * bb[0] + invA[7] * bb[1] + invA[8] * bb[2];
        return true;
    } else {
        x[0] = x[1] = x[2] = 0;
        return false;
    }
}

// x^2 + y^2 + A*x + B*y + C = 0;
// Least square fit for A, B, C;(참조: kipl.tistory.com/207)
int CircleFit_LS(std::vector<double>& xp, std::vector<double>& yp, double cparams[4]) ; 

더보기

int circleFit_LS(std::vector<double>& xp, std::vector<double>& yp, double cparams[4]) {  
    double sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0;
    double sxxx = 0, sxxy = 0, sxyy = 0, syyy = 0 ;
    for (int i = xp.size(); i-->0;) {
        double x = xp[i], y = yp[i];
        double xx = x * x, yy = y * y;
        sx += x;        sy += y;
        sxx += xx;      sxy += x*y ;       syy += yy ;
        sxxx += xx * x; sxxy += xx * y;
        sxyy += x * yy; syyy += yy * y; 
    } ;
    double A[9], b[3], sol[3];
    A[0] = sxx, A[1] = sxy, A[2] = sx,
    A[3] = sxy, A[4] = syy, A[5] = sy,
    A[6] = sx,  A[7] = sy,  A[8] = N,
    b[0] = -sxxx - sxyy,
    b[1] = -sxxy - syyy,
    b[2] = -sxx - syy;
    if (!SolveLinearEQ3x3(A, b, sol)) return 0;
    double det = sol[0] * sol[0] + sol[1] * sol[1] - 4 * sol[2];
    if (det <= 0) return 0;
    cparams[0] = -sol[0] / 2 ; //cx;
    cparams[1] = -sol[1] / 2 ; //cy;   
    cparams[2] = sqrt(det)/2.; //radius;
    return 1;
};

// 주어진 원에서 임계 거리 이내의 데이터만 골라냄;

double findInlier(std::vector<double>& xp, std::vector<double>& yp,
                double cparams[4], double dist_th,
                std::vector<double>& xinliers, std::vector<double>& yinliers) {
    xinliers.clear(); yinliers.clear();    
    double rms = 0;           // root-mean-square;
    double cx = cparams[0], cy = cparams[1], rad = cparams[2];
    for (int k = xp.size(); k-->0;) {
        double dist_dev = fabs(hypot(xp[k]-cx, yp[k]-cy)-rad)/rad ;
        if (dist_dev <= dist_th) { // collect maybe_inliers;
            xinliers.push_back(xp[k]);
            yinliers.push_back(yp[k]);
            rms += SQR(dist - rad);
        }
    }
    return sqrt(rms / xinliers.size());
}

int RansacCircleFit(std::vector<double>& xp, std::vector<double>& yp,
                    int sample_th,      // # of inliers; >= 50% of data(N), 66%; 
                    double dist_th,     // 25% of radius;   |dist-rad|/rad< dist_th
                    int max_iter,
                    double cparams[4]) {    
    double pr = double(sample_th) / double(xp.size());
    double trials99 = log(1. - 0.99)/log(1. - pow(pr, 3));
    int iter = min(max_iter, trials99);
    bool found = false;   
    double min_rms = 1.e+10;
    double cx, cy, rad;
    if (sample_th < 3) sample_th = 3;
    while (iter--) { 
        int tri[3]; 
        double tx[3], ty[3];
        getRandTriple(xp.size()-1, tri);
        for (int i = 0; i < 3; i++)
            tx[i] = xp[tri[i]], ty[i] = yp[tri[i]];
            
        if (!circumCircle(tx, ty, cparams)  
            // if three points are degenerate or on a line, discard them!
            continue ;
        std::vector<double> consensusx, consensusy;
        sdev = findInlier(xp, yp, cparams, dist_th, consensusx, consensusy);
        if (consensusx.size() >= sample_th) {          
            // estimate model params using the maybe-inliers;
            if (!circleFit_LS(consensusx, consensusy, cparams)) 
                continue; // least-square fitting fails;
            // collect real-inliers;
            double rms = findInlier(xp, yp, cparams, dist_th/2, consensusx, consensusy); 
            if (consensusx.size() < sample_th) continue;
            // estimate model params using inliers again;
            if (!circleFit_LS(consensusx, consensusy, cparams)) continue;            
            TRACE("cx = %f, cy = %f, rad=%f\n", cparams[0], cparams[1], cparams[2]);
            if (rms < min_rms) {
                min_rms = rms; found = true;
                // we need more elaborate termination conditions;
#if _DEBUG
#endif
            }
        }
    }
    return found;
};

• sample_th = 2 * N / 3;
• dist_deviate_th = 0.25;
• 파란색: 최소자승법을 이용한 피팅 (outlier의 영향을 많이 받음);
• 붉은색: RANSAC 피팅 결과 (2017.01.04 수정)

see also http://blog.naver.com/helloktk/80029898273