코너 검출에서 가장 많이 쓰이는 캐니 에지 디텍션(canny edge detection)에 대해 알아보자
캐니 에디 검출이란?
The Canny edge detector is an edge detection operator that uses a multi-stage algorithm to detect a wide range of edges in images. It was developed by John F. Canny in 1986.
대부분 에지 추출 마스크는 잡음에 대해 민감하므로, 작은 잡음이라도 그것을 에지로 간주하여 추출하는 경우가 많다. 이러한 단점을 보완하는 캐니 마스크를 이용한 에지 추출 기법이 있는데, 실제로 잡음에 민감하지 않게 하며 강한 에지를 추출하는 것에 목적을 둔다.
Jogn Canny는 에지 추출기에 있어서 세 가지 기준을 만족한다는 제안을 내놓았다.
첫째, 탐지성(good detection) : 모든 실제 에지를 탐지하는 능력
둘째, 국부성(good localization) : 실제 에지와 탐지된 에지의 차이를 최소화
셋째, 응답성(clear response) : 각 에지에 대하여 단일한 응답
세 가지 조건을 만족하는 캐니 마스크를 이용한 기법은 여러 가지가 있으며,
기본적인 알고리즘은.
Canny edge detection is a four step process.
- A Gaussian blur is applied to clear any speckles and free the image of noise.
- A gradient operator is applied for obtaining the gradients' intensity and direction.
- Non-maximum suppression determines if the pixel is a better candidate for an edge than its neighbours.
- Hysteresis thresholding finds where edges begin and end.
1. 가우시안 스무딩 필터링
2. X, Y축으로 기울기 계산(gradient), 기울기의 강도를 계산
3. Non-maximum suppression 수행 :
( ->그레이 스케일(단색)에서 모든 픽셀에 대해서 각각픽셀의 세로/가로 소벨미분 각의 방향으로 자기자신과 각(degree)방향에 있는 픽셀이 일정한 threshold값을 넘으면 그 픽셀을 살리고 아니면 다 죽인다.결과적으로 얇은 선이 나온다. )
4. Hystersis 수행.
cannyEdge.exe
lena.raw
///////////////////////////////////////////////////////////////////////////////
// Canny Edge Detection (Procedural Thinning Edges Algorithm)
// 1) smooth image using Gaussian filter
// 2) compute gradient X & Y
// 3) compute magnitude of gradient X & Y
// 4) non-maximal suppression (thinning ridge)
// 5) hysteresis thresholding
//
// in: the pointer to input image data (8-bit gray-scale image)
// x: image width
// y: image height
// sigma: std deviation, it determines the kernel size
// tLow: low threhold value for hysteresis
// tHigh: ligh treshold value for hysteresis
///////////////////////////////////////////////////////////////////////////////
void cannyEdge(unsigned char *in, int x, int y, float sigma, float tLow, float tHigh, unsigned char *out)
{
// check params
if(!in || !out) return;
if(x <= 0 && y <= 0) return;
short *deltaX, *deltaY; // store horizontal and vertical gradient
unsigned short *mag; // store magnitude of gradient
unsigned char *nms; // store non-maximal suppression result
unsigned char *buf; // store temporal result
// allocate memory space
deltaX = new short[x*y];
deltaY = new short[x*y];
mag = new unsigned short[x*y];
nms = new unsigned char[x*y];
buf = new unsigned char[x*y];
if(!deltaX || !deltaY || !mag || !nms || !buf)
{
delete [] deltaX;
delete [] deltaY;
delete [] mag;
delete [] nms;
delete [] buf;
cout << "ERROR: Memory Allocation Failed.\n";
return;
}
// 1) smoothing
if(sigma != 0) // skip smoothing when sigma = 0.
gaussian(in, x, y, sigma, buf);
// 2) compute derivative X & Y
gradient(buf, x, y, deltaX, deltaY);
// 3) compute magnitude of dx & dy
magnitude(deltaX, deltaY, x, y, mag);
// 4) non-maximal suppression
nonMaxSupp(mag, deltaX, deltaY, x, y, nms);
// 5) hysteresis thresholding
// RECOMMENDATION: 0.8 <= tHigh <= 0.9, 0.3 <= tLow <= 0.5
hysteresis(nms, x, y, tLow, tHigh, out);
// free up memory
delete [] deltaX;
delete [] deltaY;
delete [] mag;
delete [] nms;
delete [] buf;
}
///////////////////////////////////////////////////////////////////////////////
// smooth image with Gaussian filter
///////////////////////////////////////////////////////////////////////////////
void gaussian(unsigned char *inBuf, int x, int y, float sigma, unsigned char *outBuf)
{
int kernelSize;
float *kernel;
// determine size of kernel (odd #)
// 0.0 <= sigma < 0.5 : 3
// 0.5 <= sigma < 1.0 : 5
// 1.0 <= sigma < 1.5 : 7
// 1.5 <= sigma < 2.0 : 9
// 2.0 <= sigma < 2.5 : 11
// 2.5 <= sigma < 3.0 : 13 ...
kernelSize = 2 * int(2*sigma) + 3;
// allocate memory space for kernel
kernel = new float[kernelSize];
if(!kernel)
{
cout << "ERROR: Memory Allocation Failedfor kernel.\n";
return;
}
// generate 1-D gaussian kernel
makeGaussianKernel(sigma, kernel, kernelSize);
// smoothing image with Gaussian filter
// it is convolution of image array and kernel
// Use 2D separable convolution because Gaussian Filter is seperable
// Horizontal and vertical kernel are same for Gaussian kernel(symmetry).
convolve2DSeparable(inBuf, outBuf, x, y, kernel, kernelSize, kernel, kernelSize);
// free up memory space for kernel
delete [] kernel;
// DEBUG //
cout << "sigma: " << sigma << ", Kernel Size: " << kernelSize << endl;
}
///////////////////////////////////////////////////////////////////////////////
// generate 1D Gaussian kernel
// kernel size should be odd number (3, 5, 7, 9, ...)
///////////////////////////////////////////////////////////////////////////////
void makeGaussianKernel(float sigma, float *kernel, int kernelSize)
{
//const double PI = 3.14159265; // PI
int i, center;
float sum = 0; // used for normalization
double result = 0; // result of gaussian func
// compute kernel elements normal distribution equation(Gaussian)
// do only half(positive area) and mirror to negative side
// because Gaussian is even function, symmetric to Y-axis.
center = kernelSize / 2; // center value of n-array(0 ~ n-1)
if(sigma == 0)
{
for(i = 0; i <= center; ++i)
kernel[center+i] = kernel[center-i] = 0;
kernel[center] = 1.0;
}
else
{
for(i = 0; i <= center; ++i)
{
//result = exp(-(i*i)/(double)(2*sigma*sigma)) / (sqrt(2*PI)*sigma);
// dividing (sqrt(2*PI)*sigma) is not needed because normalizing result later
result = exp(-(i*i)/(double)(2*sigma*sigma));
kernel[center+i] = kernel[center-i] = (float)result;
sum += (float)result;
if(i != 0) sum += (float)result;
}
// normalize kernel
// make sum of all elements in kernel to 1
for(i = 0; i <= center; ++i)
kernel[center+i] = kernel[center-i] /= sum;
}
// DEBUG //
#if 0
cout << "1-D Gaussian Kernel\n";
cout << "===================\n";
sum = 0;
for(i = 0; i < kernelSize; ++i)
{
cout << i << ": " << kernel[i] << endl;
sum += kernel[i];
}
cout << "Kernel Sum: " << sum << endl;
#endif
}
///////////////////////////////////////////////////////////////////////////////
// Separable 2D Convolution (unsigned char version)
// If the MxN kernel can be separable to (Mx1) and (1xN) matrices, the
// multiplication can be reduced to M+N comapred to MxN in normal convolution.
// It does not check the output is excceded max for performance reason. And we
// assume the kernel contains good(valid) data, therefore, the result cannot be
// larger than max.
///////////////////////////////////////////////////////////////////////////////
bool convolve2DSeparable(unsigned char* in, unsigned char* out, int dataSizeX, int dataSizeY,
float* kernelX, int kSizeX, float* kernelY, int kSizeY)
{
int i, j, k, m, n;
float *tmp, *sum; // intermediate data buffer
unsigned char *inPtr, *outPtr; // working pointers
float *tmpPtr, *tmpPtr2; // working pointers
int kCenter, kOffset, endIndex; // kernel indice
// check validity of params
if(!in || !out || !kernelX || !kernelY) return false;
if(dataSizeX <= 0 || kSizeX <= 0) return false;
// allocate temp storage to keep intermediate result
tmp = new float[dataSizeX * dataSizeY];
if(!tmp) return false; // memory allocation error
// store accumulated sum
sum = new float[dataSizeX];
if(!sum) return false; // memory allocation error
// covolve horizontal direction ///////////////////////
// find center position of kernel (half of kernel size)
kCenter = kSizeX >> 1; // center index of kernel array
endIndex = dataSizeX - kCenter; // index for full kernel convolution
// init working pointers
inPtr = in;
tmpPtr = tmp; // store intermediate results from 1D horizontal convolution
// start horizontal convolution (x-direction)
for(i=0; i < dataSizeY; ++i) // number of rows
{
kOffset = 0; // starting index of partial kernel varies for each sample
// COLUMN FROM index=0 TO index=kCenter-1
for(j=0; j < kCenter; ++j)
{
*tmpPtr = 0; // init to 0 before accumulation
for(k = kCenter + kOffset, m = 0; k >= 0; --k, ++m) // convolve with partial of kernel
{
*tmpPtr += *(inPtr + m) * kernelX[k];
}
++tmpPtr; // next output
++kOffset; // increase starting index of kernel
}
// COLUMN FROM index=kCenter TO index=(dataSizeX-kCenter-1)
for(j = kCenter; j < endIndex; ++j)
{
*tmpPtr = 0; // init to 0 before accumulate
for(k = kSizeX-1, m = 0; k >= 0; --k, ++m) // full kernel
{
*tmpPtr += *(inPtr + m) * kernelX[k];
}
++inPtr; // next input
++tmpPtr; // next output
}
kOffset = 1; // ending index of partial kernel varies for each sample
// COLUMN FROM index=(dataSizeX-kCenter) TO index=(dataSizeX-1)
for(j = endIndex; j < dataSizeX; ++j)
{
*tmpPtr = 0; // init to 0 before accumulation
for(k = kSizeX-1, m=0; k >= kOffset; --k, ++m) // convolve with partial of kernel
{
*tmpPtr += *(inPtr + m) * kernelX[k];
}
++inPtr; // next input
++tmpPtr; // next output
++kOffset; // increase ending index of partial kernel
}
inPtr += kCenter; // next row
}
// END OF HORIZONTAL CONVOLUTION //////////////////////
// start vertical direction ///////////////////////////
// find center position of kernel (half of kernel size)
kCenter = kSizeY >> 1; // center index of vertical kernel
endIndex = dataSizeY - kCenter; // index where full kernel convolution should stop
// set working pointers
tmpPtr = tmpPtr2 = tmp;
outPtr = out;
// clear out array before accumulation
for(i = 0; i < dataSizeX; ++i)
sum[i] = 0;
// start to convolve vertical direction (y-direction)
// ROW FROM index=0 TO index=(kCenter-1)
kOffset = 0; // starting index of partial kernel varies for each sample
for(i=0; i < kCenter; ++i)
{
for(k = kCenter + kOffset; k >= 0; --k) // convolve with partial kernel
{
for(j=0; j < dataSizeX; ++j)
{
sum[j] += *tmpPtr * kernelY[k];
++tmpPtr;
}
}
for(n = 0; n < dataSizeX; ++n) // convert and copy from sum to out
{
// covert negative to positive
*outPtr = (unsigned char)((float)fabs(sum[n]) + 0.5f);
sum[n] = 0; // reset to zero for next summing
++outPtr; // next element of output
}
tmpPtr = tmpPtr2; // reset input pointer
++kOffset; // increase starting index of kernel
}
// ROW FROM index=kCenter TO index=(dataSizeY-kCenter-1)
for(i = kCenter; i < endIndex; ++i)
{
for(k = kSizeY -1; k >= 0; --k) // convolve with full kernel
{
for(j = 0; j < dataSizeX; ++j)
{
sum[j] += *tmpPtr * kernelY[k];
++tmpPtr;
}
}
for(n = 0; n < dataSizeX; ++n) // convert and copy from sum to out
{
// covert negative to positive
*outPtr = (unsigned char)((float)fabs(sum[n]) + 0.5f);
sum[n] = 0; // reset for next summing
++outPtr; // next output
}
// move to next row
tmpPtr2 += dataSizeX;
tmpPtr = tmpPtr2;
}
// ROW FROM index=(dataSizeY-kCenter) TO index=(dataSizeY-1)
kOffset = 1; // ending index of partial kernel varies for each sample
for(i=endIndex; i < dataSizeY; ++i)
{
for(k = kSizeY-1; k >= kOffset; --k) // convolve with partial kernel
{
for(j=0; j < dataSizeX; ++j)
{
sum[j] += *tmpPtr * kernelY[k];
++tmpPtr;
}
}
for(n = 0; n < dataSizeX; ++n) // convert and copy from sum to out
{
// covert negative to positive
*outPtr = (unsigned char)((float)fabs(sum[n]) + 0.5f);
sum[n] = 0; // reset for next summing
++outPtr; // next output
}
// move to next row
tmpPtr2 += dataSizeX;
tmpPtr = tmpPtr2; // next input
++kOffset; // increase ending index of kernel
}
// END OF VERTICAL CONVOLUTION ////////////////////////
// deallocate temp buffers
delete [] tmp;
delete [] sum;
return true;
}
///////////////////////////////////////////////////////////////////////////////
// compute gradient (first order derivative x and y)
///////////////////////////////////////////////////////////////////////////////
void gradient(unsigned char *in, int x, int y, short *deltaX, short *deltaY)
{
int i, j;
int t;
// compute delta X ***************************
// deltaX = f(x+1) - f(x-1)
for(i = 0; i < y; ++i)
{
t = i * x; // current position X per line
// gradient at the first pixel of each line
// note: the edge,pix[t-1] is NOT exsit
deltaX[t] = in[t+1] - in[t];
// gradients where NOT edge
for(j = 1; j < x-1; ++j)
{
t++;
deltaX[t] = in[t+1] - in[t-1];
}
// gradient at the last pixel of each line
t++;
deltaX[t] = in[t] - in[t-1];
}
// compute delta Y ***************************
// deltaY = f(y+1) - f(y-1)
for(j = 0; j < x; ++j)
{
t = j; // current Y position per column
// gradient at the first pixel
deltaY[t] = in[t+x] - in[t];
// gradients for NOT edge pixels
for(i = 1; i < y-1; ++i)
{
t += x;
deltaY[t] = in[t+x] - in[t-x];
}
// gradient at the last pixel of each column
t += x;
deltaY[t] = in[t] - in[t-x];
}
}
///////////////////////////////////////////////////////////////////////////////
// compute magnitude of gradient(deltaX & deltaY) per pixel
///////////////////////////////////////////////////////////////////////////////
void magnitude(short *deltaX, short *deltaY, int x, int y, unsigned short *mag)
{
int i, j, t;
#if 1 // accurate computation (SLOW)
t = 0;
for(i = 0; i < y; ++i)
for(j = 0; j < x; ++j, ++t)
mag[t] = (unsigned short)(sqrt((double)deltaX[t]*deltaX[t] + (double)deltaY[t]*deltaY[t]) + 0.5);
#endif
#if 0 // faster routine (use absolute value)
t = 0;
for(i = 0; i < y; ++i)
for(j = 0; j < x; ++j, ++t)
mag[t] = (unsigned short)(abs(deltaX[t]) + abs(deltaY[t]));
#endif
}
///////////////////////////////////////////////////////////////////////////////
// compute direction of edges for each pixel (angle of 1st derivative of image)
// quantize the normal directions into 16 plus 0(dx=dy=0)
///////////////////////////////////////////////////////////////////////////////
void direction(short *deltaX, short *deltaY, int x, int y, unsigned char *orient)
{
int i, j, t;
#if 1
// compute direction into 16 ways plus 0(dx=dy=0)
// it is symmetry on origin where theta > 180
// 0 : no direction, when dx = dy = 0
// 1 : theta = 0, when dy = 0
// 2 : 0 < theta < 45, when dx > dy > 0
// 3 : theta = 45, when dx = dy > 0
// 4 : 45 < theta < 90, when dy > dx > 0
// 5 : theta = 90, when dx = 0, dy > 0
// 6 : 90 < theta < 135, when dy > -dx > 0
// 7 : theta = 135, when dy = -dx > 0
// 8 : 135 < theta < 180, when -dx > dy > 0
// 9 : theta = 180, when dy = 0, dx > 0
//10 : 180 < theta < 225, when -dx > -dy > 0
//11 : theta = 225, when dx = dy < 0
//12 : 225 < theta < 270, when dy < dx < 0
//13 : theta = 270, when dx = 0, dy < 0
//14 : 270 < theta < 315, when -dy > dx > 0
//15 : theta = 315, when dx = -dy > 0
//16 : 315 < theta < 360, when dx > -dy > 0
t = 0;
for(j = 0; j < y; j++)
{
for(i = 0; i < x; i++)
{
if(deltaX[t] == 0) // all axis directions
{
if(deltaY[t] == 0) orient[t] = 0;
else if(deltaY[t] > 0) orient[t] = 5;
else orient[t] = 13;
}
else if(deltaX[t] > 0)
{
if(deltaY[t] == 0) orient[t] = 1;
else if(deltaY[t] > 0)
{
if(deltaX[t] - deltaY[t] == 0) orient[t] = 3;
else if(deltaX[t] - deltaY[t] > 0) orient[t] = 2;
else orient[t] = 4;
}
else
{
if(deltaX[t] + deltaY[t] == 0) orient[t] = 15;
else if(deltaX[t] + deltaY[t] > 0) orient[t] = 16;
else orient[t] = 14;
}
}
else
{
if(deltaY[t] == 0) orient[t] = 9;
else if(deltaY[t] > 0)
{
if(deltaY[t] + deltaX[t] == 0) orient[t] = 7;
else if(deltaY[t] + deltaX[t] > 0) orient[t] = 6;
else orient[t] = 8;
}
else
{
if(deltaY[t] - deltaX[t] == 0) orient[t] = 11;
else if(deltaY[t] - deltaX[t] > 0) orient[t] = 10;
else orient[t] = 12;
}
}
//if(orient[t] == 16) printf("%d,%d, %d ", deltaX[t], deltaY[t],orient[t]);
t++;
}
}
#endif
#if (0)
t = 0;
for(j = 0; j < y; j++)
{
for(i = 0; i < x; i++)
{
if(deltaX[t] == 0 && deltaY[t] == 0)
orient[t] = 0;
else if(deltaX[t] == 0
else if(deltaX[t] >= 0 && deltaY[t] >= 0) // quadrant I
deltaX[t] - deltaY[t] >= 0 ? orient[t] = 1 : orient[t] = 2;
else if(deltaX[t] < 0 && deltaY[t] >= 0) // quadrant II
deltaY[t] + deltaX[t] >= 0 ? orient[t] = 3 : orient[t] = 4;
else if(deltaX[t] < 0 && deltaY[t] < 0) // quadrant III
deltaY[t] - deltaX[t] >= 0 ? orient[t] = 1 : orient[t] = 2;
else if(deltaX[t] >= 0 && deltaY[t] < 0) // quadrant IV
deltaX[t] + deltaY[t] >= 0 ? orient[t] = 4 : orient[t] = 3;
//printf("%d, %d %d ", deltaX[t], deltaY[t],orient[t]);
t++; // next pixel
}
}
#endif
}
///////////////////////////////////////////////////////////////////////////////
// Non Maximal Suppression
// If the centre pixel is not greater than neighboured pixels in the direction,
// then the center pixel is set to zero.
// This process results in one pixel wide ridges.
///////////////////////////////////////////////////////////////////////////////
void nonMaxSupp(unsigned short *mag, short *deltaX, short *deltaY,
int x, int y, unsigned char *nms)
{
int i, j, t;
float alpha;
float mag1, mag2;
const unsigned char SUPPRESSED = 0;
// put zero all boundaries of image
// TOP edge line of the image
for(j = 0; j < x; ++j)
nms[j] = 0;
// BOTTOM edge line of image
t = (y-1)*x;
for(j = 0; j < x; ++j, ++t)
nms[t] = 0;
// LEFT & RIGHT edge line
t = x;
for(i = 1; i < y; ++i, t+=x)
{
nms[t] = 0;
nms[t+x-1] = 0;
}
t = x + 1; // skip boundaries of image
// start and stop 1 pixel inner pixels from boundaries
for(i = 1; i < y-1; i++, t+=2)
{
for(j = 1; j < x-1; j++, t++)
{
// if magnitude = 0, no edge
if(mag[t] == 0) nms[t] = SUPPRESSED;
else{
// looking for 8 directions (clockwise)
// NOTE: positive dy means down(south) direction on screen(image)
// because dy is defined as f(x,y+1)-f(x,y-1).
// 1: dy/dx <= 1, dx > 0, dy > 0
// 2: dy/dx > 1, dx > 0, dy > 0
// 3: dy/dx >= -1, dx < 0, dy > 0
// 4: dy/dx < -1, dx < 0, dy > 0
// 5: dy/dx <= 1, dx < 0, dy < 0
// 6: dy/dx > 1, dx < 0, dy < 0
// 7: dy/dx <= -1, dx > 0, dy < 0
// 8: dy/dx > -1, dx > 0, dy < 0
// After determine direction,
// compute magnitudes of neighbour pixels using linear interpolation
if(deltaX[t] >= 0)
{
if(deltaY[t] >= 0) // dx >= 0, dy >= 0
{
if((deltaX - deltaY[t]) >= 0) // direction 1 (SEE, South-East-East)
{
alpha = (float)deltaY[t] / deltaX[t];
mag1 = (1-alpha)*mag[t+1] + alpha*mag[t+x+1];
mag2 = (1-alpha)*mag[t-1] + alpha*mag[t-x-1];
}
else // direction 2 (SSE)
{
alpha = (float)deltaX[t] / deltaY[t];
mag1 = (1-alpha)*mag[t+x] + alpha*mag[t+x+1];
mag2 = (1-alpha)*mag[t-x] + alpha*mag[t-x-1];
}
}
else // dx >= 0, dy < 0
{
if((deltaX[t] + deltaY[t]) >= 0) // direction 8 (NEE)
{
alpha = (float)-deltaY[t] / deltaX[t];
mag1 = (1-alpha)*mag[t+1] + alpha*mag[t-x+1];
mag2 = (1-alpha)*mag[t-1] + alpha*mag[t+x-1];
}
else // direction 7 (NNE)
{
alpha = (float)deltaX[t] / -deltaY[t];
mag1 = (1-alpha)*mag[t+x] + alpha*mag[t+x-1];
mag2 = (1-alpha)*mag[t-x] + alpha*mag[t-x+1];
}
}
}
else
{
if(deltaY[t] >= 0) // dx < 0, dy >= 0
{
if((deltaX[t] + deltaY[t]) >= 0) // direction 3 (SSW)
{
alpha = (float)-deltaX[t] / deltaY[t];
mag1 = (1-alpha)*mag[t+x] + alpha*mag[t+x-1];
mag2 = (1-alpha)*mag[t-x] + alpha*mag[t-x+1];
}
else // direction 4 (SWW)
{
alpha = (float)deltaY[t] / -deltaX[t];
mag1 = (1-alpha)*mag[t-1] + alpha*mag[t+x-1];
mag2 = (1-alpha)*mag[t+1] + alpha*mag[t-x+1];
}
}
else // dx < 0, dy < 0
{
if((-deltaX[t] + deltaY[t]) >= 0) // direction 5 (NWW)
{
alpha = (float)deltaY[t] / deltaX[t];
mag1 = (1-alpha)*mag[t-1] + alpha*mag[t-x-1];
mag2 = (1-alpha)*mag[t+1] + alpha*mag[t+x+1];
}
else // direction 6 (NNW)
{
alpha = (float)deltaX[t] / deltaY[t];
mag1 = (1-alpha)*mag[t-x] + alpha*mag[t-x-1];
mag2 = (1-alpha)*mag[t+x] + alpha*mag[t+x+1];
}
}
}
// non-maximal suppression
// compare mag1, mag2 and mag[t]
// if mag[t] is smaller than one of the neighbours then suppress it
if((mag[t] < mag1) || (mag[t] < mag2))
nms[t] = SUPPRESSED;
else
{
if(mag[t] > 255) nms[t] = 255;
else nms[t] = (unsigned char)mag[t];
}
} // END OF ELSE (mag != 0)
} // END OF FOR(j)
} // END OF FOR(i)
}
///////////////////////////////////////////////////////////////////////////////
// hysteresis thresholding
// Remove weak edges and connect splitted edges
// To connect edges, starting at a pixel which is greater than and equal to
// tHigh and search all 8 neighbours. If the neighbour is greater than and
// equal to tLow, then it will also become an edge.
// The range of threshold is 0 < tLow < tHigh < 1.
// The input data for this routine should be the result of NMS, which is thin
// edges of the original images.
///////////////////////////////////////////////////////////////////////////////
void hysteresis(unsigned char *in, int width, int height, float tLow, float tHigh,
unsigned char *out)
{
// check NULL pointer
if(!in || !out) return;
int i, j, t;
int ii, jj, tt;
int iMax, jMax;
unsigned short histo[256] = {0};
int sum, highCount, lowCount;
unsigned char lowValue, highValue;
const unsigned char EDGE = 255;
// clear out buffer(set all to 0)
memset(out, (char)0, size_t(width*height));
// get histogram to compute high & low level trigger values
histogram(in, width, height, histo);
// find the number of pixels where tHigh(percentage) of total pixels except zeros
highCount = (int)((width*height - histo[0]) * tHigh + 0.5f);
// compute high level trigger value using histogram distribution
i = 255; // max value in unsigned char
sum = width * height - histo[0]; // total sum of histogram array except zeros
while(sum > highCount)
sum -= histo[i--];
highValue = (unsigned char)++i;
// find the number of pixels where tLow(percentage) of total pixels except zeros
lowCount = (int)((width*height - histo[0]) * tLow + 0.5f);
// compute low trigger value using histogram distribution
i = 0;
sum = 0;
while(sum < lowCount)
sum += histo[i++];
lowValue = (unsigned char)i;
lowValue = (int)(highValue * tLow + 0.5f);
// search a pixel which is greater than high level trigger value
// then mark it as an edge and start to trace edges using low trigger value
// If the neighbors are greater than low trigger, they become edges
// NOTE: Use general thresholding method for faster computation, but
// hysteresis thresholding eliminates noisy pixels and connects
// the diconnected edges.
t = width + 1; // skip the first line
iMax = height - 1; // max Y value to traverse, ignore bottom border
jMax = width - 1; // max X value to traverse, ignore right border
for(i = 1; i < iMax; ++i)
{
for(j = 1; j < jMax; ++j, ++t)
{
if(in[t] >= highValue && out[t] != EDGE) // trace only if not checked before
{
out[t] = EDGE; // set itself as an edge
traceEdge(in, t, width, lowValue, out); // start tracing from this point
}
}
}
// DEBUG //
#if 1
cout << "HighTreshold: " << tHigh << ", LowTreshold: " << tLow << endl;
cout << "HighTrigger: " << (int)highValue << ", LowTrigger: " << (int)lowValue << endl;
#endif
}
///////////////////////////////////////////////////////////////////////////////
// search 8 neighbour pixels and trace edges if they are greater than threshold
// "in" is NMS data as input.
// "t" is the current position index of the input array.
// "width" is the width of the image, used to compute the index of next line.
// "threshold" is the low trigger value to stop tracing edges.
// "out" is the array to store edges.
///////////////////////////////////////////////////////////////////////////////
void traceEdge(unsigned char *in, int t, int width, int threshold,
unsigned char *out)
{
// check NULL pointers
if(!in || !out) return;
const unsigned char EDGE = 255;
list<int> edges; // list of edges to be checked
int nw, no, ne, we, ea, sw, so, se; // indice of 8 neighbours
// get indice of 8 neighbours
nw = t - width -1; // north-west
no = nw + 1; // north
ne = no + 1; // north-east
we = t - 1; // west
ea = t + 1; // east
sw = t + width - 1; // south-west
so = sw + 1; // south
se = so + 1; // south-east
// initially, test 8 neighbours and add them into edge list only if they are also edges
if(in[nw] >= threshold && out[nw] != EDGE)
{
out[nw] = EDGE;
edges.push_back(nw);
}
if(in[no] >= threshold && out[no] != EDGE)
{
out[no] = EDGE;
edges.push_back(no);
}
if(in[ne] >= threshold && out[ne] != EDGE)
{
out[ne] = EDGE;
edges.push_back(ne);
}
if(in[we] >= threshold && out[we] != EDGE)
{
out[we] = EDGE;
edges.push_back(we);
}
if(in[ea] >= threshold && out[ea] != EDGE)
{
out[ea] = EDGE;
edges.push_back(ea);
}
if(in[sw] >= threshold && out[sw] != EDGE)
{
out[sw] = EDGE;
edges.push_back(sw);
}
if(in[so] >= threshold && out[so] != EDGE)
{
out[so] = EDGE;
edges.push_back(so);
}
if(in[se] >= threshold && out[se] != EDGE)
{
out[se] = EDGE;
edges.push_back(se);
}
// loop until all edge candiates are tested
while(!edges.empty())
{
t = edges.back();
edges.pop_back(); // remove the last after read
nw = t - width -1; // north-west
no = nw + 1; // north
ne = no + 1; // north-east
we = t - 1; // west
ea = t + 1; // east
sw = t + width - 1; // south-west
so = sw + 1; // south
se = so + 1; // south-east
// test 8 neighbours and add it to list so we can trace from it at next loop
if(in[nw] >= threshold && out[nw] != EDGE) // north-west
{
out[nw] = EDGE;
edges.push_back(nw);
}
if(in[no] >= threshold && out[no] != EDGE) // north
{
out[no] = EDGE;
edges.push_back(no);
}
if(in[ne] >= threshold && out[ne] != EDGE) // north-east
{
out[ne] = EDGE;
edges.push_back(ne);
}
if(in[we] >= threshold && out[we] != EDGE) // west
{
out[we] = EDGE;
edges.push_back(we);
}
if(in[ea] >= threshold && out[ea] != EDGE) // east
{
out[ea] = EDGE;
edges.push_back(ea);
}
if(in[sw] >= threshold && out[sw] != EDGE) // south-west
{
out[sw] = EDGE;
edges.push_back(sw);
}
if(in[so] >= threshold && out[so] != EDGE) // south
{
out[so] = EDGE;
edges.push_back(so);
}
if(in[se] >= threshold && out[se] != EDGE) // south-east
{
out[se] = EDGE;
edges.push_back(se);
}
}
}
'Computer Vision' 카테고리의 다른 글
| 머신 비전(Machine Vision)이란 (0) | 2014.10.27 |
|---|---|
| 그림자 영역 히스토그램(histogram) 테스트 (0) | 2014.09.02 |
| 가우시안 혼합 모델(Gaussian Mixture Models) GMM (0) | 2014.07.23 |
| 라인 허프 변환(hough transform for line detection) (7) | 2014.07.14 |
| 소벨 윤곽선 검출(Sobel Edge Detector) 소벨 검출기 소벨 에지 디텍터 (7) | 2014.07.10 |