遥感影像多光谱与全色IHS融合

目前常用的颜色模型一种是RGB三原色模型,另外一种广泛采用的颜色模型是亮度、色调、饱和度((IHS)颜色模型。亮度表示光谱的整体亮度大小,对应于图像的空间信息属性,色调描述纯色的属性,决定于光谱的主波长,是光谱在质的方面的区别,饱和度表征光谱的主波长在强度中的比例,色调和饱和度代表图像的光谱分辨率。
IHS变换图像融合就是建立在IHS空间模型的基础上,其基本思想就是在IHS空间中,将低空间分辨率的多光谱图像的亮度分量用高空间分辨率的灰度图象替换。
IHS变换融合法基本流程:
这里写图片描述
基本步骤:
(1) 将多光谱影像重采样到与全色影像具有相同的分辨率。
(2) 将多光谱图像的R、G、B三个波段转换到IHS空间,得到I、H、S三个分量;
(3) 以I分量为参考,对全色影像进行直方图匹配
(4) 用直方图匹配后的全色影像替代I分量,然后同H、S分量一起逆变换到RGB空间,从而得到融合图象,

RGB与IHS变换的基本公式:
RGB->HIS:
这里写图片描述
HIS->RGB:
这里写图片描述这里写图片描述这里写图片描述

代码:

 void RGBtoHIS(uchar r,uchar g,uchar b,float* HIS)
 {
    float H;
    float I;
    float S;

    float angle;
    float R,G,B;
    R=((float)r)/255.0;         //归一化
    G=((float)g)/255.0;
    B=((float)b)/255.0;


    I=(R+G+B)/3.0;
    S=1-min(min(R,G),B)/I;
    if (S ==0)
    {
        H = 0;
    }
    if (r == g && r ==b)
    {
        b = 0;
    }
    else
    {
        angle=acos(0.5*(R-G+R-B)/(sqrt((R-G)*(R-G)+(R-B)*(G-B))));
        angle=angle*180.0/PI;
    }
    if (B>G)
    {
        H=360-angle;        //H量计算
    }
    else
    {
        H=angle;
    }


    H=H*255.0/360.0;
    I=I*255.0;
    S=S*255.0;

    HIS[0]=H;
    HIS[1]=I;
    HIS[2]=S;
 }
void HIStoRGB(float H,float I,float S,float* m)
{
    float s=(S)/255;
    float i=(I)/255;
    float h=(H)*360/255;
    if (h<0)
    {h=h+360;
    }
    if(h<120&&h>=0)
    {
        m[0]=i*(1.0+((s*cos(h*PI/180))/cos((60-h)*PI/180)));
        m[2]=i*(1.0-s);
        m[1]=3.0*i-(m[0]+m[2]);
    }

    if(h<240&&h>=120)
    {
        h=h-120;
        m[1]=i*(1.0+((s*cos(h*PI/180))/cos((60-h)*PI/180)));
        m[0]=i*(1.0-s);
        m[2]=3.0*i-(m[0]+m[1]);
    }
    if(h>=240&&h<=360)
    {
        h=h-240;
        m[2]=i*(1.0+((s*cos(h*PI/180))/cos((60-h)*PI/180)));
        m[1]=i*(1.0-s);
        m[0]=3*i-(m[1]+m[2]);
    }
    m[0]=(m[0]*255);
    m[1]=(m[1]*255);
    m[2]=(m[2]*255);
}
void imageFusion()
{
    dlg=new InPutDialog(this);
    dlg->exec();

    QString pan = "D:/img/GF1_pan_1_1.tif";
    QByteArray ba = pan.toUtf8();
    const char *str = ba.data();
    GDALDataset *picPan = (GDALDataset*)GDALOpen(str,GA_ReadOnly);

    QByteArray ba1 =dlg->stringpicR.toUtf8();
    const char *str1=ba1.data();
    GDALDataset *picR1=(GDALDataset*)GDALOpen(str1,GA_ReadOnly);

    const char* path1 = "D:\\Data\\matchR.tif";
    HistogramMatch(picR1,picPan,path1);

    QByteArray ba2 =dlg->stringpicG.toUtf8();
    const char *str2=ba2.data();
    GDALDataset *picG1=(GDALDataset*)GDALOpen(str2,GA_ReadOnly);

    const char* path2 = "D:\\Data\\matchG.tif";
    HistogramMatch(picG1,picPan,path2);

    QByteArray ba3 =dlg->stringpicB.toUtf8();
    const char *str3=ba3.data();
    GDALDataset *picB1=(GDALDataset*)GDALOpen(str3,GA_ReadOnly);

    const char* path3 = "D:\\Data\\matchB.tif";
    HistogramMatch(picB1,picPan,path3);

    GDALDataset *picR=(GDALDataset *)GDALOpen("D:\\Data\\pan.tif",GA_ReadOnly); 
    GDALDataset *picG=(GDALDataset *)GDALOpen("D:\\Data\\pan.tif",GA_ReadOnly); 
    GDALDataset *picB=(GDALDataset *)GDALOpen("D:\\Data\\pan.tif",GA_ReadOnly); 

    QByteArray ba4 =dlg->stringpicTotal.toUtf8();
    const char *str4=ba4.data();
    GDALDataset *picT=(GDALDataset*)GDALOpen(str4,GA_ReadOnly);

    int isrcwidth=picR->GetRasterXSize();
    int isrcheight=picR->GetRasterYSize();

    float *RDataBuff=new float[isrcwidth];
    float *GDataBuff=new float[isrcwidth];
    float *BDataBuff=new float[isrcwidth];
    float *TotalBuffer=new float[isrcwidth];

    GDALDriver *poDriver1=(GDALDriver*)GDALGetDriverByName("GTiff");
    GDALDataset *pDstH=poDriver1->Create("D:\\Data\\H.tif",isrcwidth,isrcheight,1,GDT_Byte,NULL);
    GDALDriver *poDriver2=(GDALDriver*)GDALGetDriverByName("GTiff");
    GDALDataset *pDstI=poDriver2->Create("D:\\Data\\I.tif",isrcwidth,isrcheight,1,GDT_Byte,NULL);
    GDALDriver *poDriver3=(GDALDriver*)GDALGetDriverByName("GTiff");
    GDALDataset *pDstS=poDriver3->Create("D:\\Data\\S.tif",isrcwidth,isrcheight,1,GDT_Byte,NULL);
    GDALDriver *poDriver4=(GDALDriver*)GDALGetDriverByName("GTiff");
    //GDALDataset *HighPic=poDriver4->Create("D:\\Data\\result.tif",isrcwidth,isrcheight,3,GDT_Float32,NULL);
    GDALDataset *HighPic=poDriver4->Create("D:\\Data\\result.tif",isrcwidth,isrcheight,3,GDT_Byte,NULL);

    float *HBuffer=new float[isrcwidth];
    float *IBuffer=new float[isrcwidth];
    float *SBuffer=new float[isrcwidth];
    uchar m1,m2,m3;
    float m[3];
    QProgressDialog *dlg1=new QProgressDialog("RGB转HIS","OK",0,isrcheight,this);
    QProgressDialog *dlg2=new QProgressDialog("HIS转RGB","OK",0,isrcheight,this);
    for(int i=0;i<isrcheight;i++)
    {
        dlg1->setWindowModality(Qt::WindowModal);
        dlg1->show();
        picR->GetRasterBand(1)->RasterIO(GF_Read,0,i,isrcwidth,1,RDataBuff,isrcwidth,1,GDT_Float32,0,0);
        picG->GetRasterBand(1)->RasterIO(GF_Read,0,i,isrcwidth,1,GDataBuff,isrcwidth,1,GDT_Float32,0,0);
        picB->GetRasterBand(1)->RasterIO(GF_Read,0,i,isrcwidth,1,BDataBuff,isrcwidth,1,GDT_Float32,0,0);
        for(int j=0;j<isrcwidth;j++)
        {
            m1=RDataBuff[j];
            m2=GDataBuff[j];
            m3=BDataBuff[j];
            RGBtoHIS(m1,m2,m3,m);
            HBuffer[j]=m[0];
            IBuffer[j]=m[1];
            SBuffer[j]=m[2];
        }
        dlg1->setValue(i);
        QCoreApplication::processEvents();
        if (dlg1->wasCanceled())
        {break;
        }
        pDstH->GetRasterBand(1)->RasterIO(GF_Write,0,i,isrcwidth,1,HBuffer,isrcwidth,1,GDT_Float32,0,0);
        pDstI->GetRasterBand(1)->RasterIO(GF_Write,0,i,isrcwidth,1,IBuffer,isrcwidth,1,GDT_Float32,0,0);
        pDstS->GetRasterBand(1)->RasterIO(GF_Write,0,i,isrcwidth,1,SBuffer,isrcwidth,1,GDT_Float32,0,0);
    }
    dlg1->setValue(isrcheight);

    const char* path = "D:\\Data\\pan.tif";

    HistogramMatch(picT,pDstI,path);

    float *RBuff=new float[isrcwidth];
    float *GBuff=new float[isrcwidth];
    float *BBuff=new float[isrcwidth];
    //DT_8U *RBuff=new DT_8U[isrcwidth];
    //DT_8U *GBuff=new DT_8U[isrcwidth];
    //DT_8U *BBuff=new DT_8U[isrcwidth];

    float n[3];
    float n1,n2,n3;
    GDALDataset *trans_I=(GDALDataset *)GDALOpen("D:\\Data\\pan.tif",GA_ReadOnly);

    for(int i=0;i<isrcheight;i++)
    {
        dlg2->setWindowModality(Qt::WindowModal);
        dlg2->show();
        pDstH->GetRasterBand(1)->RasterIO(GF_Read,0,i,isrcwidth,1,HBuffer,isrcwidth,1,GDT_Float32,0,0);
        pDstS->GetRasterBand(1)->RasterIO(GF_Read,0,i,isrcwidth,1,SBuffer,isrcwidth,1,GDT_Float32,0,0);
        trans_I->GetRasterBand(1)->RasterIO(GF_Read,0,i,isrcwidth,1,TotalBuffer,isrcwidth,1,GDT_Float32,0,0);
        for (int j=0;j<isrcwidth;j++)
        {
            n1=HBuffer[j];
            n2=TotalBuffer[j];
            n3=SBuffer[j];
            HIStoRGB(n1,n2,n3,n);
            RBuff[j]=n[0];
            GBuff[j]=n[1];
            BBuff[j]=n[2];
        }
        dlg2->setValue(i);
        QCoreApplication::processEvents();
        if (dlg2->wasCanceled())
        {break;
        }
        HighPic->GetRasterBand(1)->RasterIO(GF_Write,0,i,isrcwidth,1,RBuff,isrcwidth,1,GDT_Float32,0,0);
        HighPic->GetRasterBand(2)->RasterIO(GF_Write,0,i,isrcwidth,1,GBuff,isrcwidth,1,GDT_Float32,0,0);
        HighPic->GetRasterBand(3)->RasterIO(GF_Write,0,i,isrcwidth,1,BBuff,isrcwidth,1,GDT_Float32,0,0);
    }
    dlg2->setValue(isrcheight);

    delete[] RBuff;
    delete[] GBuff;
    delete[] BBuff;
    delete[] RDataBuff;
    delete[] GDataBuff;
    delete[] BDataBuff; 
    delete[] TotalBuffer;
    delete[] HBuffer;
    delete[] IBuffer;
    delete[] SBuffer;
    GDALClose((GDALDatasetH)trans_I);
    GDALClose((GDALDatasetH)picR);
    GDALClose((GDALDatasetH)picG);
    GDALClose((GDALDatasetH)picB);
    GDALClose((GDALDatasetH)picT);
    GDALClose((GDALDatasetH)pDstS);
    GDALClose((GDALDatasetH)pDstI);
    GDALClose((GDALDatasetH)pDstH);
    GDALClose((GDALDatasetH)HighPic);
}
void HistogramMatch(GDALDataset * poDSPan, GDALDataset * poDSI, const char* str)
{
    int iWidth = poDSI->GetRasterXSize();
    int iHeight = poDSI->GetRasterYSize();
    int Size = iWidth * iHeight;
    //计算I分量的累计直方图
    float fpPro[256];
    memset( fpPro, 0, sizeof(float) * 256 );
    GDALRasterBand *poBand_I = poDSI->GetRasterBand(1);
    int anHistogram_I[256];
    poBand_I->GetHistogram( -0.5, 255.5, 256, anHistogram_I, FALSE, FALSE, GDALDummyProgress, NULL );
    //直方图归一化          
    for (int i = 0; i < 256; i++)
    {
        fpPro[i] = (float) anHistogram_I[i] / Size;
    }
    //计算累计直方图  
    float temp1[256];
    for (int i = 0; i < 256; i++)
    {
        if( i == 0 )
            temp1[0] = fpPro[0];
        else
            temp1[i] = temp1[i-1] + fpPro[i];
        fpPro[i] = temp1[i];
    } 

    //计算Pan波段的累计直方图
    float pPro[256];
    memset( pPro, 0, sizeof(float) * 256 );
    GDALRasterBand *poBand_Pan = poDSPan->GetRasterBand(1);
    int anHistogram_Pan[256];
    poBand_Pan->GetHistogram( -0.5, 255.5, 256, anHistogram_Pan, FALSE, FALSE, GDALDummyProgress, NULL );
    //直方图归一化          
    for (int i = 0; i < 256; i++)
    {
        pPro[i] = (float) anHistogram_Pan[i] / Size;
    }
    //计算累计直方图
    float temp2[256];
    for (int i = 0; i < 256; i++)
    {
        if( i == 0 )
            temp2[0] = pPro[0];
        else
            temp2[i] = temp2[i-1] + pPro[i];
        pPro[i] = temp2[i];
    } 

    //确定映射关系//
    unsigned char Map[256];
    for ( int i = 0; i < 256; i++ )
    {
        int c = 0;                           //最接近的规定直方图灰度级变量//
        float min_value = 1.0;              //最小差值变量//
        //枚举规定直方图各个灰度//
        for ( int j = 0; j < 256; j++ )
        {
            float now_value = 0.0;         //当前插值变量//
            //差值计算//
            if( pPro[i]-fpPro[j] >= 0.0f )
                now_value = pPro[i] - fpPro[j];
            else
                now_value = fpPro[j] - pPro[i];
            //寻找最接近的规定直方图灰度级//
            if( now_value < min_value )
            {
                c = j;                      //最接近的灰度级//
                min_value = now_value;      //最小差值//
            }
        }
        //建立灰度映射表//
        Map[i] = unsigned char( c );
    }
    //一一进行映射
    unsigned char * poDataPan = new unsigned char[Size];
    unsigned char * poNewDataPan = new unsigned char[Size];
    poBand_Pan->RasterIO(GF_Read, 0, 0, iWidth, iHeight, poDataPan, iWidth, iHeight, GDT_Byte,0,0);
    int m = 0;
    for ( int ny = 0; ny < iHeight; ny++ )
    {
        for ( int nx = 0; nx < iWidth; nx++ )
        {   
            unsigned char x = poDataPan[m];
            poNewDataPan[m] = (unsigned char)( Map[x] + 0.5f );
            m++;
        }
    }


    GDALDriver * poDriver = (GDALDriver *)GDALGetDriverByName("GTiff");
    GDALDataset * poDSResult = poDriver->Create(str,iWidth,iHeight,1,GDT_Float32,NULL);
//  GDALDataset * poDSResult = poDriver->Create("D:\\Data\\pan.tif",iWidth,iHeight,1,GDT_Float32,NULL);
    poDSResult->GetRasterBand(1)->RasterIO(GF_Write,0,0,iWidth,iHeight,poNewDataPan,iWidth,iHeight,GDT_Byte,0,0);
    GDALClose((GDALDatasetH) poDSResult); 
}

代码主要包含了三个函数,顾名思义为从RGB->HIS、HIS->RGB和直方图匹配三个函数,前两个比较好搞定,直方图匹配费了较大力气,其中用了Qt和gdal, Qt主要是用来显示。

“`

  • 7
    点赞
  • 100
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值