• 【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

    在这里插入图片描述

    
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 <
		}
	}
	
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49

    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);
			
		}
	}
}


    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48

    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;
}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38

    4. 完整程序

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

    #include
#include
#include
#include 
#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 <
		}
	}
	
}

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");

}

    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49
    • 50
    • 51
    • 52
    • 53
    • 54
    • 55
    • 56
    • 57
    • 58
    • 59
    • 60
    • 61
    • 62
    • 63
    • 64
    • 65
    • 66
    • 67
    • 68
    • 69
    • 70
    • 71
    • 72
    • 73
    • 74
    • 75
    • 76
    • 77
    • 78
    • 79
    • 80
    • 81
    • 82
    • 83
    • 84
    • 85
    • 86
    • 87
    • 88
    • 89
    • 90
    • 91
    • 92
    • 93
    • 94
    • 95
    • 96
    • 97
    • 98
    • 99
    • 100
    • 101
    • 102
    • 103
    • 104
    • 105
    • 106
    • 107
    • 108
    • 109
    • 110
    • 111
    • 112
    • 113
    • 114
    • 115
    • 116
    • 117
    • 118
    • 119
    • 120
    • 121
    • 122
    • 123
    • 124
    • 125
    • 126
    • 127
    • 128
    • 129
    • 130
    • 131
    • 132
    • 133
    • 134
    • 135
    • 136
    • 137
    • 138
    • 139
    • 140
    • 141
    • 142
    • 143
    • 144
    • 145
    • 146
    • 147
    • 148
    • 149
    • 150
    • 151
    • 152
    • 153
    • 154
    • 155
    • 156
    • 157
    • 158
    • 159
    • 160
    • 161
    • 162
    • 163
    • 164
    • 165
    • 166
    • 167
    • 168
    • 169
    • 170
    • 171
    • 172
    • 173
    • 174
    • 175
    • 176
    • 177
    • 178
    • 179
    • 180
    • 181
    • 182
    • 183
    • 184
    • 185
    • 186
    • 187
    • 188
    • 189
    • 190
    • 191
    • 192
    • 193
    • 194
    • 195
    • 196
    • 197
    • 198
    • 199
    • 200
    • 201
    • 202
    • 203
    • 204
    • 205
    • 206
    • 207
    • 208
    • 209
    • 210
    • 211
    • 212
    • 213
    • 214
    • 215
    • 216
    • 217
    • 218
    • 219
    • 220
    • 221
    • 222
    • 223
    • 224
    • 225
    • 226
    • 227
    • 228
    • 229
    • 230
    • 231
    • 232
    • 233
    • 234
    • 235
    • 236
    • 237
    • 238
    • 239
    • 240
    • 241
    • 242
    • 243
    • 244
    • 245
  • 相关阅读:
    【nlp】2.4 GRU模型
    自动填充字段值
    appium+python自动化测试
    OTN光层保护
    FRP原理 实战
    Web前端开发者的福音，这款APP让你更上一层楼
    【ue5】滑铲系统蓝图笔记
    vulnhub之Ripper
    HTML+CSS+JS——动漫风二次元论坛(2页) HTML5网页设计成品_学生DW静态网页设计代做_web课程设计网页制作
    【虚幻引擎UE】UE5 实现相机录制视频并导出（C++调用外部exe）
  • 原文地址：https://blog.csdn.net/gsgs1234/article/details/132945952