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

}
  • 2
    点赞
  • 20
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值