目前常用的顏色模型一種是RGB三原色模型，另外一種廣泛採用的顏色模型是亮度、色調、飽和度（(IHS)顏色模型。亮度表示光譜的整體亮度大小，對應於影象的空間資訊屬性，色調描述純色的屬性，決定於光譜的主波長，是光譜在質的方面的區別，飽和度表徵光譜的主波長在強度中的比例，色調和飽和度代表影象的光譜解析度。
IHS變換影象融合就是建立在IHS空間模型的基礎上，其基本思想就是在IHS空間中，將低空間解析度的多光譜影象的亮度分量用高空間解析度的灰度圖象替換。
IHS變換融合法基本流程：

基本步驟：
(1) 將多光譜影像重取樣到與全色影像具有相同的解析度。
(2) 將多光譜影象的Ｒ、Ｇ、Ｂ三個波段轉換到IHS空間，得到Ｉ、Ｈ、Ｓ三個分量；
(3) 以Ｉ分量為參考，對全色影像進行直方圖匹配
(4) 用直方圖匹配後的全色影像替代Ｉ分量，然後同Ｈ、Ｓ分量一起逆變換到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主要是用來顯示。

“`