【c++GDAL】IHS融合

news/2024/7/21 5:25:22 标签: c++, 图像处理

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

}

http://www.niftyadmin.cn/n/5032748.html

相关文章

Json-Jackson和FastJson

Json-Jackson和FastJson

狂神&#xff1a; 测试Jackson 纯Java解决日期格式化 设置ObjectMapper FastJson&#xff1a; 知乎&#xff1a;Jackson使用指南 1、常见配置 方式一&#xff1a;yml配置 spring.jackson.date-format指定日期格式&#xff0c;比如yyyy-MM-dd HH:mm:ss&#xff0c;或者具体的…
阅读更多...
R语言绘制PCA双标图、碎石图、变量载荷图和变量贡献图

R语言绘制PCA双标图、碎石图、变量载荷图和变量贡献图

1、原论文数据双标图 代码&#xff1a; setwd("D:/Desktop/0000/R") #更改路径#导入数据 df <- read.table("Input data.csv", header T, sep ",")# ----------------------------------- #所需的包: packages <- c("ggplot2&quo…
阅读更多...
Kubernetes学习大纲

Kubernetes学习大纲

Kubernetes&#xff08;通常缩写为 K8s&#xff09;是一个用于自动化部署、扩展和管理容器化&#xff08;通常是 Docker 容器&#xff09;应用程序的开源平台。以下是一些 Kubernetes 的核心概念和相关的命令。 1. 集群&#xff08;Cluster&#xff09;和节点&#xff08;Node…
阅读更多...
win部署CRM

win部署CRM

win部署crm&#xff09; 1.phpstudy2.composer3.代码4.其他配置 周末锴哥让我帮他部署了一个CRM&#xff0c;写个教程&#xff0c;方便之后他用。锴哥用的是 NxCrm&#xff0c;先把代码下下来。 1.phpstudy 1.首先是下载小皮面板&#xff0c;配置php的环境。这里面下载了php8…
阅读更多...
线性代数基础-行列式

线性代数基础-行列式

一、行列式之前的概念 1.全排列&#xff1a; 把n个不同的元素排成一列&#xff0c;称为n个元素的全排列&#xff0c;简称排列 &#xff08;实际上就是我们所说的排列组合&#xff0c;符号是A&#xff0c;arrange&#xff09; 2.标准序列&#xff1a; 前一项均小于后一项的序列…
阅读更多...
华为云香港S3云服务器性能测评_99元一年租用价格

华为云香港S3云服务器性能测评_99元一年租用价格

华为云香港S3云服务器1核2G1M带宽99元一年性能测评&#xff0c;配置为S3云服务器1核2G1M带宽&#xff0c;S系列热卖机型&#xff0c;适用于个人建站、普通web应用等负载较低场景&#xff0c;系统盘为高IO40G系统盘&#xff0c;华为云百科分享华为云香港云服务器配置费用&#x…
阅读更多...
SQL Server对象类型（5）——4.5. 同义词（Synonym）

SQL Server对象类型（5）——4.5. 同义词（Synonym）

4.5. 同义词（Synonym） 4.5.1. 同义词概念 与Oracle中相同，SQL Server中的同义词是虚的、被定义的模式对象，其本身并不存储任何数据。其用途之一就是为其他类型基础对象提供一个别名；用途之二就是为应用提供一个抽象层，以方便后期应用相关的基础对象的更改和维护。用户可…
阅读更多...
What does rail-to-rail mean (Rail-to-Rail Op amp) ?

What does rail-to-rail mean (Rail-to-Rail Op amp) ?

What does rail-to-rail mean (Rail-to-Rail Op amp) ?
阅读更多...
最新文章

Copyright @ 2022~2023