Canny边缘检测

float InvSqrt (float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i >> 1);        // 计算第一个近似根
x = *(float*)&i;
x = x*(1.5f - xhalf*x*x);       // 牛顿迭代法
return 1/x;
}




void cannyFollowEdge(int y, int x, int nThredLow, unsigned char * cannyImageData, int * magnitudeImage, int width, int height)  
{  
    //在强边缘周围寻找弱边缘
int k = 0, index = 0;
    int xNum[8] = {1, 1, 0, -1, -1, -1, 0, 1};  
    int yNum[8] = {0, 1, 1, 1, 0, -1 ,-1, -1};  
    int yy = 0, xx = 0;  
    for(k = 0; k < 8; k++)  
    {  
        yy = y + yNum[k];  
        xx = x + xNum[k];
index = yy * width + xx;
        if(cannyImageData[index] == 128 && magnitudeImage[index] >= nThredLow )  
        {  
            //该点设为边界点  
            cannyImageData[index] = 255;  
            //以该点为中心再搜索  
            cannyFollowEdge(yy,xx,nThredLow,cannyImageData,magnitudeImage,width,height);  
        }  
    }  
}  


void NonMaxValueControl(int * MagnitudeGradt, unsigned char * cannyImageData, \
short *xDerivative, short *yDerivative, double *Theta,  int width, int height)
{  
int g1 = 0, g2 = 0, g3 = 0, g4 = 0;//用于进行插值，得到亚像素点坐标值  
double dTmp1 = 0.0, dTmp2 = 0.0;//保存两个亚像素点插值得到的灰度数据  
double dWeight=0.0; 
int i = 0, j = 0;


//边缘初始化
for(i = 0; i < width; i++)  
{  
cannyImageData[i] = 0;  
cannyImageData[(height - 1) * width + i] = 0;  
}  
for(j = 0; j < height; j++)  
{  
cannyImageData[j * width] = 0;  
cannyImageData[j * width + (width - 1)] = 0;  
}


//非极大值抑制
for( i = 1; i < (width - 1); i++)  
{  
for(j = 1; j < ( height - 1); j++)  
{  
int nPointIdx = i + j * width;


//幅度值为0的点直接赋值0
if(MagnitudeGradt[nPointIdx] == 0)
cannyImageData[nPointIdx] = 0;  
else  
{  
/* 1,2,8,9*/  
if( ((Theta[nPointIdx] >= 90) && (Theta[nPointIdx] < 135)) ||\
((Theta[nPointIdx] >= 270) && (Theta[nPointIdx] < 315)))
{
//根据斜率和四个中间值进行插值求解
g1 = MagnitudeGradt[nPointIdx - width - 1];
g2 = MagnitudeGradt[nPointIdx - width];
g3 = MagnitudeGradt[nPointIdx + width];
g4 = MagnitudeGradt[nPointIdx + width + 1];


dWeight = abs(xDerivative[nPointIdx]) * 1.0/ abs(yDerivative[nPointIdx]);//反正切  
dTmp1 = g1 * dWeight + g2 * ( 1 - dWeight);
dTmp2 = g4 * dWeight + g3 * ( 1 - dWeight);
}  


/* 1,4,6,9*/
else if( ((Theta[nPointIdx] >= 135)&&(Theta[nPointIdx] < 180)) ||\
((Theta[nPointIdx]>=315)&&(Theta[nPointIdx]<360)))  
{  
g1 = MagnitudeGradt[nPointIdx - width - 1];  
g2 = MagnitudeGradt[nPointIdx - 1];  
g3 = MagnitudeGradt[nPointIdx + 1];  
g4 = MagnitudeGradt[nPointIdx + width + 1];  
dWeight = abs(yDerivative[nPointIdx]) * 1.0 / abs(xDerivative[nPointIdx]);   //正切  
dTmp1 = g2 * dWeight + g1 * (1 - dWeight);  
dTmp2 = g4 * dWeight + g3 * (1 - dWeight);
}  
/* 2,3,7,8 */
else if( ((Theta[nPointIdx]>=45)&&(Theta[nPointIdx]<90)) ||\
((Theta[nPointIdx]>=225)&&(Theta[nPointIdx]<270)))  
{  
g1 = MagnitudeGradt[nPointIdx-width];  
g2 = MagnitudeGradt[nPointIdx-width+1];  
g3 = MagnitudeGradt[nPointIdx+width];  
g4 = MagnitudeGradt[nPointIdx+width-1];  
dWeight = abs(xDerivative[nPointIdx]) * 1.0 / abs(yDerivative[nPointIdx]);   //反正切  
dTmp1 = g2*dWeight+g1*(1-dWeight);  
dTmp2 = g3*dWeight+g4*(1-dWeight);
}  
/* 3,4,6,7 */
else if( ((Theta[nPointIdx]>=0)&&(Theta[nPointIdx]<45)) ||\
((Theta[nPointIdx]>=180)&&(Theta[nPointIdx]<225)))  
{  
g1 = MagnitudeGradt[nPointIdx - width + 1];  
g2 = MagnitudeGradt[nPointIdx + 1];  
g3 = MagnitudeGradt[nPointIdx + width - 1];  
g4 = MagnitudeGradt[nPointIdx - 1];  
dWeight = abs(yDerivative[nPointIdx]) * 1.0 / abs(xDerivative[nPointIdx]);   //正切  
dTmp1 = g1*dWeight+g2*(1-dWeight);  
dTmp2 = g3*dWeight+g4*(1-dWeight);  
}  
}         
/* 局部最大值判断，写入检测结果 */  
if((MagnitudeGradt[nPointIdx]>=dTmp1) && (MagnitudeGradt[nPointIdx]>=dTmp2))  
cannyImageData[nPointIdx] = 128;
else
cannyImageData[nPointIdx] = 0;
}  
}  
}


/*===========================================================================*\
函数名：CannyEdge
函数功能：Canny Edge Detection
参数：[in]imageData:输入单通道图像数据
 [out]CannyImageData
 [in]width:输入图像帧宽度
 [in]height:输入图像帧高度
 [in]channel:输入图像通道数
返回值：
作者：ivan.zou
日期：2016.04.27
\*===========================================================================*/
int CannyEdge(unsigned char * grayImageData, unsigned char* cannyImageData, int width, int height, int lowThred, int highThred)
{
int i = 0, j = 0, index = 0;
int imageSize = height * width;
unsigned char gray_var[9] = {0};
short * xDerivative = (short * )malloc(imageSize * sizeof(short));//x向偏导数  
short * yDerivative = (short * )malloc(imageSize * sizeof(short));//y向偏导数  
int * MagititudeGradt = (int *)malloc(imageSize * sizeof(int));//梯度幅值  
double * Theta = (double * )malloc(imageSize * sizeof(double));//梯度方向
if(xDerivative == NULL || yDerivative == NULL || MagititudeGradt == NULL || Theta == NULL)
return -1;


//计算sobel梯度
index = 0;
for(i = 1; i < height - 1; ++i)
{
for(j = 1; j < width - 1; ++j)
{
index = i * width + j;
gray_var[0] = grayImageData[(i - 1) * width + (j - 1)];
gray_var[1] = grayImageData[i * width + j - 1];
gray_var[2] = grayImageData[(i + 1) * width + (j - 1)];
gray_var[3] = grayImageData[(i - 1) * width + j];
gray_var[4] = grayImageData[i * width + j];
gray_var[5] = grayImageData[(i + 1) * width + j];
gray_var[6] = grayImageData[(i - 1) * width + (j + 1)];
gray_var[7] = grayImageData[i * width + j + 1];
gray_var[8] = grayImageData[(i + 1) * width + j + 1];


xDerivative[index] = (gray_var[6] + 2 * gray_var[7] + gray_var[8]) - \
(gray_var[0] + 2 * gray_var[1] + gray_var[2]);
yDerivative[index] = (gray_var[0] + 2 * gray_var[3] + gray_var[6]) - \
(gray_var[2] + 2 * gray_var[5] + gray_var[8]);
}
}


//计算梯度幅值和梯度的方向
index = 0;
for( i = 0; i < height; i++)
{
for( j = 0; j < width; j++)
{
MagititudeGradt[index] = (int)(sqrt(xDerivative[index] * xDerivative[index] + \
yDerivative[index] * yDerivative[index]) + 0.5);
Theta[index] = atan2(yDerivative[index], xDerivative[index]) * 57.3;
if(Theta[index] < 0)
Theta[index] += 360;//角度转换到0~360范围
index++;
}
}


//非极大值抑制
NonMaxValueControl(MagititudeGradt, cannyImageData, xDerivative, yDerivative, Theta, width,height);


//双阈值检测
index = 0;
for(i = 0; i < height; i++)
{
for(j = 0; j < width; j++)
{
if((MagititudeGradt[index] >= highThred) && (128 == cannyImageData[index]))
{
cannyImageData[index] = 255;
cannyFollowEdge(i, j, lowThred, cannyImageData, MagititudeGradt, width, height);
}
index++;
}
}


//将还没有设置为边界的点设置为非边界点  
index = 0;
for(i = 0; i < height; i++)  
{  
for( j = 0; j < width; j++)  
{  
if(cannyImageData[index] != 255)  
{  
cannyImageData[index]  = 0 ;   // 设置为非边界点
}
index++;
}
}


if(xDerivative != NULL)
free(xDerivative);
if(yDerivative != NULL)
free(yDerivative);
if(MagititudeGradt != NULL)
free(MagititudeGradt);
if(Theta != NULL)
free(Theta);
return 0;
}


/*===========================================================================*\
函数名：GaussFilterGray1D
函数功能：高斯一维滤波
参数：[in]imageData:输入单通道图像数据
 [in]width:输入图像帧宽度
 [in]height:输入图像帧高度
 [in]channel:输入图像通道数
返回值：
作者：ivan.zou
日期：2016.04.27
\*===========================================================================*/
int GaussFilterGray1D(unsigned char *imageData, int width, int height, int channel)
{
int i = 0, j = 0;
int m = 0, index = 0;
unsigned int swapVar = 0;
unsigned char *pCurF = NULL;
unsigned char *pCurS = NULL;
unsigned int imageSize = width * height;
unsigned int GaussKernel[3] = {1122,1851,1122};//4096
unsigned char * swapImage = (unsigned char *)malloc(imageSize * sizeof(unsigned char));

if(NULL == swapImage || NULL == imageData)
return -1;
if(channel != 1 || width < 0 || height < 0)
return -1;


//纵向滤波
for(i = 0; i < width; i++)
{
for(j = 1; j < height - 1; j++)
{
index = j * width + i;
swapVar = imageData[index - width] * GaussKernel[0] + \
imageData[index] * GaussKernel[1] + imageData[index + width] * GaussKernel[2];
swapImage[index] = (unsigned char)(swapVar >> 12);  
}  
}


//横向滤波
for(i = 0; i < height; i++)
{  
for(j = 1; j < width - 1; j++)  
{
index = i * width + j;
swapVar = swapImage[index - 1] * GaussKernel[0] + \
swapImage[index] * GaussKernel[1] + swapImage[index + 1] * GaussKernel[2];  
imageData[index] = (unsigned char)(swapVar >> 12);  
}  
}  


//边缘处理
pCurF = imageData;
pCurS = imageData + width;
for( i = 0; i < width; ++i,++pCurF,++pCurS)
pCurF[0] = pCurS[0];


pCurF = imageData + (height - 1) * width;
pCurS = imageData + (height - 2) * width;
for( i = 0; i < width; ++i,++pCurF,++pCurS)
pCurF[0] = pCurS[0];


if( swapImage != NULL)
free(swapImage);
return 0;
}