Loading [MathJax]/jax/output/CommonHTML/jax.js

주어진 점집합을 기술하는 직선을 얻는 방법 중에 하나로 각각의 점들이 직선에서 벗어난 거리의 제곱을 더한 값(square of residual)을 최소화시키는 기울기와 절편을 찾는 최소자승법(least square method: linear regression)이 있다. 점집합이 {(xi,yi}를 직선 y=ax+b로 피팅을 하는 경우 직선에서 벗어난 정도(residual)는 직선까지의 거리를 사용할 수도 있고 또는 y 값의 차이를 이용할 수도 있다. 먼저 해를 closed form으로 쓸 수 있는 y 값의 offset을 residual로 사용하자. 그러면 fitting error는 각 점에서 residual의 제곱의 합으로 주어진다. 

R2(a,b)=i|yi(axi+b)|2

여기서 주어진 점집합의 moment를 Sxx=x2i, Syy=y2i, Sxy=xiyi, Sx=xi, Sy=yi로 놓으면 

R2(a,b)=Syy+a2Sxx+nb22aSxy+2abSx2bSy

이다. 주어진 점집합을 피팅하는 직선의 파라미터는 타원의 방정식으로 주어짐을 알 수 있고, 주워진 타원상의 임의의 파라미터에 해당하는 직선의 fittign error는 항상 일정한 값을 가짐을 알 수 있다.

좌표의 원점을 각 성분의 평균값만큼 이동하면 Sx=Sy=0이 되어 더 식이 단순해진다. 분산
σ2x=1nSxxˉx2  σ2x=1nxx

σ2y=1nSyyˉy2  σ2y=1nyy 공분산

cov(x,y)=1nSxyˉxˉycov(x,y)=1nxy

을 써서 표현하면,

R2=nσ2xa22ncov(x,y)a+nb2+nσ2y=nσ2x(acov(x,y)σ2x)2+nb2+ncov(x,y)2σ2x+nσ2yncov(x,y)2σ2x+nσ2y 

따라서 residual을 최소로 만들어 주는 직선의 기울기 ay절편 b

a=cov(x,y)σ2xb=0

로 주어지는데, 이 직선은 원점을 통과하는 직선이 된다. 원래의 좌표계로 돌아가면 기울기는 원점의 이동에 무관하므로 변화가 없고 직선이 (ˉx,ˉy)을 통과해야 하므로 절편 b값은

b=ˉyaˉx

로 주어진다. 상관계수 r(x,y)를 이용하면 Fitting이 잘된 정도를 정량적으로 표현이 된다.

r=cov(x,y)σ2xσ2y

즉,

R=nσ2y(r2+1)

로 주어진다.

728x90

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

Ellipse Fitting  (0) 2024.03.02
Bilateral Filter  (0) 2024.02.18
영상에 Impulse Noise 넣기  (2) 2023.02.09
Canny Edge: Non-maximal suppression  (0) 2023.01.11
Statistical Region Merging  (0) 2022.06.11
,

p=0.2, imax=255, imin=0

int AddImpulseNoise(CRaster& raster, double p, int imax, int imin, CRaster& noised) {
    srand(unsigned(time(NULL)));
    CSize sz = raster.GetSize();
    p = max(min(p, 1),0); //noise probability;
    int np = int(p * sz.cx * sz.cy);
    // clone;
    noised = raster;
    imax = max(min(255, imax), 0); //maximum impulse value;
    imin = min(max(0, imin), 255); //minimum impulse value;
    for (int count = 0, turn = 0; count < np; count++) {
        //x in [0, sz.cx - 1];
        int x = int((sz.cx - 1) * double(rand()) / RAND_MAX); 
        //y in [0, sz.cy - 1];
        int y = int((sz.cy - 1) * double(rand()) / RAND_MAX); 
        if (turn == 0) {
            turn = 1;
            noised.SetPixel(x, y, imax);
        } else {
            turn = 0;
            noised.SetPixel(x, y, imin);
        }
    }
    return 1;
}
728x90

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

Bilateral Filter  (0) 2024.02.18
파라미터 공간에서 본 최소자승 Fitting  (0) 2023.05.21
Canny Edge: Non-maximal suppression  (0) 2023.01.11
Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
,

Canny 알고리즘은 영상에서 edge을 찾을 수 있게 설계된 알고리즘으로 여러 단계의 영상 처리 과정으로 구성되어 있다. 일반적으로 edge 검출기는 노이즈에 매우 민감하기 때문에 low pass filter인 Gaussian filter를 적용하여 영상의 노이즈를 제거한다. edge에 해당하는 위치에서는 edge에 수직인 방향으로 gray 값의 변동이 크다는 사실을 이용해서 edge를 찾을 수 있다. 이를 위해서 영상의 gradient field g=(gx,gy) 를 구하고, 이를 이용해서 gradient의 크기(|g|=g2x+g2y)와 방향 정보(θ=tan1(gy/gx))를 얻는다. 대부분의 경우 edge 근방에서 gradient의 크기는 연속적으로 변하므로 edge가 두꺼운 형태가 된다. edge를 한 픽셀 두께로 줄이기 위해서는 gradient 크기의 극대점들을 찾는 추가적인 과정을 거쳐야 한다. 이 추가 처리를 non-maximal suppression이라 하는데, 이는 극대점에서 gradient 방향 또는 그 반대 방향으로 움직이는 경우 gradient 크기가 감소한다는 점을 이용한다. 아래 그림에서 1번 지점이 edge 위에 있으면 gradient 방향(edge 방향에 수직임)의 앞/뒤 지점인 2와 3에서 크기보다 커야 한다.

edge 조건:  |g|1>|g|2 and |g|1>|g|3

보통 gradient field의 방향을 수평(0, 180도), 수직(90도, 270도), 45도(225도), 135도(325도)로 나누고, 해당 방향의 인접 픽셀에서 gradient 크기를 비교하여 극대점을 판별한다. 주어진 픽셀의 8 방향 근방에서의 gradient  세기는 각각 nn, nw, ww, sw, ss, se, ee, ne로 정의했다. 예를 들면 gradient field의 방향이 거의 수평인 경우(방향이 -22.5~22.5, -180~-157.5, 157.5~180 범위에 들어가는 경우다. 이 경우 edge 방향은 거의 수직) 극대가 되려면  수평 방향의 최근방에서 크기 ee나 ww 보다 더 커야 한다.

void NonmaxSuppress(int radius, int width, int height, 
    std::vector<double> &gradx, std::vector<double> &grady,
    std::vector<int> &edgemap) {
    const double rad2deg = 45.0 / atan(1.0);
    std::vector<double> mag(width * height, 0);
    std::vector<double> dir(width * height, 0);
    int startx = radius, endx = width - radius;
    int starty = width * radius, endy = width * (height - radius); 
    for (int x = startx; x < endx; x++) 
        for (int y = starty; y < endy; y += width) {
            int pos = x + y;
            double gx = gradx[pos, gy = grady[pos];
            mag[pos] = hypot(gx, gy);
            dir[pos] = rad2deg * atan2(gy, gx);
        }
    for (int x = startx; x < endx; x++) { // non-maximal supression
        for (int y = starty; y < endy; y += width) {
            int pos = x + y;            
            double gmag = mag[pos];
            double ang = dir[pos];
            double ww = mag[pos-1];
            double ee = mag[pos+1];
            double nn = mag[pos-width];
            double ss = mag[pos+width];
            double ne = mag[pos-width+1];
            double se = mag[pos+width+1];
            double sw = mag[pos+width-1];
            double nw = mag[pos-width-1];          
            bool bedge = false;
            if ((0 <= ang && ang < 22.5) || (157 <= ang && ang <= 180)
                || (-22.5 <= ang && ang < 0) || ( -180 <= ang < -157.5))
                bedge = (gmag >= ee) && (gmag >= ww); // 수평;
            else if ((22.5 <= ang && ang < 67.5) || (-157.5 <= ang && ang > -112.5))
                bedge = (gmag >= se) && (gmag >= nw); // 45도
            else if ((67.5 <= ang && ang < 112.5) || (-112.5 <= ang && ang > -67.5))
                bedge = (gmag >= nn) && (gmag >= ss); // 수직
            else 
                bedge = (gmag >= sw) && (gmag >= ne); // 135도
            edgemap[pos] = bedge ? int(MAG_SCALE * gmag): 0;
        }
    }
}

gradient field의 방향 정보를 이용한 non-maximal suppression은 atan() 함수의 호출로 인해서 연산 비용이 많이 든다. 이를 해소하기 위해서 gradient 방향을 양자화하지 말고 직접 gradient  방향으로 한 픽셀 근방에서 gradient의 크기를 주변 8개 픽셀에서의 gradient 값(nn, nw, ww, sw, ss, se, ee, nw)을 이용해서 구하는 방법을 사용하자. 이 방법을 사용하면 atan()를 호출할 필요가 없어진다. 예를 들면 gradient의 방향이 0에서 45도 사이에 있다고 하자. 그러면 A 점이 극대가 되려면 A에서  gradient 방향의 연장선 상에 있는 한 픽셀 정도 떨어진 지점 B와 C에서 gradient 크기보다 더 커야 한다. B와 C에서 gradient 크기는 주변 8 방향 중에서 가장 가까이 있는 픽셀에서의 크기 ee와 se를 선형보간을 사용해서 얻을 수 있다.

|g|Bsetanθ+ee(1tanθ)   (linear interpolation)=segygx+ee(1gygx)|g|Cnwtanθ+ww(1tanθ)=nwgygx+ww(1gygx)

극대점 판별을 다음과 같이 쓰면 나눗셈 연산도 피할 수 있다.

|g|A|g|B    gx|g|Asegy+ee(gxgy)

|g|A|g|C    gx|g|Anwgy+ww(gxgy) 

void NonmaxSuppress(int radius, int width, int height, 
    std::vector<double> &gradx, std::vector<double> &grady,
    std::vector<int> &edgemap) {
    std::vector<double> mag(width * height, 0);
    int startx = radius, endx = width - radius;
    int starty = width * radius, endy = width * (height - radius);
    for (int x = starx; x < endx; x++) 
        for (int y = 0; y < endy; y += width) {
            int pos = x + y;
            mag[pos] = hypot(gradx[pos], grady[pos]);
        }
    for (int x = startx; x < endx; x++) { // non-maximal supression
        for (int y = starty; y < endy; y += width) {
            int pos = x + y;            
            double gx = gradx[pos];
            double gy = grady[pos];
            double gmag = mag[pos];
            double ww = mag[pos-1];
            double ee = mag[pos+1];
            double nn = mag[pos-width];
            double ss = mag[pos+width];
            double ne = mag[pos-width+1];
            double se = mag[pos+width+1];
            double sw = mag[pos+width-1];
            double nw = mag[pos-width-1];
            bool bedge = false;
            if (gx * gy > = 0) {//1,3 분면
                if (abs(gx) >= abs(gy)) {//0~45,180~225;
                    double tmp = abs(gx*gmag);
                    bedge = (tmp >= abs(gy*se + (gx-gy)*ee)) && (tmp > abs(gy*nw + (gx-gy)*ww));
                } else {//45~90,225~270;
                    double tmp = abs(gy*gmag);
                    bedge = (tmp >= abs(gx*se + (gy-gx)*ss)) && (tmp > abs(gx*nw + (gy-gx)*nn));
                }
            } else {//2,4 분면
                if (abs(gx) >= abs(gy)) {//135~180, 315~360;
                    double tmp = abs(gx*gmag);
                    bedge = (tmp >= abs(gy*ne - (gx+gy)*ee)) && (tmp > abs(gy*sw - (gx+gy)*ww));
                } else {//90~135, 270~315;
                    double tmp = abs(gy*gmag);
                    bedge = (tmp >= abs(gx*ne - (gy+gx)*nn)) && (tmp > abs(gx*sw - (gy+gx)*ss)); 
                }   
            }
            edgemap[pos] = bedge ? int(MAG_SCALE * gmag): 0;
        }
    }
}

728x90
,