欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识！
您现在的位置是： 首页

ENVI%2线性拉伸算法实现

程序员文章站 2024-03-24 22:57:46
...

一、在ENVI里面有Linear和Linear2%的线性拉伸的方法，当然还有其它各种各样的拉伸方式，用的最多的就是Linear2%

二、Linear方法较为简单，原理如下所示：我们需要增强的图像范围一般[0,255]，因此下面公式中的c=0，d=255，得到一般公式为g(x,y)=255/(b-a)*(f(x,y)-a)，然后将小于a的灰度值赋值为0，大于b的灰度值赋值为255即可。

ENVI%2线性拉伸算法实现

三、Linear算法实现（本算法是在Qt下编译运行的，需要调用GDAL库，同时设置了输出文件方便查看效果）

    GDALDataset *dataSet = (GDALDataset*)GDALOpen(pszSrcFile,GA_ReadOnly);//原始图像，路径自己改

    int width = dataSet->GetRasterXSize();
    int height = dataSet->GetRasterYSize();
    GDALDataType dataType = dataSet->GetRasterBand(1)->GetRasterDataType();
    GDALDriver *poDriver = (GDALDriver*)GDALGetDriverByName("GTiff");
    GDALDataset *poDstDS = poDriver->Create(pszDstFile,width,height,3,dataType,NULL);//创建输出路径，路径自己改

    GDALRasterBand *poBand1,*poBand2,*poBand3;
    GDALRasterBand *pDstband1 = poDstDS->GetRasterBand(1);//输出图像R波段
    GDALRasterBand *pDstband2 = poDstDS->GetRasterBand(2);//输出图像G波段
    GDALRasterBand *pDstband3 = poDstDS->GetRasterBand(3);//输出图像B波段

    poBand1 = dataSet->GetRasterBand(1);//读取R波段

    unsigned char *pBuf1 = new unsigned char[width*height];//存储原始图像灰度值
    unsigned char *BufR = new unsigned char[width*height];//存储增强图像灰度值
    poBand1->RasterIO(GF_Read,0,0,width,height,pBuf1,width,height,GDT_Byte,0,0);
    double adfMinMax[2];//获取图像灰度最大最小值
    GDALComputeRasterMinMax((GDALRasterBandH)poBand1,FALSE,adfMinMax);
    int min = adfMinMax[0];
    int max = adfMinMax[1];
    for(int i=0;i<width*height;i++)
    {
        if(pBuf1[i]<min)
            BufR[i] = 0;
        else if(pBuf1[i]>max)
            BufR[i] = 255;
        else
            BufR[i] = 255/(max - min)*(pBuf1[i] - min);
    }
    pDstband1->RasterIO(GF_Write,0,0,width,height,BufR,width,height,GDT_Byte,0,0);//写入R波段

    poBand2 = dataSet->GetRasterBand(2);//读取G波段

    unsigned char *pBuf2 = new unsigned char[width*height];
    unsigned char *BufG = new unsigned char[width*height];
    poBand2->RasterIO(GF_Read,0,0,width,height,pBuf2,width,height,GDT_Byte,0,0);
    double adfMinMax2[2];
    GDALComputeRasterMinMax((GDALRasterBandH)poBand2,FALSE,adfMinMax2);
    int min2 = adfMinMax2[0];
    int max2 = adfMinMax2[1];
    for(int i=0;i<width*height;i++)
    {
        if(pBuf2[i]<min2)
            BufG[i] = 0;
        else if(pBuf2[i]>max2)
            BufG[i] = 255;
        else
            BufG[i] = 255/(max2 - min2)*(pBuf2[i] - min2);
    }
    pDstband2->RasterIO(GF_Write,0,0,width,height,BufG,width,height,GDT_Byte,0,0);//写入G波段


    poBand3 = dataSet->GetRasterBand(3);//读取B波段

    unsigned char *pBuf3 = new unsigned char[width*height];
    unsigned char *BufB = new unsigned char[width*height];
    poBand3->RasterIO(GF_Read,0,0,width,height,pBuf3,width,height,GDT_Byte,0,0);
    double adfMinMax3[2];
    GDALComputeRasterMinMax((GDALRasterBandH)poBand3,FALSE,adfMinMax3);
    int min3 = adfMinMax3[0];
    int max3 = adfMinMax3[1];
    for(int i=0;i<width*height;i++)
    {
        if(pBuf3[i]<min3)
            BufB[i] = 0;
        else if(pBuf3[i]>max3)
            BufB[i] = 255;
        else
            BufB[i] = 255/(max3 - min3)*(pBuf3[i] - min3);
    }
    pDstband3->RasterIO(GF_Write,0,0,width,height,BufB,width,height,GDT_Byte,0,0);//写入B波段
    //关闭数据集，删除内存
    GDALClose((GDALDatasetH)dataSet);
    GDALClose((GDALDatasetH)poDstDS);
    delete [] pBuf1;
    delete [] BufR;
    delete [] pBuf2;
    delete [] BufG;
    delete [] pBuf3;
    delete [] BufB;

四、Linear2%方法复杂一点，原理如下：我们依然假设增强的图像范围为[0,255]，读取原始图像后首先进行直方图统计，得到0~255每个灰度值对应的个数（如下面的HistBand1），然后用对应的灰度个数除以图像总的像元数（width*height），将得到的结果进行累计（比如0对应的值是0.001,1对应的值是0.002，那么0的累计值就死0.001，而1的累计值就是0对应的值再加上1对应的值，及时0.003），放入HistR、HistG、HistB，这样就得到了累计直方图，然后找到2%对应的灰度值minR和98%对应的灰度值maxR（其他波段类似），将这两个值作为原始图像灰度的最大最小值，运用上述Linear方法即可得到最终给结果。

五、Linear2%算法实现（本算法是在Qt下编译运行的，需要调用GDAL库，同时设置了输出文件方便查看效果）

    GDALDataset *dataSet = (GDALDataset*)GDALOpen(pszSrcFile,GA_ReadOnly);//输入文件路径自己改

    int width = dataSet->GetRasterXSize();
    int height = dataSet->GetRasterYSize();
    GDALDataType dataType = dataSet->GetRasterBand(1)->GetRasterDataType();
    GDALDriver *poDriver = (GDALDriver*)GDALGetDriverByName("GTiff");
    GDALDataset *poDstDS = poDriver->Create(pszDstFile,width,height,3,dataType,NULL);//输出文件路径自己改

    //统计直方图
    unsigned long long HistBand1[256] = {0};
    unsigned long long HistBand2[256] = {0};
    unsigned long long HistBand3[256] = {0};

    //累计直方图
    float HistR[256] = {0};
    float HistG[256] = {0};
    float HistB[256] = {0};

    GDALRasterBand *poBand1,*poBand2,*poBand3;
    GDALRasterBand *pDstband1 = poDstDS->GetRasterBand(1);
    GDALRasterBand *pDstband2 = poDstDS->GetRasterBand(2);
    GDALRasterBand *pDstband3 = poDstDS->GetRasterBand(3);

    //R波段累计直方图
    poBand1 = dataSet->GetRasterBand(1);
    poBand1->GetHistogram(-0.5,255.5,256,HistBand1,TRUE,FALSE,NULL,NULL);
    HistR[0] = (float)HistBand1[0]/(width*height);
    for(int i=1;i<=255;i++)
    {
        HistR[i] = HistR[i-1] + (float)HistBand1[i]/(width*height);
    }

    unsigned char *pBuf1 = new unsigned char[width*height];
    unsigned char *BufR = new unsigned char[width*height];
    poBand1->RasterIO(GF_Read,0,0,width,height,pBuf1,width,height,GDT_Byte,0,0);
    //R波段赋值
    int minR,maxR;//2%处的灰度值和98%处的灰度值
    for(int i=0;i<=255;i++)
    {
        if(HistR[i]<=0.02)
            minR = i;
        if(HistR[i]<=0.98)
            maxR = i;
    }

    for(int i=0;i<width*height;i++)
    {
        if(pBuf1[i]<=minR)
            BufR[i] = 0;
        else if(pBuf1[i]>=maxR)
            BufR[i] = 255;
        else
            BufR[i] = 255/(maxR - minR)*(pBuf1[i] - minR);
    }
    pDstband1->RasterIO(GF_Write,0,0,width,height,BufR,width,height,GDT_Byte,0,0);

    //G波段累计直方图
    poBand2 = dataSet->GetRasterBand(2);
    poBand2->GetHistogram(-0.5,255.5,256,HistBand2,TRUE,FALSE,NULL,NULL);
    HistG[0] = (float)HistBand2[0]/(width*height);
    for(int i=1;i<=255;i++)
    {
        HistG[i] = HistG[i-1] + (float)HistBand2[i]/(width*height);
    }

    unsigned char *pBuf2 = new unsigned char[width*height];
    unsigned char *BufG = new unsigned char[width*height];
    poBand2->RasterIO(GF_Read,0,0,width,height,pBuf2,width,height,GDT_Byte,0,0);
    //G波段赋值
    int minG,maxG;
    for(int i=0;i<=255;i++)
    {
        if(HistG[i]<=0.02)
            minG = i;
        if(HistG[i]<=0.98)
            maxG = i;
    }

    for(int i=0;i<width*height;i++)
    {
        if(pBuf2[i]<=minG)
            BufG[i] = 0;
        else if(pBuf2[i]>=maxG)
            BufG[i] = 255;
        else
            BufG[i] = 255/(maxG - minG)*(pBuf2[i] - minG);
    }
    pDstband2->RasterIO(GF_Write,0,0,width,height,BufG,width,height,GDT_Byte,0,0);

    //B波段累计直方图
    poBand3 = dataSet->GetRasterBand(3);
    poBand3->GetHistogram(-0.5,255.5,256,HistBand3,TRUE,FALSE,NULL,NULL);
    HistB[0] = (float)HistBand3[0]/(width*height);
    for(int i=1;i<=255;i++)
    {
        HistB[i] = HistB[i-1] + (float)HistBand3[i]/(width*height);
    }

    unsigned char *pBuf3 = new unsigned char[width*height];
    unsigned char *BufB = new unsigned char[width*height];
    poBand3->RasterIO(GF_Read,0,0,width,height,pBuf3,width,height,GDT_Byte,0,0);
    //B波段赋值
    int minB,maxB;
    for(int i=0;i<=255;i++)
    {
        if(HistB[i]<=0.02)
            minB = i;
        if(HistB[i]<=0.98)
            maxB = i;
    }

    for(int i=0;i<width*height;i++)
    {
        if(pBuf3[i]<=minB)
            BufB[i] = 0;
        else if(pBuf3[i]>=maxB)
            BufB[i] = 255;
        else
            BufB[i] = 255/(maxB - minB)*(pBuf3[i] - minB);
    }
    pDstband3->RasterIO(GF_Write,0,0,width,height,BufB,width,height,GDT_Byte,0,0);

    GDALClose((GDALDatasetH)dataSet);
    GDALClose((GDALDatasetH)poDstDS);
    delete [] pBuf1;
    delete [] BufR;
    delete [] pBuf2;
    delete [] BufG;
    delete [] pBuf3;
    delete [] BufB;

六、效果对比，和ENVI效果一样，没有ENVI中的红色框框。

原始图像：

ENVI%2线性拉伸算法实现

Linear图像：

ENVI%2线性拉伸算法实现

Linear2%图像：

ENVI%2线性拉伸算法实现

 

 

相关标签： GDAL C

上一篇： SSL简要介绍以及Tomcat的SSL配置

下一篇： 【android】开发干货-图片合成之生成《飞行员驾驶证》

推荐阅读