Geometry & Recognition

Edge를 보존하는 low-pass filter로 주어진 픽셀 값을 그 위치에서 edge 방향을 따라 구한 평균값으로 대체한다. 방향은 8방향 (45도) 기준이다. Edge의 방향은 그림처럼 중심 픽셀을 기준으로 주변의 8개를 픽셀을 이용해서 45도씩 회전하면서  차이(평균 변화율)를 계산하여 가장 작은 값을 주는 방향으로 잡는다. 중심 픽셀이 line[1][1] 일 때

     0도 edge: abs(line[1][0] - line[1][2])가 최소

   45도 edge: abs(line[0][0] - line[2][2])가 최소

   90도 edge: abs(|ine[0][1] - line[2][1])가 최소

 135도 edge: abs(line[0][2] - line[2][0])가 최소 

Edge 방향의 3 픽셀 평균값으로 중심 픽셀 값을 대체하면 edge는 보존하면서 mean filter를 적용한 효과를 얻을 수 있다.

void smooth_ep_3x3(BYTE *src, int width, int height, BYTE* dst) {
    const int xend = width - 2, yend = height - 2;
    BYTE *line[3];
    line[0] = src; line[1] = line[0] + width; line[2] = line[1] + width;
    BYTE *dst_line = dst + width;        // starting dst row;
    for (int y = 0; y < yend; ++y) {
        for (int x = 0; x < xend; ++x) {
            int diff1 = line[1][x] - line[1][x + 2];
            if (diff1 < 0) diff1 = - diff1;
            int diff2 = line[0][x] - line[2][x + 2];
            if (diff2 < 0) diff2 = - diff2;
            int diff3 = line[0][x + 1] - line[2][x + 1];
            if (diff3 < 0) diff3 = - diff3;
            int diff4 = line[0][x + 2] - line[2][x];
            if (diff4 < 0) diff4 = - diff4;
            
            if (diff1 <= diff2 && diff1 <= diff3 && diff1 <= diff4)         //0-도
                dst_line[x + 1] = (line[1][x] + line[1][x + 1] + line[1][x + 2]) / 3;
            else if (diff2 <= diff3 && diff2 <= diff4)                      //45-도
                dst_line[x + 1] = (line[0][x] + line[1][x + 1] + line[2][x + 2]) / 3;
            else if (diff3 <= diff4)                                        //90-도;
                dst_line[x + 1] = (line[0][x + 1] + line[1][x + 1] + line[2][x + 1]) / 3;
            else                                                            //135-도;
                dst_line[x + 1] = (line[0][x + 2] + line[1][x + 1] + line[2][x]) / 3;
        }   
        dst_line += width; //move to next dst line;
        // increases src line ptr;
        BYTE *tptr = line[2] + width;
        line[0] = line[1]; line[1] = line[2]; line[2] = tptr;
    }
};

9번 반복 적용 결과: 지문처럼 규칙적인 패턴이 나오는 경우는 Gabor 필터를 이용하면 보다 좋은 결과를 얻을 수 있다.

iteration에 따른 correlation의 계산:

iter=0, corr w.r.t original=0.925144, corr w.r.t previous=0.925144
iter=1, corr w.r.t original=0.903661, corr w.r.t previous=0.985120
iter=2, corr w.r.t original=0.888224, corr w.r.t previous=0.994821
iter=3, corr w.r.t original=0.876582, corr w.r.t previous=0.997403
iter=4, corr w.r.t original=0.866254, corr w.r.t previous=0.998183
iter=5, corr w.r.t original=0.856620, corr w.r.t previous=0.998596
iter=6, corr w.r.t original=0.847365, corr w.r.t previous=0.998841
iter=7, corr w.r.t original=0.838400, corr w.r.t previous=0.998981
iter=8, corr w.r.t original=0.829703, corr w.r.t previous=0.999088

방 안에 뜨거운 물체가 놓이면 그 열은 시간이 지나면서 사방으로 퍼져나간다. 열원에서 주위로 열이 시간에 따라 전달이 되는 프로세스를 기술하는 방정식이 열(전도)방정식이다. 이것을 이미지 필터링에서도 적용이 가능하다. 이미지에 나타난 각각의 노이즈를 열원으로 보면 주변으로 열의 전달은 열원의 온도를 낮추게 되는데 이를 이미지에서 노이즈를 제거하는 filtering 과정으로 생각할 수 있다. 그러나 열방정식은 방향성이 없는 등방성 확산(isotropic diffusion process)을 기술한다. 이미지에 열전도 방정식을 바로 적용하면 이미지의 블러링과 유사한 효과를 만든다. filtering를 하더라도 에지를 보존해야 이미지의 정보 손실이 적게 된다. 에지 영역에서는 gradient의 공분산 행렬($\Sigma$)의 고윳값의 차이가 크고, 노이즈나 corner 근방에서는 두 고윳값의 크기가 비슷하다. 따라서, 에지를 보존하기 위해서는 공간으로 확산하는 역할을 하는 이미지의 Laplace 항 ($\nabla^2 I = I_{xx} + I_{yy}$)을 gradient의 공분산 행렬의 고유 벡터와 고유치, 그리고 이미지의 hessian 행렬 등을 조합하여서 불변인 형태를 만들어야 한다. 여러 가능한 후보군 중에 하나의 예를 들면

$$ I(t+dt) - I(t) = f_1(s) \text{Tr}({\bf  U U^T H}) + f_2(s) \text{Tr}({\bf V V^T H}) \\ {\bf H}=\begin{pmatrix} I_{xx} & I_{xy } \\ I_{xy} & I_{yy} \end{pmatrix}, \quad \nabla^2 I = \text{Tr} ({\bf H}) $$

여기서, $\bf U, V$는 그라디엔트 공분산 행렬의 두 고유 벡터이고, $s=\lambda_1 + \lambda_2$는 고유값의 합이다. $f_1$과 $f_2$의 함수 형태를 고유치가 큰 방향에서는 억제되고, 고유치가 작은 방향에서는 덜 억제가 되도록 적절히 선택하는 경우에 에지를 보존하는 필터링의 효과를 얻을 수 있다.

아래의 예는 이 식을 써서 300번 iteration을 한 결과이다. 원본 이미지는 첨부한 논문에서 추출하였다.

https://kipl.tistory.com/563

이미지를 획득할 때 조명 조건이 균일하지 않으면 이미지 처리 과정에서 복잡성이 증가하게 된다. 이런 경우에 배경의 불균일성을 제거하는 전처리 과정을 수행할 필요가 생기게 된다. 일반적으로 조명의 불균일성은 공간적으로 넓은 범위에서 천천히 일어나게 되는데, 이를 이용하면 이미지에서 배경을 좀 더 수월하게 제거할 수 있다. 보통 큰 윈도를 (전경의 물체의 크기에 비해서) 가지는 median filter를 연속적으로 적용하여 배경을 추출할 수 있지만, median filter는 일반적으로 계산 비용이 상당히 비싸다.

좀 더 저렴한 연산비용이 드는 방법으로는 이미지가 각 픽셀 좌표에서 하나의 값을 할당하므로 3차원 공간에서 곡면으로 생각할 수 있다는 점을 이용하면 된다. 일정하게 변하는 조명에 의한 배경은 3차원에서 평면을 형성하는 것으로 근사할 수 있다. 주어진 이미지 곡면에서 배경 평면을 추출하는 쉬운 방법은 최소자승법을 사용하여 이미지 곡면을 가장 잘 근사하는 평면을 찾는 것이다. 그러나 이 방법은 전경의 물체가 fitting에 있어 일종의 outlier 역할을 하게 된다. 전경에 해당하는 물체가 많은 이미지에서도 잘 동작하려면 보다 robust 한 fitting 알고리즘을 고안해야 한다.

3차원 공간에서 이미지 곡면을 고려할 때, 전경(밝은 색)이 이미지 영역에서 차지하는 비중이 작은 경우에는 그림처럼(A) 충분히 큰 공을 곡면 아래쪽에서 접하게 굴리면 공의 표면이 형성하는 곡면은 이미지의 작은 요철(전경 물체)에 영향을 거의 받지 않는다. 따라서 구르는 공에 의해 형성된 접곡면은 이미지의 배경으로 생각할 수 있다 (B는 전경이 어두운 경우).

이미지 곡면이 요철이 심하더라도 반지름이 충분히 큰 공을 사용하게 되면 부드러운 배경 표면을 얻을 수 있다. 공 표면과 이미지 표면과의 접촉 관계는 공의 표면에 놓인 patch의 좌표점에서 (미리 계산) 높이를 이미지의 픽셀 값과 비교하여 얻는다. 공의 반경이 큰 경우에는 연산을 줄이기 위해서 원본 이미지를 축소해서 작업을 하고, 다시 간단한 보간을 써서 원본 크기로 복원하면 충분하다.

다음 예는 밝은 배경에 어두운 전경인 이미지에 대해서 이 알고리즘을 적용한 경우다. 먼저, 원본 이미지의 negative 이미지를  만들어서 (어두운 배경에 밝은 전경) rolling ball transformation을 적용해서 배경을 얻었다.

 negative image:

radius=16;

radius=5:

배경을 뺀 결과(radius=16):

#define SQR(x) ((x)*(x))
struct RollingBall {
    std::vector<BYTE> data;
    int patchwidth;
    RollingBall(int radius) {
        int artrimper;
        if (radius <= 30)
            artrimper = 12; // trim 24% in x and y
        else if (radius <= 100) 
            artrimper = 16; // trim 32% in x and y
        else 
            artrimper = 20; // trim 40% in x and y
        
        Build(radius, artrimper);
    }
    void Build(int ballradius, int artrimper) {
        if (ballradius < 1) ballradius = 1;
        const int diam = ballradius << 1;
        const int xtrim = (artrimper * diam) / 100;
        const int hwpatch = ballradius - xtrim;
        patchwidth = (hwpatch<<1) + 1;			// make odd number;
        data.resize(patchwidth * patchwidth, 0);

        const int rsquare = SQR(ballradius);
        for (int k = data.size(); k-->0;) {
            int x = k % (patchwidth) - hwpatch; //relative to patch center-x.
            int y = k / (patchwidth) - hwpatch; //relative to patch center-y;
            int tmp = rsquare - SQR(x) - SQR(y);
            if (tmp > 0) data[k] = int(sqrt(double(tmp)));
        }
    }
};

// processing을 빠르게 하기 위해서, ball의 반경이 큰 경우에는 원본이미지를 
// 축소하여서(2,4,8,..) 처리한 후에, 처리 안된 점들은 interpolation을 하여서
// bacground image을 생성하면 된다.
void RollBall(const RollingBall& ball,                  //rolling ball object;
              BYTE **image, int width, int height,      //source image;
              BYTE **background)                        //background; 
{
    const int hwpatch = ball.patchwidth >> 1;         
    const BYTE *patch = &ball.data[0] ;         //pre-calculated patch z-value;
    int zcenter = 0;                            //start z-center in the xy-plane
    for (int y = 0; y < height; y++) {
        int top = y - hwpatch;
        int bot = y + hwpatch;
        for (int x = 0; x < width; x++) {
            // 현패치 내에서(중심(x-hwpatch, y-hwpatch)) 픽셀의 z값(그레이 레벨)과 
            // zcenter만큼 shift한 공의 윗표면에 놓인 패치의 z와 차이를 구함. 차이가
            // 가장 작을 때 값만큼 공의 중심을 높이거나 내리면 공과 픽셀 surface가 밑에서 
            // 접하게 된다. 
            int zmin   = 255; 
            int left  = x - hwpatch;
            int right = x + hwpatch;
            for (int iy = top, ball_pos = 0; iy <= bot; iy++) 
                for (int ix = left; ix <= right; ix++, ball_pos++) 
                    if (ix >= 0 && ix < width && iy >= 0 && iy < height) {
                        int zdif = image[iy][ix] - (zcenter + (patch[ball_pos]));
                        if (zdif < zmin) zmin = zdif;
                    }

            if (zmin != 0) zcenter += zmin; 
            for (int iy = top, ball_pos = 0; iy <= bot; iy++) 
                for (int ix = left; ix <= right; ix++, ball_pos++)
                    if (ix >= 0 && ix < width && iy >= 0 && iy < height) {
                       int zadd = zcenter + patch[ball_pos];  
                       if (zadd > background[iy][ix])           
                           background[iy][ix] = (BYTE)zadd;
                    } 
        } 
    }
}

Mean-shift는 확률 밀도 함수의 극대점(mode)을 연속적으로 찾는 방법이다. 주어진 데이터(특징점)가 샘플링된 확률 밀도 함수의 평균과 분산을 알 수 없을 때 적용할 수 있는 non-parametric 알고리즘이다. Clustering, tracking, 또는 filtering 등에 사용될 수 있다.
Mean-shift filter는 각 픽셀을 중심으로 하는 윈도 내의 모든 픽셀에 대해서 특징점 공간(= 픽셀 좌표 공간 + 픽셀 컬러 공간)에서 밀도가 극대인 특징점(특징점 공간에서 질량중심)을 찾고, 이 특징점에 해당하는 픽셀의 컬러로 윈도 중심 픽셀의 컬러를 대체한다. 

아래 예제는 컬러 공간을 RGB에서 YIQ로 바꾸어서 알고리즘을 적용한 후에 그 결과를 RGB로 다시 변환하였다. 윈도를 크게 하면 알고리즘의 수행 속도가 매우 느려진다. 공간 윈도를 원에서 정사각형으로 바꾸어도 된다.

Clustering 알고리즘으로 볼 때 Mean-Shift 알고리즘의 특징은

1. 단순한 알고리즘 구조를 갖지만 연산 비용이 매우 크다. 고차원 특징에 대해서는 부적합
2. 클러스터 수에 대한 제한이 없다.(k-means와 비교)
3. 클러스터 초기 위치 설정이 필요없다.(k-means와 비교)
4. outliers의 영향이 적다.(k-means와 비교)
5. 유일한 매개변수는 윈도 크기다. 원도 크기를 미리 설정해야 한다.

// 구현 code

void RGB2YIQ(BYTE *color, int w, int h, float *ycomp, float *icomp, float *qcomp);

더보기

void RGB2YIQ(BYTE *color, int w, int h, float *ycomp, float *icomp, float *qcomp) {
    for (int i = w * h; i-- > 0;) {
        int b = color[3*i+0], g = color[3*i+1], r = color[3*i+2];
        ycomp[i] = 0.299f  *r + 0.587f *g + 0.114f  *b; // Y
        icomp[i] = 0.5957f *r - 0.2744f*g - 0.3212f *b; // I
        qcomp[i] = 0.2114f *r - 0.5226f*g + 0.3111f *b; // Q
    }
}

// YIQ2RGB color conversion;

void YIQ2RGB(float Yc, float Ic, float Qc, int *b, int *g, int *r) ;

더보기

void YIQ2RGB(float Yc, float Ic, float Qc, int *b, int *g, int *r) {
     *r = int(Yc + 0.9563f*Ic + 0.6210f*Qc);
     *g = int(Yc - 0.2721f*Ic - 0.6473f*Qc);
     *b = int(Yc - 1.1070f*Ic + 1.7046f*Qc);
}

#define SQR(x) ((x) * (x))
void MeanShiftFilterRGB(BYTE *image/*color*/, int width, int height,
                    int rad/*=spatial*/, int radcol/*=color space(YIQ)*/,
                    BYTE *out) {
    std::vector<float> ycomp(width*height) ;
    std::vector<float> icomp(width*height) ;
    std::vector<float> qcomp(width*height) ;
    int rad2    = rad * rad ;
    int radcol2 = radcol*radcol;
    // RGB2YIQ;
    RGB2YIQ(image, width, height, &ycomp[0], &icomp[0], &qcomp[0]);
    CRect rect(0, 0, width, height);
    rect.DeflateRect(0, 0, 1, 1);
    for (int y=0, pos=0; y<height; y++) {
        for (int x=0; x<width; x++, pos++) {
            int xc = x, yc = y;
            float Yc = ycomp[pos], Ic = icomp[pos], Qc = qcomp[pos];
            int iters = 0;
            while (1) {
                int oldxc = xc, oldyc = yc;
                float YcOld = Yc, IcOld = Ic, QcOld = Qc;
                float xsum = 0, ysum = 0; 
                float sumY = 0, sumI = 0, sumQ = 0;
                int count = 0;
                CRect rc = CRect(xc-rad, yc-rad, xc+rad, yc+rad) & rect;
                for (int y2 = rc.top; y2 <= rc.bottom; y2++) {
                    int ry = y2 - yc;
                    for (int x2 = rc.left; x2 <= rc.right; x2++) {
                        int rx = x2 - xc;
                        if ((SQR(ry) + SQR(rx)) <= rad2) {    
                            int curpos = y2*width + x2;
                            float Y2 = ycomp[curpos];
                            float I2 = icomp[curpos];
                            float Q2 = qcomp[curpos];
                            float dY = Y2 - Yc, dI = I2 - Ic, dQ = Q2 - Qc;
                            if ((SQR(dY) + SQR(dI) + SQR(dQ)) <= radcol2) {
                                xsum += x2; ysum += y2;
                                sumY += Y2; sumI += I2; sumQ += Q2;
                                count++;
                            }
                        }
                    }
                }
                if (count==0) count = 1;
                //average(YIQ);
                Yc = sumY / count; Ic = sumI / count; Qc = sumQ / count;
                // average(x,y);
                xc = int(xsum / count + 0.5); 
                yc = int(ysum / count + 0.5);
                int dx = xc - oldxc, dy = yc - oldyc;
                float dY = Yc - YcOld, dI = Ic - IcOld, dQ = Qc - QcOld;
                // shift;
                float shift = SQR(dx) + SQR(dy) + SQR(dY) + SQR(dI) + SQR(dQ); 
                if (shift <= 3 || ++iters >= 100) break;
            }
            int r, g, b;
            YIQ2RGB(Yc, Ic, Qc, &b, &g, &r);
            out[3*pos+0] = b > 255 ? 255: b < 0 ? 0: b;
            out[3*pos+1] = g > 255 ? 255: g < 0 ? 0: g;
            out[3*pos+2] = r > 255 ? 255: r < 0 ? 0: r;
        }       
    }
};

(rad=6, radcol=50)인 예.

YIQ space에서 컬러분포:

(rad=6, radcol=20) 인 예.

아래 이미지의 검정색 부분은 unstable point이다. 즉, 이 점들은 항상 다른 점으로 이동한다. 에지점들은 unstable 함이 극명하게 보인다.