Geometry & Recognition

영상에는 보통 우리가 관심이 있는 물체와 배경을 동시에 담고 있다. 그런데 어떤 영상의 경우 전체 픽셀 값의 다이내믹 레인지(dynamic range)가 매우 좁아서, 달리 이야기하면 영상의 히스토그램을 살펴볼 때 픽셀 값의 분포가 매우 좁은 영역에 한정이 되어있으면 물체와 배경의 구별이 잘 안되는 영상이 되게 된다. 이러한 영상의 경우 몇 가지 영상 처리 방법을 이용하면 좋은 contrast를 갖는 영상을 만들 수 있다.

Histogram Equalization(HE)은 이러한 기법 중의 하나로 좁은 영역에 분포하는 픽셀 값을 가능한 모든 범위에 골고루 분포하도록 재배치하여서 영상의 contrast를 증대시키는 방법이다. Histogram Equalization에서 사용하는 픽셀 값의 mapping은 

(1) 출력 영상의 픽셀 값 범위는 전 그레이 레벨이어야 하고 

(2) 각각의 그레이 레벨에 해당하는 픽셀의 수가 거의 일정하여야 한다: 픽셀 값이 이산적이어서 완벽하게 상수가 될 수는 없다, 좀 더 정확히는 출력 영상의 히스토그램 누적합이 그레이 레벨에 선형 비례해야 한다는 조건으로 결정된다.

입력 영상에서 값이 k 인 픽셀은 출력 영상에서 어떤 그레이 값을 가져야 할까? 조건에 따르면 입력 영상의 k 그레이 레벨이 출력 영상의 j 그레이 레벨로 매핑이 된다면, 출력 영상에서는 0부터 j까지 그레이 레벨에 해당하는 픽셀 수는 모든 그레이 레벨에서 픽셀 분포가 균일하다는 조건 때문에 

로 주어진다. 이 값과 입력 영상에서 0부터 k까지 그레이 레벨을 갖는 픽셀의 수가 같으면 된다.

입력 영상에서 그레이 값이 0부터 k까지인 픽셀 수는 

총 픽셀 수는 cumulative histogram[255]로 주어지므로 원하는 매핑 j = T [k]는

이 매핑에 의하면 입력 영상의 픽셀이 가지는 가장 큰 그레이 값은 출력 영상에는 255에 해당하게 된다. 출력 영상은 정의에 의해서 cumulative histogram이 직선적으로 증가하는 형태를 가지게 된다. 큰 영상의 경우 이 매핑을 Lookup 테이블로 저장하면 연산의 횟수를 줄일 수 있다.

note)
1. 0차 cumulative histogram을 이용하면 0부터 k까지의 histogram의 합은 한 번의 호출로 얻어진다.
2. 0차 cumulative histogram의 마지막 bin의 값은 전체 픽셀 수의 합이 된다.

// histogram의 합은 cumulative histogram을 구성하면 된다;
void HistogramEqualize(BYTE *src, int width, int height, BYTE* dst) {
    int histogram[256] = {0}, cumhist[256];
    for (int k = width * height; k-->0;) ++histogram[src[k]];
    
    // make 0-st order cumulative histogram;
    cumhist[0] = histogram[0];
    for (int k = 1; k < 256; k++) cumhist[k] = cumhist[k-1] + histogram[k] ;
    
    int map[256];
    double v = 255.0 / cumhist[255]; 	// cumhist[255] = number of pixels;
    for (int k = 0; k < 256; k++) {
    	int v = int(v * cumhist[k]);
        map[k] = v > 255 ? 255: v;
    }
    
    for (int k = width * height; k-->0;) dst[k] = map[src[k]];
};

평탄화 전의 히스토그램(붉은색) 및 누적 히스토그램(연두색)

평탄화 후의 히스토그램(붉은색) 및 누적 히스토그램(연두색)

• Forward transformation:
  $$X_{k} =  \sum_{n = 0} ^{N-1} x_n e^{ -2i  \pi k n /N }, ~\quad k=0,1,..., N-1;$$
• Backward transformation:
  $$x_n = \frac{1}{N}\sum_{k = 0} ^{N-1} X_k  e^{ + 2i  \pi k n /N }, ~\quad n=0,1,..., N-1;$$
• Butterworth Filter:
  $$H(x, y) = \frac{1}{1+\left( \frac{x^2 +y^2}{H_c^2 }\right)^N} ,\quad N=\text{order},~H_c =\text{cut-off};$$
• 주어진 수가 $2^n$인가?(kipl.tistory.com/84). 입력 데이터 개수가 2의 거듭제곱이 아니면 부족한 부분은 0으로 채워서 $2^n \ge \text{데이터 개수}$ 크기가 되도록 만든다. 

Cooley–Tukey FFT algorithm:  

/* 1차원 FFT:: dir=1 for forward, -1 for backward;
** x = real-component, and y = imaginary componet.
** nn = length of data: 2의 거듭제곱이어야 한다.
*/
int fft1d(int dir, int nn, double x[], double y[]) {
    /* Calculate the number of points */
    int m = 0;
    while (nn > 1) {
        nn >>= 1;  m++;
    }
    nn = 1 << m;
    
    /* Do the bit reversal */
    int i2 = nn >> 1;
    int j = 0;
    for (int i = 0; i < nn - 1; i++) {
        if (i < j) { // swap(x[i], x[j]), swap(y[i], y[j]);
            double tx = x[i], ty = y[i];
            x[i] = x[j];  y[i] = y[j];
            x[j] = tx;    y[j] = ty;
        }
        int k = i2;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
    
    /* Compute the FFT */
    double c1 = -1.0;
    double c2 = 0.0;
    int l2 = 1;
    for (l = 0; l < m; l++) {
        int l1 = l2;
        l2 <<= 1;
        double u1 = 1.0;
        double u2 = 0.0;
        for (int j = 0; j < l1; j++) {
            for (int i = j; i < nn; i += l2) {
                int i1 = i + l1;
                double t1 = u1 * x[i1] - u2 * y[i1];
                double t2 = u1 * y[i1] + u2 * x[i1];
                x[i1] = x[i] - t1;
                y[i1] = y[i] - t2;
                x[i] += t1;
                y[i] += t2;
            }
            double z = u1 * c1 - u2 * c2;
            u2 = u1 * c2 + u2 * c1;
            u1 = z;
        }
        c2 = sqrt((1.0 - c1) / 2.0);
        if (dir == 1) c2 = -c2;
        c1 = sqrt((1.0 + c1) / 2.0);
    }
    
    /* Scaling for forward transform */
    if (dir == 1) {
        for (i = 0; i < nn; i++) {
            x[i] /= double(nn);
            y[i] /= double(nn);
        }
    }
    return 1;
};

/* 2차원 FFT. data = (2 * nn) x mm 크기의 버퍼;
** 각 행은 nn개의 실수 성분 다음에 nn개의 허수 성분이 온다.
** nn, mm은 각각 행과 열은 나타내며, 2의 거듭제곱이어야 한다.
** isign = 1 for forward, -1 for backward. 
** copy = buffer of length 2*mm ;
*/
void fft2d(double* data, int nn, int mm, int isign, double* copy) { 
    int stride = nn << 1; //row_stride;
    /* Transform by ROWS for forward transform */
    if (isign == 1) {
        int index1 = 0;
        for (int i = 0; i < mm; i++) {
            // real = &data[index1], imag = &data[index1 + nn];
            fft1d(isign, nn, &data[index1], &data[index1 + nn]);
            index1 += stride;
        }
    }
    
    /* Transform by COLUMNS */
    for (int j = 0; j < nn; j++) {
        /* Copy pixels into temp array */
        int index1 = j ;
        int index2 = 0;
        for (i = 0; i < mm; i++) {
            copy[index2] = data[index1];
            copy[index2 + mm] = data[index1 + nn];
            index2++;
            index1 += stride;
        }
        
        /* Perform transform */
        fft1d(isign, mm, &copy[0], &copy[mm]);
        
        /* Copy pixels back into data array */
        index1 = j;
        index2 = 0;
        for (i = 0; i < mm; i++) {
            data[index1] = copy[index2];
            data[index1 + mm] = copy[index2 + mm];
            index2++;
            index1 += stride;
        }
    }
    
    /* Transform by ROWS for inverse transform */
    if (isign == -1) {
        int index1 = 0;
        for (i = 0; i < mm; i++) {
            fft1d(isign, nn, &data[index1], &data[index1 + nn]);
            index1 += stride;
        }
    }
}

void butterworth( double * data, int Xdim, int Ydim, 
                 int Homomorph, int LowPass,
                 double Power, double CutOff, double Boost);

// H * F;
void butterworth(double * data, int Xdim, int Ydim, 
                 int Homomorph, int LowPass,
                 double Power, double CutOff, double Boost) 
{
    double CutOff2 = CutOff * CutOff;
    /* Prepare to filter */
    int stride = Xdim << 1;  // width of a sigle row(real + imag);    
    int halfx = Xdim >> 1;
    int halfy = Ydim >> 1;
    double *rdata = &data[0];   //real part;
    double *idata = &data[Xdim];//imag part;
    /* Loop over Y axis */
    for (int y = 0; y < Ydim; y++) {
        int y1 = (y < halfy) ? y: y - Ydim;
        /* Loop over X axis */
        for (int x = 0; x < Xdim; x++) {
            int x1 = (x < halfx) ? x: x - Xdim;
            /* Calculate value of Butterworth filter */
            double filter;
            if (LowPass)
                Filter = (1 / (1 + pow((x1 * x1 + y1 * y1) / CutOff2, Power)));
            else if ((x1 != 0) || (y1 != 0))
                Filter = (1 / (1 + pow( CutOff2 / (x1 * x1 + y1 * y1), Power)));
            else
                Filter = 0.0;
            if (Homomorph)
                Filter = Boost + (1 - Boost) * Filter;
            /* Do pointwise multiplication */
            rdata[x] *= Filter;
            idata[x] *= Filter;
        };
        rdata += stride; //go to next-line;
        idata += stride; 
    }
};

Butterworth Low Pass filter 적용의 예(order=2, cutoff=20)

PowerSpectrum 변화

주어진 N개의 관측값을 k개의 가우시안 mixture로 모델링한 후 EM 알고리즘을 적용하여 모델을 결정하는 파라미터를 추정하고자 한다. 추정해야 할 모델 파라미터는 k개의 가우시안을 기술하는 평균값(μ)과 공분산 행렬(Covariance Matrix: Σ)과 이들 가우시안 분포의 mixing(α) 정도를 나타내는 파라미터이다.

다음은 이 과정을 수행하는 C++ class와 사용 예제, 및 결과를 보여준다.

//가우시안 kernel Class;
class GaussKernel2D {
    double mx, my;
    double sdx, sdy, sdxy;
    double weight;
public:
    std::vector<double> posterior;
    std::vector<double> wgauss;

    void init(std::vector<POINT> &pts, int sx, int sy, int w) {
        posterior.resize(pts.size());
        wgauss.resize(pts.size());
        weight = w ;
        // initialize model parameters(random하게 선택함-->try another method!!::
        // 주어진 데이터로 전체 분포의 범위와 중심을 알 수 있고, 이를 주어진 클래스 수만큼 임의로 
        // 분배하는 방식으로 초기조건을 설정하면 보다 안정적으로 동작한다)
        mx = sx * double(rand()) / RAND_MAX;
        my = sy * double(rand()) / RAND_MAX;
        sdx = sx / 4 + 1;
        sdy = sy / 4 + 1;
        mx += rand() % 100 ;
        my += rand() % 100 ;
        sdxy = 0;
    }
    double gauss2d(double x, double y) { 
        double varx = sdx * sdx;
        double vary = sdy * sdy;
        double det = varx * vary - sdxy * sdxy;
        double idet = 1.0 / det ;
        double dxx = (x - mx) * (x - mx);
        double dyy = (y - my) * (y - my);
        double dxy = (x - mx) * (y - my);
        return (1.0 / sqrt(det) / 6.28319) * exp(-0.5 * (dxx * vary \
                               + dyy * varx - 2. * dxy * dxy) * idet); 
    }
    void getParams(std::vector<POINT> &pts, double prior = 0) {
        double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
        for (int j = pts.size(); j-->0;) {
            double x = pts[j].x, y = pts[j].y;
            sx  += posterior[j] * x;
            sy  += posterior[j] * y;
            sxx += posterior[j] * x * x;
            syy += posterior[j] * y * y;
            sxy += posterior[j] * x * y;
        }
        double denorm = np * weight;
        mx = sx/denorm; my= sy / denorm;
        double devx = sxx / denorm - mx * mx ;
        if (devx <= 0) devx = 0.001;
        sdx = sqrt(devx);
        double devy=syy / denorm - my * my;
        if (devy <= 0) devy = 0.001;
        sdy = sqrt(devy);
        sdxy = sxy / denorm - mx * my;
        // if prior = non-zero -> weight = weight*(1-alpha)+alpha*prior; alpha=0.1?
    };
    // weight; // posterior;
    void estimate(std::vector<double> &px) {
        weight = 0;
        for (int j = px.size(); j-->0;) {
            posterior[j] = wgauss[j] / px[j];    
            weight += posterior[j];
        }
        weight /= px.size(); 
    } 
    //P(x|thetal) * prior;
    void setProb(std::vector<POINT> &pts) {
        for (int i = pts.size(); i-->0;)
            wgauss[i] = weight * gauss2d(pts[i].x, pts[i].y);
    }
    void Draw(CDC* pDC, DWORD color = RGB(0xFF, 0, 0)) {
        CPen pen0(PS_SOLID, 1, color);
        CPen* pOld = pDC->SelectObject(&pen0);
        drawCon(pDC, mx, my, sdx, sdy, sdxy);        // draw ellipses;
        pDC->Ellipse(mx - 2, my - 2, mx + 2, my + 2);// draw center of ellipse;
        pDC->SelectObject(pOld);            
    }
};

void em_main(std::vector<POINT> &pts) {
    GaussKernel2D kernel[NKERNELS];
    int nclass = NKERNELS;
    double weights[20] = {1};
    std::vector<double> px(pts.size());
    double wsum = 0 ;

    for (int i = 0; i < nclass; i++) {
        kernel[i].init(pts, 400, 400, weights[i]);
        wsum += weights[i];
    };    
#define MAX_ITER 50
    for (int iter = 0; iter < MAX_ITER; ++iter) {   
        for (i = px.size(); i-->0;) px[i] = 0;
        for (int k = 0; k < nclass; k++){
            GaussKernel2D &gker = kernel[k];
            gker.setProb(pts); 
            for (int i = px.size(); i-->0;)
                px[i] += gker.wgauss[i] ;
        }        
        for (k = 0; k < nclass; k++) {
            kernel[k].estimate(px);
            kernel[k].getParams(pts);
        }
        //또는 log-likelihood를 계산하여서 그 변화가 적으면 loop-끝내면 된다..
    }
}

//참고 : 아래의 데이터는 사전에 라벨링이 된 것이 아니다. 컬러링은 한번 계산한 후에 분포에 맞게 컬러를 조절하여서 다시 계산한 것이다.

 

f(y|θ);