图像边缘提取算法

Canny边缘检测算法

经典的Canny边缘检测算法通常都是从高斯模糊开始，到基于双阈值实现边缘连接结束。但是在实际工程应用中，考虑到输入图像都是彩色图像，最终边缘连接之后的图像要二值化输出显示，所以完整的Canny边缘检测算法实现步骤如下：

1.      彩色图像转换为灰度图像

2.      对图像进行高斯模糊

3.      计算图像梯度，根据梯度计算图像边缘幅值与角度

4.      非最大信号压制处理（边缘细化）

5.      双阈值边缘连接处理

6.      二值化图像输出结果

0、彩色图像转换为灰度图像

根据彩色图像RGB转灰度公式：gray  =  R * 0.299 + G * 0.587 + B * 0.114

1、对图像进行高斯模糊，实现图形的平滑，去除尖锐噪声

定义：高斯滤波是一种线性平滑滤波，适用于消除高斯噪声，广泛应用于图像处理的减噪过程。通俗的讲，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。高斯滤波的具体操作是：用一个模板（或称卷积、掩模）扫描图像中的每一个像素，用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值。

高斯滤波的相关解释：http://blog.csdn.net/hhygcy/article/details/43290

void gaussianFilter2 (int width, int height) 　　　　　　　//高斯滤波器  
{  
    int templates[25] = { 1, 4, 7, 4, 1,              　　　　//高斯滤波窗口
        4, 16, 26, 16, 4,   
        7, 26, 41, 26, 7,  
        4, 16, 26, 16, 4,   
        1, 4, 7, 4, 1 };        
    memcpy(smooth, gray, width*height*sizeof(int) );         //给smooth赋值gray
    for (int j=2;j<height-2;j++)  
    {  
        for (int i=2;i<width-2;i++)                         //i,j用于循环gray中每一个像素（除去图像的边缘）
        {  
            int sum = 0;  
            int index = 0;  
            for ( int m=j-2; m<j+3; m++)  
            {  
                for (int n=i-2; n<i+3; n++)  
                {  
                    sum += gray[ m*width + n] * templates[index++] ;  
                }  
            }  
            sum /= 273;  　　　　　　　　　　　　　　　　　　　　//sum等于使用高斯滤波平滑后的像素值　　　　　　　　　　　　　　　　　　　　　　　　　　
            if (sum > 255)  
                sum = 255;  
            smooth[ j*width+i ] = sum;  
        }  
    }  
}

2、梯度求解函数

void gradient(int width, int height)
{
    for (int row = 0; row < height-1; row++) 
    {
        for (int col = 0; col < width-1; col++) 
        {
            int index = row * width + col;
            // 计算X方向梯度
            //float xg = 0.5*(smooth[(row+1)*width+col]- smooth[row*width+col]+ smooth[(row+1)*width+col+1]- smooth[row*width+col]); //此处算法有点错误修改如下
            float xg = 0.5*(smooth[(row+1)*width+col]- smooth[row*width+col]+ smooth[(row+1)*width+col+1]- smooth[row*width+col+1]);
            float yg =  0.5*(smooth[row*width+col]- smooth[row*width+col+1]+ smooth[(row+1)*width+col]- smooth[(row+1)*width+col+1]);
            // 计算振幅与角度
            data[index] =sqrt(xg*xg+yg*yg);
            if(xg == 0)
            {
                if(yg > 0)
                {
                    magnitudes[index]=90;                      
                }
                if(yg < 0)
                {
                    magnitudes[index]=-90;
                }
            }
            else if(yg == 0)
            {
                magnitudes[index]=0;
            }
            else
            {
                magnitudes[index] = (float)((atan(yg/xg) * 180)/3.1415926);                 
            }
            // make it 0 ~ 180
            magnitudes[index] += 90;
        }
    }
}

3、非最大信号压制（边缘细化）

信号压制本来是数字信号处理中经常用的，这里的非最大信号压制主要目的是实现边缘细化，通过该步处理边缘像素进一步减少。非最大信号压制主要思想是假设3x3的像素区域，中心像素P(x,y) 根据上一步中计算得到边缘角度值angle，可以将角度分为四个离散值0、45、90、135分类依据如下：

其中黄色区域取值范围为0~22.5 与157.5～180

绿色区域取值范围为22.5 ~ 67.5

蓝色区域取值范围为67.5~112.5

红色区域取值范围为112.5~157.5

分别表示上述四个离散角度的取值范围。得到角度之后，比较中心像素角度上相邻

两个像素，如果中心像素小于其中任意一个，则舍弃该边缘像素点，否则保留。一

个简单的例子如下：

void  byxh(int width, int height)                                     // 非最大信号压制
{
    for (int row = 1; row < height-1; row++) {
        for (int col = 1; col < width-1; col++) {
            int index = row * width + col;
            float angle = magnitudes[index];
            float m0 = data[index];
            magnitudes[index] = m0;
            if(angle >=0 && angle < 22.5)                        // angle 0
            {
                float m1 =data[row*width+col-1];
                float m2 = data[row*width+col+1];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
            else if(angle >= 22.5 && angle < 67.5)                        // angle +45
            {
                float m1 = data[(row-1)*width+col+1];
                float m2 = data[(row+1)*width+col-1];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
            else if(angle >= 67.5 && angle < 112.5)                // angle 90
            {
                float m1 = data[(row+1)*width+col];
                float m2 = data[(row-1)*width+col];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
            else if(angle >=112.5 && angle < 157.5)                // angle 135 / -45
            {
                float m1 = data[(row-1)*width+col-1];
                float m2 = data[(row+1)*width+col+11];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
            else if(angle >=157.5) // angle 0
            {
                float m1 = data[(row+1)*width+col];
                float m2 =data[(row-1)*width+col];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
        }
    }
}

4、 模糊阈值和边缘连接

非最大信号压制以后，输出的幅值如果直接显示结果可能会少量的非边缘像素被包含到结果中，所以要通过选取阈值进行取舍，传统的基于一个阈值的方法如果选择的阈值较小起不到过滤非边缘的作用，如果选择的阈值过大容易丢失真正的图像边

缘，Canny提出基于双阈值(Fuzzy threshold)方法很好的实现了边缘选取，在实际应用中双阈值还有边缘连接的作用。双阈值选择与边缘连接方法通过假设两个阈值其中一个为高阈值TH另外一个为低阈值TL则有：

a.      对于任意边缘像素低于TL的则丢弃

b.      对于任意边缘像素高于TH的则保留

c.      对于任意边缘像素值在TL与TH之间的，如果能通过边缘连接到一个像素大于

TH而且边缘所有像素大于最小阈值TL的则保留，否则丢弃。代码实现如下：

void FuzzyThreshold(int width, int height)
{
    float lowThreshold =1.5;
    float highThreshold =3.75;    //通过改变lowThreshold、highThreshold两个值来修正提取效果
    int offset = 0;
    for(int i=0;i<width*height;i++)
        data[i]=0;

    for (int row = 0; row < height; row++) {
        for (int col = 0; col < width; col++) {
            if(magnitudes[offset] >= highThreshold && data[offset] == 0)
            {
                edgeLink(width, height,col, row, offset, lowThreshold);
            }
            offset++;
        }
    }
}


void edgeLink(int width, int height,int x1, int y1, int index, float threshold) {
    int x0 = (x1 == 0) ? x1 : x1 - 1;
    int x2 = (x1 == width - 1) ? x1 : x1 + 1;
    int y0 = y1 == 0 ? y1 : y1 - 1;
    int y2 = y1 == height -1 ? y1 : y1 + 1;
    data[index] = magnitudes[index];
    for (int x = x0; x <= x2; x++) {
        for (int y = y0; y <= y2; y++) {
            int i2 = x + y * width;
            if ((y != y1 || x != x1)
                && data[i2] == 0 
                && magnitudes[i2] >= threshold)
            {
                edgeLink(width, height,x, y, i2,threshold);
                return;
            }
        }
    }
}

5、图像二值化

void arrayToImg(int width, int height,int bitcount)//二值化
{
    DWORD dwLineBytes=GetLineBytes(width,bitcount);
    int k=0;
    BYTE s;
    for(int i=0;i<height;i++)
    {
        for(int j=0;j<width*3;j++)
        {
            s=data[k];
            *(imgData+dwLineBytes*(height-1-i)+j)=data[k];
            j++;
            *(imgData+dwLineBytes*(height-1-i)+j)=data[k];;
            j++;
            *(imgData+dwLineBytes*(height-1-i)+j)=data[k];;  
            k++;
        }
    }
}

// laoma.cpp : 定义控制台应用程序的入口点。
//所有代码

#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <malloc.h>
#include <ctype.h>
#include <process.h>
#include <math.h>
#include "BMP.h"
BITMAPFILEHEADER bmfh;
BITMAPINFOHEADER bmih;
BYTE *imgData;//定义一个字节型的指针变量
int * R;
int * G;
int * B;
int * gray;
float * smooth;
float *data;
float *magnitudes ;
//检查路径是否合法：文件能打开；以bmp为后缀名
int CheckFilePath(char *filepath);
//读入位图的文件头
int ReadFileHeader(char *filepath,BITMAPFILEHEADER *bmfh);
//打印位图的文件头
void PrintFileHeader(BITMAPFILEHEADER *bmfh);
//读入位图的信息头
int ReadInfoHeader(char *filepath,BITMAPINFOHEADER *bmih);
//打印位图的信息头
void PrintInfoHeader(BITMAPINFOHEADER *bmih);
//创建8位位图的调色板
int CreatePalette(RGBQUAD pal[]);
//读入位图的像素数据
int ReadPixelData(char *filepath,BYTE *imgData);
//计算每行像素所占的字节数
LONG GetLineBytes(int imgWidth,int bitCount);
//打印位图的像素数据
void PrintPixelData(BYTE *imgData,int width,int height,int bitCount);
//打印菜单选项
void PrintMenu();
//另存为位图
int SaveAsImage(char *filepath);
//显示位图
void ShowImage(char * filepath);
void byxh(int width, int height);
void FuzzyThreshold(int width, int height);
void gaussianFilter2 (int width, int height);
void gradient(int width, int height);
void arrayToImg(int width, int height,int bitcount);
int main()
{
    char filepath[256]="1.bmp";//开辟了一个256大小的字符型一维数组
    char saveasfilepath[256];
    int i;
    int width;
    int height;
    int bitCount;
    DWORD dwLineBytes;

    i=CheckFilePath(filepath);//调用函数
    if(i==-1)
    {
        return -1;
    }
    CheckFilePath(filepath); 
    ReadFileHeader(filepath,&bmfh);
    ReadInfoHeader(filepath,&bmih);
    height=bmih.biHeight;
    width=bmih.biWidth;
    bitCount=bmih.biBitCount;
    dwLineBytes=GetLineBytes(width,bitCount);
    imgData=(BYTE*)malloc(dwLineBytes*height*sizeof(BYTE));
    R=(int*)malloc(width*height*sizeof(int));
    G=(int*)malloc(width*height*sizeof(int));
    B=(int*)malloc(width*height*sizeof(int));
    gray=(int*)malloc(width*height*sizeof(int));
    data=(float*)malloc(width*height*sizeof(float));;
    magnitudes=(float*)malloc(width*height*sizeof(float)); ;
    smooth=(float *)malloc(sizeof(float)* width*height);
    ReadPixelData(filepath,imgData);
    //PrintFileHeader(&bmfh);
    //PrintInfoHeader(&bmih);
    //ShowImage(filepath);
    PrintPixelData(imgData,width,height,bitCount);
    //printf("Input the path(ex. d://poon.bmp) you want to save:
");
    //scanf("%s",saveasfilepath);
    //i=SaveAsImage(saveasfilepath);
    gaussianFilter2(width,height);
    gradient(width,height);
    byxh(width,height);
    FuzzyThreshold(width,height);
    //FILE *streamgray= fopen("zsggray.dat", "wb");
    for(int i=0;i<width*height;i++)
    {
        int s=0;
        if(data[i]>2.5)
            data[i]=255;
        else data[i]=0;
        // fwrite(&s, sizeof(int), 1, streamgray); 
    }
    //fclose(streamgray);
    arrayToImg(width,height,bitCount);
    SaveAsImage("zsgResult.bmp");
    ShowImage("zsgResult.bmp");
}
int ReadFileHeader(char *filepath,BITMAPFILEHEADER *bmfh)
{
    FILE *fp;
    fp=fopen(filepath,"rb");
    if(!fp)
    {
        printf("Can not open the file:%s
",filepath);
        return -1;
    }
    if(fread(&bmfh->bfType,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not read bfType in the file header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmfh->bfSize,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not read bfSize in the file header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmfh->bfReserved1,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not read bfReserved1 in the file header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmfh->bfReserved2,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not read bfReserved2 in the file header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmfh->bfOffBits,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not read bfOffBits in the file header.
");
        fclose(fp);
        return -1;
    }
    fclose(fp);
    return 0;
}
int ReadInfoHeader(char *filepath,BITMAPINFOHEADER *bmih)
{
    FILE *fp;
    fp=fopen(filepath,"rb");
    if(!fp)
    {
        printf("Can not open the file:%s
",filepath);
        return -1;
    }
    fseek(fp,14,SEEK_SET);
    if(fread(&bmih->biSize,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not read biSize in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biWidth,sizeof(LONG),1,fp)!=1)
    {
        printf("Can not read biWidth in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biHeight,sizeof(LONG),1,fp)!=1)
    {
        printf("Can not read biHeight in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biPlanes,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not read biPlanes in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biBitCount,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not read biBitCount in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biCompression,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not read biCompression in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biSizeImage,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not read biSizeImage in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biXPelsPerMeter,sizeof(LONG),1,fp)!=1)
    {
        printf("Can not read biXPelsPerMeter in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biYPelsPerMeter,sizeof(LONG),1,fp)!=1)
    {
        printf("Can not read biYPelsPerMeter in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biClrUsed,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not read biClrUsed in the info header.
");
        fclose(fp);
        return -1;
    }
    if(fread(&bmih->biClrImportant,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not read biClrImportant in the info header.
");
        fclose(fp);
        return -1;
    }
    fclose(fp);
    return 0;
}
int CreatePalette(RGBQUAD pal[])
{
    int i;
    if(sizeof(pal)/sizeof(RGBQUAD)!=256)
    {
        printf("The size of the palette must be 256.
");
        return -1;
    }
    for(i=0;i<256;i++)
    {
        pal[i].rgbBlue=i;
        pal[i].rgbGreen=i;
        pal[i].rgbRed=i;
        pal[i].rgbReserved=0;
    }
    return 0;
}
int ReadPixelData(char *filepath,BYTE *imgData)
{
    BITMAPINFOHEADER bmih;
    BITMAPFILEHEADER bmfh;
    BYTE *data;
    FILE *fp;
    int n;
    int width;
    int height;
    int bitCount;
    DWORD dwLineBytes;//扫描行的Windows规定一个扫描行所占的字节数必须是4的倍数（即以long为单位），不足的以0填充
    ReadFileHeader(filepath,&bmfh);
    ReadInfoHeader(filepath,&bmih);
    width=bmih.biWidth;
    height=bmih.biHeight;
    bitCount=bmih.biBitCount;
    dwLineBytes=GetLineBytes(width,bitCount);
    if(_msize(imgData)!=(dwLineBytes*height))
    {
        printf("The size you allocate for the pixel data is not right.
");
        printf("Fittable size: %ld bytes.
",(dwLineBytes*height));
        printf("Your size: %ld bytes.
",sizeof(imgData));
        return -1;
    }
    data=(BYTE*)malloc(dwLineBytes*height*sizeof(BYTE));
    if(!data)
    {
        printf("Can not allocate memory for the pixel data.
");
        return -1;
    }
    fp=fopen(filepath,"rb");
    if(!fp)
    {
        printf("Can not open the file: %s
",filepath);
        free(data);
        return -1;
    }
    if(bitCount==8)
    {
        fseek(fp,bmfh.bfOffBits,SEEK_SET);
    }
    else if(bitCount==24)
    {
        fseek(fp,bmfh.bfOffBits,SEEK_SET);
    }
    else
    {
        printf("Only Support: 8 or 24 bits.
");
        free(data);
        fclose(fp);
        return -1;
    }
    if(fread(data,dwLineBytes*height*sizeof(BYTE),1,fp)!=1)
    {
        printf("Can not read the pixel data.
");
        free(data);
        fclose(fp);
        return -1;
    }
    memcpy(imgData,data,dwLineBytes*height*sizeof(BYTE));
    free(data);
    fclose(fp);
    return 0;
}
void PrintFileHeader(BITMAPFILEHEADER *bmfh)
{
    printf("The contents in the file header of the BMP file:
");
    printf("bfOffBits: %ld
",bmfh->bfOffBits);
    printf("bfReserved1: %ld
",bmfh->bfReserved1);
    printf("bfReserved2: %ld
",bmfh->bfReserved2);
    printf("bfSize: %ld
",bmfh->bfSize);
    printf("bfType: %ld
",bmfh->bfType);
}
void PrintInfoHeader(BITMAPINFOHEADER *bmih)
{
    printf("The content in the info header of the BMP file:
");
    printf("biBitCount: %ld
",bmih->biBitCount);
    printf("biClrImportant: %ld
",bmih->biClrImportant);
    printf("biClrUsed: %ld
",bmih->biClrUsed);
    printf("biCompression: %ld
",bmih->biCompression);
    printf("biHeight: %ld
",bmih->biHeight);
    printf("biPlanes: %ld
",bmih->biPlanes);
    printf("biSize: %ld
",bmih->biSize);
    printf("biSizeImage: %ld
",bmih->biSizeImage);
    printf("biWidth: %ld
",bmih->biWidth);
    printf("biXPelsPerMeter: %ld
",bmih->biXPelsPerMeter);
    printf("biYPelsPerMeter: %ld
",bmih->biYPelsPerMeter);
}
LONG GetLineBytes(int imgWidth,int bitCount)
{
    return (imgWidth*bitCount+31)/32*4;
}
void PrintPixelData(BYTE *imgData,int width,int height,int bitCount)
{
    int i;
    int j;
    int p;
    DWORD dwLineBytes=GetLineBytes(width,bitCount);
    if(bitCount==8)
    {
        for(i=0;i<height;i++)
        {
            for(j=0;j<width;j++)
            {
                p=*(imgData+dwLineBytes*(height-1-i)+j);
                printf("%d,",p);
            }
            printf("
");
        }
    }
    else if(bitCount==24)
    {
        int k=0;
        //FILE *streamR= fopen("zsgRR.dat", "wb");
        //FILE *streamG= fopen("zsgGG.dat", "wb");
        //FILE *streamB= fopen("zsgBB.dat", "wb");
        for(i=0;i<height;i++)
        {
            for(j=0;j<width*3;j++)
            {
                B[k]=*(imgData+dwLineBytes*(height-1-i)+j);
                j++;
                G[k] =*(imgData+dwLineBytes*(height-1-i)+j);
                j++;
                R[k]=*(imgData+dwLineBytes*(height-1-i)+j);  
                //fwrite(&R[k], sizeof(R), 1, streamR); 
                //fwrite(&G[k], sizeof(R), 1, streamG); 
                //fwrite(&B[k], sizeof(R), 1, streamB);
                gray[k]= (R[k]*299 + G[k]*587 + B[k]*114 + 500) / 1000;
                k++;
            }
        }
        //fclose(streamR);fclose(streamG);
        //fclose(streamB);

    }
    else
    {
        printf("Only supported: 8 or 24 bits.
");
    }
}
int CheckFilePath(char *filepath)
{
    FILE *fp;
    int len=strlen(filepath)/sizeof(char);
    char ext[3];
    if(filepath[0]!=int('"'))
    {
        strncpy(ext,&filepath[len-3],3);
        if(!(ext[0]=='b' && ext[1]=='m' && ext[2]=='p'))
        {
            printf("Error: The file is not a BMP file.
");
            printf("Error: The extention of the filename must be 'bmp',not 'BMP'
");
            return -1;
        }
        fp=fopen(filepath,"r");
        if(!fp)
        {
            printf("Error: The path is not correct.
");
            return -1;
        }
        fclose(fp);
    }
    else
    {
        printf("Error: The path must not include blank space.
");
        return -1;
    }
    return 0;
}

int SaveAsImage(char *filepath)
{
    FILE *fp;
    RGBQUAD pal[256];
    int i;
    int height;
    DWORD dwLineBytes;

    if(bmih.biBitCount!=8 && bmih.biBitCount!=24)
    {
        printf("Error: only supported 8 or 24 bits.
");
        return -1;
    }
    height=bmih.biHeight;
    dwLineBytes=GetLineBytes(bmih.biWidth,bmih.biBitCount);
    fp=fopen(filepath,"wb");
    if(!fp)
    {
        printf("Error: Can not open the file:%s
",filepath);
        return -1;
    }

    for(i=0;i<256;i++)
    {
        pal[i].rgbReserved=0;
        pal[i].rgbBlue=i;
        pal[i].rgbGreen=i;
        pal[i].rgbRed=i;
    }

    if(fwrite(&bmfh.bfType,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not write bfType in the file header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmfh.bfSize,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not write bfSize in the file header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmfh.bfReserved1,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not write bfReserved1 in the file header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmfh.bfReserved2,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not write bfReserved2 in the file header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmfh.bfOffBits,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not write bfOffBits in the file header.
");
        fclose(fp);
        return -1;
    }
    if(fwrite(&bmih.biSize,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not write biSize in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biWidth,sizeof(LONG),1,fp)!=1)
    {
        printf("Can not write biWidth in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biHeight,sizeof(LONG),1,fp)!=1)
    {
        printf("Can not write biHeight in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biPlanes,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not write biPlanes in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biBitCount,sizeof(WORD),1,fp)!=1)
    {
        printf("Can not write biBitCount in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biCompression,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not write biCompression in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biSizeImage,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not write biSizeImage in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biXPelsPerMeter,sizeof(LONG),1,fp)!=1)
    {
        printf("Can not write biXPelsPerMeter in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biYPelsPerMeter,sizeof(LONG),1,fp)!=1)
    {
        printf("Can not write biYPelsPerMeter in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biClrUsed,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not write biClrUsed in the info header.
");
        fclose(fp);
        return -1;
    }

    if(fwrite(&bmih.biClrImportant,sizeof(DWORD),1,fp)!=1)
    {
        printf("Can not write biClrImportant in the info header.
");
        fclose(fp);
        return -1;
    }

    if(bmih.biBitCount==8)
    {
        if(fwrite(pal,sizeof(RGBQUAD),256,fp)!=256)
        {
            printf("Error: can not write the color palette.
");
            fclose(fp);
            return -1;
        }
    }

    if(fwrite(imgData,height*dwLineBytes,1,fp)!=1)
    {
        printf("Error: can not write the pixel data.
");
        fclose(fp);
        return -1;
    }

    fclose(fp);
    printf("Save As the image successfully.
");
    return 0;
}
void ShowImage(char * filepath)
{
    char cmd[266];
    strcpy(cmd,"start ");
    strcat(cmd,filepath);
    printf("%s
",cmd);
    system(cmd);
}

void gaussianFilter2 (int width, int height)  
{  
    int templates[25] = { 1, 4, 7, 4, 1,   
        4, 16, 26, 16, 4,   
        7, 26, 41, 26, 7,  
        4, 16, 26, 16, 4,   
        1, 4, 7, 4, 1 };        
    memcpy(smooth, gray, width*height*sizeof(int) );  
    for (int j=2;j<height-2;j++)  
    {  
        for (int i=2;i<width-2;i++)  
        {  
            int sum = 0;  
            int index = 0;  
            for ( int m=j-2; m<j+3; m++)  
            {  
                for (int n=i-2; n<i+3; n++)  
                {  
                    sum += gray[ m*width + n] * templates[index++] ;  
                }  
            }  
            sum /= 273;  
            if (sum > 255)  
                sum = 255;  
            smooth[ j*width+i ] = sum;  
        }  
    }  
}

void gradient(int width, int height)
{
    for (int row = 0; row < height-1; row++) 
    {
        for (int col = 0; col < width-1; col++) 
        {
            int index = row * width + col;

            // 计算X方向梯度
            //float xg = 0.5*(smooth[(row+1)*width+col]- smooth[row*width+col]+ smooth[(row+1)*width+col+1]- smooth[row*width+col]); //此处算法有点错误修改如下
            float xg = 0.5*(smooth[(row+1)*width+col]- smooth[row*width+col]+ smooth[(row+1)*width+col+1]- smooth[row*width+col+1]);
            float yg =  0.5*(smooth[row*width+col]- smooth[row*width+col+1]+ smooth[(row+1)*width+col]- smooth[(row+1)*width+col+1]);
            // 计算振幅与角度
            data[index] =sqrt(xg*xg+yg*yg);
            if(xg == 0)
            {
                if(yg > 0)
                {
                    magnitudes[index]=90;                      
                }
                if(yg < 0)
                {
                    magnitudes[index]=-90;
                }
            }
            else if(yg == 0)
            {
                magnitudes[index]=0;
            }
            else
            {
                magnitudes[index] = (float)((atan(yg/xg) * 180)/3.1415926);                 
            }
            // make it 0 ~ 180
            magnitudes[index] += 90;
        }
    }
}

void  byxh(int width, int height)// 非最大信号压制
{
    for (int row = 1; row < height-1; row++) {
        for (int col = 1; col < width-1; col++) {
            int index = row * width + col;
            float angle = magnitudes[index];
            float m0 = data[index];
            magnitudes[index] = m0;
            if(angle >=0 && angle < 22.5) // angle 0
            {
                float m1 =data[row*width+col-1];
                float m2 = data[row*width+col+1];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
            else if(angle >= 22.5 && angle < 67.5) // angle +45
            {
                float m1 = data[(row-1)*width+col+1];
                float m2 = data[(row+1)*width+col-1];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
            else if(angle >= 67.5 && angle < 112.5) // angle 90
            {
                float m1 = data[(row+1)*width+col];
                float m2 = data[(row-1)*width+col];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
            else if(angle >=112.5 && angle < 157.5) // angle 135 / -45
            {
                float m1 = data[(row-1)*width+col-1];
                float m2 = data[(row+1)*width+col+11];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
            else if(angle >=157.5) // angle 0
            {
                float m1 = data[(row+1)*width+col];
                float m2 =data[(row-1)*width+col];
                if(m0 < m1 || m0 < m2)
                {
                    magnitudes[index] = 0;
                }
            }
        }
    }

}
void edgeLink(int width, int height,int x1, int y1, int index, float threshold) {
    int x0 = (x1 == 0) ? x1 : (x1 - 1);
    int x2 = (x1 == width - 1) ? x1 : (x1 + 1);
    int y0 = (y1 == 0) ? y1 : (y1 - 1);
    int y2 = (y1 == height -1) ? y1 : (y1 + 1);
    data[index] = magnitudes[index];
    for (int x = x0; x <= x2; x++) {
        for (int y = y0; y <= y2; y++) {
            int i2 = x + y * width;
            if ((y != y1 || x != x1)
                && data[i2] == 0 
                && magnitudes[i2] >= threshold)
            {
                edgeLink(width, height,x, y, i2,threshold);
                return;
            }
        }
    }
}
void FuzzyThreshold(int width, int height)
{
    float lowThreshold =2.5;
    float highThreshold =7.5;
    int offset = 0;
    for(int i=0;i<width*height;i++)
        data[i]=0;

    for (int row = 0; row < height; row++) {
        for (int col = 0; col < width; col++) {
            if(magnitudes[offset] >= highThreshold && data[offset] == 0)
            {
                edgeLink(width, height,col, row, offset, lowThreshold);
            }
            offset++;
        }
    }
}
void arrayToImg(int width, int height,int bitcount)//二值化
{
         DWORD dwLineBytes=GetLineBytes(width,bitcount);
         int k=0;
         BYTE s;
         for(int i=0;i<height;i++)
         {
             for(int j=0;j<width*3;j++)
             {
                s=data[k];
                 *(imgData+dwLineBytes*(height-1-i)+j)=data[k]> 0 ? -1 : 0xff000000;
                  j++;
                  *(imgData+dwLineBytes*(height-1-i)+j)=data[k]> 0 ? -1 : 0xff000000;
                 j++;
                  *(imgData+dwLineBytes*(height-1-i)+j)=data[k]> 0 ? -1 : 0xff000000;  
            k++;
             }
         }
//     for(int i=0; i<width*height; i++)  
//     {
//         int gray = data[i]; 
//         imgData[i] = gray > 0 ? -1 : 0xff000000;
//     }
}

 参考：

http://blog.csdn.net/jia20003/article/details/41173767

http://blog.csdn.net/jia20003/article/details/7234741

http://blog.csdn.net/jia20003/article/details/7664777