（C++）GDAL学习笔记——5全⾊影像与多光谱影像的融合任务将全⾊影像与多光谱影像融合。
原理使⽤的⽅法是HSI变换法，即先将多光谱影像选择RGB三个波段进⾏HSI变化得到H（⾊调）、S（饱和度）、I（亮度）三个波段。
之后按照全⾊波段的尺⼨将H和S两个波段进⾏缩放得到新的H和S波段，⽤全⾊波段替代I波段，将新的H和S波段与全⾊波段组合后进⾏HSI逆变换得到新的RGB波段影像。
HSI变换法的实现这⾥不详细介绍，可以参考这篇帖⼦。
代码/*实验5 将全⾊影像与多光谱影像融合时间 20210712这⾥使⽤的多光谱影像为BGR3波段*/#include<iostream>#include"gdal.h"#include"gdal_priv.h"#include"ogr_spatialref.h"#include<cmath>#include<algorithm>#include<Eigen/Dense>#define PI 3.1415using namespace std;using namespace Eigen;// ⽤来对图像进⾏缩放bool ImResize(MatrixXf srcM,MatrixXf &dstM);int main(){GDALAllRegister();char MulPath[]="D:\\Practice\\org\\exp5_mul.tif";char PanPath[]="D:\\Practice\\org\\exp5_pan.tif";char SavePath[]="D:\\Practice\\res\\exp5_res.tif";GDALDataset* mMulSet =(GDALDataset*)GDALOpen(MulPath, GA_ReadOnly);// 多光谱cout <<" E"<< endl;GDALDataset* mPanSet =(GDALDataset*)GDALOpen(PanPath, GA_ReadOnly);// 全⾊// 获取图像参数int mmX = mMulSet->GetRasterXSize();int mmY = mMulSet->GetRasterYSize();int mpX = mPanSet->GetRasterXSize();int mpY = mPanSet->GetRasterYSize();// 拉伸变换的参数float k1 =float(mpX)/ mmX;float k2 =float(mpY)/ mmY;// 取出波段GDALRasterBand* mMBandB = mMulSet->GetRasterBand(1);GDALRasterBand* mMBandG = mMulSet->GetRasterBand(2);GDALRasterBand* mMBandR = mMulSet->GetRasterBand(3);GDALRasterBand* mPBand = mPanSet->GetRasterBand(1);// 申请缓冲区unsigned short* mMBandB_ =new unsigned short[mmX*mmY];unsigned short* mMBandG_ =new unsigned short[mmX*mmY];unsigned short* mMBandR_ =new unsigned short[mmX*mmY];unsigned char* mPBand_ =new unsigned char[mpX*mpY];// 把数据写⼊缓冲区mMBandR->RasterIO(GF_Read,0,0, mmX, mmY, mMBandR_, mmX, mmY, GDT_Int16,0,0); mMBandG->RasterIO(GF_Read,0,0, mmX, mmY, mMBandG_, mmX, mmY, GDT_Int16,0,0); mMBandB->RasterIO(GF_Read,0,0, mmX, mmY, mMBandB_, mmX, mmY, GDT_Int16,0,0);mPBand->RasterIO(GF_Read,0,0, mpX, mpY, mPBand_, mpX, mpY, GDT_Byte,0,0);// 申请缓冲区IHSMatrixXf I1(mmX, mmY);MatrixXf H1(mmX, mmY);MatrixXf S1(mmX, mmY);// RGB 转 IHSfor(int m =0; m < mmY; m++)for(int n =0; n < mmX; n++){float R = mMBandR_[m*mmX + n];float G = mMBandG_[m*mmX + n];float B = mMBandB_[m*mmX + n];float r = R /(R + G + B);float g = G /(R + G + B);float b = B /(R + G + B);float h, s, i;float theta =acos(0.5*(r - g + r - b)/sqrt(pow(r - g,2)+(r - b)*(g - b)));if(b <= g){ h = theta;}else{ h =2* PI - theta;}s =1-3*min(r,min(g, b));i =(R + G + B)/(3*65535);I1(m,n)= i *65535;S1(m,n)= s *100;H1(m,n)= h *180/ PI;}// 将S1 和 H1 拉伸MatrixXf S2(mpX, mpY);MatrixXf H2(mpX, mpY);MatrixXf I2(mpX, mpY);ImResize(H1, H2);ImResize(S1, S2);//for (int ii = 0; ii < 20; ii++) { cout << H1(ii, ii) <<" "<< H2(ii, ii) << endl; }for(int i =0; i < mpY; i++)for(int j =0; j < mpX; j++){I2(i, j)=float(mPBand_[i*mpX + j]);}cout <<"RGB 转 HSI完成"<< endl;/******************HSI 转 RGB**********************/MatrixXf RR(mpX, mpY);MatrixXf GG(mpX, mpY);MatrixXf BB(mpX, mpY);for(int k =0; k < mpX*mpY; k++){float H =H2(k), S =S2(k), I =I2(k);float h = H*PI /180;float s = S /100;float i = I /255;float r, g, b;if(h >=0&& h <2* PI /3){b = i*(1- s);r = i*(1+ s*cos(h)/cos(PI /3- h));g =3* i -(b + r);}else if(h <4* PI /3){r = i*(1- s);g = i*(1+ s*cos(h -2* PI /3)/cos(PI - h));b =3* i -(g + r);}else{g = i*(1- s);b = i*(1+ s*cos(h -4* PI /3)/cos(5* PI /3- h));r =3* i -(g + b);}RR(k)= r *255;GG(k)= g *255;BB(k)= b *255;}// 将计算得到的RGB申请缓冲区并放⼊unsigned short* mRR =new unsigned short[mpX*mpY];unsigned short* mGG =new unsigned short[mpX*mpY];unsigned short* mBB =new unsigned short[mpX*mpY];for(int y =0; y < mpY; y++){for(int x =0; x < mpX; x++){mRR[y*mpX + x]=(unsigned short)RR(y, x);mGG[y*mpX + x]=(unsigned short)GG(y, x);mBB[y*mpX + x]=(unsigned short)BB(y, x);}}cout <<"HSI 转 RGB 完成"<< endl;// 创建输出数据集GDALDriver* hDriver =GetGDALDriverManager()->GetDriverByName("GTiff");GDALDataset* mSaveSet = hDriver->Create(SavePath, mpX, mpY,3, GDT_UInt16,NULL); // 获取影像的地理信息和投影信息double geoInformation[6];mPanSet->GetGeoTransform(geoInformation);const char* gdalProjection = mPanSet->GetProjectionRef();// 将影像的地理信息写到纠正后影像数据集中// 将影像的地理信息写到纠正后影像数据集中mSaveSet->SetGeoTransform(geoInformation);mSaveSet->SetProjection(gdalProjection);// 将波段数据写⼊mSaveSet->GetRasterBand(1)->RasterIO(GF_Write,0,0, mpX, mpY, mRR, mpX, mpY, GDT_UInt16,0,0); mSaveSet->GetRasterBand(2)->RasterIO(GF_Write,0,0, mpX, mpY, mGG, mpX, mpY, GDT_UInt16,0,0); mSaveSet->GetRasterBand(3)->RasterIO(GF_Write,0,0, mpX, mpY, mBB, mpX, mpY, GDT_UInt16,0,0);delete[] mMBandB_, mMBandG_, mMBandR_;delete[] mRR, mGG, mBB;GDALClose(mSaveSet);GDALClose(mMulSet);GDALClose(mPanSet);GDALDestroyDriverManager();cout <<"程序运⾏完毕"<< endl;getchar();return0;}bool ImResize(MatrixXf srcM, MatrixXf & dstM){int msW = srcM.cols();int msH = srcM.rows();int mdW = dstM.cols();int mdH = dstM.rows();// 计算缩放参数float k1 =float(mdW)/ msW;float k2 =float(mdH)/ msH;for(int h =0; h < mdH; h++){for(int w =0; w < mdW; w++){float tmpX = w / k1;float tmpY = h / k2;dstM(h, w)=srcM(int(tmpY),int(tmpX));}}return false;}结果展⽰多光谱影像：全⾊影像：融合影像：ENVI的Brovey融合结果：结语对⽐了⼀下使⽤ENVI得到的融合结果，是差不多的，不过在使⽤MATLAB进⾏的时候，不知为何融合影像的⾊调有些不同，暂时没有深究，留待以后吧。
