ENVI%2线性拉伸算法实现

落日映苍穹つ 2022-05-15 04:27 322阅读 0赞

一、在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，a与b之间的按照公式计算**即可。

![70][]

三、**Linear算法实现**（本算法是在Qt下编译运行的，C++语言，调用了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中的红色框框。

**原始图像：**

![70 1][]

**Linear图像：**

![70 2][]

**Linear2%图像：**

![70 3][]

[70]: /images/20220515/458c4d81866344bdb9d8a4eb0bf740fb.png
[70 1]: /images/20220515/86e3cab2c5b047df82fe45d1a549493f.png
[70 2]: /images/20220515/a2c16e010ba149afba4a19704f7146ac.png
[70 3]: /images/20220515/fac21650755c4dc58f2d9fb1d1fa8316.png