Hough变换检测直线

直线的霍夫变换：

霍夫空间极坐标与图像空间的转换公式：
p = y * sin(theta) + x * cos(theta);

之后遍历图像的每个坐标点，每个坐标点以一度为增量，求取对应的p值，存入数组中，查找数组中数目大于一定阈值的p和theta，再在图像空间中把 直线 恢复出来

霍夫变换就是把图像左边空间上的线段转换到霍夫空间一个点，然后通过点的数目多少来确定是否为一条线段（但是画出的结果为一条直线）

实现代码如下：

#include <iostream>
#include "gdal_priv.h"
#include <string>
using namespace std;


//图像的膨胀
//算法的实现依然有重复的地方，会造成空间和时间上的浪费，但是效果还可以，暂定如此，如有好的算法再进行改进
void Expand(unsigned char* date,unsigned char* ExpandImage,int Width,int Height)
{
    int x, y;
    int Direction[8][2] = {-1, -1, 0, -1, 1, -1, 1, 0, 1, 1, 0, 1, -1, 1, -1, 0};
    int Index;


    for(int i = 1;i < Height - 1;++i)
    {
        for(int j = 1;j < Width - 1;++j)
        {
            if(date[i * Width + j] == 255)
            {
                for(int k = 0;k < 8;++k)
                {
                    x = i + Direction[k][0];
                    y = j + Direction[k][1];
                    Index = x * Width + y;

                    ExpandImage[Index] = 255;

                }
            }
        }
    }


    return;
}

//图像的腐蚀
void Erosion(unsigned char* date,unsigned char* ErosionImage,int Width,int Height)
{
    int x, y;
    int Direction[8][2] = {-1, -1, 0, -1, 1, -1, 1, 0, 1, 1, 0, 1, -1, 1, -1, 0};
    int Index;


    for(int i = 1;i < Height - 1;++i)
    {
        for(int j = 1;j < Width - 1;++j)
        {
            if(date[i * Width + j] == 0)
            {
                for(int k = 0;k < 8;++k)
                {
                    x = i + Direction[k][0];
                    y = j + Direction[k][1];
                    Index = x * Width + y;

                    ErosionImage[Index] = 0;

                }
            }
        }
    }



    return;
}

//图像相减求边界
void Subtraction(unsigned char* ResultImage,unsigned char* leftImage,unsigned char* rightImage,int Width,int Height)
{
    int Index;
    for(int i = 0;i < Height;++i)
        for(int j = 0;j < Width;++j)
        {
            Index = i * Width + j;
            ResultImage[Index] = leftImage[Index] - rightImage[Index];
        }

    return ;
}

void FindBoundary(unsigned char* image,unsigned char* tempImage,int Width,int Height)
{
    //生成对于图像膨胀的图像
    unsigned char* ExpandImage = new unsigned char[Width * Height];
    memset(ExpandImage,0,sizeof(unsigned char) * Width * Height);

    Expand(image ,ExpandImage,Width ,Height);


    Subtraction(tempImage,ExpandImage,image,Width,Height);
}


/*************************************************************************
* 直线的Hough检测
* 参数：image0为原图形，image1为边缘检测结果，w、h为图像的宽和高
* 由于得到的Hough变换结果图像与原图像大小不同，为了得到新的宽高信息
* w、h使用引用类型
*************************************************************************/
unsigned char** HoughLine(unsigned char* image0, unsigned char* &tempImage, int &Width, int &Height,int scale=1)
{
    //定义三角函数表
    double sinValue[360];
    double cosValue[360];
    int i,x,y;
    int k = 100;
    int p = (int)(sqrt((double)(Width * Width + Height * Height) + 1));   //计算对角线长度
    //申请临时存储空间 用来保存边缘检测结果
//  tempImage = new unsigned char[Width * Height];
    memset(tempImage,0,sizeof(unsigned char) * Width * Height);
    //边缘检测
//  SideGrandiant(image0, tempImage, Width, Height);
    FindBoundary(image0, tempImage, Width, Height);





//  //根据Hough变换结果图的大小 重新为输出图象分配空间
//  if(image1 != NULL) 
//      delete image1;
//
    image1 = (BYTE*)malloc(sizeof(BYTE)*p*360*4);
//
//  image1 = new unsigned char[p * 360];
//  memset(image1,0,p * 360);

    //将图像转换为矩阵形式
//  BYTE** HoughBuf =CreatImage(image1,360,p);

    unsigned char** HoughBuf = new unsigned char* [p];
    for(int i = 0;i < p;++i)
    {
        HoughBuf[i] = new unsigned char[360];
        memset(HoughBuf[i],0,sizeof(unsigned char) * 360);
        //HoughBuf[i] = image1 + i * 360;
    }

    //for(int i = 0;i < p;++i)
    //{
    //  for(int j = 0;j < 360;++j)
    //      HoughBuf[i][j] = image1[i * 360 + j];
    //}
    //计算三角函数表
    for(i=0; i<360 ; i++)
    {
        sinValue[i] = sin(i*3.1415926/180);
        cosValue[i] = cos(i*3.1415926/180);
    }
    int tp;
    //遍历原图象中的每个像素
    for(y = 0; y < Height; y++)
        for(x = 0; x < Width; x++)
        {
            //对经过当前像素的任何直线区域进行检测
            for(i = 0; i < 360; i++)
            {
                if( tempImage[(y * Width + x)] > k )
                {
                    tp = (int)( x * sinValue[i] + y * cosValue[i]);
                    //忽略负数同时防止计数器溢出
                    if (tp < 0 || HoughBuf[tp][i] == 255) continue;
                    HoughBuf[tp][i] += scale;
                }
            }
        }
    //重新设定图象大小
    //Width = 360;
    //Height = p;

//  delete tempImage;

    return HoughBuf;
}


//画检测到的直线
void DrawLine(unsigned char* LineIamge,int Width,int Height,int p,int theta)
{
    double k,b;
    int x,y;
    if(theta != 90)    //如果斜率存在
    {
        //计算直线方程的参数
        b = p / cos(theta * 3.1415926535 / 180);
        k = -sin(theta * 3.1415926535 / 180) / cos(theta * 3.1415926535 / 180);
        y=0;
        x=0;
        //斜率小于1的情况
        if(abs(k) <= 1)
        {
            for(x = 0;x < Width;x++)
            {
                y=(int)(k * x + b);
                if(y >= 0 && y < Height)
                {
                    LineIamge[y * Width + x] = 255;
                }
            }
        }
        //斜率大于1的情况
        else
        {
            for(y = 0;y < Height;y++)
            {
                x = (int)(y / k - b / k);
                if(x >= 0 && x < Width)
                {
                    /*imageBuf[y][x*4]=255;
                    imageBuf[y][x*4+1]=0;
                    imageBuf[y][x*4+2]=0;
                    imageBuf[y][x*4+3]=255;*/
                    LineIamge[y * Width + x] = 255;
                }
            }
        }
    }
    //斜率不存在的情况
    else
    {
        for(y = 0; y < Height;y++)
        {
            /*imageBuf[y][p*4]=255;
            imageBuf[y][p*4+1]=0;
            imageBuf[y][p*4+2]=0;
            imageBuf[y][p*4+3]=255;*/
            LineIamge[y * Width + p] = 255;
        }
    }
}



int main()
{
    GDALAllRegister();
    CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");

    string str1="1.1.bmp";
    string str2="1.2.tif";
    string str3 = "1.3.tif";

    GDALDataset* pInDataset=(GDALDataset*)GDALOpen(str1.data() ,GA_ReadOnly);

    if(pInDataset == nullptr)
    {
        cout<<"未找到输入图像"<<endl;
        getchar();
        return 1;
    }

    int Width=pInDataset->GetRasterXSize();
    int Height = pInDataset->GetRasterYSize();
    int band = pInDataset->GetRasterCount();

    double dGeoTrans[6]={};
    pInDataset->GetGeoTransform(dGeoTrans);
    const char* cProjectionRef = pInDataset->GetProjectionRef();

    unsigned char* Image = new unsigned char[Width * Height];
    GDALRasterBand* pRasterBand = pInDataset->GetRasterBand(1);

    CPLErr err=pRasterBand->RasterIO(GF_Read ,0 ,0 ,Width ,Height ,Image ,Width ,Height ,GDT_Byte ,0,0);

    unsigned char* ResultImage = new unsigned char[Width * Height];//保存图像的边缘信息，通过膨胀后的图像减去原始图像得到的边缘图像，用来检测是否得到了边缘
    //ResultImage = nullptr;
    //unsigned char** HoughBuf = nullptr;
    unsigned char** HoughBuf;
    HoughBuf = HoughLine(Image,ResultImage,Width,Height);



    GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");

    GDALDataset* pOutDataset = pDriver->Create("1.2.tif" ,Width ,Height ,1 ,GDT_Byte ,NULL);
    GDALRasterBand* pOutRasterband=pOutDataset->GetRasterBand(1);
    pOutRasterband->RasterIO(GF_Write ,0 ,0 ,Width ,Height ,ResultImage ,Width ,Height ,GDT_Byte ,0 ,0);


    int p = (int)(sqrt((double)(Width * Width + Height * Height) + 1));   //计算对角线长度

    unsigned char* LineImage = new unsigned char[Width * Height];
    memset(LineImage,0,sizeof(unsigned char) * Width * Height);

    if(HoughBuf != nullptr)
        for(int i = 0;i < p;++i)
        {
            for(int j = 0;j < 360;++j)
            {
                if(HoughBuf[i][j] > 200)        //设置阈值为200，可自行设置
                {
                    cout<<(int)HoughBuf[i][j]<<" ";
                    DrawLine(LineImage,Width,Height,i,j);//得到的结果为一条直线
                }
            }
            cout<<"***************"<<endl;
        }

    GDALDriver* pLineDriver = GetGDALDriverManager()->GetDriverByName("GTiff");

    GDALDataset* pOutLineDataset = pLineDriver->Create("1.3.tif" ,Width ,Height ,1 ,GDT_Byte ,NULL);
    GDALRasterBand* pOutLineRasterband=pOutLineDataset->GetRasterBand(1);
    pOutLineRasterband->RasterIO(GF_Write ,0 ,0 ,Width ,Height ,LineImage ,Width ,Height ,GDT_Byte ,0 ,0);

    delete Image;
    delete ResultImage;
    delete LineImage;
    for(int i = 0;i < p;++i)
        delete HoughBuf[i];
    delete HoughBuf;

    GDALClose(pOutDataset);
    GDALClose(pInDataset);
    GDALClose(pOutLineDataset);
    GetGDALDriverManager()->DeregisterDriver(pDriver);
    GetGDALDriverManager()->DeregisterDriver(pLineDriver);

    system("pause");
    return 0;
}