ENVI%2线性拉伸算法实现

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

三、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中的红色框框。

原始图像：

Linear图像：

Linear2%图像：