【c++&GDAL】IHS融合

2023-09-17 17:13:41

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

}
更多推荐

记一次 .NET 某仪器测量系统 CPU爆高分析

一：背景1.讲故事最近也挺奇怪，看到了两起CPU爆高的案例，且诱因也是一致的，觉得有一些代表性，合并分享出来帮助大家来避坑吧，闲话不多说，直接上windbg分析。二：WinDbg分析1.CPU真的爆高吗这里要提醒一下，别人说爆高不一定真的就是爆高，我们一定要拿数据说话，可以用!tp观察下。0:000>!tplogSta

关于ITSS认证资质整改和降级

最近来我公司咨询ITSS年审换证的企业比较多，很多企业伙伴的ITSS信息技术服务运维符合性证书2023年到期将面临换证，很多企业觉得拿证三年了都没有问题，换证随便糊弄一下就行了。但是在年审换证再评估中会遇到很多问题，就一些问题我们武汉好地科技小编整理出了一套年审换证再评估相关注意事项。ITSS认证再评估相关事项：1．持

ELK日志分析系统

日志服务器提高安全性集中存放日志缺陷：对日志的分析困难ELK日志分析系统ElasticsearchLogstashKibana介绍ELK日志分析系统是一套完整的日志集中处理解决方案，基于Elasticsearch、Logstash、Kibana三种开源工具进行日志收集、存储和可视化elk可以帮助用户快速定位和分析应用程

Vue3函数式编程

文章目录前言一、三种编程风格1.template2.jsx/tsx3.函数式编写风格二、函数式编程1.使用场景2.参数3.例子3.render渲染函数总结前言本文主要记录vue3中的函数式编程以及其他编程风格的简介一、三种编程风格1.templateVue使用一种基于HTML的模板语法，使我们能够声明式地将其组件实例的

SPA首屏加载速度慢

什么是首屏加载首屏时间（FirstContentfulPaint），指的是浏览器从响应用户输入网址地址，到首屏内容渲染完成的时间，此时整个网页不一定要全部渲染完成，但需要展示当前视窗需要的内容首屏加载可以说是用户体验中最重要的环节关于计算首屏时间通过DOMContentLoad或者performance来计算出首屏时间

亚马逊登山扣CPC认证ASTMF1774测试和UIAA121测试报告申请

一.什么是登山扣答:登山扣是扣子的一种，顾名思义其就是用来在登山的时候配合绳子起到一个承重悬挂的作用.采用铝吕合金、铁或者是不锈钢等材料制作而成的一种登山工具之一。其形状多样，比较常见的是椭圆形和圆形的，除此之外还有长方形、三角形等样式的登山扣。铝合金登山扣由于质地较轻所以重量也比较交轻，所以携带方便，其耐腐蚀性和防锈

算法通关村-----LRU的设计与实现

LRU缓存问题描述请你设计并实现一个满足LRU(最近最少使用)缓存约束的数据结构。实现LRUCache类：LRUCache(intcapacity)以正整数作为容量capacity初始化LRU缓存。intget(intkey)如果关键字key存在于缓存中，则返回关键字的值，否则返回-1。voidput(intkey,i

YashanDB第三期YCA认证培训圆满结束

9月11日，由YashanDB举办的“第三期YCA认证培训”圆满结束。本次培训吸引了超过110名学员报名，华润数字科技有限公司、北京中亦安图科技股份有限公司、迪思杰(北京)数据管理技术有限公司等多家合作伙伴积极参与。经过7天的学习，在9月6日举行的YCA认证考试中，77名学员获得“YashanDB数据库V22.2认证管

Linux命令行操作：使用“more“命令进行分页显示

文章目录1.引言1.1介绍Linux操作系统和命令行界面什么是Linux操作系统？为什么命令行界面在Linux中如此重要？1.2介绍Linux中的分页显示命令分页显示命令的作用与意义不同分页显示命令的比较2."more"命令的基本用法2.1安装和启动"more"命令如何安装"more"命令？如何从命令行中启动"more

golang使用高阶函数优化业务功能

业务描述两个接口（新增Tag和更新Tag），在业务层均需要添加两个校验，校验Tag信息是否重复和Tag的数据中的编码是否重复。基本实现方式对应的增加两个校验的函数/方法，在接口实现中依次调用两个函数/方法进行校验。优缺点实现简单；但重复代码多，后期再增加其他校验，扩展性较差。高阶函数方式一方式因为业务方法参数相同，业务

FOXBORO FBM233 P0926GX控制脉冲模块

FOXBOROFBM233P0926GX是一种控制脉冲模块，通常用于工业自动化和控制系统中。这个模块的主要功能是生成和控制脉冲信号，以用于执行特定的操作或控制过程。以下是可能适用于FOXBOROFBM233P0926GX控制脉冲模块的一些常见特点：脉冲生成：FBM233P0926GX模块通常能够生成可控的脉冲信号，包括

热文推荐
  1. WebGL 选中物体

    目录前言如何实现选中物体示例程序（PickObject.js）代码详解gl.readPixels（）函数规范示例效果前言有些三维应用程序需要允许用户能够交互地操纵三维物体，要这样做首先就得允许用户选中某个物体。对物体进行选中操作的用处很广泛。比如，让用户选中三维用户界面上的一个按钮，或者让用户选中三维场景中的多张照片中

  2. 提升前端开发效率：基于vue的van-radio-group组件封装指南

    前言vant作为一款流行的ui框架，其中，van-radio-group组件是一个常用的单选框组件，但有时我们需要根据项目需求进行定制化封装。本文将介绍如何基于vue框架封装van-radio-group组件，让我们一起来探索吧！封装文件在这个组件中，使用了vant框架提供的van-radio-group和van-ra

  3. Linux-网卡和网络配置

    链接一篇大佬的博客：Linux之手把手教会修改网卡名称文章目录修改网卡名称步骤1：修改“/etc/default/grub”步骤2：修改“/etc/sysconfig/network-scripts”下的文件步骤3：修改“ifcfg-eth0”配置步骤4：判断操作系统的引导模式步骤5：根据不同的引导模式重新读取配置文件

  4. el-table 列背景色渐变

    最初的想法是，给每一行添加背景色，逐行递减透明度，发现结果比较突兀，效果如下：如果有需要这种样式的，代码如下：<template><div><el-table:data="tableData":header-cell-style="headerCellStyle":cell-style="cellStyle"style

  5. Redis Part1

    单体架构：一台Web服务器、一台数据库服务器。回顾，关系型数据库：基于二维表来存储数据的数据库就是关系型数据库。MySQL跟Redis的区别：MySQL是关系型数据库，它是基于表来存储数据的，MySQL数据是写在磁盘的，它是跟磁盘进行交互的；Redis是非关系型数据库，它是把数据存储在内存当中的，是跟内存进行交互的。M

  6. Go基础语法：数组

    6数组6.1数组的定义数组在定义时就需要声明其元素数量和类型：//T即元素类型var数组变量名[元素数量]T如：vara[5]int.数组的长度必须是常量，并且长度是数组类型的一部分，一旦定义，长度不能变。所以，[5]int和[10]int是不同的类型。packagemainimport"fmt"funcmain(){

  7. 功能基础篇3——Python中的输入输出、文件读写、序列化

    IO文件读写openopen()为内置函数，用于读写文件mode读写x，create，创建，文件存在报错，可写不可读r，read，读入，默认读写方式，文件不存在报错，可读不可写w，write，写入，文件不存在会创建文件，存在清空文件内容，可写不可读a，append，追加，文件不存在会创建文件，光标移至文件末尾，可写不可

  8. CentOS常用命令

    一、yum命令1.1列表源服务器软件包yumlist//列出yum源服务器上所有的软件包，太多了，一般不单独使用，大多数要配合管道grep使用yumlistopenssl*//列出yum源服务器上所有以openssl开头的软件包yumlist|grepopenssl//列出yum源服务器上所有包含openssl的软件包

  9. 实时可追溯性和事后可追溯性

    在当今的世界里，成功的组织都将其创新和复杂产品开发流程建立在跨多个工作流的互联数据上。为此，他们需要收集利益相关者的想法以构建系统架构，管理高级需求以创建详细的用户故事，并实施验证程序以检查问题。同时，在满足标准规定并在当前的市场环境中保持竞争力的过程中，与各利益相关者的协作需要深度的追踪。许多行业标准都要求有需求追踪