【影像配准】遥感影像配准结果输出tif影像（附有完整代码）

        遥感影像经过配准后会输出两幅影像：配准后的参考影像和待配准影像；

        注：png 转 tif

        因为两幅多时像影像在拍摄时即便是同源也无法保证拍摄到的景象完全一致，换句话说：配准时为了使两幅影像的同名点能够对应上，两幅影像（参考影像和待配准影像）的尺寸难免会发生变化（多出的部分使用0填充），甚至待配准影像还会发生些微旋转；

        因此在配准结果输出 .tif 影像时绝对不是简单的 png 转 tif ，下面对遥感影像配准结果输出 .tif 做出了具体实现；

        具体原理：根据投影坐标和像素坐标的换算关系，计算出新 .tif 影像的投影信息，然后新建一个新 .tif 影像，把新 .tif 影像的投影信息写进去，并把配准后的待配准影像（png）的相关信息写入新 .tif 影像，最后生成配准后的待配准影像（.tif）,配准后的参考影像（.tif）的生成以此类推；

说明：初始遥感影像经过配准后输出两幅影像为 png；

下面具体给出代码上的实现和具体的编译环境！！！

一、代码实现

1、pngTotif.h

#pragma once
#include <iostream>
#include"gdal_priv.h"
#include <gdal.h>
#include <gdal_priv.h>
#include <streambuf>
#include <atlstr.h>

#include <opencv2\opencv.hpp>
#include<opencv2\imgproc\types_c.h>
#include<opencv2\imgproc\imgproc.hpp>
#include<opencv2\highgui\highgui.hpp>
#include<opencv2\features2d\features2d.hpp>		//特征提取

using namespace std;


class pngTotif
{
public:
	pngTotif();

	//获取tif影像的投影信息
	void getProj(double* geo);

	//png转tif
	void pngtotif(string srcTifPath, string srcPngPath, string dstTifPath);

	~pngTotif();
};

2、pngTotif.cpp

#include "pngTotif.h"


pngTotif::pngTotif() {}


//获取tif影像的投影信息
void pngTotif::getProj(double* geo)
{
	GDALAllRegister();

	GDALDataset* podataset = (GDALDataset*)GDALOpen("..\\test_data\\02_041.tif", GA_ReadOnly);

	//保存 .tif 文件的投影信息' x=geo[0]、y=geo[3]、像素大小=geo[1]
	//double geo[6];

	podataset->GetGeoTransform(geo);
}


//png转tif
void pngTotif::pngtotif(string srcTifPath, string srcPngPath, string dstTifPath)
{
	//cv::Mat image_1 = cv::imread("..\\test_data\\02_04.tif", -1);
	//cv::Mat image_2 = cv::imread("..\\test_data\\1.png", -1);

	cv::Mat image_1 = cv::imread(srcPngPath, -1);
	cv::Mat image_2 = cv::imread(srcPngPath, -1);

	/*************** 获取配准后影像的左上角顶点 **************/

	cv::Point point;

	//获取配准后影像(出去黑边后)的左上角地理坐标
	for (int i = 0; i < 30; i++)
	{
		for (int j = 0; j < 30; j++)
		{
			//注意行和列与i，j的对应关系
			//左上角在原点
			if ((image_2.at<cv::Vec3b>(i, j)[0] > 30 || image_2.at<cv::Vec3b>(i, j)[1] > 30 || image_2.at<cv::Vec3b>(i, j)[2] > 30) && i == 0 && j == 0)
			{
				point.x = j;
				point.y = i;
				break;
			}
			//左上角在左边界
			else if ((image_2.at<cv::Vec3b>(i, j)[0] > 30 || image_2.at<cv::Vec3b>(i, j)[1] > 30 || image_2.at<cv::Vec3b>(i, j)[2] > 30) && j == 0 && i > 0 &&
				(image_2.at<cv::Vec3b>(i - 1, j)[0] < 30 && image_2.at<cv::Vec3b>(i - 1, j)[1] < 30 && image_2.at<cv::Vec3b>(i - 1, j)[2] == 0))
			{
				point.x = j;
				point.y = i;
				break;
			}
			//左上角在上边界
			else if ((image_2.at<cv::Vec3b>(i, j)[0] > 30 || image_2.at<cv::Vec3b>(i, j)[1] > 30 || image_2.at<cv::Vec3b>(i, j)[2] > 30) && i == 0 && j > 0 &&
				(image_2.at<cv::Vec3b>(i, j - 1)[0] < 30 && image_2.at<cv::Vec3b>(i, j - 1)[1] < 30 && image_2.at<cv::Vec3b>(i, j - 1)[2] < 30))
			{
				point.x = j;
				point.y = i;
				break;
			}
			//左上角在中间情况
			else if ((image_2.at<cv::Vec3b>(i, j)[0] > 30 || image_2.at<cv::Vec3b>(i, j)[1] > 30 || image_2.at<cv::Vec3b>(i, j)[2] > 30) && i > 0 && j > 0 &&
				(image_2.at<cv::Vec3b>(i - 1, j)[0] < 30 && image_2.at<cv::Vec3b>(i - 1, j)[1] < 30 && image_2.at<cv::Vec3b>(i - 1, j)[2] < 30) &&
				(image_2.at<cv::Vec3b>(i, j - 1)[0] < 30 && image_2.at<cv::Vec3b>(i, j - 1)[1] < 30 && image_2.at<cv::Vec3b>(i, j - 1)[2] < 30))
			{
				point.x = j;
				point.y = i;
				break;
			}
		}
	}

	cout << "point.x=" << point.x << "   " << "point.y=" << point.y << endl;

	//获取tif影像的投影信息
	//保存 .tif 文件的投影信息' x=geo[0]、y=geo[3]、像素大小=geo[1]
	double geo[6];

	getProj(geo);

	//x=geo[0]、y=geo[3]是左上角像素坐标对应的投影坐标

	//获取新tif影像的左上角投影坐标
	cv::Point2d point0;

	point0.x = geo[0] - geo[1] * point.x;
	point0.y = geo[3] + geo[1] * point.y;

	cout << "point0.x=" << point0.x << "   " << "point0.y" << point0.y << endl;

	/***************** 创建新的tif影像 ******************/

	//注意图像名称不能是中文
	//string srcPngPath = "../test_data/1.png";						//读取图像png
	//string dstTifPath = "./image_save/2020_02_04.tif";					//新tif保存路径

	char* path1 = (char*)srcPngPath.c_str();
	char* path2 = (char*)dstTifPath.c_str();

	GDALDataset* poDataset1; //待创建的GDAL数据集 
	GDALDataset* poDataset2; //png图信息
	GDALDriver* poDriver;    //驱动，用于创建新的文件

	//gdal打开png图像
	GDALAllRegister();
	poDataset2 = (GDALDataset*)GDALOpen(path1, GA_ReadOnly);

	//判断是否打开文件
	if (poDataset2 == NULL) { printf("文件打开失败"); }
	else printf("文件打开成功");

	//获取待加入坐标的图像的通道长宽
	int nBandCount = poDataset2->GetRasterCount();
	int nImgSizeX = poDataset2->GetRasterXSize();
	int nImgSizeY = poDataset2->GetRasterYSize();

	//图像的仿射信息见geo[6]，需要修改仿射信息，以适应新的tif，即修改坐上角坐标信息（重点!!!）
	geo[0] = point0.x;
	geo[3] = point0.y;

	//创建新文件
	char fomat[] = "GTiff";
	poDriver = GetGDALDriverManager()->GetDriverByName(fomat);

	//获取格式类型
	char** papszMetadata = poDriver->GetMetadata();

	//根据文件路径文件名，图像宽，高，波段数，数据类型，文件类型，创建新的数据集
	poDataset1 = poDriver->Create(path2, nImgSizeX, nImgSizeY, nBandCount, GDT_Byte, papszMetadata);

	//坐标赋值
	poDataset1->SetGeoTransform(geo);

	//将原图像数据读出，进行相应处理后，写入新文件
	BYTE* ppafScan = new BYTE[nImgSizeX * nImgSizeY * nBandCount];

	poDataset2->RasterIO(GF_Read, 0, 0, nImgSizeX, nImgSizeY, ppafScan, nImgSizeX, nImgSizeY, GDT_Byte, nBandCount, 0, 0, 0, 0);

	//将图像数据写入新图像中
	poDataset1->RasterIO(GF_Write, 0, 0, nImgSizeX, nImgSizeY, ppafScan, nImgSizeX, nImgSizeY, GDT_Byte, nBandCount, 0, 0, 0, 0);

	delete poDataset1;
	delete poDataset2;
}


pngTotif::~pngTotif() {}

3、main.cpp

#include<iostream>
#include "pngTotif.h"

using namespace std;


int main()
{
	//对象实例化
	pngTotif getTif;

	string srcTifpath = "..\\test_data\\02_04.tif";
	string srcPngpath = "..\\test_data\\1.png";
	string dstTifpath = ".\\image_save\\2020_02_04.tif";

	getTif.pngtotif(srcTifpath, srcPngpath, dstTifpath);
	
	return 0;
}

二、编译环境 

vs2019+gdal2.3.2+opencv4.4.0