BaseLibraryCPP/ImageOperatorBase.h

#pragma once
/**
 * 影像数据操作项
 * 所有数据相关的操作，都是基于ENVI的数据格式
 * 影像读写操作时基于 GDAL 库
 * **/



#ifndef IMAGEOPERATORBASE_H
#define IMAGEOPERATORBASE_H
#include "BaseConstVariable.h"
#include "GeoOperator.h"
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <complex>
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <string>
#include <memory>
#include <QString>
#include <cpl_conv.h> // for CPLMalloc()
#include "LogInfoCls.h"
#include <QObject>

enum ProjectStripDelta {
	Strip_6, // 6度带
	Strip_3
};

enum CoordinateSystemType { // 坐标系类型
	GeoCoordinateSystem,
	ProjectCoordinateSystem,
	UNKNOW
};

struct PointRaster { // 影像坐标点
	double x;
	double y;
	double z;
};


struct PointXYZ {
	double x, y, z;
};

struct  PointGeo {
	double lon, lat, ati;
};

struct PointImage {
	double pixel_x, pixel_y;
};



struct ImageGEOINFO {
	int width;
	int height;
	int bandnum;
};


enum GDALREADARRCOPYMETHOD {
	MEMCPYMETHOD,		// 直接拷贝
	VARIABLEMETHOD		// 变量赋值

};
 


class ShowProessAbstract{

public:
	virtual void showProcess(double precent,QString tip);
	virtual void showToolInfo(  QString tip)  ;
};




/// 根据经纬度获取
/// EPSG代码，根据经纬度返回对应投影坐标系统，其中如果在中华人民共和国境内，默认使用
/// CGCS2000坐标系统(EPSG 4502 ~ 4512 6度带,EPSG 4534 ~ 4554 3度带)，其余地方使用WGS坐标系统，
/// 投影方法 高斯克吕格（国内）， 高斯克吕格
/// \param long 经度
/// \param lat 纬度
/// \return 对应投影坐标系统的 EPSG编码,-1 表示计算错误
long getProjectEPSGCodeByLon_Lat(double long, double lat,ProjectStripDelta stripState );

long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat);

long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat);


QString GetProjectionNameFromEPSG(long epsgCode);


long GetEPSGFromRasterFile(QString filepath);

std::shared_ptr<PointRaster> GetCenterPointInRaster(QString filepath);

CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long EPSGCODE);
 
 
// 文件打开 // 当指令销毁时，调用GDALClose 销毁类型
std::shared_ptr<GDALDataset>   OpenDataset(const QString& in_path, GDALAccess rwmode= GA_ReadOnly); 
void   CloseDataset(GDALDataset* ptr);

// 数据格式转换
int   TIFF2ENVI(QString in_tiff_path,QString out_envi_path);
int   ENVI2TIFF(QString in_envi_path,QString out_tiff_path);

// 保存影像数据 --直接保存 ENVI 文件

int   CreateDataset(QString new_file_path, int height, int width, int band_num, double* gt, QString projection, GDALDataType gdal_dtype, bool need_gt); // 创建文件

int   saveDataset(QString new_file_path, int start_line, int start_cols, int band_ids, int datacols, int datarows, void* databuffer);

// 根据限制条件估算分块大小
int   block_num_pre_memory(int width, int height, GDALDataType gdal_dtype,double memey_size);

// 将结果转换为复数 或者 实数
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>   ReadComplexMatrixData(int start_line,int width, int line_num, std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype);

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>    ReadMatrixDoubleData(int start_line, int width, int line_num, std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype,int band_idx);

Eigen::MatrixXd   getGeoTranslationArray(QString in_path);
ImageGEOINFO   getImageINFO(QString in_path);

GDALDataType   getGDALDataType(QString fileptah);


struct RasterExtend {
	double min_x; //纬度
	double min_y;//经度
	double max_x;//纬度
	double max_y;//经度
};

 



/// <summary>
/// gdalImage图像操作类
/// </summary>

class   gdalImage
{

public: // 方法
	gdalImage();
	gdalImage(const QString& raster_path);
	~gdalImage();
	virtual void setHeight(int);
	virtual void setWidth(int);
	virtual void setTranslationMatrix(Eigen::MatrixXd gt);
	virtual void setData(Eigen::MatrixXd,int data_band_ids=1);
	virtual Eigen::MatrixXd getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
	virtual Eigen::MatrixXi getDatai(int start_row, int start_col, int rows_count, int cols_count, int band_ids);

	virtual Eigen::MatrixXd getGeoTranslation();
	virtual GDALDataType getDataType();
	virtual void saveImage(Eigen::MatrixXd, int start_row, int start_col, int band_ids);
	virtual void saveImage(Eigen::MatrixXi, int start_row, int start_col, int band_ids);
	virtual void saveImage();
	virtual void setNoDataValue(double nodatavalue, int band_ids);
	virtual void setNoDataValuei(int nodatavalue, int band_ids);
	virtual double getNoDataValue(int band_ids);
	virtual int getNoDataValuei(int band_ids);
	virtual int InitInv_gt();
	virtual Landpoint getRow_Col(double lon, double lat);
	virtual Landpoint getLandPoint(double i, double j, double ati);

	virtual void getLandPoint(double  i, double  j, double  ati, Landpoint& Lp);

	virtual double mean(int bandids = 1);
	double BandmaxValue(int bandids = 1);
	double BandminValue(int bandids = 1);
	virtual GDALRPCInfo getRPC();
	virtual Eigen::MatrixXd getLandPoint(Eigen::MatrixXd points);
 
	virtual Eigen::MatrixXd getHist(int bandids);

	virtual RasterExtend getExtend();

public:
	QString img_path; // 图像文件
	int height; // 高
	int width; // 宽
	int band_num;// 波段数
	int start_row;// 
	int start_col;//
	int data_band_ids;
	Eigen::MatrixXd gt; // 变换矩阵
	Eigen::MatrixXd inv_gt; // 逆变换矩阵
	Eigen::MatrixXd data;
	QString projection;
};



/// <summary>
/// gdalImage图像操作类
/// </summary>
class   gdalImageComplex:public gdalImage
{

public: // 方法
	gdalImageComplex(const QString& raster_path);
	~gdalImageComplex();
	void setData(Eigen::MatrixXcd);
	void saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids);
	Eigen::MatrixXcd getDataComplex(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
 
	void saveImage() override;
	void savePreViewImage();
public:
	Eigen::MatrixXcd data;
};

// 创建影像
gdalImage   CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false, bool isEnvi = false);
gdalImage   CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long espgcode, GDALDataType eType=GDT_Float32, bool need_gt = true, bool overwrite = false,bool isENVI=false);

gdalImageComplex   CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false);

gdalImageComplex   CreateEchoComplex(const QString& img_path, int height, int width, int band_num);

ErrorCode  DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QString demsloperPath);

int   ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample);

void resampleRaster(const char* inputRaster, const char* outputRaster, double targetPixelSizeX, double targetPixelSizeY);

void cropRasterByLatLon(const char* inputFile, const char* outputFile, double minLon, double maxLon, double minLat, double maxLat);

int   ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample = GRIORA_Bilinear);

void transformRaster(const char* inputFile, const char* outputFile, int sourceEPSG, int targetEPSG);

ErrorCode transformCoordinate(double x, double y, int sourceEPSG, int targetEPSG, Point2& p);

int alignRaster(QString  inputPath, QString referencePath, QString outputPath, GDALResampleAlg eResample);

//--------------------- 保存文博 -------------------------------

int   saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path);
 
//----------------------------------------------------

//--------------------- 图像合并流程    ------------------------------
enum MERGEMODE
{
	MERGE_GEOCODING,
};


ErrorCode MergeRasterProcess(QVector<QString> filepath, QString outfileptah, QString mainString, MERGEMODE  mergecode = MERGEMODE::MERGE_GEOCODING, bool isENVI = false, ShowProessAbstract* dia=nullptr);


ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> inimgs, gdalImage resultimg,gdalImage maskimg, ShowProessAbstract* dia = nullptr);
 

// 保存矩阵转换为envi文件，默认数据格式为double
bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath);



template<typename T>
std::shared_ptr<T> readDataArr(gdalImage& imgds, int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method)
{
	std::shared_ptr<T> result = nullptr;

	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(imgds.img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷

	GDALDataType	gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
	GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids);

	rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row;
	cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col;

	//Eigen::MatrixXd datamatrix(rows_count, cols_count);

	if (gdal_datatype == GDT_Byte) {
		char* temp = new char[rows_count * cols_count];
		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
			std::memcpy(result.get(), temp, rows_count * cols_count);
		}
		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
			long count = rows_count * cols_count;
			for (long i = 0; i < count; i++) {
				result.get()[i] = T(temp[i]);
			}
		}

		delete[] temp;
	}
	else if (gdal_datatype == GDT_UInt16) {
		unsigned short* temp = new unsigned short[rows_count * cols_count];
		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
			std::memcpy(result.get(), temp, rows_count * cols_count);
		}
		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
			long count = rows_count * cols_count;
			for (long i = 0; i < count; i++) {
				result.get()[i] = T(temp[i]);
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == GDT_Int16) {
		short* temp = new short[rows_count * cols_count];
		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
			std::memcpy(result.get(), temp, rows_count * cols_count);
		}
		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
			long count = rows_count * cols_count;
			for (long i = 0; i < count; i++) {
				result.get()[i] = T(temp[i]);
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == GDT_UInt32) {
		unsigned int* temp = new unsigned int[rows_count * cols_count];
		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
			std::memcpy(result.get(), temp, rows_count * cols_count);
		}
		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
			long count = rows_count * cols_count;
			for (long i = 0; i < count; i++) {
				result.get()[i] = T(temp[i]);
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == GDT_Int32) {
		int* temp = new int[rows_count * cols_count];
		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
			std::memcpy(result.get(), temp, rows_count * cols_count);
		}
		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
			long count = rows_count * cols_count;
			for (long i = 0; i < count; i++) {
				result.get()[i] = T(temp[i]);
			}
		}
		delete[] temp;
	}
	// else if (gdal_datatype == GDT_UInt64) {
	//	unsigned long* temp = new unsigned long[rows_count * cols_count];
	//	demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
	//rows_count, gdal_datatype, 0, 0); 	for (int i = 0; i < rows_count; i++) { 		for (int j = 0; j <
	//cols_count; j++) { 			datamatrix(i, j) = temp[i * cols_count + j];
	//		}
	//	}
	//	delete[] temp;
	// }
	// else if (gdal_datatype == GDT_Int64) {
	//	long* temp = new long[rows_count * cols_count];
	//	demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
	//rows_count, gdal_datatype, 0, 0); 	for (int i = 0; i < rows_count; i++) { 		for (int j = 0; j <
	//cols_count; j++) { 			datamatrix(i, j) = temp[i * cols_count + j];
	//		}
	//	}
	//	delete[] temp;
	// }
	else if (gdal_datatype == GDT_Float32) {
		float* temp = new float[rows_count * cols_count];
		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);

		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
			std::memcpy(result.get(), temp, rows_count * cols_count);
		}
		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
			long count = rows_count * cols_count;
			for (long i = 0; i < count; i++) {
				result.get()[i] = T(temp[i]);
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == GDT_Float64) {
		double* temp = new double[rows_count * cols_count];
		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);

		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
			std::memcpy(result.get(), temp, rows_count * cols_count);
		}
		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
			long count = rows_count * cols_count;
			for (long i = 0; i < count; i++) {
				result.get()[i] = T(temp[i]);
			}
		}
		delete[] temp;
	}
	//else if (gdal_datatype == GDT_CFloat32) {
	//	if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
	//		float* temp = new float[rows_count * cols_count * 2];
	//		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);

	//		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
	//		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
	//			std::memcpy(result.get(), temp, rows_count * cols_count);
	//		}
	//		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
	//			long count = rows_count * cols_count;
	//			for (long i = 0; i < count; i++) {
	//				result.get()[i] = T(temp[i * 2],
	//					temp[i * 2 + 1]);
	//			}
	//		}
	//		delete[] temp;
	//	}
	//	else {
	//		result = nullptr;
	//	}
	//}
	//else if (gdal_datatype == GDT_CFloat64 ) {
	//	if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
	//		double* temp = new double[rows_count * cols_count * 2];
	//		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);

	//		result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
	//		if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
	//			std::memcpy(result.get(), temp, rows_count * cols_count);
	//		}
	//		else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
	//			long count = rows_count * cols_count;
	//			for (long i = 0; i < count; i++) {
	//				result.get()[i] = T(temp[i * 2],
	//					temp[i * 2 + 1]);
	//			}
	//		}
	//		delete[] temp;
	//	}
	//	else {
	//		result = nullptr;
	//	}
	//}
	else {
	}
	GDALClose((GDALDatasetH)rasterDataset);
	omp_unset_lock(&lock);	 // 锟酵放伙拷斤拷
	omp_destroy_lock(&lock); // 劫伙拷斤拷
	// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	return result;
}

template<typename T>
std::shared_ptr<T> readDataArrComplex(gdalImageComplex& imgds, int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method)
{
	std::shared_ptr<T> result = nullptr;

	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(imgds.img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷

	GDALDataType	gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
	GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids);

	rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row;
	cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col;

	//Eigen::MatrixXd datamatrix(rows_count, cols_count);

	if (gdal_datatype == GDT_CFloat32) {
		if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
			float* temp = new float[rows_count * cols_count * 2];
			demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);

			result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
			if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
				std::memcpy(result.get(), temp, rows_count * cols_count);
			}
			else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
				long count = rows_count * cols_count;
				for (long i = 0; i < count; i++) {
					result.get()[i] = T(temp[i * 2],
						temp[i * 2 + 1]);
				}
			}
			delete[] temp;
		}
		else {
			result = nullptr;
		}
	}
	else if (gdal_datatype == GDT_CFloat64) {
		if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
			double* temp = new double[rows_count * cols_count * 2];
			demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);

			result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
			if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
				std::memcpy(result.get(), temp, rows_count * cols_count);
			}
			else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
				long count = rows_count * cols_count;
				for (long i = 0; i < count; i++) {
					result.get()[i] = T(temp[i * 2],
						temp[i * 2 + 1]);
				}
			}
			delete[] temp;
		}
		else {
			result = nullptr;
		}
	}
	else {
	}
	GDALClose((GDALDatasetH)rasterDataset);
	omp_unset_lock(&lock);	 // 锟酵放伙拷斤拷
	omp_destroy_lock(&lock); // 劫伙拷斤拷
	// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	return result;
}


//---------------------  图像分块   ------------------------------







#endif