열방정식은 매질에서 열이 전달되는 과정을 기술한다. 열은 온도차이가 있을 때 전달되는 에너지로 온도 분포 $u(\vec{r},t)$가 주어진 경우 주변으로 나가는 열에너지는 온도의 gradient $\nabla u$에 비례한다. 

$$ \text{heat flux: }\quad  \vec{j}= \kappa (u)\nabla u$$

$\kappa$는 매질에서 열이 얼마나 잘 전도는지를 표현하는 일반적으로 $u$에 의존할 수 있다. 균일하고 등방적인 매질인 경우는 $\kappa$가 일정한 상수에 해당한다. 국소적으로 단위부피당 단위시간당 빠져나가는 열에너지는 thermal flux의 divergence에 비례하고 이 값이 0이 아니면 그 지점의 온도 변화를 만든다. 따라서 온도는 다음과 같은 방정식을 만족한다.

$$ u_t = \nabla \cdot \vec{j}  = \nabla\cdot(\kappa \nabla u)$$

이 방정식은 물질의 확산 현상에도 적용이 가능한데 이 경우 $u$는 물질의 밀도분포에 해당할 것이다.

$\kappa$가 상수인 경우 이 방정식은 선형방정식으로 초기 온도분포 $u(x,t=0)=f(x)$가 주어지면 해는 사이즈가 $\sigma = \sqrt{2\kappa t}$인 가우시안 커널(heat kernel 또는 Greeen function)과 초기 온도분포의 convolution으로 표현할 수 있다. 1차원 선형 열방정식의 해는 구체적으로 다음과 같이 표현된다.

\begin{align} u(x,t) =  \int \frac{1}{\sqrt{4\pi \kappa t}} e^{- (x-y)^2 / 4\kappa t } f(y) dy = ( G_{\sigma} * f)(x,t)\\G_\sigma (x) = \frac{1}{\sqrt{2\pi\sigma} }  e^{- x^2/2\sigma^2} \end{align}

열방정식은 뜨거운 커피가 담긴 잔을 방에 놓았을 때 방 안의 온도가 어떻게 변할지도 알 수 있게 해 준다. 물론 경험적으로 커피의 온도는 내려가고 방안은 잔의 주변부터 점차로 데워져서 나중에서는 모든 부분이 일정한 온도를 가지는 평형상태에 도달한다. 이는  가우시안 커널의 크기가 $\sqrt{t}$에 비례하므로 시간이 흐를수록 모든 지점에서 온도는 일정한 값으로 수렴할 것임을 쉽게 예측할 수 있다.

이 열방정식을 이미지에 대해서도 적용할 수 있다. 이 경우 초기 온도분포는 주어진 이미지의 명암값에 해당한다. 열방정식의 우변을 이미지에 적용하면 Gaussian 블러링에 해당하고, 시간은 블러링의 순차적인 적용 횟수를 의미한다. 이미지가 블러링 되는 경우 영상이 담고 있는 정보는 점점 잃게 되는데 이를 방지하기 위해서는 영상의 에지 정보를 보존하는 방법으로 열방정식을 변형시켜 보자. 영상에서 에지영역에서는 gradient 값이 커지므로 이 지점에서는 커널사이즈를 작게 (열방정식 관점에서는 열전도도를 줄임) 만들면 블러링 효과가 줄어들게 된다. 이를 위해서 다음과 같은 열전도도 함수를 도입하자.

$$ \kappa (u) = \frac {1}{1+ |\nabla u|^2/\lambda^2 }$$

여기서 $\lambda$은 contrast 파라미터로 에지의 세기가 $|\nabla| \gg \lambda$ 영역에서는 블러링의 영향을 거의 받지 않게 되어 열방정식을 적용하더라도 에지에 대한 정보 손실은 거의 일어나지 않는다. 그러나 $\kappa $가 상수가 아닌 경우는 열방정식은 비선형이므로 해를 명시적으로 쓸 수 없고 오직 수치해석적으로만 구할 수 있다.

먼저 미분연산자(Sobel operator: $\lambda$의 설정은 미분연산자의 normalization을 고려해야 함)을 이용해서 이미지의 gradient을 구하면 flux $\vec{j}$을 구할 수 있고, 다시 $j_x$에 대해서 $x-$방향 미분, $j_y$에 대해서 $y-$방향 미분을 해서 $\vec{j}$의 divergence, 즉 이미지의 일반화된 Laplace 연산 결과를 얻을 수 있다. 이 결과에 시간간격을 곱한 양을 이전 이미지에 더하여 새로운 이미지를 얻는다. 이 과정을 반복적으로 수행하면 된다.

$$ u(x,y, t+\Delta t) = u(x, y, t) + \Delta t \times \nabla \cdot (\kappa (|\nabla u| ) \vec{j})$$

$\kappa$가 상수일 때 $t \to \kappa t$로 시간변수를 바꾸면 열방정식은 차원이 없는 형태로 되고, $\kappa <1$이므로 시간간격은 $\Delta t \ll 1$를 만족하도록 선택해야 한다.

 
 
 

 

$\lambda=50,~\Delta t=0.01$일 때 5 스텝 간격으로의 변화

 
double conductivity(double gmagsq, double lambdasq) {
    return 1.0/ (1.0 + gmagsq/lambdasq);
}
int LaplaceImage(double *img, int w, int h, double lambda, double *laplace) {
    const double lambdasq = lambda * lambda;
    // Sobel gradient 3x3
    const int nn[] = { -w - 1, -w, -w + 1, -1, 0, 1, w - 1, w, w + 1};
    const int sobelX[] = { -1, 0, 1, -2, 0, 2, -1, 0, 1};
    const int sobelY[] = { -1, -2, -1, 0, 0, 0, 1, 2, 1};
    std::vector<double> Gx(w * h, 0);
    std::vector<double> Gy(w * h, 0);
    std::vector<double> Gmag2(w * h, 0);
    const int xmax = w - 1, ymax = h - 1;
    // gradient-x, gradient-y, and mag of gradient.
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                double v = img[pos + nn[k]];
                sx += sobelX[k] * v;
                sy += sobelY[k] * v;
            }
            Gx[pos] = sx;
            Gy[pos] = sy;
            // gx^2 + gy^2;
            Gmag2[pos] = sx * sx + sy * sy;
        }
        pos++; // skip x=xmax;
    }
    // flux;
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double kappa = conductivity(Gmag2[pos], lambdasq);
            // multiply conductivity;
            Gx[pos] *= kappa;
            Gy[pos] *= kappa;
        }
        pos++; // skip x=xmax;
    }
    // divergence -> laplace;
    for (int y = 1, pos = w; y < ymax; y++) { // pos = w; starting address;
        pos++; //skip x=0;
        for (int x = 1; x < xmax; x++, pos++) {
            double sx = 0, sy = 0;
            for (int k = 0; k < 9; k++) {
                sx += sobelX[k] * Gx[pos + nn[k]];
                sy += sobelY[k] * Gy[pos + nn[k]];
            }
            laplace[pos] = sx + sy;
        }
        pos++; // skip x=xmax;
    }
    return 1;
}
void EpSmooth(BYTE *image, int w, int h, BYTE *out) {
    const double lambda = 50.0;
    const double dt = 0.01; //time step;
    const int maxiter = 100;
    std::vector<double> fimage(w * h, 0);
    std::vector<double> laplace(w * h, 0);
    for (int k = fimage.size(); k-- > 0;) fimage[k] = image[k];
    for (int iter = maxiter; iter-- > 0;) {
        LaplaceImage(&fimage[0], w, h, lambda, &laplace[0]);
        // update image;
        for (int k = fimage.size(); k-- > 0; ) 
            fimage[k] += dt * laplace[k];   
    }
    // preparing output;
    for (int k = fimage.size(); k-- > 0; ) {
        int a = int(fimage[k]);
        out[k] = a > 255 ? 255: a < 0 ? 0: a;
    }
}

 

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
728x90

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

Anisotropic Diffusion Filter  (0) 2024.02.23
Watershed Segmentation  (0) 2021.02.27
Local Ridge Orientation  (0) 2021.02.21
Contrast Limited Adaptive Histogram Equalization (CLAHE)  (3) 2021.02.15
Fixed-Point Bicubic Interpolation  (1) 2021.01.19
Posted by helloktk
,

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' 카테고리의 다른 글

Ellipse Fitting  (0) 2024.03.02
Bilateral Filter  (0) 2024.02.18
Canny Edge: Non-maximal suppression  (0) 2023.01.11
Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
Posted by helloktk
,

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

$$ \text{edge 조건}:~~|\nabla g|_1 > |\nabla g|_2  ~\text{and}~ |\nabla g|_1 > |\nabla 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를 선형보간을 사용해서 얻을 수 있다.

\begin{align} \tt |\nabla g|_B &\approx \tt se* \tan \theta + ee* (1-\tan \theta)~~~(\text{linear interpolation}) \\  &=  \tt se*\frac{gy}{gx} + ee* \left( 1-\frac{gy}{gx} \right)  \\  \tt |\nabla g|_C &\approx \tt nw* \tan \theta + ww* (1-\tan \theta) \\  &= \tt nw*\frac{gy}{gx} + ww* \left( 1-\frac{gy}{gx} \right) \end{align}

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

$$ \tt |\nabla g|_A \ge  |\nabla g|_B ~~\Leftrightarrow~~  gx * |\nabla g|_A \ge se* gy  + ee*(gx - gy) \\ \tt |\nabla g|_A \ge  |\nabla g|_C ~~\Leftrightarrow~~  gx * |\nabla g|_A \ge nw* gy  + ww*(gx - gy) $$ 

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

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

Bilateral Filter  (0) 2024.02.18
영상에서 Impulse Noise 넣기  (2) 2023.02.09
Statistical Region Merging  (0) 2022.06.11
Moment-preserving Thresholding  (0) 2022.05.29
Minimum Cross Entropy Thresholding  (0) 2022.05.29
Posted by helloktk
,

Statistical region merging은 이미지의 픽셀을 일정한 기준에 따라 더 큰 영역으로 합병하는 bottom-up 방식의 과정이다. 두 영역 $R_1$과 $R_2$가 하나의 영역으로 합병이 되기 위해서는 두 영역의 평균 픽셀 값의 차이가 

$$ g \sqrt{ \frac {\ln(2/\delta)}{2Q}  \Big( \frac {1}{|R_1|}+ \frac {1}{|R_2|}\Big) }$$

를 넘지 않아야 한다. $g=256$으로  gray 레벨의 갯수를 의미하고, $|R_i|$는 $R_i$ 영역에 포함된 픽셀 수를 나타낸다. $\delta$는 작은 수로 이미지의 픽셀 수의 제곱에 반비례한다. 보통 $\delta = 1/(6\times \text {width}\times \text {height})^2$로 선택한다. $Q$는 이미지의 통계적인 복잡성을 정량화하는 양으로 이 알고리즘에서는 외부에서 설정이 되는 값이다. 낮은 $Q$값을 선택하면 분할된 영역의 수가 작아지고(undersegmentation), 반대로 높은 $Q$ 값을 입력하면 분할된 영상에 너무 많은 영역이 나타나게 된다(oversegmentation).

Ref:

https://en.wikipedia.org/wiki/Statistical_region_merging

http://www.lix.polytechnique.fr/~nielsen/Srmjava.java

graph-based segmentation

 

Srm::Srm(int width, int height, BYTE *image) {
    this->width  = width;
    this->height = height;
    int n        = width * height;
    this->count = new int [n];
    this->Ravg  = new float [n];
    this->Gavg  = new float [n];
    this->Bavg  = new float [n];
    this->image = image;
    // disjoint sets with n pixels;
    this->UF = new Universe(n);
    // initialize to each pixel a leaf region;
    for (int i = 0, pos = 0; i < n; i++, pos += 3) {
        count[i] = 1;
        Bavg[i] = image[pos    ];
        Gavg[i] = image[pos + 1];
        Ravg[i] = image[pos + 2];
    }
    this->Q = 32;		// adjustable.					    
    this->g = 256.0;
    this->logDelta = 2. * log(6.0 * n);
}
bool Srm::Predicate(int reg1, int reg2) {
    double dR = (Ravg[reg1] - Ravg[reg2]); dR *= dR;
    double dG = (Gavg[reg1] - Gavg[reg2]); dG *= dG;
    double dB = (Bavg[reg1] - Bavg[reg2]); dB *= dB;
    double logreg1 = min(g, count[reg1]) * log(1.0 + count[reg1]);
    double logreg2 = min(g, count[reg2]) * log(1.0 + count[reg2]);
    double factor = g * g / (2.0 * Q);
    double dev1 = factor * (logreg1 + logDelta) / count[reg1] ;
    double dev2 = factor * (logreg2 + logDelta) / count[reg2] ;
    double dev = dev1 + dev2;
    return ( (dR < dev) && (dG < dev) && (dB < dev) );
}
void Srm::Merge(int root1, int root2) {
    if (root1 == root2) return;
    int w1 = count[root1], w2 = count[root2];
    int root = UF->Union(root1, root2);
    //update the merged region;
    count[root] = w1 + w2;
    double count_sum = w1 + w2;
    Ravg[root] = (w1 * Ravg[root1] + w2 * Ravg[root2]) / count_sum;
    Gavg[root] = (w1 * Gavg[root1] + w2 * Gavg[root2]) / count_sum;
    Bavg[root] = (w1 * Bavg[root1] + w2 * Bavg[root2]) / count_sum;
}
Edge* Srm::Pairs(int nedge) {
    // 4-connectivity;
    int ymax = height - 1, xmax = width - 1;
    Edge* edgeList = new Edge[nedge];
    int cnt = 0;
    for (int y = 0; y < ymax; y++) {
        for (int x = 0; x < xmax; x++) {
            int pos = y * width + x;
            int b1 = image[3 * pos + 0];
            int g1 = image[3 * pos + 1];
            int r1 = image[3 * pos + 2];
            //right: x--x
            edgeList[cnt].r1 = pos;     //current
            edgeList[cnt].r2 = pos + 1; //right
            int bdiff = abs(b1 - image[3 * (pos + 1) + 0]);
            int gdiff = abs(g1 - image[3 * (pos + 1) + 1]);
            int rdiff = abs(r1 - image[3 * (pos + 1) + 2]);
            edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff) ;
            //below: x
            //       |
            //       x
            edgeList[cnt].r1 = pos;
            edgeList[cnt].r2 = pos + width;
            bdiff = abs(b1 - image[3 * (pos + width) + 0]);
            gdiff = abs(g1 - image[3 * (pos + width) + 1]);
            rdiff = abs(r1 - image[3 * (pos + width) + 2]);
            edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
        }
    }
    //x=width-1;
    for (int y = 0; y < ymax; y++) {
        int pos = y * width + (width - 1); // (x,y) = (width-1, y)
        // x
        // |
        // x
        edgeList[cnt].r1 = pos;
        edgeList[cnt].r2 = pos + width;
        int bdiff = abs((int)image[3 * pos + 0] - image[3 * (pos + width) + 0]);
        int gdiff = abs((int)image[3 * pos + 1] - image[3 * (pos + width) + 1]);
        int rdiff = abs((int)image[3 * pos + 2] - image[3 * (pos + width) + 2]);
        edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
    }
    //y=height-1;
    for (int x = 0; x < xmax; x++) {
        int pos = (height - 1) * width + x;      //(x,y)=(x, height-1);
        //right; x--x
        edgeList[cnt].r1 = pos;
        edgeList[cnt].r2 = pos + 1;
        int bdiff = abs((int)image[3 * pos + 0] - image[3 * (pos + 1) + 0]);
        int gdiff = abs((int)image[3 * pos + 1] - image[3 * (pos + 1) + 1]);
        int rdiff = abs((int)image[3 * pos + 2] - image[3 * (pos + 1) + 2]);
        edgeList[cnt++].diff = max3(bdiff, gdiff, rdiff);
    }
    return edgeList;
}
int Srm::Segment() {
    // 4-connectivity 
    int nedge = 2 * (width - 1) * (height - 1) + (height - 1) + (width - 1);
    Edge* edgeList = Pairs(nedge);
    BucketSort(edgeList, nedge);
    for (int i = 0; i < nedge; i++) {
        int root1 = UF->Find(edgeList[i].r1);
        int root2 = UF->Find(edgeList[i].r2);
        if ((root1 != root2) && (Predicate(root1, root2)))
            Merge(root1, root2);
    }
    delete [] edgeList;
    int rgn_count = 0;
    for (int node = width * height; node-- > 0;)
        if (UF->IsRoot(node)) rgn_count++;
    return rgn_count;
}
// sorting with buckets; returns an ordered edgeList;
void BucketSort(Edge* &edgeList, int n) {
    int hist[256] = {0}, chist[256];
    for (int i = 0; i < n; i++) hist[edgeList[i].diff]++;
    // cumulative histogram
    chist[0] = 0;  // Note, chist[0] ne hist[0];
    for (int i = 1; i < 256; i++)
        chist[i] = chist[i - 1] + hist[i - 1];

    Edge *ordered = new Edge [n];
    for (int i = 0; i < n; i++)
        ordered[chist[pair[i].diff]++] = pair[i];        
    delete[] edgeList;
    edgeList = ordered;
}

 

 

728x90

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

영상에서 Impulse Noise 넣기  (2) 2023.02.09
Canny Edge: Non-maximal suppression  (0) 2023.01.11
Moment-preserving Thresholding  (0) 2022.05.29
Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Posted by helloktk
,

영상의 히스토그램($h[z]$)이 bimodal로 주어지는 경우 적절한 threshold 값을 선택해서 전경과 배경을 분리할 수 있다. 전경을 대표하는 픽셀 값을 $z_f$, 배경을 대표하는 픽셀 값을 $z_b$라면 이진화 후 정규화된 히스토그램은 

$$ \tilde{h}[z] =p_b \delta_{z, z_{b}} + p_f \delta_{z, z_{f}}$$

로 표현된다. $p_b$은 배경에 해당하는 픽셀 비율이고, $p_f$은 전경에 해당하는 픽셀 비율이다.

threshold 값을 어떻게 선택하면 이진화된 영상의 히스토그램이 원 영상의 히스토그램의 특성을 최대한 담게 할 수 있을까? 이에 대한 기준으로 두 히스토그램의 $n$차 moment가 같은 값을 갖도록 하자. 주어진 미지수가 $p_b$, $p_f$, $z_b$, $z_f$이 있으므로 최소한 4개의 moment가 같도록 만들어야 한다. 가장 낮은 찾수의 moment로부터 시작해서 3차까지 4개의 moments가 같다는 조건에서 아래의 식들을 얻을 수 있다.

$$ \text{0-차 moment: }  m_0 \equiv \sum_{z=0}^{255} h[z] = p_b  + p_f = 1$$

$$ \text{1-차 moment: }  m_1 \equiv \sum_{z=0}^{255}z h[z] = p_b z_b + p_f z_f$$

$$ \text{2-차 moment: }  m_2 \equiv \sum_{z=0}^{255}z^2 h[z] = p_b z_b^2 + p_f z_f^2$$

$$ \text{3-차 moment: }  m_3 \equiv \sum_{z=0}^{255}z^3 h[z] = p_b z_b^3+ p_f z_f^3$$

원 영상의 moment $m_0$, $m_1$, $m_2$, $m_3$을 계산해서 풀면

\begin{gather} c_0 = \frac{m_3 m_1 - m_2^2 }{ m_0 m_2 - m_1^2 } ,\quad  c_1 = \frac{m_1 m_2 - m_0 m_3}{ m_0 m_2 - m_1^2 } \\ z_b = \frac{1}{2} \left(-c_1 - \sqrt{ c_1^2 - 4 c_0} \right) \\  z_f = \frac{1}{2}\left( -c_1 + \sqrt{ c_1^2 -4 c_0} \right) \\ p_b = \frac{z_f - m_1}{z_f - z_b} \\ p_f = 1 - p_b\end{gather}

따라서 threshold 값

$$ \sum_{z=0}^{T-1} h[z] = p_b  $$

을 만족하는 $T$을 선택하면 된다.

 

Ref: W. Tsai, "Moment-preserving thresholding: a new approach," Computer Vision, Graphics, and Image Processing, vol. 29, pp. 377-393, 1985.

int MomentsPreseving_threshold(int histogram[256]) {
    int tot = 0;
    for (int i = 0; i < 256; i++)
        tot += histogram[i];
    //normalised histogram
    double hist[256];
    for (int i = 0; i < 256; i++)
        hist[i] = double(histogram[i]) / tot;
    /* moments calculation: zero moment is 1 by defintion*/
    double m0 = 1, m1 = 0, m2 = 0, m3 = 0;
    for (int i = 0; i < 256; i++ ) {
        double h = hist[i];
        m1 += i * h;
        m2 += i * i * h;
        m3 += i * i * i * h;
    }
    double det = m0 * m2 - m1 * m1;
    double c0 = (m1 * m3 - m2 * m2) / det;
    double c1 = (m2 * m1 - m3 * m0) / det;
    double zb = 0.5 * (-c1 - sqrt (c1 * c1 - 4.0 * c0));
    double zf = 0.5 * (-c1 + sqrt (c1 * c1 - 4.0 * c0));
    double pb = (zf - m1) / (zf - zb);  
    double s = 0;
    for (int i = 0; i < 256; i++) {
        s += hist[i];
        if (s > pb)
            return i; // threshold
    }
    return 0;
}
728x90

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

Canny Edge: Non-maximal suppression  (0) 2023.01.11
Statistical Region Merging  (0) 2022.06.11
Minimum Cross Entropy Thresholding  (0) 2022.05.29
Quadtree Segmentation  (0) 2022.05.21
Harris Corner Detector  (0) 2022.04.07
Posted by helloktk
,