【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.1415926
using 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");
}
