【c++GDAL】IHS融合

【c++&GDAL】IHS融合
基于IHS变换融合，实现多光谱和全色影像之间的融合。IHS分别指亮度(I)、色度(H)、饱和度(S)。IHS变换融合基于亮度I进行变换，色度和饱和度空间保持不变。
IHS融合步骤：
（1）将多光谱RGB影像变换到IHS空间；
（2）基于一定融合规则使用亮度分量I与全色影像进行变换，得到新的全色I’，
（3）将I’HS逆变换到RGB空间，得到融合影像。

文章目录

- - 1.RGB2IHS
  - 2.IHS2RGB
  - 3. IHS融合
  - 4. 完整程序

1.RGB2IHS

在这里插入图片描述


void RGBtoHIS(double* R, double* G, double* B, double* pan, int w, int h,double* H,double* I,double* S)
{int sum = w * h * sizeof(double);   //图像所占容量memcpy((double *)H, (double *)R, sum);memcpy((double *)I, (double *)R, sum);memcpy((double *)S, (double *)R, sum);int i, j;double theta = 0,n;for (j = 0; j < h; j++){for (i = 0; i < w; i++){int m = j * w + i;//HIS转换公式中的RGB均需要归一化至[0,1]区间内，matlab的im2double()转换后已然至该区间内R[m] = R[m] / 255;G[m] = G[m] / 255;B[m] = B[m] / 255;//I,S，H分量转弧度，分量范围[0,1],I[m] = (R[m] + G[m] + B[m]) / 3;S[m] = 1 - min(min(R[m], G[m]), B[m]) / I[m];//acos()返回以弧度表示的 x 的反余弦，弧度区间为 [0, pi]theta = acos(0.5*((R[m] - G[m]) + (R[m] - B[m])) / sqrt((R[m] - G[m])*(R[m] - G[m]) + (R[m] - B[m]) * (G[m] - B[m])));theta = theta * 180 / pi;   //转角度if (B[m] <= G[m]){H[m] = theta;}else{H[m] = 360 - theta;}if (S[m] == 0 )    //H的非法值  && S[m]==NULL{H[m] = 0;S[m] = 0;}H[m] = H[m] * 255 /360;S[m] = S[m] * 255;I[m] = I[m] * 255;//cout <<I[m] <<"   ";  //为什么S会出现非法值}}}

2.IHS2RGB

在这里插入图片描述

void HIStoRGB(double* H, double* I, double* S, double* R, double* G, double* B, int w, int h)
{int sum = w * h * sizeof(double);   //图像所占容量memcpy((double *)R, (double *)H, sum);memcpy((double *)G, (double *)S, sum);memcpy((double *)B, (double *)I, sum);int i, j,m;for (j = 0; j < h; j++){for (i = 0; i < w; i++){m = j * w + i;H[m] = H[m] * 360 / 255;   //区间[0,360]S[m] = S[m] / 255;    //S,I的范围都在区间[0,1]上，计算得出R,G,B范围也在区间[0,1]上I[m] = I[m] / 255;if (H[m] >= 0 && H[m] < 120){B[m] = I[m] * (1 - S[m]);R[m] = I[m] * (1 + (S[m] * cos(H[m] * pi / 180)) / cos((60 - H[m])* pi / 180));G[m] = 3 * I[m] - (R[m] + B[m]);}else if (H[m] >= 120 && H[m] < 240){H[m] = H[m] - 120;R[m]= I[m] * (1 - S[m]);G[m] = I[m] * (1 + (S[m] * cos(H[m] * pi / 180)) / cos((60 - H[m])* pi / 180));B[m] = 3 * I[m] - (R[m] + G[m]);}else //(H[m] >= 240 && H[m] < 360){H[m] = H[m] - 240;G[m] = I[m] * (1 - S[m]);B[m] = I[m] * (1 + (S[m] * cos(H[m] * pi / 180)) / cos((60 - H[m])* pi / 180));R[m] = 3 * I[m] - (G[m] + B[m]);}R[m] = max(min(1.0, R[m]), 0.0);G[m] = max(min(1.0, G[m]), 0.0);B[m] = max(min(1.0, B[m]), 0.0);}}
}

3. IHS融合

一般而言融合规则为使用I和pan进行直方图匹配，下列代码使用的融合规则为线性拉伸。融合的步骤即将高分辨率影像进行线性拉伸，使之与多光谱影像亮度分量灰度范围一致，拉伸后的作为新的亮度分量newI。
线性拉伸公式：

void HIS_fusion(double* H, double* I, double* S,double * pan,double *newI,int w,int h)
{int sum = w * h * sizeof(double);   //图像所占容量memcpy((double *)newI, (double *)pan, sum);int i, j;//全色波段与I的直方图匹配double max1, min1, max2, min2;//将高分辨率影像拉伸与亮度分量一致，变换前范围[min1,max1],后[min2,max2]//获取全色影像范围[min1,max1],和多光谱I分量范围[min2,max2]max1 = pan[0]; min1 = pan[0];max2 = I[0]; min2 = I[0];for (i = 0; i < w*h; i++){max1 = max(pan[i], max1);min1 = min(pan[i], min1);max2 = max(I[i], max1);min2 = min(I[i], min1);}double A, B;A = (max2 - min2) / (max1 - min1);B = (max1*min2 - min1 * max2) / (max1 - min1);for (i = 0; i < w*h; i++){	newI[i] = pan[i] * A + B;newI[i] = newI[i] / 255;}GDALDriver* imgDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); const char* outFilename = "Inew.tif";   GDALDataset* o = imgDriver->Create(outFilename,w, h, 1, GDT_Float64, NULL);o->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, w, h, newI, w, h, GDT_Float64, 0, 0);cout << "基于HIS变换的融合完成" << endl;
}

4. 完整程序

在进行匹配前，一般要将多光谱影像采样至全色影像范围内，直接设置RasterIO()的第七八个参数(nBufXSize,nBufYSize)为全色影像的大小，来进行多光谱影像的缩放，GDAL默认最邻近采样。

#include<iostream>
#include<math.h>
#include<iomanip>
#include <algorithm>
#include "gdal_priv.h"
#include "gdalwarper.h"
#define pi 3.1415926using namespace std;void RGBtoHIS(double* R, double* G, double* B, double* pan, int w, int h,double* H,double* I,double* S)
{int sum = w * h * sizeof(double);   //图像所占容量memcpy((double *)H, (double *)R, sum);memcpy((double *)I, (double *)R, sum);memcpy((double *)S, (double *)R, sum);int i, j;double theta = 0,n;for (j = 0; j < h; j++){for (i = 0; i < w; i++){int m = j * w + i;//HIS转换公式中的RGB均需要归一化至[0,1]区间内，matlab的im2double()转换后已然至该区间内R[m] = R[m] / 255;G[m] = G[m] / 255;B[m] = B[m] / 255;//I,S，H分量转弧度，分量范围[0,1],I[m] = (R[m] + G[m] + B[m]) / 3;S[m] = 1 - min(min(R[m], G[m]), B[m]) / I[m];//acos()返回以弧度表示的 x 的反余弦，弧度区间为 [0, pi]theta = acos(0.5*((R[m] - G[m]) + (R[m] - B[m])) / sqrt((R[m] - G[m])*(R[m] - G[m]) + (R[m] - B[m]) * (G[m] - B[m])));theta = theta * 180 / pi;   //转角度if (B[m] <= G[m]){H[m] = theta;}else{H[m] = 360 - theta;}if (S[m] == 0 )    //H的非法值  && S[m]==NULL{H[m] = 0;S[m] = 0;}H[m] = H[m] * 255 /360;S[m] = S[m] * 255;I[m] = I[m] * 255;//cout <<I[m] <<"   ";  //为什么S会出现非法值}}}void HIStoRGB(double* H, double* I, double* S, double* R, double* G, double* B, int w, int h)
{int sum = w * h * sizeof(double);   //图像所占容量memcpy((double *)R, (double *)H, sum);memcpy((double *)G, (double *)S, sum);memcpy((double *)B, (double *)I, sum);int i, j,m;for (j = 0; j < h; j++){for (i = 0; i < w; i++){m = j * w + i;H[m] = H[m] * 360 / 255;   //区间[0,360]S[m] = S[m] / 255;    //S,I的范围都在区间[0,1]上，计算得出R,G,B范围也在区间[0,1]上I[m] = I[m] / 255;if (H[m] >= 0 && H[m] < 120){B[m] = I[m] * (1 - S[m]);R[m] = I[m] * (1 + (S[m] * cos(H[m] * pi / 180)) / cos((60 - H[m])* pi / 180));G[m] = 3 * I[m] - (R[m] + B[m]);}else if (H[m] >= 120 && H[m] < 240){H[m] = H[m] - 120;R[m]= I[m] * (1 - S[m]);G[m] = I[m] * (1 + (S[m] * cos(H[m] * pi / 180)) / cos((60 - H[m])* pi / 180));B[m] = 3 * I[m] - (R[m] + G[m]);}else //(H[m] >= 240 && H[m] < 360){H[m] = H[m] - 240;G[m] = I[m] * (1 - S[m]);B[m] = I[m] * (1 + (S[m] * cos(H[m] * pi / 180)) / cos((60 - H[m])* pi / 180));R[m] = 3 * I[m] - (G[m] + B[m]);}R[m] = max(min(1.0, R[m]), 0.0);G[m] = max(min(1.0, G[m]), 0.0);B[m] = max(min(1.0, B[m]), 0.0);}}
}void HIS_fusion(double* H, double* I, double* S,double * pan,double *newI,int w,int h)
{int sum = w * h * sizeof(double);   //图像所占容量memcpy((double *)newI, (double *)pan, sum);int i, j;//全色波段与I的直方图匹配double max1, min1, max2, min2;//将高分辨率影像拉伸与亮度分量一致，变换前范围[min1,max1],后[min2,max2]max1 = pan[0]; min1 = pan[0];max2 = I[0]; min2 = I[0];for (i = 0; i < w*h; i++){max1 = max(pan[i], max1);min1 = min(pan[i], min1);max2 = max(I[i], max1);min2 = min(I[i], min1);}double A, B;A = (max2 - min2) / (max1 - min1);B = (max1*min2 - min1 * max2) / (max1 - min1);for (i = 0; i < w*h; i++){	newI[i] = pan[i] * A + B;newI[i] = newI[i] / 255;}GDALDriver* imgDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); const char* outFilename = "Inew.tif";   GDALDataset* o = imgDriver->Create(outFilename,w, h, 1, GDT_Float64, NULL);o->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, w, h, newI, w, h, GDT_Float64, 0, 0);cout << "基于HIS变换的融合完成" << endl;
}
void main()
{GDALAllRegister();CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");const char* file1 = "多光谱.tif"; const char* file2 = "全色.tif";  GDALDataset* Mul = (GDALDataset*)GDALOpen(file1, GA_ReadOnly);GDALDataset* Pan = (GDALDataset*)GDALOpen(file2, GA_ReadOnly);if (Mul == NULL || Pan == NULL){cout << "读取图像失败" << endl;}else {cout << "读取成功" << endl;}int MulW = Mul->GetRasterXSize();int MulH = Mul->GetRasterYSize();int MulC = Mul->GetRasterCount();int PanW = Pan->GetRasterXSize();int PanH = Pan->GetRasterYSize();int PanC = Pan->GetRasterCount();GDALDataType Mtype = Mul->GetRasterBand(1)->GetRasterDataType();GDALDataType Ptype = Pan->GetRasterBand(1)->GetRasterDataType();GDALRasterBand* MulR = Mul->GetRasterBand(1);GDALRasterBand* MulG = Mul->GetRasterBand(2);GDALRasterBand* MulB = Mul->GetRasterBand(3);GDALRasterBand* P = Pan->GetRasterBand(1);//Uint16 --多光谱   Uint8 --全色unsigned short* r = new unsigned short[PanW*PanH*Mtype];unsigned short* g= new unsigned short[PanW*PanH*Mtype];unsigned short* b = new unsigned short[PanW*PanH*Mtype];unsigned char* p = new unsigned char[PanW*PanH*Ptype];P->RasterIO(GF_Read, 0, 0, PanW, PanH, p, PanW, PanH, Ptype, 0, 0);//注：设置RasterIO()的第七八个参数(nBufXSize,nBufYSize)为全色影像的大小，来进行多光谱影像的缩放,GDAL默认最邻近采样MulR->RasterIO(GF_Read, 0, 0, MulW, MulH, r , PanW, PanH, Mtype, 0, 0);MulG->RasterIO(GF_Read, 0, 0, MulW, MulH, g, PanW, PanH, Mtype, 0, 0);MulB->RasterIO(GF_Read, 0, 0, MulW, MulH, b, PanW, PanH, Mtype, 0, 0);//类型转换------------------------------------------double* R = new double[PanW*PanH];double* G = new double[PanW*PanH];double* B = new double[PanW*PanH];double* pan = new double[PanW*PanH];int i;for (i = 0; i < PanW*PanH; i++){R[i] = double(r[i]);G[i] = double(g[i]);B[i] = double(b[i]);pan[i] = double(p[i]);}GDALDriver* imgDriver = GetGDALDriverManager()->GetDriverByName("GTiff");  const char* outFilename = "Result.tif"; GDALDataset* out = imgDriver->Create(outFilename, PanW, PanH ,MulC, GDT_Float64, NULL);double* H = new double[PanW*PanH];double* I = new double[PanW*PanH];double* S = new double[PanW*PanH];RGBtoHIS(R,G,B,pan, PanW, PanH, H, I, S);double* newI = new double[PanW*PanH];HIS_fusion(H, I, S, pan, newI, PanW, PanH);   //全色波段拉伸替代I分量//最后融合的结果以RGB的形式显示double* newr = new double[PanW*PanH];double* newg = new double[PanW*PanH];double* newb = new double[PanW*PanH];HIStoRGB(H, newI, S, newr, newg, newb, PanW, PanH);out->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, PanW, PanH, newr, PanW, PanH, GDT_Float64, 0, 0);out->GetRasterBand(2)->RasterIO(GF_Write, 0, 0, PanW, PanH, newg, PanW, PanH, GDT_Float64, 0, 0);out->GetRasterBand(3)->RasterIO(GF_Write, 0, 0, PanW, PanH, newb, PanW, PanH, GDT_Float64, 0, 0);/*计算得出R,G,B范围也在区间[0,1]上则以GDT_Float64存储，若转换到区间[0,255]上，若是char类型的以GDT_Byte存储*/GDALClose(Mul);GDALClose(Pan);GDALClose(out);delete R, G, B, P;delete r,g,b,pan;delete H,I,S,newI;delete newr, newg, newb;system("pause");}