diff --git a/ALLRelease/ALLRelease.vcxproj b/ALLRelease/ALLRelease.vcxproj
index 0552d3d..fd90192 100644
--- a/ALLRelease/ALLRelease.vcxproj
+++ b/ALLRelease/ALLRelease.vcxproj
@@ -196,9 +196,6 @@
{b8b40c54-f7fe-4809-b6fb-8bc014570d7b}
-
- {8c8ca066-a93a-4098-9a46-b855efeaadf2}
-
{4e6e79a3-048c-4fb4-bbb0-43c518a3e6d4}
@@ -211,9 +208,6 @@
{070c157e-3c30-4e2b-a80c-cbc7b74df03f}
-
- {d603a623-132d-4304-ab03-638fc438f084}
-
{ed06dfcd-4b9f-41f7-8f25-1823c2398142}
diff --git a/BaseCommonLibrary/BaseCommonLibrary.vcxproj b/BaseCommonLibrary/BaseCommonLibrary.vcxproj
index e6163eb..002db1e 100644
--- a/BaseCommonLibrary/BaseCommonLibrary.vcxproj
+++ b/BaseCommonLibrary/BaseCommonLibrary.vcxproj
@@ -220,10 +220,13 @@
pch.h
true
true
- StreamingSIMDExtensions2
+ NoExtensions
true
stdcpp14
stdc11
+ true
+ false
+ MaxSpeed
Windows
@@ -274,11 +277,12 @@
-
+
+
+
-
@@ -287,15 +291,20 @@
+
+
-
+
+
+
+
Create
Create
diff --git a/BaseCommonLibrary/BaseCommonLibrary.vcxproj.filters b/BaseCommonLibrary/BaseCommonLibrary.vcxproj.filters
index 81cb853..74f915b 100644
--- a/BaseCommonLibrary/BaseCommonLibrary.vcxproj.filters
+++ b/BaseCommonLibrary/BaseCommonLibrary.vcxproj.filters
@@ -33,33 +33,36 @@
BaseTool
-
- BaseTool
-
BaseTool
BaseTool
-
- BaseTool
-
BaseTool
BaseTool
-
- BaseTool
-
BaseTool
BaseTool
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
@@ -71,30 +74,18 @@
BaseTool
-
- BaseTool
-
BaseTool
BaseTool
-
- BaseTool
-
BaseTool
-
- BaseTool
-
BaseTool
-
- BaseTool
-
BaseTool
@@ -104,16 +95,45 @@
BaseTool
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
+
+ BaseTool
+
-
- 头文件
-
ToolAbstract
+
+ BaseTool
+
-
+
+ BaseTool
+
\ No newline at end of file
diff --git a/BaseCommonLibrary/BaseTool/BaseConstVariable.h b/BaseCommonLibrary/BaseTool/BaseConstVariable.h
index 979b12f..67114db 100644
--- a/BaseCommonLibrary/BaseTool/BaseConstVariable.h
+++ b/BaseCommonLibrary/BaseTool/BaseConstVariable.h
@@ -15,12 +15,15 @@
+/** 定义常见文件格式*********/
+#define ENVI_FILE_FORMAT_FILTER u8"ALL File(*.*);;ENVI Bin(*.bin);;ENVI Data(*.dat);;ENVI Data(*.data);;tiff影像(*.tif);;tiff影像(*.tiff)"
+#define XML_FILE_FORMAT_FILTER u8"ALL File(*.*);;XML File(*.xml);;tiff影像(*.tiff)"
-
-#define EIGEN_USE_BLAS
-#define EIGEN_USE_LAPACK
-#define EIGEN_VECTORIZE_SSE2
+//
+//#define EIGEN_USE_BLAS
+//#define EIGEN_USE_LAPACK
+//#define EIGEN_VECTORIZE_SSE2
//#define DEBUGSHOWDIALOG
@@ -49,7 +52,7 @@ inline char* get_cur_time() {
return s;
}
-#define PRINT(fmt, ...) printf("%s " fmt, get_cur_time(), ##__VA_ARGS__)
+#define PRINT(fmt, ...) printf("%s " fmt, get_cur_time(), ##__VA_ARGS__);
@@ -136,6 +139,14 @@ struct Point3 {
void setZ(double iz) { z = iz; }
};
+struct Point_3d {
+ double x;
+ double y;
+ double z;
+};
+
+
+
struct DemBox {
double min_lon;
double max_lon;
@@ -154,7 +165,7 @@ struct Vector3_single {
/*********************************************** FEKO仿真参数 ********************************************************************/
-struct SatellitePos {
+extern "C" struct SatellitePos {
double time;
double Px ;
double Py ;
@@ -165,7 +176,7 @@ struct SatellitePos {
};
-struct SatelliteAntPos {
+extern "C" struct SatelliteAntPos {
double time; // 0
double Px;
double Py;
@@ -188,7 +199,7 @@ struct SatelliteAntPos {
};
-struct PatternImageDesc {
+extern "C" struct PatternImageDesc {
long phinum;
long thetanum;
double startTheta;
@@ -307,7 +318,7 @@ inline void delArrPtr(void* p)
}
delete[] p;
p = nullptr;
-}
+};
inline void delPointer(void* p)
{
@@ -316,28 +327,34 @@ inline void delPointer(void* p)
}
delete p;
p = nullptr;
-}
+};
inline void PrintTime() {
time_t current_time;
time(¤t_time);
printf("Current timestamp in seconds: %ld\n", (long)current_time);
-}
+};
/** 计算分块 ******************************************************************/
-inline long getBlockRows(long sizeMB, long cols,long sizeMeta,long maxRows) {
- long rownum= (round(Memory1MB * 1.0 / sizeMeta / cols * sizeMB) + cols - 1);
+inline long getBlockRows(long sizeMB, long cols, long sizeMeta, long maxRows) {
+ long rownum = (round(Memory1MB * 1.0 / sizeMeta / cols * sizeMB) + cols - 1);
rownum = rownum < 0 ? 1 : rownum;
- rownum =rownum < maxRows ? rownum : maxRows;
+ rownum = rownum < maxRows ? rownum : maxRows;
return rownum;
-}
+};
inline long nextpow2(long n) {
long en = ceil(log2(n));
- return pow(2,en);
-}
+ return pow(2, en);
+};
+
+inline void releaseVoidArray(void* a)
+{
+ free(a);
+ a = NULL;
+};
diff --git a/BaseCommonLibrary/BaseTool/BaseTool.cpp b/BaseCommonLibrary/BaseTool/BaseTool.cpp
index 5cd20e5..5ea4ed0 100644
--- a/BaseCommonLibrary/BaseTool/BaseTool.cpp
+++ b/BaseCommonLibrary/BaseTool/BaseTool.cpp
@@ -1,5 +1,5 @@
#include "stdafx.h"
-#define EIGEN_USE_BLAS
+//#define EIGEN_USE_BLAS
#define EIGEN_VECTORIZE_SSE4_2
//#include
@@ -38,7 +38,10 @@
#include // 包含SSE指令集头文件
#include // 包含SSE2指令集头文件
#include // 包含OpenMP头文件
-
+#include
+#include
+#include
+#include
QString longDoubleToQStringScientific(long double value) {
std::ostringstream stream;
@@ -679,4 +682,94 @@ void initializeMatrixWithSSE2(Eigen::MatrixXf& mat, const float* data, long rowc
}
}
+Eigen::MatrixXd BASECONSTVARIABLEAPI MuhlemanSigmaArray(Eigen::MatrixXd& eta_deg)
+{
+ Eigen::MatrixXd sigma = Eigen::MatrixXd::Zero(eta_deg.rows(), eta_deg.cols());
+ for (long i = 0; i < sigma.rows(); i++) {
+ for (long j = 0; j < sigma.cols(); j++) {
+ sigma(i,j) = calculate_MuhlemanSigma(eta_deg(i, j));
+ }
+ }
+ return sigma;
+}
+
+Eigen::MatrixXd BASECONSTVARIABLEAPI dB2Amp(Eigen::MatrixXd& sigma0)
+{
+ Eigen::MatrixXd sigma = Eigen::MatrixXd::Zero(sigma0.rows(), sigma0.cols());
+ for (long i = 0; i < sigma.rows(); i++) {
+ for (long j = 0; j < sigma.cols(); j++) {
+ sigma(i, j) =std::pow(10.0, sigma0(i,j)/20.0);
+ }
+ }
+ return sigma;
+}
+
+
+
+QDateTime parseCustomDateTime(const QString& dateTimeStr) {
+ // 手动解析日期时间字符串
+ int year = dateTimeStr.left(4).toInt();
+ int month = dateTimeStr.mid(5, 2).toInt();
+ int day = dateTimeStr.mid(8, 2).toInt();
+ int hour = dateTimeStr.mid(11, 2).toInt();
+ int minute = dateTimeStr.mid(14, 2).toInt();
+ int second = dateTimeStr.mid(17, 2).toInt();
+ int msec = dateTimeStr.mid(20, 6).toInt(); // 只取毫秒的前三位,因为QDateTime支持到毫秒
+
+ // 创建 QDate 和 QTime 对象
+ QDate date(year, month, day);
+ QTime time(hour, minute, second, msec ); // 转换为微秒,但QTime只支持毫秒精度
+
+ // 构造 QDateTime 对象
+ QDateTime result(date, time);
+
+ return result;
+}
+
+
+bool isLeapYear(int year) {
+ return (year % 4 == 0 && (year % 100 != 0 || year % 400 == 0));
+}
+
+int daysInMonth(int year, int month) {
+ if (month == 2) return isLeapYear(year) ? 29 : 28;
+ else if (month == 4 || month == 6 || month == 9 || month == 11) return 30;
+ else return 31;
+}
+
+
+TimestampMicroseconds parseAndConvert( std::string dateTimeStr) {
+ // 解析年、月、日、时、分、秒和微秒
+ int year = std::stoi(dateTimeStr.substr(0, 4));
+ int month = std::stoi(dateTimeStr.substr(5, 2));
+ int day = std::stoi(dateTimeStr.substr(8, 2));
+ int hour = std::stoi(dateTimeStr.substr(11, 2));
+ int minute = std::stoi(dateTimeStr.substr(14, 2));
+ int second = std::stoi(dateTimeStr.substr(17, 2));
+ int microsec = std::stoi(dateTimeStr.substr(20, 6));
+
+ // 计算从1970年至目标年份前一年的总天数
+ long long totalDays = 0;
+ for (int y = 1970; y < year; ++y) {
+ totalDays += isLeapYear(y) ? 366 : 365;
+ }
+
+ // 加上目标年份从1月到上一个月的天数
+ for (int m = 1; m < month; ++m) {
+ totalDays += daysInMonth(year, m);
+ }
+
+ // 加上本月的天数
+ totalDays += day - 1;
+
+ // 转换为总秒数,再加上小时、分钟、秒
+ long long totalSeconds = totalDays * 24 * 60 * 60 + hour * 60 * 60 + minute * 60 + second;
+
+ // 转换为毫秒和微秒
+ long long msecsSinceEpoch = totalSeconds * 1000 + microsec / 1000;
+ int microseconds = microsec % 1000;
+
+ return { msecsSinceEpoch, microseconds };
+}
+
diff --git a/BaseCommonLibrary/BaseTool/BaseTool.h b/BaseCommonLibrary/BaseTool/BaseTool.h
index bbed5a7..6ad57ff 100644
--- a/BaseCommonLibrary/BaseTool/BaseTool.h
+++ b/BaseCommonLibrary/BaseTool/BaseTool.h
@@ -31,7 +31,8 @@
#include
#include
#include
-
+#include
+#include
///////////////////////////////////// 基础数学函数 /////////////////////////////////////////////////////////////
@@ -101,7 +102,7 @@ Eigen::Matrix3d BASECONSTVARIABLEAPI rotationMatrix(const Eigen::Vector3d& axis
long double BASECONSTVARIABLEAPI convertToMilliseconds(const std::string& dateTimeStr);
-
+QDateTime BASECONSTVARIABLEAPI parseCustomDateTime(const QString& dateTimeStr);
///
/// list 应该是按照从小到大的顺序排好
///
@@ -129,28 +130,55 @@ Eigen::VectorXd BASECONSTVARIABLEAPI linspace(double start, double stop, int nu
void initializeMatrixWithSSE2(Eigen::MatrixXd& mat, const double* data, long rowcount, long colcount);
void initializeMatrixWithSSE2(Eigen::MatrixXf& mat, const float* data, long rowcount, long colcount);
+Eigen::MatrixXd BASECONSTVARIABLEAPI MuhlemanSigmaArray(Eigen::MatrixXd& eta_deg);
+Eigen::MatrixXd BASECONSTVARIABLEAPI dB2Amp(Eigen::MatrixXd& sigma0);
+struct TimestampMicroseconds {
+ boost::int64_t msecsSinceEpoch; // 自1970-01-01T00:00:00 UTC以来的毫秒数
+ int microseconds; // 额外的微秒(精确到微秒)
+};
+
+
+bool BASECONSTVARIABLEAPI isLeapYear(int year);
+int BASECONSTVARIABLEAPI daysInMonth(int year, int month);
+
+TimestampMicroseconds BASECONSTVARIABLEAPI parseAndConvert( std::string dateTimeStr);
+
/** 模板函数类 ***********************************************************************************************************/
+
+inline double calculate_MuhlemanSigma(double eta_deg) {
+ const double eta_rad = eta_deg * M_PI / 180.0; // 角度转弧度
+
+ const double cos_eta = std::cos(eta_rad);
+ const double sin_eta = std::sin(eta_rad);
+
+ const double denominator = sin_eta + 0.1 * cos_eta;
+
+ return (0.0133 * cos_eta) / std::pow(denominator, 3);
+};
+
+
+
template
-inline void memsetInitArray(std::shared_ptr ptr, long arrcount,T ti) {
+inline void memsetInitArray(std::shared_ptr ptr, long arrcount, T ti) {
for (long i = 0; i < arrcount; i++) {
ptr.get()[i] = ti;
}
-}
+};
template
inline void memcpyArray(std::shared_ptr srct, std::shared_ptr dest, long arrcount) {
for (long i = 0; i < arrcount; i++) {
dest.get()[i] = srct.get()[i];
}
-}
+};
template
-inline void minValueInArr(T* ptr, long arrcount, T& minvalue) {
-
+inline void minValueInArr(T* ptr, long arrcount, T& minvalue) {
+
if (arrcount == 0)return;
minvalue = ptr[0];
@@ -159,7 +187,7 @@ inline void minValueInArr(T* ptr, long arrcount, T& minvalue) {
minvalue = ptr[i];
}
}
-}
+};
template
inline void maxValueInArr(T* ptr, long arrcount, T& maxvalue) {
@@ -172,7 +200,7 @@ inline void maxValueInArr(T* ptr, long arrcount, T& maxvalue) {
maxvalue = ptr[i];
}
}
-}
+};
@@ -182,7 +210,7 @@ inline void maxValueInArr(T* ptr, long arrcount, T& maxvalue) {
template
inline T complexAbs(std::complex ccdata) {
return T(sqrt(pow(ccdata.real(), 2) + pow(ccdata.imag(), 2)));
-}
+};
template
inline void complex2dB(std::complex* ccdata, T* outdata, long long count) {
@@ -191,7 +219,7 @@ inline void complex2dB(std::complex* ccdata, T* outdata, long long count) {
outdata[i] = 20 * log10(complexAbs(ccdata[i]));
}
-}
+};
diff --git a/BaseCommonLibrary/BaseTool/EchoDataFormat.cpp b/BaseCommonLibrary/BaseTool/EchoDataFormat.cpp
index e78b9d8..bc50ed5 100644
--- a/BaseCommonLibrary/BaseTool/EchoDataFormat.cpp
+++ b/BaseCommonLibrary/BaseTool/EchoDataFormat.cpp
@@ -234,6 +234,7 @@ ErrorCode EchoL0Dataset::Open(QString folder, QString filename)
{
return ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR;
}
+ return ErrorCode::SUCCESS;
}
QString EchoL0Dataset::getxmlName()
@@ -511,6 +512,19 @@ ErrorCode EchoL0Dataset::loadFromXml() {
return ErrorCode::SUCCESS;
}
+std::shared_ptr EchoL0Dataset::getAntPosVelc()
+{
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ long prfcount = this->PluseCount;
+ std::shared_ptr antposlist= SatelliteAntPosOperator::readAntPosFile(this->GPSPointFilePath, prfcount);
+
+ omp_unset_lock(&lock); //
+ omp_destroy_lock(&lock); //
+ return antposlist;
+}
+
std::shared_ptr EchoL0Dataset::getAntPos()
{
omp_lock_t lock;
@@ -592,6 +606,66 @@ std::shared_ptr> EchoL0Dataset::getEchoArr(long startPRF, l
return temp;
}
+
+std::vector> EchoL0Dataset::getEchoArrVector(long startPRF, long& PRFLen) {
+ if (!(startPRF < this->PluseCount)) {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE)) << startPRF << " " << this->PluseCount;
+ return std::vector>(0);
+ }
+ else {}
+
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ std::shared_ptr rasterDataset = OpenDataset(this->echoDataFilePath, GDALAccess::GA_ReadOnly);
+
+
+ GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
+ GDALRasterBand* poBand = rasterDataset->GetRasterBand(1);
+ if (NULL == poBand || nullptr == poBand) {
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
+ return std::vector>(0);
+ }
+ else {}
+ long width = rasterDataset->GetRasterXSize();
+ long height = rasterDataset->GetRasterYSize();
+ long band_num = rasterDataset->GetRasterCount();
+
+ PRFLen = (PRFLen + startPRF) < height ? PRFLen : height - startPRF;
+
+ std::vector> tempArr(size_t(PRFLen) * width);
+ if (height != this->PluseCount || width != this->PlusePoints) {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
+ }
+ else {
+ if (gdal_datatype == GDT_CFloat64) {
+ std::shared_ptr> temp(new std::complex[PRFLen * width], delArrPtr);
+ poBand->RasterIO(GF_Read, 0, startPRF, width, PRFLen, temp.get(), width, PRFLen, GDT_CFloat64, 0, 0);
+ GDALFlushCache((GDALDatasetH)rasterDataset.get());
+
+
+ for (long i = 0; i < PRFLen; i++){
+ for (long j = 0; j < width; j++){
+ tempArr[i * width + j] = temp.get()[i * width + j];
+ }
+ }
+
+ }
+ else {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
+ }
+ }
+
+ rasterDataset.reset();
+ omp_unset_lock(&lock); //
+ omp_destroy_lock(&lock); //
+ return tempArr;
+}
+
+
std::shared_ptr> EchoL0Dataset::getEchoArr()
{
return this->getEchoArr(0,this->PluseCount);
@@ -674,8 +748,10 @@ ErrorCode EchoL0Dataset::saveEchoArr(std::shared_ptr> echoP
}
GDALClose(rasterDataset);
rasterDataset = nullptr;
+ GDALDestroy();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
+
return ErrorCode::SUCCESS;
}
diff --git a/BaseCommonLibrary/BaseTool/EchoDataFormat.h b/BaseCommonLibrary/BaseTool/EchoDataFormat.h
index ca0e5ce..ec2f685 100644
--- a/BaseCommonLibrary/BaseTool/EchoDataFormat.h
+++ b/BaseCommonLibrary/BaseTool/EchoDataFormat.h
@@ -209,9 +209,11 @@ public: //
public:
// ȡļ
+ std::shared_ptr< SatelliteAntPos> getAntPosVelc();
std::shared_ptr getAntPos();
std::shared_ptr> getEchoArr(long startPRF, long& PRFLen);
std::shared_ptr> getEchoArr();
+ std::vector> getEchoArrVector(long startPRF, long& PRFLen);
//ļ
ErrorCode saveAntPos(std::shared_ptr ptr); // עΣգдǰǷȷ
ErrorCode saveEchoArr(std::shared_ptr> echoPtr, long startPRF, long PRFLen);
diff --git a/BaseCommonLibrary/BaseTool/FileOperator.cpp b/BaseCommonLibrary/BaseTool/FileOperator.cpp
index d2b1343..6d6ba6c 100644
--- a/BaseCommonLibrary/BaseTool/FileOperator.cpp
+++ b/BaseCommonLibrary/BaseTool/FileOperator.cpp
@@ -93,7 +93,8 @@ QString getFileNameFromPath(const QString& path)
QString getFileNameWidthoutExtend(QString path)
{
QFileInfo fileInfo(path);
- QString fileNameWithoutExtension = fileInfo.baseName(); // 获取无后缀文件名
+ QString fileNameWithoutExtension = fileInfo.completeBaseName(); // 获取无后缀文件名
+ qDebug() << u8"File name without extension: " << fileNameWithoutExtension;
return fileNameWithoutExtension;
}
diff --git a/BaseCommonLibrary/BaseTool/FileOperator.h b/BaseCommonLibrary/BaseTool/FileOperator.h
index 72d542a..b016ede 100644
--- a/BaseCommonLibrary/BaseTool/FileOperator.h
+++ b/BaseCommonLibrary/BaseTool/FileOperator.h
@@ -16,7 +16,7 @@
#include
#include
#include
-
+
bool BASECONSTVARIABLEAPI isDirectory(const QString& path);
bool BASECONSTVARIABLEAPI isExists(const QString& path);
@@ -29,7 +29,7 @@ unsigned long BASECONSTVARIABLEAPI convertToULong(const QString& input);
///
///
///
-std::vector BASECONSTVARIABLEAPI getFilelist(const QString& folderpath, const QString& FilenameExtension = ".*",int (*logfun)(QString logtext,int value)=nullptr);
+std::vector BASECONSTVARIABLEAPI getFilelist(const QString& folderpath, const QString& FilenameExtension = ".*", int (*logfun)(QString logtext, int value) = nullptr);
QString BASECONSTVARIABLEAPI getParantFolderNameFromPath(const QString& path);
@@ -41,11 +41,11 @@ QString BASECONSTVARIABLEAPI getFileExtension(QString path);
int BASECONSTVARIABLEAPI write_binfile(char* filepath, char* data, size_t data_len);
-char* read_textfile(char* text_path, int* length);
+char* read_textfile(char* text_path, int* length);
-bool BASECONSTVARIABLEAPI exists_test(const QString& name);
+bool BASECONSTVARIABLEAPI exists_test(const QString& name);
-size_t BASECONSTVARIABLEAPI fsize(FILE* fp);
+size_t BASECONSTVARIABLEAPI fsize(FILE* fp);
QString BASECONSTVARIABLEAPI getParantFromPath(const QString& path);
void BASECONSTVARIABLEAPI copyFile(const QString& sourcePath, const QString& destinationPath);
diff --git a/BaseCommonLibrary/BaseTool/GeoOperator.cpp b/BaseCommonLibrary/BaseTool/GeoOperator.cpp
index 9784c3a..ce49a3e 100644
--- a/BaseCommonLibrary/BaseTool/GeoOperator.cpp
+++ b/BaseCommonLibrary/BaseTool/GeoOperator.cpp
@@ -18,7 +18,9 @@
#include
#include
#include "GeoOperator.h"
-
+#include // OGRSpatialReference ڿռοת
+#include // GDALWarp
+#include
Landpoint operator +(const Landpoint& p1, const Landpoint& p2)
@@ -135,6 +137,44 @@ Landpoint XYZ2LLA(const Landpoint& XYZ) {
+void XYZ2BLH_FixedHeight(double x, double y, double z,double ati, Landpoint& point) {
+ const double a = 6378137.0; // WGS84
+ const double f = 1.0 / 298.257223563;
+ const double e2 = 2 * f - f * f; // һƫƽ
+
+ // 㾭L ()
+ const double L_rad = std::atan2(y, x);
+ point.lon = L_rad * 180.0 / M_PI; // תΪ
+
+ const double p = std::sqrt(x * x + y * y);
+ const double H = ati; // ʹ֪߶
+
+ // ʼγȹ㣨֪߶H
+ double B_rad = std::atan2(z, p * (1 - e2 * (a / (a + H))));
+
+ // γB̶H
+ for (int i = 0; i < 10; ++i) { // ֪Hʱ
+ const double sin_B = std::sin(B_rad);
+ const double N = a / std::sqrt(1 - e2 * sin_B * sin_B);
+ const double delta = e2 * N / (N + H); // ߶ȹ̶ʱ
+
+ const double B_new = std::atan2(z, p * (1 - delta));
+
+ if (std::abs(B_new - B_rad) < 1e-9) {
+ B_rad = B_new;
+ break;
+ }
+ B_rad = B_new;
+ }
+
+ point.lat = B_rad * 180.0 / M_PI; // ת
+
+ // ȷΧ [-180, 180]
+ point.lon = std::fmod(point.lon + 360.0, 360.0);
+ if (point.lon > 180.0) point.lon -= 360.0;
+ point.ati = ati;
+}
+
double getAngle(const Landpoint& a, const Landpoint& b)
{
@@ -408,3 +448,53 @@ double getlocalIncAngle(Vector3D& satepoint, Vector3D& landpoint, Vector3D& slop
return angle;
}
}
+
+bool BASECONSTVARIABLEAPI ConvertResolutionToDegrees(int sourceEPSG, double resolutionMeters, double refLon, double refLat, double& degreePerPixelX, double& degreePerPixelY){
+ // ʼԴϵƽͶӰĿϵWGS84 γȣ
+ OGRSpatialReference sourceSRS, targetSRS;
+ sourceSRS.importFromEPSG(sourceEPSG); // Դϵȷ EPSG
+ targetSRS.importFromEPSG(4326); // ĿΪ WGS84 γ
+
+ // תγ -> ƽ
+ OGRCoordinateTransformation* toPlane = OGRCreateCoordinateTransformation(&targetSRS, &sourceSRS);
+ if (!toPlane) return false;
+
+ // ο㾭γתΪƽ
+ double x = refLon, y = refLat;
+ if (!toPlane->Transform(1, &x, &y)) {
+ OGRCoordinateTransformation::DestroyCT(toPlane);
+ return false;
+ }
+
+ // תƽ -> γ
+ OGRCoordinateTransformation* toGeo = OGRCreateCoordinateTransformation(&sourceSRS, &targetSRS);
+ if (!toGeo) {
+ OGRCoordinateTransformation::DestroyCT(toPlane);
+ return false;
+ }
+
+ // 㶫ֱʣȱ仯
+ double eastX = x + resolutionMeters, eastY = y;
+ double eastLon = eastX, eastLat = eastY;
+ if (!toGeo->Transform(1, &eastLon, &eastLat)) {
+ OGRCoordinateTransformation::DestroyCT(toPlane);
+ OGRCoordinateTransformation::DestroyCT(toGeo);
+ return false;
+ }
+ degreePerPixelX = (eastLon - refLon) / resolutionMeters; // ȷÿӦ
+
+ // 㱱ֱʣγȱ仯
+ double northX = x, northY = y + resolutionMeters;
+ double northLon = northX, northLat = northY;
+ if (!toGeo->Transform(1, &northLon, &northLat)) {
+ OGRCoordinateTransformation::DestroyCT(toPlane);
+ OGRCoordinateTransformation::DestroyCT(toGeo);
+ return false;
+ }
+ degreePerPixelY = (northLat - refLat) / resolutionMeters; // γȷÿӦ
+
+ // ͷԴ
+ OGRCoordinateTransformation::DestroyCT(toPlane);
+ OGRCoordinateTransformation::DestroyCT(toGeo);
+ return true;
+}
diff --git a/BaseCommonLibrary/BaseTool/GeoOperator.h b/BaseCommonLibrary/BaseTool/GeoOperator.h
index 6d4eb8a..66ce525 100644
--- a/BaseCommonLibrary/BaseTool/GeoOperator.h
+++ b/BaseCommonLibrary/BaseTool/GeoOperator.h
@@ -30,6 +30,7 @@ Eigen::MatrixXd BASECONSTVARIABLEAPI LLA2XYZ(Eigen::MatrixXd landpoint);
/// 经纬度--degree
Landpoint BASECONSTVARIABLEAPI XYZ2LLA(const Landpoint& XYZ);
+void BASECONSTVARIABLEAPI XYZ2BLH_FixedHeight(double x, double y, double z, double ati, Landpoint& point);
Landpoint BASECONSTVARIABLEAPI operator +(const Landpoint& p1, const Landpoint& p2);
@@ -123,4 +124,24 @@ CartesianCoordinates BASECONSTVARIABLEAPI sphericalToCartesian(const Spheric
double BASECONSTVARIABLEAPI getlocalIncAngle(Vector3D& satepoint, Vector3D& landpoint, Vector3D& slopeVector);
+
+/**
+ * @brief 将平面坐标分辨率(米)转换为经纬度分辨率(度)
+ * @param sourceEPSG 平面坐标系 EPSG 码(如 UTM Zone 50N 对应 32650)
+ * @param resolutionMeters 输入分辨率(单位:米)
+ * @param refLon 参考点经度(十进制度,用于计算局部转换系数)
+ * @param refLat 参考点纬度(十进制度,用于计算局部转换系数)
+ * @param[out] degreePerPixelX 经度方向分辨率(度/像素)
+ * @param[out] degreePerPixelY 纬度方向分辨率(度/像素)
+ * @return 是否转换成功
+ */
+bool BASECONSTVARIABLEAPI ConvertResolutionToDegrees(
+ int sourceEPSG,
+ double resolutionMeters,
+ double refLon,
+ double refLat,
+ double& degreePerPixelX,
+ double& degreePerPixelY
+);
+
#endif
\ No newline at end of file
diff --git a/BaseCommonLibrary/BaseTool/ImageOperatorBase.cpp b/BaseCommonLibrary/BaseTool/ImageOperatorBase.cpp
index b3c419e..0edf2f0 100644
--- a/BaseCommonLibrary/BaseTool/ImageOperatorBase.cpp
+++ b/BaseCommonLibrary/BaseTool/ImageOperatorBase.cpp
@@ -431,1467 +431,6 @@ GDALDataType getGDALDataType(QString fileptah)
return gdal_datatype;
}
-gdalImage::gdalImage()
-{
- this->height = 0;
- this->width = 0;
- this->data_band_ids = 1;
- this->start_row = 0;
- this->start_col = 0;
-}
-
-///
-/// 斤拷图饺∮帮拷锟?1锟?7
-///
-///
-gdalImage::gdalImage(const QString& raster_path)
-{
- omp_lock_t lock;
- omp_init_lock(&lock); // 锟斤拷始斤拷斤拷
- omp_set_lock(&lock); // 锟斤拷没斤拷锟?1锟?7
- this->img_path = raster_path;
-
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- // 锟斤拷DEM影锟斤拷
- GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
- raster_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
- this->width = rasterDataset->GetRasterXSize();
- this->height = rasterDataset->GetRasterYSize();
- this->band_num = rasterDataset->GetRasterCount();
-
- double* gt = new double[6];
- // 锟斤拷梅斤拷锟斤拷
- rasterDataset->GetGeoTransform(gt);
- this->gt = Eigen::MatrixXd(2, 3);
- this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5];
-
- this->projection = rasterDataset->GetProjectionRef();
- // 斤拷斤拷
- // double* inv_gt = new double[6];;
- // GDALInvGeoTransform(gt, inv_gt); // 斤拷斤拷
- // 斤拷投影
- GDALFlushCache((GDALDatasetH)rasterDataset);
- GDALClose((GDALDatasetH)rasterDataset);
- rasterDataset = NULL; // 指矫匡拷
- this->InitInv_gt();
- delete[] gt;
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
-}
-
-gdalImage::~gdalImage() {}
-
-void gdalImage::setHeight(int height)
-{
- this->height = height;
-}
-
-void gdalImage::setWidth(int width)
-{
- this->width = width;
-}
-
-void gdalImage::setTranslationMatrix(Eigen::MatrixXd gt)
-{
- this->gt = gt;
-}
-
-void gdalImage::setData(Eigen::MatrixXd, int data_band_ids)
-{
- this->data = data;
- this->data_band_ids = data_band_ids;
-}
-
-Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, int cols_count,
- int band_ids = 1)
-{
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
- this->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 <= this->height ? rows_count : this->height - start_row;
- cols_count = start_col + cols_count <= this->width ? cols_count : this->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);
- 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_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);
- 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_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);
- 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_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);
- 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_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);
- 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_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);
-
- 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_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);
-
- 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 {
- }
- GDALClose((GDALDatasetH)rasterDataset);
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- return datamatrix;
-}
-
-Eigen::MatrixXf gdalImage::getDataf(int start_row, int start_col, int rows_count, int cols_count,
- int band_ids = 1)
-{
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
- this->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 <= this->height ? rows_count : this->height - start_row;
- cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
-
- Eigen::MatrixXf 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);
- 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_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);
- 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_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);
- 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_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);
- 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_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);
- 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_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);
-
- 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_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);
-
- 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 {
- }
- GDALClose((GDALDatasetH)rasterDataset);
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- return datamatrix;
-}
-
-
-
-Eigen::MatrixXi gdalImage::getDatai(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
-{
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
- this->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 <= this->height ? rows_count : this->height - start_row;
- cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
-
- Eigen::MatrixXi 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);
- 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_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);
- 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_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);
- 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_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);
- 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_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);
- 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_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);
-
- 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_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);
-
- 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 {
- }
- GDALClose((GDALDatasetH)rasterDataset);
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- return datamatrix;
-}
-
-ErrorCode gdalImage::getData(double* data, int start_row, int start_col, int rows_count, int cols_count, int band_ids)
-{
- ErrorCode state =ErrorCode::SUCCESS;
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->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 <= this->height ? rows_count : this->height - start_row;
- cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
-
-
-
- if (gdal_datatype == GDT_Float64) {
- demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, data, cols_count,rows_count, gdal_datatype, 0, 0);
- }
- else {
- state = ErrorCode::FAIL;
- }
- GDALClose((GDALDatasetH)rasterDataset);
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
- return state;
-}
-
-ErrorCode gdalImage::getData(long* data, int start_row, int start_col, int rows_count, int cols_count, int band_ids)
-{
- ErrorCode state = ErrorCode::SUCCESS;
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->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 <= this->height ? rows_count : this->height - start_row;
- cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
-
- if (gdal_datatype == GDT_Int32) {
- demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, data, cols_count,rows_count, gdal_datatype, 0, 0);
- }
- else {
- state = ErrorCode::FAIL;
- }
- GDALClose((GDALDatasetH)rasterDataset);
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
- return state;
-}
-
-Eigen::MatrixXd gdalImage::getGeoTranslation()
-{
- return this->gt;
-}
-
-GDALDataType gdalImage::getDataType()
-{
- GDALDataset* rasterDataset =
- (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));
- GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
- return gdal_datatype;
-}
-
-///
-///
-///
-///
-///
-///
-///
-void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0,
- int band_ids = 1)
-{
- GDALDataType datetype = this->getDataType();
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if(start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
- QString tip = u8"file path: " + this->img_path+" image size :( "+QString::number(this->height)+" , "+ QString::number(this->width)+" ) "+" input size ("+ QString::number(start_row + data.rows())+", "+ QString::number(start_col + data.cols())+") ";
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
-
- QString filesuffer = getFileExtension(this->img_path).toLower();
- bool isTiff = filesuffer.contains("tif");
- GDALDriver* poDriver = isTiff? GetGDALDriverManager()->GetDriverByName("GTiff"): GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if(exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- } else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, datetype, NULL); // 斤拷锟斤拷
-
- if (nullptr == poDstDS) {
- QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- return;
- }
-
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for(int i = 0; i < this->gt.rows(); i++) {
- for(int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- //delete gt_ptr;
- }
-
-
- int datarows = data.rows();
- int datacols = data.cols();
- void* databuffer = nullptr;
- if (datetype == GDT_Float32) {
- databuffer = new float[datarows * datacols];
- for (int i = 0; i < datarows; i++) {
- for (int j = 0; j < datacols; j++) {
- ((float*)databuffer)[i * datacols + j] = float(data(i, j));
- }
- }
-
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, datetype, 0, 0);
- }
- else if (datetype == GDT_Float64) {
- databuffer = new double[datarows * datacols];
- for (int i = 0; i < datarows; i++) {
- for (int j = 0; j < datacols; j++) {
- ((double*)databuffer)[i * datacols + j] = double(data(i, j));
- }
- }
-
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, datetype, 0, 0);
- }
- else {
-
- }
- GDALFlushCache(poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- delete[] databuffer;
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
-}
-
-void gdalImage::saveImage(Eigen::MatrixXf data, int start_row = 0, int start_col = 0,
- int band_ids = 1)
-{
- GDALDataType datetype = this->getDataType();
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
- QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- QString filesuffer = getFileExtension(this->img_path).toLower();
- bool isTiff = filesuffer.contains("tif");
- GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if (exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- }
- else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, datetype, NULL); // 斤拷锟斤拷
-
- if (nullptr == poDstDS) {
- QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- return;
- }
-
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for (int i = 0; i < this->gt.rows(); i++) {
- for (int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- //delete gt_ptr;
- }
-
-
- int datarows = data.rows();
- int datacols = data.cols();
- void* databuffer = nullptr;
- if (datetype == GDT_Float32) {
- databuffer = new float[datarows * datacols];
- for (int i = 0; i < datarows; i++) {
- for (int j = 0; j < datacols; j++) {
- ((float*)databuffer)[i * datacols + j] = float(data(i, j));
- }
- }
-
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, datetype, 0, 0);
- }
- else if (datetype == GDT_Float64) {
- databuffer = new double[datarows * datacols];
- for (int i = 0; i < datarows; i++) {
- for (int j = 0; j < datacols; j++) {
- ((double*)databuffer)[i * datacols + j] = double(data(i, j));
- }
- }
-
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, datetype, 0, 0);
- }
- else {
-
- }
- GDALFlushCache(poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- delete[] databuffer;
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
-}
-
-
-void gdalImage::saveImage(Eigen::MatrixXi data, int start_row, int start_col, int band_ids)
-{
- GDALDataType datetype=this->getDataType();
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
- QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- QString filesuffer = getFileExtension(this->img_path).toLower();
- bool isTiff = filesuffer.contains("tif");
- GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if (exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- }
- else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, GDT_Float32, NULL); // 斤拷锟斤拷
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for (int i = 0; i < this->gt.rows(); i++) {
- for (int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- //delete gt_ptr;
- }
-
- long datarows = data.rows();
- long datacols = data.cols();
-
- long* databuffer = new long[datarows * datacols]; // (float*)malloc(datarows * datacols * sizeof(float));
-
- for (long i = 0; i < datarows; i++) {
- for (long j = 0; j < datacols; j++) {
- databuffer[i * datacols + j] = data(i, j);
- }
- }
- // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
- // datarows, GDT_Float32,band_ids, num,0,0,0);
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, datetype, 0, 0);
-
- GDALFlushCache(poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- delete[] databuffer;
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
-}
-
-void gdalImage::saveImage(std::shared_ptr data, int start_row, int start_col, int rowcount, int colcount, int band_ids)
-{
- GDALDataType datetype = this->getDataType();
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if (start_row + rowcount > this->height || start_col + colcount > this->width) {
- QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + rowcount) + ", " + QString::number(start_col + colcount) + ") ";
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- QString filesuffer = getFileExtension(this->img_path).toLower();
- bool isTiff = filesuffer.contains("tif");
- GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if (exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- }
- else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, GDT_Float64, NULL); // 斤拷锟斤拷
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for (int i = 0; i < this->gt.rows(); i++) {
- for (int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- }
-
- long datarows = rowcount;
- long datacols = colcount;
- double* databuffer = new double[datarows * datacols];
- if (datetype == GDT_Float64) {
- memcpy(databuffer, data.get(), sizeof(double) * datarows * datacols);
- }
- else {
- for (long i = 0; i < datarows; i++) {
- for (long j = 0; j < datacols; j++) {
- databuffer[i * datacols + j] = data.get()[i * datacols + j];
- }
- }
- }
-
- // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
- // datarows, GDT_Float32,band_ids, num,0,0,0);
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, datetype, 0, 0);
-
- GDALFlushCache(poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- delete[] databuffer;
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
-}
-
-void gdalImage::saveImage(std::shared_ptr data, int start_row, int start_col, int rowcount, int colcount, int band_ids)
-{
- GDALDataType datetype = this->getDataType();
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if (start_row + rowcount > this->height || start_col + colcount > this->width) {
- QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + rowcount) + ", " + QString::number(start_col + colcount) + ") ";
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- QString filesuffer = getFileExtension(this->img_path).toLower();
- bool isTiff = filesuffer.contains("tif");
- GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if (exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- }
- else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, GDT_Float32, NULL); // 斤拷锟斤拷
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for (int i = 0; i < this->gt.rows(); i++) {
- for (int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- }
-
- long datarows = rowcount;
- long datacols = colcount;
- float* databuffer = new float[datarows * datacols];
- if (datetype == GDT_Float32) {
- memcpy(databuffer, data.get(), sizeof(float) * datarows * datacols);
- }
- else {
- for (long i = 0; i < datarows; i++) {
- for (long j = 0; j < datacols; j++) {
- databuffer[i * datacols + j] = data.get()[i * datacols + j];
- }
- }
- }
-
- // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
- // datarows, GDT_Float32,band_ids, num,0,0,0);
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, datetype, 0, 0);
-
- GDALFlushCache(poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- //delete poDstDS;
- //poDstDS = nullptr;
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- delete[] databuffer;
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
-}
-
-
-void gdalImage::saveImage(std::shared_ptr data, int start_row, int start_col, int rowcount, int colcount, int band_ids)
-{
- GDALDataType datetype = this->getDataType();
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if (start_row + rowcount > this->height || start_col + colcount > this->width) {
- QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + rowcount) + ", " + QString::number(start_col + colcount) + ") ";
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- QString filesuffer = getFileExtension(this->img_path).toLower();
- bool isTiff = filesuffer.contains("tif");
- GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if (exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- }
- else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, GDT_Float32, NULL); // 斤拷锟斤拷
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for (int i = 0; i < this->gt.rows(); i++) {
- for (int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- }
-
- long datarows = rowcount;
- long datacols = colcount;
- int* databuffer = new int[datarows * datacols];
- if (datetype == GDT_Int32) {
- memcpy(databuffer, data.get(), sizeof(int) * datarows * datacols);
- }
- else {
- for (long i = 0; i < datarows; i++) {
- for (long j = 0; j < datacols; j++) {
- databuffer[i * datacols + j] = data.get()[i * datacols + j];
- }
- }
- }
-
- // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
- // datarows, GDT_Float32,band_ids, num,0,0,0);
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, datetype, 0, 0);
-
- GDALFlushCache(poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- delete[] databuffer;
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
-}
-
-
-void gdalImage::saveImage()
-{
- this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
-}
-
-void gdalImage::setNoDataValue(double nodatavalue = -9999, int band_ids = 1)
-{
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
- // GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
- GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
- poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue);
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
-}
-
-void gdalImage::setNoDataValuei(int nodatavalue, int band_ids)
-{
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
- // GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
- GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
- poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue);
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
-}
-
-double gdalImage::getNoDataValue(int band_ids)
-{
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
- // GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
- GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
- double v= poDstDS->GetRasterBand(band_ids)->GetNoDataValue( );
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- return v;
-}
-
-int gdalImage::getNoDataValuei(int band_ids)
-{
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
- // GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
- GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
- int v= poDstDS->GetRasterBand(band_ids)->GetNoDataValue( );
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- return v;
-}
-
-
-int gdalImage::InitInv_gt()
-{
- // 1 lon lat = x
- // 1 lon lat = y
- Eigen::MatrixXd temp = Eigen::MatrixXd::Zero(2, 3);
- this->inv_gt = temp;
- double a = this->gt(0, 0);
- double b = this->gt(0, 1);
- double c = this->gt(0, 2);
- double d = this->gt(1, 0);
- double e = this->gt(1, 1);
- double f = this->gt(1, 2);
- double g = 1;
- double det_gt = b * f - c * e;
- if(det_gt == 0) {
- return 0;
- }
- this->inv_gt(0, 0) = (c * d - a * f) / det_gt; // 2
- this->inv_gt(0, 1) = f / det_gt; // lon
- this->inv_gt(0, 2) = -c / det_gt; // lat
- this->inv_gt(1, 0) = (a * e - b * d) / det_gt; // 1
- this->inv_gt(1, 1) = -e / det_gt; // lon
- this->inv_gt(1, 2) = b / det_gt; // lat
- return 1;
-}
-
-Landpoint gdalImage::getRow_Col(double lon, double lat)
-{
- Landpoint p{ 0, 0, 0 };
- p.lon = this->inv_gt(0, 0) + lon * this->inv_gt(0, 1) + lat * this->inv_gt(0, 2); // x
- p.lat = this->inv_gt(1, 0) + lon * this->inv_gt(1, 1) + lat * this->inv_gt(1, 2); // y
- return p;
-}
-
-Landpoint gdalImage::getLandPoint(double row, double col, double ati = 0)
-{
- Landpoint p{ 0, 0, 0 };
- p.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); // x
- p.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); // y
- p.ati = ati;
- return p;
-}
-
-void gdalImage::getLandPoint(double row, double col, double ati, Landpoint& Lp)
-{
- Lp.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); // x
- Lp.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); // y
- Lp.ati = ati;
-
-}
-
-
-
-
-
-
-double gdalImage::mean(int bandids)
-{
- double mean_value = 0;
- double count = this->height * this->width;
- int line_invert = 100;
- int start_ids = 0;
- do {
- Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
- mean_value = mean_value + sar_a.sum() / count;
- start_ids = start_ids + line_invert;
- } while(start_ids < this->height);
- return mean_value;
-}
-
-double gdalImage::BandmaxValue(int bandids)
-{
- double max_value = 0;
- bool state_max = true;
- int line_invert = 100;
- int start_ids = 0;
- double temp_max = 0;
- do {
- Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
- if(state_max) {
- state_max = false;
- max_value = sar_a.maxCoeff();
- } else {
- temp_max = sar_a.maxCoeff();
- if(max_value < temp_max) {
- max_value = temp_max;
- }
- }
- start_ids = start_ids + line_invert;
- } while(start_ids < this->height);
- return max_value;
-}
-
-double gdalImage::BandminValue(int bandids)
-{
- double min_value = 0;
- bool state_min = true;
- int line_invert = 100;
- int start_ids = 0;
- double temp_min = 0;
- do {
- Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
- if(state_min) {
- state_min = false;
- min_value = sar_a.minCoeff();
- } else {
- temp_min = sar_a.minCoeff();
- if(min_value < temp_min) {
- min_value = temp_min;
- }
- }
- start_ids = start_ids + line_invert;
- } while(start_ids < this->height);
- return min_value;
-}
-
-GDALRPCInfo gdalImage::getRPC()
-{
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
- CPLSetConfigOption("GDAL_DATA", "./data");
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注斤拷锟斤拷
- // 斤拷锟斤拷
- GDALDataset* pDS = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
- // 锟斤拷元斤拷锟叫伙拷取RPC锟斤拷息
- char** papszRPC = pDS->GetMetadata("RPC");
-
- // 斤拷取锟斤拷RPC锟斤拷息斤拷山峁癸拷锟?1锟?7
- GDALRPCInfo oInfo;
- GDALExtractRPCInfo(papszRPC, &oInfo);
-
- GDALClose((GDALDatasetH)pDS);
-
- return oInfo;
-}
-
-Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points)
-{
- if(points.cols() != 3) {
- throw new std::exception("the size of points is equit 3!!!");
- }
-
- Eigen::MatrixXd result(points.rows(), 3);
- result.col(2) = points.col(2); // 锟竭筹拷
- points.col(2) = points.col(2).array() * 0 + 1;
- points.col(0).swap(points.col(2)); // 斤拷
- Eigen::MatrixXd gts(3, 2);
- gts.col(0) = this->gt.row(0);
- gts.col(1) = this->gt.row(1);
-
- result.block(0, 0, points.rows(), 2) = points * gts;
- return result;
-}
-
-Eigen::MatrixXd gdalImage::getHist(int bandids)
-{
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
- // 锟斤拷DEM影锟斤拷
- GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
- this->img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
-
- GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
- GDALRasterBand* xBand = rasterDataset->GetRasterBand(bandids);
-
- double dfMin = this->BandminValue(bandids);
- double dfMax = this->BandmaxValue(bandids);
- int count = int((dfMax - dfMin) / 0.01);
- count = count > 255 ? count : 255;
- GUIntBig* panHistogram = new GUIntBig[count];
- xBand->GetHistogram(dfMin, dfMax, count, panHistogram, TRUE, FALSE, NULL, NULL);
- Eigen::MatrixXd result(count, 2);
- double delta = (dfMax - dfMin) / count;
- for(int i = 0; i < count; i++) {
- result(i, 0) = dfMin + i * delta;
- result(i, 1) = double(panHistogram[i]);
- }
- delete[] panHistogram;
- GDALClose((GDALDatasetH)rasterDataset);
- return result;
-}
-
-RasterExtend gdalImage::getExtend()
-{
- RasterExtend extend{ 0,0,0,0 };
- double x1 = this->gt(0, 0);
- double y1 = this->gt(1, 0);
-
- double x2 = this->gt(0, 0) + (this->width - 1) * gt(0, 1) + (0) * gt(0, 2);
- double y2 = this->gt(1, 0) + (this->width - 1) * gt(1, 1) + (0) * gt(1, 2);
-
- double x3 = this->gt(0, 0) + (0) * gt(0, 1) + (this->height - 1) * gt(0, 2);
- double y3 = this->gt(1, 0) + (0) * gt(1, 1) + (this->height - 1) * gt(1, 2);
-
- double x4 = this->gt(0, 0) + (this->width - 1) * gt(0, 1) + (this->height - 1) * gt(0, 2);
- double y4 = this->gt(1, 0) + (this->width - 1) * gt(1, 1) + (this->height - 1) * gt(1, 2);
-
-
- extend.min_x = x1 < x2 ? x1 : x2;
- extend.max_x = x1 < x2 ? x2 : x1;
- extend.min_y = y1 < y2 ? y1 : y2;
- extend.max_y = y1 < y2 ? y2 : y1;
-
-
- extend.min_x = extend.min_x < x3 ? extend.min_x : x3;
- extend.max_x = extend.max_x > x3 ? extend.max_x : x3;
- extend.min_y = extend.min_y < y3 ? extend.min_y : y3;
- extend.max_y = extend.max_y > y3 ? extend.max_y : y3;
-
-
- extend.min_x = extend.min_x < x4 ? extend.min_x : x4;
- extend.max_x = extend.max_x > x4 ? extend.max_x : x4;
- extend.min_y = extend.min_y < y4 ? extend.min_y : y4;
- extend.max_y = extend.max_y > y4 ? extend.max_y : y4;
-
- return extend;
-}
-
-gdalImage BASECONSTVARIABLEAPI CreategdalImageDouble(QString& img_path, int height, int width, int band_num, bool overwrite, bool isEnvi)
-{
-
- if (exists_test(img_path.toUtf8().constData())) {
- if (overwrite) {
- gdalImage result_img(img_path);
- return result_img;
- }
- else {
- throw "file has exist!!!";
- exit(1);
- }
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
- GDALDriver* poDriver = nullptr;
- if (isEnvi) {
- poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- }
- else {
- poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
- }
-
-
- GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,GDT_Float64, NULL); // 锟斤拷锟斤拷锟斤拷
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- gdalImage result_img(img_path);
- return result_img;
-
-}
-
-gdalImage CreategdalImageDouble(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite, bool isEnvi)
-{
- if (exists_test(img_path.toUtf8().constData())) {
- if (overwrite) {
- gdalImage result_img(img_path);
- return result_img;
- }
- else {
- throw "file has exist!!!";
- exit(1);
- }
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
- GDALDriver* poDriver = nullptr;
- if (isEnvi) {
- poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- }
- else {
- poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
- }
-
-
- GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
- GDT_Float64, NULL); // 锟斤拷锟斤拷锟斤拷
- if (need_gt) {
- if (!projection.isEmpty()) {
- poDstDS->SetProjection(projection.toUtf8().constData());
- }
- double gt_ptr[6] = { 0 };
- for (int i = 0; i < gt.rows(); i++) {
- for (int j = 0; j < gt.cols(); j++) {
- gt_ptr[i * 3 + j] = gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- }
- for (int i = 1; i <= band_num; i++) {
- poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
- }
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- gdalImage result_img(img_path);
- return result_img;
-}
-
-gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num,
- Eigen::MatrixXd gt, QString projection,bool need_gt, bool overwrite, bool isEnvi, GDALDataType datetype)
-{
- if(exists_test(img_path.toUtf8().constData())) {
- if(overwrite) {
- gdalImage result_img(img_path);
- return result_img;
- } else {
- throw "file has exist!!!";
- exit(1);
- }
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
- GDALDriver* poDriver = nullptr;
- if (isEnvi) {
- poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- }
- else {
- poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
- }
-
-
- GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
- datetype, NULL); // 锟斤拷锟斤拷锟斤拷
- if(need_gt) {
- if (!projection.isEmpty()) {
- poDstDS->SetProjection(projection.toUtf8().constData());
- }
- double gt_ptr[6] = { 0 };
- for(int i = 0; i < gt.rows(); i++) {
- for(int j = 0; j < gt.cols(); j++) {
- gt_ptr[i * 3 + j] = gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- }
- for(int i = 1; i <= band_num; i++) {
- poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
- }
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- gdalImage result_img(img_path);
- return result_img;
-}
-
-gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long epsgCode, GDALDataType eType, bool need_gt, bool overwrite, bool isENVI)
-{
- if (exists_test(img_path.toUtf8().constData())) {
- if (overwrite) {
- gdalImage result_img(img_path);
- return result_img;
- }
- else {
- throw "file has exist!!!";
- exit(1);
- }
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
- GDALDriver* poDriver = isENVI? GetGDALDriverManager()->GetDriverByName("ENVI"): GetGDALDriverManager()->GetDriverByName("GTiff");
- GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, eType, NULL); // 锟斤拷锟斤拷锟斤拷
- if (need_gt) {
- OGRSpatialReference oSRS;
-
- if (oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) {
- qDebug() << "Failed to import EPSG code " << epsgCode ;
- throw "Failed to import EPSG code ";
- exit(1);
- }
- char* pszWKT = NULL;
- oSRS.exportToWkt(&pszWKT);
- qDebug() << "WKT of EPSG:"<< epsgCode <<" :\n" << pszWKT ;
- poDstDS->SetProjection(pszWKT);
- double gt_ptr[6] = { 0 };
- for (int i = 0; i < gt.rows(); i++) {
- for (int j = 0; j < gt.cols(); j++) {
- gt_ptr[i * 3 + j] = gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- }
- for (int i = 1; i <= band_num; i++) {
- poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
- }
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- gdalImage result_img(img_path);
- return result_img;
-}
-
-gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width,
- int band_num, Eigen::MatrixXd gt, QString projection,
- bool need_gt, bool overwrite)
-{
- if(exists_test(img_path.toUtf8().constData())) {
- if(overwrite) {
- gdalImageComplex result_img(img_path);
- return result_img;
- } else {
- throw "file has exist!!!";
- exit(1);
- }
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
- GDT_CFloat64, NULL);
- if(need_gt) {
- poDstDS->SetProjection(projection.toUtf8().constData());
-
- // 锟斤拷锟斤拷转锟斤拷锟斤拷锟斤拷
- double gt_ptr[6] = { 0 };
- for(int i = 0; i < gt.rows(); i++) {
- for(int j = 0; j < gt.cols(); j++) {
- gt_ptr[i * 3 + j] = gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- }
- //for(int i = 1; i <= band_num; i++) {
- // poDstDS->GetRasterBand(i)->SetNoDataValue(0); // 回波部分
- //}
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- gdalImageComplex result_img(img_path);
- return result_img;
-}
-
-gdalImageComplex BASECONSTVARIABLEAPI CreategdalImageComplexNoProj(const QString& img_path, int height, int width, int band_num, bool overwrite)
-{
-
- // 创建新文件
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
-
- GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = (poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_CFloat64, NULL));
- GDALFlushCache((GDALDatasetH)poDstDS);
- GDALClose((GDALDatasetH)poDstDS);
- //poDstDS.reset();
- omp_unset_lock(&lock); //
- omp_destroy_lock(&lock); //
-
- gdalImageComplex result_img(img_path);
- return result_img;
-}
-
-gdalImageComplex CreateEchoComplex(const QString& img_path, int height, int width, int band_num)
-{
- // 创建图像
- Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
- //Xgeo = GeoTransform[0] + Xpixel * GeoTransform[1] + Ypixel * GeoTransform[2]
- //Ygeo = GeoTransform[3] + Xpixel * GeoTransform[4] + Ypixel * GeoTransform[5]
- // X
- gt(0, 0) = 0; gt(0, 2) = 1; gt(0, 2) = 0;
- gt(1, 0) = 0; gt(1, 1) = 0; gt(1, 2) = 1;
- // Y
- QString projection = "";
- gdalImageComplex echodata = CreategdalImageComplex(img_path, height, width, 1, gt, projection, false, true);
- return echodata;
-
-}
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width,
int new_height, GDALResampleAlg eResample)
@@ -2163,8 +702,6 @@ int alignRaster(QString inputPath, QString referencePath, QString outputPath, GD
}
-
-
int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path)
{
@@ -2208,723 +745,6 @@ void clipRaster(QString inRasterPath, QString outRasterPath, long minRow, long m
}
-ErrorCode MergeRasterProcess(QVector filepaths, QString outfileptah, QString mainString, MERGEMODE mergecode, bool isENVI, ShowProessAbstract* dia )
-{
- // 参数检查
- if (!isExists(mainString)) {
- qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND) )<< "\t" << mainString;
- return ErrorCode::FILENOFOUND;
- }
- else {}
- gdalImage mainimg(mainString);
- QVector imgdslist(filepaths.count());
- for (long i = 0; i < filepaths.count(); i++) {
- if (!isExists(filepaths[i])) {
- qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND)) << "\t" << filepaths[i];
- return ErrorCode::FILENOFOUND;
- }
- else {
- imgdslist[i] = gdalImage(filepaths[i]);
- if (imgdslist[i].band_num != mainimg.band_num) {
- qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTERBAND_NOTEQUAL)) << "\t" << imgdslist[i].band_num <<" != "<< mainimg.band_num;
- return ErrorCode::RASTERBAND_NOTEQUAL;
- }
- }
- }
-
- // 检查坐标系是否统一
- long EPSGCode = GetEPSGFromRasterFile(mainString);
- long tempCode = 0;
- for (long i = 0; i < filepaths.count(); i++) {
- tempCode = GetEPSGFromRasterFile(filepaths[i]);
- if (EPSGCode != tempCode) {
- qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSAME)) << "\t" << EPSGCode <<"!="<< tempCode;
- return ErrorCode::EPSGCODE_NOTSAME;
- }
- }
-
- // 检查影像类型是否统一
- GDALDataType mainType = mainimg.getDataType();
- for (long i = 0; i < imgdslist.count(); i++) {
- if (mainType != imgdslist[i].getDataType()) {
- qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTER_DATETYPE_NOTSAME)) << "\t" << mainType << "!=" << imgdslist[i].getDataType();
- return ErrorCode::RASTER_DATETYPE_NOTSAME;
- }
- }
-
- Eigen::MatrixXd maingt = mainimg.getGeoTranslation();
- Eigen::MatrixXd rgt = Eigen::MatrixXd::Zero(2,3);
- RasterExtend mainExtend = mainimg.getExtend();
- rgt(0, 1) = (mainExtend.max_x - mainExtend.min_x) / (mainimg.width - 1); //dx
- rgt(1, 2) = -1*std::abs(( (mainExtend.max_y - mainExtend.min_y) / (mainimg.height - 1)));//dy
- QVector extendlist(imgdslist.count());
- for (long i = 0; i < imgdslist.count(); i++) {
- extendlist[i] = imgdslist[i].getExtend();
- mainExtend.min_x = mainExtend.min_x < extendlist[i].min_x ? mainExtend.min_x : extendlist[i].min_x;
- mainExtend.max_x = mainExtend.max_x > extendlist[i].max_x ? mainExtend.max_x : extendlist[i].max_x;
- mainExtend.min_y = mainExtend.min_y < extendlist[i].min_y ? mainExtend.min_y : extendlist[i].min_y;
- mainExtend.max_y = mainExtend.max_y > extendlist[i].max_y ? mainExtend.max_y : extendlist[i].max_y;
- }
-
- rgt(0, 0) = mainExtend.min_x;
- rgt(1, 0) = mainExtend.max_y;
-
- // 计算数量
-
- long width = std::ceil((mainExtend.max_x - mainExtend.min_x) / rgt(0, 1) + 1);
- long height = std::ceil(std::abs((mainExtend.min_y - mainExtend.max_y) / rgt(1, 2)) + 1);
- OGRSpatialReference oSRS;
- if (oSRS.importFromEPSG(EPSGCode) != OGRERR_NONE) {
- qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSUPPORT)) << "\t" << EPSGCode;
- return ErrorCode::EPSGCODE_NOTSUPPORT;
- }
-
- gdalImage resultImage = CreategdalImage(outfileptah, height, width, mainimg.band_num, rgt, EPSGCode, mainType, true, true, isENVI);
-
- QString resultMaskString = addMaskToFileName(outfileptah, QString("_MASK"));
- gdalImage maskImage = CreategdalImage(resultMaskString, height, width,1 , rgt, EPSGCode, GDT_Int32, true, true, isENVI);
-
- // 初始化
- long resultline = Memory1MB * 500 / 8 / resultImage.width;
- resultline = resultline < 10000 ? resultline : 10000; // 最多100行
- resultline = resultline > 0 ? resultline : 2;
- long bandnum = resultImage.band_num + 1;
- long starti = 0;
- long rasterCount = imgdslist.count();
-
-
-
- QProgressDialog progressDialog(u8"初始化影像", u8"终止", 0, resultImage.height);
- progressDialog.setWindowTitle(u8"初始化影像");
- progressDialog.setWindowModality(Qt::WindowModal);
- progressDialog.setAutoClose(true);
- progressDialog.setValue(0);
- progressDialog.setMaximum(resultImage.height);
- progressDialog.setMinimum(0);
- progressDialog.show();
-
-
- for (starti = 0; starti < resultImage.height; starti = starti + resultline) {
- long blocklines = resultline;
- blocklines = starti + blocklines < resultImage.height ? blocklines : resultImage.height - starti;
- for (long b = 1; b < bandnum; b++) {
- Eigen::MatrixXd data = resultImage.getData(starti, 0, blocklines, resultImage.width, b);
- Eigen::MatrixXi maskdata = maskImage.getDatai(starti, 0, blocklines, resultImage.width, b);
- data = data.array() * 0;
- maskdata = maskdata.array() * 0;
- resultImage.saveImage(data, starti, 0, b);
- maskImage.saveImage(maskdata, starti, 0, b);
- }
- if (nullptr != dia) {
- dia->showProcess(starti * 1.0 / resultImage.height, u8"初始化影像数据");
- }
- progressDialog.setValue(starti+ blocklines);
- }
- progressDialog.close();
-
-
-
- switch (mergecode)
- {
- case MERGE_GEOCODING:
- return MergeRasterInGeoCoding(imgdslist, resultImage, maskImage, dia);
- default:
- break;
- }
-
-
- return ErrorCode::SUCCESS;
-}
-
-ErrorCode MergeRasterInGeoCoding(QVector imgdslist, gdalImage resultimg, gdalImage maskimg, ShowProessAbstract* dia)
-{
- omp_set_num_threads(Paral_num_thread);
- // 逐点合并计算
- QVector extendlist(imgdslist.count());
- for (long i = 0; i < imgdslist.count(); i++) {
- extendlist[i] = imgdslist[i].getExtend();
- imgdslist[i].InitInv_gt();
- }
-
- // 分块计算
- long resultline = Memory1MB * 1000 / 8 / resultimg.width;
- resultline = resultline < 300 ? resultline : 300; // 最多100行
-
- long bandnum = resultimg.band_num+1;
- long starti = 0;
- long rasterCount = imgdslist.count();
-
- long processNumber = 0;
- QProgressDialog progressDialog(u8"合并影像", u8"终止", 0, resultimg.height);
- progressDialog.setWindowTitle(u8"合并影像");
- progressDialog.setWindowModality(Qt::WindowModal);
- progressDialog.setAutoClose(true);
- progressDialog.setValue(0);
- progressDialog.setMaximum(resultimg.height);
- progressDialog.setMinimum(0);
- progressDialog.show();
- omp_lock_t lock;
- omp_init_lock(&lock);
-
-#pragma omp parallel for
- for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
- long blocklines = resultline;
- blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
-
- long rid = starti;
- long cid = 0;
-
- Landpoint pp = {0,0,0};
- Landpoint lp = { 0,0,0 };
-
- for (long ir = 0; ir < rasterCount; ir++) {// 影像
- long minRid = imgdslist[ir].height;
- long maxRid = 0;
-
- Eigen::MatrixXd ridlist = resultimg.getData(starti, 0, blocklines, resultimg.width, 1);
- ridlist = ridlist.array() * 0;
- Eigen::MatrixXd cidlist = ridlist.array() * 0;
-
- for (long i = 0; i < blocklines; i++) {// 行号
- rid = starti + i;
- for (long j = 0; j < resultimg.width; j++) {// 列号
- cid = j;
- resultimg.getLandPoint(rid, cid,0,pp);
- lp = imgdslist[ir].getRow_Col(pp.lon, pp.lat); // 获取点坐标
- ridlist(i, j) = lp.lat;
- cidlist(i, j) = lp.lon;
- }
- }
-
- //ImageShowDialogClass* dialog = new ImageShowDialogClass;
- //dialog->show();
- //dialog->load_double_MatrixX_data(cidlist, u8"");
-
- //dialog->exec();
-
-
- if (ridlist.maxCoeff() < 0 || ridlist.minCoeff() >= imgdslist[ir].height) {
- continue;
- }
-
- if (cidlist.maxCoeff() < 0 || cidlist.minCoeff() >= imgdslist[ir].width) {
- continue;
- }
-
- minRid = std::floor(ridlist.minCoeff());
- maxRid = std::ceil(ridlist.maxCoeff());
- minRid = minRid < 0 ? 0 : minRid;
- maxRid = maxRid < imgdslist[ir].height ? maxRid : imgdslist[ir].height - 1;
-
- long rowlen = maxRid - minRid + 1;
- if(rowlen <= 0) {
- continue;
- }
- // 获取分配代码
- Landpoint p0{ 0,0,0 }, p11{ 0,0,0 }, p21{ 0,0,0 }, p12{ 0,0,0 }, p22{ 0,0,0 }, p{ 0,0,0 };
-
-
- long rowcount = 0;
- long colcount = 0;
- double ridtemp = 0, cidtemp = 0;
-
- long lastr = 0, nextr = 0;
- long lastc = 0, nextc = 0;
-
- double r0=0, c0 = 0;
-
- for (long b = 1; b < bandnum; b++) {
- Eigen::MatrixXd resultdata = resultimg.getData(starti, 0, blocklines, resultimg.width, b);
- Eigen::MatrixXi resultmask = maskimg.getDatai(starti, 0, blocklines, resultimg.width, b);
- Eigen::MatrixXd data = imgdslist[ir].getData(minRid, 0, rowlen, imgdslist[ir].width, b);
-
- double nodata = imgdslist[ir].getNoDataValue(b);
- for (long ii = 0; ii < data.rows(); ii++) {
- for (long jj = 0; jj < data.cols(); jj++) {
- if (std::abs(data(ii, jj) - nodata) < 1e-6) {
- data(ii, jj) = 0;
- }
- }
- }
- rowcount = ridlist.rows();
- colcount = ridlist.cols();
- double Bileanervalue = 0;
- for (long i = 0; i < rowcount; i++) {
- for (long j = 0; j < colcount; j++) {
- ridtemp = ridlist(i, j);
- cidtemp = cidlist(i, j);
-
- lastr = std::floor(ridtemp);
- nextr = std::ceil(ridtemp);
- lastc = std::floor(cidtemp);
- nextc = std::ceil(cidtemp);
-
- if (lastr < 0 || lastr >= imgdslist[ir].height
- || nextr < 0 || nextr >= imgdslist[ir].height
- || lastc < 0 || lastc >= imgdslist[ir].width
- || nextc <0|| nextc >=imgdslist[ir].width) {
- continue;
- }
- else {}
-
- r0 = ridtemp - std::floor(ridtemp);
- c0 = cidtemp - std::floor(cidtemp);
-
- lastr = lastr - minRid;
- nextr = nextr - minRid;
-
- p0 = Landpoint{ c0,r0,0 };
- p11 = Landpoint{ 0,0,data(lastr,lastc) };
- p21 = Landpoint{ 0,1,data(nextr,lastc) };
- p12 = Landpoint{ 1,0,data(lastr,nextc) };
- p22 = Landpoint{ 1,1,data(nextr,nextc) };
- Bileanervalue = Bilinear_interpolation(p0, p11, p21, p12, p22);
- if (std::abs(Bileanervalue) < 1e-6||resultmask(i, j)>0) {
- continue;
- }
- resultdata(i,j) = resultdata(i, j)+ Bileanervalue;
- resultmask(i, j) = resultmask(i, j) + 1;
- }
- }
- resultimg.saveImage(resultdata, starti, 0, b);
- maskimg.saveImage(resultmask, starti, 0, b);
- }
- }
-
- omp_set_lock(&lock);
- processNumber = processNumber + blocklines;
- qDebug() << "\rprocess bar:\t" << processNumber * 100.0 / resultimg.height << " % " << "\t\t\t";
- if (nullptr != dia) {
- dia->showProcess(processNumber * 1.0 / resultimg.height, u8"合并图像");
- }
- if (progressDialog.maximum() <= processNumber) {
- processNumber = progressDialog.maximum() - 1;
- }
- progressDialog.setValue(processNumber);
- omp_unset_lock(&lock);
- }
- omp_destroy_lock(&lock);
-
- progressDialog.setWindowTitle(u8"影像掩膜");
- progressDialog.setLabelText(u8"影像掩膜");
- for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
- long blocklines = resultline;
- blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
- for (long b = 1; b < bandnum; b++) {
- Eigen::MatrixXd data = resultimg.getData(starti, 0, blocklines, resultimg.width, b);
- Eigen::MatrixXd maskdata = maskimg.getData(starti, 0, blocklines, maskimg.width, b);
-
- for (long i = 0; i < data.rows(); i++) {
- for (long j = 0; j < data.cols(); j++) {
- if (maskdata(i, j) == 0) {
- data(i, j) = -9999;
- continue;
- }
- data(i, j) = data(i, j) / maskdata(i, j);
- }
- }
-
- resultimg.saveImage(data, starti, 0, b);
- maskimg.saveImage(maskdata, starti, 0, b);
- }
- if (nullptr != dia) {
- dia->showProcess((starti + blocklines) * 1.0 / resultimg.height, u8"影像掩膜");
- }
- progressDialog.setValue(starti + blocklines);
- }
- resultimg.setNoDataValue(-9999);
- progressDialog.close();
- return ErrorCode::SUCCESS;
-}
-
-bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath)
-{
-
- Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
-
- gdalImage img= CreategdalImageDouble(dataStrPath, data.rows(), data.cols(), 1, gt, "", false, true,true);
-
- img.saveImage(data, 0,0,1);
-
- return true;
-}
-
-gdalImageComplex::gdalImageComplex(const QString& raster_path)
-{
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- this->img_path = raster_path;
-
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
- raster_path.toUtf8().constData(), GA_ReadOnly));
- this->width = rasterDataset->GetRasterXSize();
- this->height = rasterDataset->GetRasterYSize();
- this->band_num = rasterDataset->GetRasterCount();
-
- double* gt = new double[6];
- rasterDataset->GetGeoTransform(gt);
- this->gt = Eigen::MatrixXd(2, 3);
- this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5];
-
- double a = this->gt(0, 0);
- double b = this->gt(0, 1);
- double c = this->gt(0, 2);
- double d = this->gt(1, 0);
- double e = this->gt(1, 1);
- double f = this->gt(1, 2);
-
- this->projection = rasterDataset->GetProjectionRef();
-
- // 释放投影
- GDALFlushCache((GDALDatasetH)rasterDataset);
- GDALClose((GDALDatasetH)rasterDataset);
- rasterDataset = NULL; // 指矫匡拷
- this->InitInv_gt();
- delete[] gt;
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- omp_unset_lock(&lock); // 锟酵放伙拷斤拷
- omp_destroy_lock(&lock); // 劫伙拷斤拷
-}
-
-gdalImageComplex::~gdalImageComplex() {}
-
-void gdalImageComplex::setData(Eigen::MatrixXcd data)
-{
- this->data = data;
-}
-
-void gdalImageComplex::saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids)
-{
-
-
-
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if(start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
- QString tip = u8"file path: " + this->img_path;
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if(exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- } else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, GDT_CFloat64, NULL); // 斤拷锟斤拷
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for(int i = 0; i < this->gt.rows(); i++) {
- for(int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- //delete[] gt_ptr;
- }
-
- int datarows = data.rows();
- int datacols = data.cols();
-
- double* databuffer = new double[data.size() * 2];
- for(int i = 0; i < data.rows(); i++) {
- for(int j = 0; j < data.cols(); j++) {
- databuffer[i * data.cols() * 2 + j * 2] = data(i, j).real();
- databuffer[i * data.cols() * 2 + j * 2 + 1] = data(i, j).imag();
- }
- }
-
- // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
- // datarows, GDT_Float32,band_ids, num,0,0,0);
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
- databuffer, datacols, datarows, GDT_CFloat64, 0, 0);
- GDALFlushCache(poDstDS);
- delete databuffer;
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- omp_unset_lock(&lock); //
- omp_destroy_lock(&lock); //
-}
-
-void gdalImageComplex::saveImage(std::shared_ptr> data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
-{
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if (start_row + rowCount > this->height || start_col + colCount > this->width) {
- QString tip = u8"file path: " + this->img_path;
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- return;
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if (exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- }
- else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, GDT_CFloat64, NULL); // 斤拷锟斤拷
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for (int i = 0; i < this->gt.rows(); i++) {
- for (int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- //delete[] gt_ptr;
- }
-
- double* databuffer = new double[rowCount * colCount * 2];
- for (long i = 0; i < rowCount; i++) {
- for (long j = 0; j < colCount; j++) {
- databuffer[i * colCount * 2 + j * 2] = data.get()[i*colCount+j].real();
- databuffer[i * colCount * 2 + j * 2 + 1] = data.get()[i * colCount + j].imag();
- }
- }
-
- // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
- // datarows, GDT_Float32,band_ids, num,0,0,0);
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, colCount, rowCount,
- databuffer, colCount, rowCount, GDT_CFloat64, 0, 0);
- GDALFlushCache(poDstDS);
- delete databuffer;
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- omp_unset_lock(&lock); //
- omp_destroy_lock(&lock); //
-
-}
-
-void gdalImageComplex::saveImage(std::complex* data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
-{
- omp_lock_t lock;
- omp_init_lock(&lock);
- omp_set_lock(&lock);
- if (start_row + rowCount > this->height || start_col + colCount > this->width) {
- QString tip = u8"file path: " + this->img_path;
- qDebug() << tip;
- throw std::exception(tip.toUtf8().constData());
- return;
- }
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
- GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- GDALDataset* poDstDS = nullptr;
- if (exists_test(this->img_path)) {
- poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
- }
- else {
- poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
- this->band_num, GDT_CFloat64, NULL); // 斤拷锟斤拷
- poDstDS->SetProjection(this->projection.toUtf8().constData());
-
- double gt_ptr[6];
- for (int i = 0; i < this->gt.rows(); i++) {
- for (int j = 0; j < this->gt.cols(); j++) {
- gt_ptr[i * 3 + j] = this->gt(i, j);
- }
- }
- poDstDS->SetGeoTransform(gt_ptr);
- //delete[] gt_ptr;
- }
-
- double* databuffer = new double[rowCount * colCount * 2];
- for (long i = 0; i < rowCount; i++) {
- for (long j = 0; j < colCount; j++) {
- databuffer[i * colCount * 2 + j * 2] = data[i * colCount + j].real();
- databuffer[i * colCount * 2 + j * 2 + 1] = data[i * colCount + j].imag();
- }
- }
-
- // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
- // datarows, GDT_Float32,band_ids, num,0,0,0);
- poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, colCount, rowCount,
- databuffer, colCount, rowCount, GDT_CFloat64, 0, 0);
- GDALFlushCache(poDstDS);
- delete databuffer;
- GDALClose((GDALDatasetH)poDstDS);
- GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
- omp_unset_lock(&lock); //
- omp_destroy_lock(&lock); //
-}
-
-Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col, int rows_count,
- int cols_count, int band_ids)
-{
- GDALDataset* poDataset;
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
-
- // 打开TIFF文件
- poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
- if(poDataset == nullptr) {
- QMessageBox::warning(nullptr, u8"错误", u8"无法打开文件:" + this->img_path);
- qDebug() << u8"无法打开文件:" + this->img_path;
- }
-
- // 获取数据集的第一个波段
- GDALRasterBand* poBand;
- poBand = poDataset->GetRasterBand(1);
-
- // 读取波段信息,假设是复数类型
- int nXSize = cols_count; poBand->GetXSize();
- int nYSize = rows_count; poBand->GetYSize();
-
- double* databuffer = new double[nXSize * nYSize * 2];
- poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count,
- rows_count, GDT_CFloat64, 0, 0);
- GDALClose((GDALDatasetH)poDataset);
- Eigen::MatrixXcd rasterData(nYSize, nXSize); // 使用Eigen的MatrixXcd
- for(size_t i = 0; i < nYSize; i++) {
- for(size_t j = 0; j < nXSize; j++) {
- rasterData(i, j) = std::complex(databuffer[i * nXSize * 2 + j * 2],
- databuffer[i * nXSize * 2 + j * 2 + 1]);
- }
- }
-
- delete[] databuffer;
- return rasterData;
-}
-
-std::shared_ptr> gdalImageComplex::getDataComplexSharePtr(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
-{
- GDALDataset* poDataset;
- GDALAllRegister();
- CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
-
- // 打开TIFF文件
- poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
- if(poDataset == nullptr) {
- QMessageBox::warning(nullptr, u8"错误", u8"无法打开文件:" + this->img_path);
- qDebug() << u8"无法打开文件:" + this->img_path;
- }
-
- // 获取数据集的第一个波段
- GDALRasterBand* poBand;
- poBand = poDataset->GetRasterBand(1);
-
- // 读取波段信息,假设是复数类型
-
-
- double* databuffer = new double[rows_count * cols_count * 2];
- poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count,
- rows_count, GDT_CFloat64, 0, 0);
- GDALClose((GDALDatasetH)poDataset);
-
- std::shared_ptr> rasterData(new std::complex[rows_count * cols_count], delArrPtr);
-
- for(size_t i = 0; i < rows_count; i++) {
- for(size_t j = 0; j < cols_count; j++) {
- rasterData.get()[i*cols_count+ j] = std::complex(databuffer[i * cols_count * 2 + j * 2],
- databuffer[i * cols_count * 2 + j * 2 + 1]);
- }
- }
-
- delete[] databuffer;
- return rasterData;
-}
-
-void gdalImageComplex::saveImage()
-{
- this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
-}
-void gdalImageComplex::savePreViewImage()
-{
- qDebug()<<"void gdalImageComplex::savePreViewImage()";
- Eigen::MatrixXd data_abs = Eigen::MatrixXd::Zero(this->height, this->width);
- data_abs = (this->data.array().real().pow(2) + this->data.array().imag().pow(2))
- .array()
- .log10()*10.0; // 计算振幅
-
- double min_abs = data_abs.minCoeff(); // 最大值
- double max_abs = data_abs.maxCoeff(); // 最小值
- double delta = (max_abs - min_abs) / 1000; // 1000分位档
- Eigen::MatrixX data_idx =
- ((data_abs.array() - min_abs).array() / delta).array().floor().cast();
-
- // 初始化
- double hist[1001];
- for(size_t i = 0; i < 1001; i++) {
- hist[i] = 0; // 初始化
- }
- for(size_t i = 0; i < this->height; i++) {
- for(size_t j = 0; j < this->width; j++) {
- hist[data_idx(i, j)]++;
- }
- }
-
- // 统计
- size_t count = this->height * this->width;
- double precent = 0;
- size_t curCount = 0;
- double pre2 = 0;
- bool findprec_2 = true;
- double pre98 = 0;
- bool findprec_98 = true;
- for(size_t i = 0; i < 1001; i++) {
- precent = precent + hist[i];
- if(findprec_2 & precent / count > 0.02) {
- pre2 = i * delta + min_abs;
- findprec_2 = false;
- }
- if(findprec_98 & precent / count > 0.98) {
- pre98 = (i-1) * delta + min_abs;
- findprec_98 = false;
- }
- }
- // 拉伸
- delta = (pre98-pre2)/200;
- data_idx=
- ((data_abs.array() - pre2).array() / delta).array().floor().cast();
-
- for(size_t i = 0; i < this->height; i++) {
- for(size_t j = 0; j < this->width; j++) {
- if(data_idx(i,j)<0){
- data_idx(i,j)=0;
- }
- else if(data_idx(i,j)>255){
- data_idx(i,j)=255;
- }else{
-
- }
- }
- }
-
- // 赋值
- QString filePath = this->img_path;
- QFile file(filePath);
- QFileInfo fileInfo(file);
-
- QString directory = fileInfo.absolutePath();
- qDebug() << "文件所在目录:" << directory;
- QString baseName = fileInfo.completeBaseName();
- qDebug() << "无后缀文件名:" << baseName;
-
- // 创建文件路径
- QString previewImagePath = JoinPath(directory, baseName+"_preview.png");
- cv::Mat img(this->height, this->width, CV_8U ,cv::Scalar(0));
-
- for(size_t i = 0; i < this->height; i++) {
- for(size_t j = 0; j < this->width; j++) {
- img.at(i,j)= (uchar)(data_idx(i,j));
- }
- }
- //std::string outimgpath=previewImagePath.toUtf8().data();
- cv::imwrite(previewImagePath.toUtf8().data(), img);
-}
-
@@ -3028,6 +848,9 @@ QString GetProjectionNameFromEPSG(long epsgCode)
// return projName;
}
+
+
+
long GetEPSGFromRasterFile(QString filepath)
{
qDebug() << "============= GetEPSGFromRasterFile ======================";
@@ -3092,6 +915,7 @@ long GetEPSGFromRasterFile(QString filepath)
}
}
}
+
std::shared_ptr GetCenterPointInRaster(QString filepath)
{
qDebug() << "============= GetCenterPointInRaster ======================";
@@ -3192,21 +1016,15 @@ CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long epsg_code)
else if (oSRS.IsProjected()) {
return CoordinateSystemType::ProjectCoordinateSystem;
}
+ else {
+ return CoordinateSystemType::UNKNOW;
+ }
}
else {
return CoordinateSystemType::UNKNOW;
}
}
-void ShowProessAbstract::showProcess(double precent, QString tip)
-{
-
-}
-
-void ShowProessAbstract::showToolInfo(QString tip)
-{
-}
-
void resampleRaster(const char* inputRaster, const char* outputRaster, double targetPixelSizeX, double targetPixelSizeY) {
// 初始化GDAL
@@ -3557,239 +1375,6 @@ ErrorCode DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QStri
return ErrorCode::SUCCESS;
}
-void testOutAmpArr(QString filename, float* amp, long rowcount, long colcount)
-{
-
-
- Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
-
- for (long hii = 0; hii < rowcount; hii++) {
- for (long hjj = 0; hjj < colcount; hjj++) {
- h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
- }
- }
- QString ampPath = getDebugDataPath(filename);
- saveEigenMatrixXd2Bin(h_amp_img, ampPath);
- qDebug() << filename.toLocal8Bit().constData() ;
- qDebug() << "max:\t" << h_amp_img.maxCoeff() ;
- qDebug() << "min:\t" << h_amp_img.minCoeff() ;
-}
-
-void testOutAmpArr(QString filename, double* amp, long rowcount, long colcount)
-{
-
-
- Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
-
- for (long hii = 0; hii < rowcount; hii++) {
- for (long hjj = 0; hjj < colcount; hjj++) {
- h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
- }
- }
- QString ampPath = getDebugDataPath(filename);
- saveEigenMatrixXd2Bin(h_amp_img, ampPath);
- qDebug() << filename.toLocal8Bit().constData() ;
- qDebug() << "max:\t" << h_amp_img.maxCoeff() ;
- qDebug() << "min:\t" << h_amp_img.minCoeff() ;
-}
-
-
-void testOutClsArr(QString filename, long* amp, long rowcount, long colcount) {
-
- Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
-
- for (long hii = 0; hii < rowcount; hii++) {
- for (long hjj = 0; hjj < colcount; hjj++) {
- h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
- }
- }
- QString ampPath = getDebugDataPath(filename);
- saveEigenMatrixXd2Bin(h_amp_img, ampPath);
- qDebug() << filename.toLocal8Bit().constData() ;
- qDebug() << "max:\t" << h_amp_img.maxCoeff() ;
- qDebug() << "min:\t" << h_amp_img.minCoeff() ;
-
-}
-
-void BASECONSTVARIABLEAPI testOutComplexDoubleArr(QString filename, std::complex* data, long rowcount, long colcount)
-{
- gdalImageComplex compleximg= CreateEchoComplex(filename, rowcount, colcount, 1);
- compleximg.saveImage( data, 0, 0, rowcount, colcount, 1);
-
- return void BASECONSTVARIABLEAPI();
-}
-
-void BASECONSTVARIABLEAPI testOutDataArr(QString filename, double* data, long rowcount, long colcount)
-{
- return testOutAmpArr(filename, data, rowcount, colcount);
-}
-
-void BASECONSTVARIABLEAPI testOutDataArr(QString filename, float* data, long rowcount, long colcount)
-{
- return testOutAmpArr(filename, data, rowcount, colcount);
-}
-
-void BASECONSTVARIABLEAPI testOutDataArr(QString filename, long* data, long rowcount, long colcount)
-{
- return testOutClsArr(filename,data,rowcount,colcount);
-}
-
-void testOutAntPatternTrans(QString antpatternfilename, double* antPatternArr,
- double starttheta, double deltetheta,
- double startphi, double deltaphi,
- long thetanum, long phinum)
-{
-
-
- Eigen::MatrixXd antPatternMatrix(thetanum, phinum);
- for (long t = 0; t < thetanum; ++t) {
- for (long p = 0; p < phinum; ++p) {
- long index = t * phinum + p;
- if (index < thetanum * phinum) {
- antPatternMatrix(t, p) = static_cast(antPatternArr[index]); // Copy to Eigen matrix
- }
- }
- }
-
- Eigen::MatrixXd gt(2, 3);
- gt(0, 0) = startphi;//x
- gt(0, 1) = deltaphi;
- gt(0, 2) = 0;
-
- gt(1, 0) = starttheta;
- gt(1, 1) = 0;
- gt(1, 2) = deltetheta;
-
- QString antpatternfilepath = getDebugDataPath(antpatternfilename);
- gdalImage ds = CreategdalImageDouble(antpatternfilepath, thetanum, phinum, 1, gt, "", true, true, true);
- ds.saveImage(antPatternMatrix, 0, 0, 1);
-}
-
-
-void MergeTiffs(QList inputFiles, QString outputFile) {
- GDALAllRegister();
-
- if (inputFiles.isEmpty()) {
- fprintf(stderr, "No input files provided.\n");
- return;
- }
-
- // Open the first file to determine the data type and coordinate system
- GDALDataset* poFirstDS = (GDALDataset*)GDALOpen(inputFiles.first().toUtf8().constData(), GA_ReadOnly);
- if (poFirstDS == nullptr) {
- fprintf(stderr, "Failed to open the first file %s\n", inputFiles.first().toUtf8().constData());
- return;
- }
-
- double adfGeoTransform[6];
- CPLErr eErr = poFirstDS->GetGeoTransform(adfGeoTransform);
- if (eErr != CE_None) {
- fprintf(stderr, "Failed to get GeoTransform for the first file %s\n", inputFiles.first().toUtf8().constData());
- GDALClose(poFirstDS);
- return;
- }
-
- int nXSize = 0;
- int nYSize = 0;
- double minX = std::numeric_limits::max();
- double minY = std::numeric_limits::max();
- double maxX = std::numeric_limits::lowest();
- double maxY = std::numeric_limits::lowest();
- double pixelWidth = adfGeoTransform[1];
- double pixelHeight = adfGeoTransform[5];
-
- // Determine the bounding box and size of the output raster
- for (const QString& inputFile : inputFiles) {
- GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile.toUtf8().constData(), GA_ReadOnly);
- if (poSrcDS == nullptr) {
- fprintf(stderr, "Failed to open %s\n", inputFile.toUtf8().constData());
- continue;
- }
-
- double adfThisTransform[6];
- eErr = poSrcDS->GetGeoTransform(adfThisTransform);
- if (eErr != CE_None) {
- fprintf(stderr, "Failed to get GeoTransform for %s\n", inputFile.toUtf8().constData());
- GDALClose(poSrcDS);
- continue;
- }
-
- minX = std::min(minX, adfThisTransform[0]);
- minY = std::min(minY, adfThisTransform[3] + adfThisTransform[5] * poSrcDS->GetRasterYSize());
- maxX = std::max(maxX, adfThisTransform[0] + adfThisTransform[1] * poSrcDS->GetRasterXSize());
- maxY = std::max(maxY, adfThisTransform[3]);
-
- GDALClose(poSrcDS);
- }
-
- nXSize = static_cast(std::ceil((maxX - minX) / pixelWidth));
- nYSize = static_cast(std::ceil((maxY - minY) / (-pixelHeight)));
-
- GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
- if (poDriver == nullptr) {
- fprintf(stderr, "GTiff driver not available.\n");
- GDALClose(poFirstDS);
- return;
- }
-
- char** papszOptions = nullptr;
- GDALDataset* poDstDS = poDriver->Create(outputFile.toUtf8().constData(), nXSize, nYSize, 1, poFirstDS->GetRasterBand(1)->GetRasterDataType(), papszOptions);
- if (poDstDS == nullptr) {
- fprintf(stderr, "Creation of output file failed.\n");
- GDALClose(poFirstDS);
- return;
- }
-
- poDstDS->SetGeoTransform(adfGeoTransform);
-
- const OGRSpatialReference* oSRS = poFirstDS->GetSpatialRef( );
- poDstDS->SetSpatialRef(oSRS);
-
- float fillValue = std::numeric_limits::quiet_NaN();
- void* pafScanline = CPLMalloc(GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * nXSize);
- memset(pafScanline, 0, GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * nXSize);
-
- // Initialize all pixels to NaN
- for (int iY = 0; iY < nYSize; ++iY) {
- GDALRasterBand* poBand = poDstDS->GetRasterBand(1);
- poBand->RasterIO(GF_Write, 0, iY, nXSize, 1, pafScanline, nXSize, 1, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
- }
-
- CPLFree(pafScanline);
-
- // Read each source image and write into the destination image
- for (const QString& inputFile : inputFiles) {
- GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile.toUtf8().constData(), GA_ReadOnly);
- if (poSrcDS == nullptr) {
- fprintf(stderr, "Failed to open %s\n", inputFile.toUtf8().constData());
- continue;
- }
-
- double adfThisTransform[6];
- poSrcDS->GetGeoTransform(adfThisTransform);
-
- int srcXSize = poSrcDS->GetRasterXSize();
- int srcYSize = poSrcDS->GetRasterYSize();
-
- int dstXOffset = static_cast(std::round((adfThisTransform[0] - minX) / pixelWidth));
- int dstYOffset = static_cast(std::round((maxY - adfThisTransform[3]) / (-pixelHeight)));
-
- GDALRasterBand* poSrcBand = poSrcDS->GetRasterBand(1);
- GDALRasterBand* poDstBand = poDstDS->GetRasterBand(1);
-
- void* pafBuffer = CPLMalloc(GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * srcXSize * srcYSize);
- poSrcBand->RasterIO(GF_Read, 0, 0, srcXSize, srcYSize, pafBuffer, srcXSize, srcYSize, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
-
- poDstBand->RasterIO(GF_Write, dstXOffset, dstYOffset, srcXSize, srcYSize, pafBuffer, srcXSize, srcYSize, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
-
- CPLFree(pafBuffer);
- GDALClose(poSrcDS);
- }
-
- GDALClose(poDstDS);
- GDALClose(poFirstDS);
-
-}
void ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode) {
// 注册所有GDAL驱动
@@ -3810,13 +1395,22 @@ void ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long o
// qDebug() << "无效的EPSG代码:" << outepsgcode;
return;
}
-
+ GDALDataType datetype = srcDataset->GetRasterBand(1)->GetRasterDataType();
// 获取目标坐标系的WKT表示
char* targetSRSWkt = nullptr;
targetSRS.exportToWkt(&targetSRSWkt);
+ bool flag = (datetype == GDT_Byte || datetype == GDT_Int8 || datetype == GDT_Int16 ||datetype == GDT_UInt16 || datetype == GDT_Int32 || datetype == GDT_UInt32 || datetype == GDT_Int64 || datetype == GDT_UInt64);
+
// 创建重投影后的虚拟数据集(Warped VRT)
- GDALDataset* warpedVRT = (GDALDataset*)GDALAutoCreateWarpedVRT(
+ GDALDataset* warpedVRT = flag? (GDALDataset*)GDALAutoCreateWarpedVRT(
+ srcDataset,
+ nullptr, // 输入坐标系(默认使用源数据)
+ targetSRSWkt, // 目标坐标系
+ GRA_NearestNeighbour, // 重采样方法:双线性插值
+ 0.0, // 最大误差(0表示自动计算)
+ nullptr // 其他选项
+ ) :(GDALDataset*)GDALAutoCreateWarpedVRT(
srcDataset,
nullptr, // 输入坐标系(默认使用源数据)
targetSRSWkt, // 目标坐标系
@@ -4148,98 +1742,52 @@ void CreateSARIntensityByLookTable(QString IntensityRasterPath,
std::function processBarShow
)
{
- GDALAllRegister();
- constexpr size_t GB4 = size_t(1) * 1024 * 1024 * 1024; // 4GB in bytes
-
- // 打开输入数据集
- GDALDataset* lookDS = (GDALDataset*)GDALOpen(LookTableRasterPath.toStdString().c_str(), GA_ReadOnly);
- GDALDataset* intensityDS = (GDALDataset*)GDALOpen(IntensityRasterPath.toStdString().c_str(), GA_ReadOnly);
-
- // 验证数据集有效性
- if (!lookDS || !intensityDS ||
- lookDS->GetRasterXSize() != intensityDS->GetRasterXSize() ||
- lookDS->GetRasterYSize() != intensityDS->GetRasterYSize())
- {
- if (lookDS) GDALClose(lookDS);
- if (intensityDS) GDALClose(intensityDS);
- return;
- }
-
- // 获取栅格参数
- const int width = lookDS->GetRasterXSize();
- const int height = lookDS->GetRasterYSize();
- const int rows_sar = max_rid - min_rid + 1;
- const int cols_sar = max_cid - min_cid + 1;
-
- // 计算分块策略
- const size_t pixelBytes =
- GDALGetDataTypeSizeBytes(lookDS->GetRasterBand(1)->GetRasterDataType()) * 2 + // 两个波段
- GDALGetDataTypeSizeBytes(intensityDS->GetRasterBand(1)->GetRasterDataType());
-
- const size_t maxPixelsPerBlock = GB4 / pixelBytes;
- int blockXSize = width;
- int blockYSize = static_cast(maxPixelsPerBlock / width);
- blockYSize = std::max(1, std::min(blockYSize, height));
-
- // 输出矩阵初始化
- std::vector sarData(rows_sar * cols_sar, 0.0);
+ gdalImage looktableds(LookTableRasterPath);
+ gdalImage geoIntensity(IntensityRasterPath);
+ gdalImage SARIntensity= CreategdalImageDouble(SARIntensityPath,max_rid-min_rid,max_cid-min_cid,1);
+ long blockYSize = Memory1GB / looktableds.width / 8 * 2;
+
+ Eigen::MatrixXd SARData = SARIntensity.getData(0, 0, SARIntensity.height, SARIntensity.width, 1);
+ SARData = SARData.array() * 0;
// 分块处理
- for (int yOff = 0; yOff < height; yOff += blockYSize)
+ for (int yOff = 0; yOff < looktableds.height; yOff += blockYSize)
{
- const int ySize = std::min(blockYSize, height - yOff);
-
- // 读取行号列号数据
- std::vector rData(width * ySize);
- std::vector cData(width * ySize);
- lookDS->GetRasterBand(1)->RasterIO(GF_Read, 0, yOff, width, ySize,
- rData.data(), width, ySize, GDT_Int32, 0, 0);
- lookDS->GetRasterBand(2)->RasterIO(GF_Read, 0, yOff, width, ySize,
- cData.data(), width, ySize, GDT_Int32, 0, 0);
+ processBarShow(yOff, looktableds.height);
+ qDebug() << "Process : [" << yOff * 100.0 / looktableds.height << " % ]";
+ Eigen::MatrixXd rowData = looktableds.getData(yOff, 0, blockYSize, looktableds.width, 1);
+ Eigen::MatrixXd colData = looktableds.getData(yOff, 0, blockYSize, looktableds.width, 2);
+ Eigen::MatrixXd geoData = geoIntensity.getData(yOff, 0, blockYSize, looktableds.width, 1);
- // 读取强度数据
- std::vector intensity(width * ySize);
- intensityDS->GetRasterBand(1)->RasterIO(GF_Read, 0, yOff, width, ySize,
- intensity.data(), width, ySize, GDT_Float32, 0, 0);
+ for (long i = 0; i < rowData.rows(); i++) {
+ for (long j = 0; j < rowData.cols(); j++) {
+ long r =round( rowData(i,j))-min_rid;
+ long c = round(colData(i, j))-min_cid;
- // 处理当前块
-#pragma omp parallel for collapse(2)
- for (int y = 0; y < ySize; ++y) {
- for (int x = 0; x < width; ++x) {
- const int idx = y * width + x;
- const long r = rData[idx];
- const long c = cData[idx];
-
- if (r >= min_rid && r <= max_rid && c >= min_cid && c <= max_cid) {
- const int row = r - min_rid;
- const int col = c - min_cid;
-#pragma omp atomic
- sarData[row * cols_sar + col] += intensity[idx];
+ if (r >= 0 && r < SARIntensity.height && c >= 0 && c < SARIntensity.width) {
+ SARData(r, c) = SARData(r, c) + geoData(i, j);
}
}
}
- if (processBarShow) {
- processBarShow(yOff, height);
- }
-
}
+ SARIntensity.saveImage(SARData, 0, 0, 1);
+ qDebug() << "Process : [ 100 % ]";
+ processBarShow(1000,1000);
+}
- // 写入输出
- GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
- GDALDataset* outDS = driver->Create(SARIntensityPath.toStdString().c_str(),
- cols_sar, rows_sar, 1, GDT_Float64, nullptr);
- outDS->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, cols_sar, rows_sar,
- sarData.data(), cols_sar, rows_sar, GDT_Float64, 0, 0);
- // 资源清理
- GDALClose(lookDS);
- GDALClose(intensityDS);
- GDALClose(outDS);
+bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath)
+{
+
+ Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
+
+ gdalImage img = CreategdalImageDouble(dataStrPath, data.rows(), data.cols(), 1, gt, "", false, true, true);
+
+ img.saveImage(data, 0, 0, 1);
+
+ return true;
}
-
-
-
diff --git a/BaseCommonLibrary/BaseTool/ImageOperatorBase.h b/BaseCommonLibrary/BaseTool/ImageOperatorBase.h
index b0f33ab..b16127d 100644
--- a/BaseCommonLibrary/BaseTool/ImageOperatorBase.h
+++ b/BaseCommonLibrary/BaseTool/ImageOperatorBase.h
@@ -28,6 +28,9 @@
#include
#include
#include
+#include "ShowProessAbstract.h"
+
+
enum ProjectStripDelta {
Strip_6, // 6度带
@@ -76,13 +79,6 @@ enum GDALREADARRCOPYMETHOD {
-class BASECONSTVARIABLEAPI ShowProessAbstract {
-
-public:
- virtual void showProcess(double precent, QString tip);
- virtual void showToolInfo(QString tip);
-};
-
@@ -93,7 +89,7 @@ public:
/// \param long 经度
/// \param lat 纬度
/// \return 对应投影坐标系统的 EPSG编码,-1 表示计算错误
-long BASECONSTVARIABLEAPI getProjectEPSGCodeByLon_Lat(double long, double lat, ProjectStripDelta stripState);
+long BASECONSTVARIABLEAPI getProjectEPSGCodeByLon_Lat(double lon, double lat, ProjectStripDelta stripState= ProjectStripDelta::Strip_6);
long BASECONSTVARIABLEAPI getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat);
@@ -139,10 +135,10 @@ GDALDataType BASECONSTVARIABLEAPI getGDALDataType(QString fileptah);
struct RasterExtend {
- double min_x; //纬度
- double min_y;//经度
- double max_x;//纬度
- double max_y;//经度
+ double min_x;
+ double min_y;
+ double max_x;
+ double max_y;
};
@@ -178,8 +174,6 @@ public: // 方法
virtual void saveImage(std::shared_ptr, int start_row, int start_col, int rowcount, int colcount, int band_ids);
virtual void saveImage(std::shared_ptr, int start_row, int start_col, int rowcount, int colcount, int band_ids);
-
-
virtual void saveImage();
virtual void setNoDataValue(double nodatavalue, int band_ids);
virtual void setNoDataValuei(int nodatavalue, int band_ids);
@@ -201,6 +195,10 @@ public: // 方法
virtual RasterExtend getExtend();
+
+
+
+
public:
QString img_path; // 图像文件
int height; // 高
@@ -234,18 +232,22 @@ public: // 方法
Eigen::MatrixXcd getDataComplex(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
std::shared_ptr> getDataComplexSharePtr(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
- void saveImage() override;
+ void saveComplexImage();//override;
void savePreViewImage();
public:
Eigen::MatrixXcd data;
};
+
+bool BASECONSTVARIABLEAPI CopyProjectTransformMatrixFromRasterAToRasterB(QString RasterAPath, QString RasterBPath);
+
+
// 创建影像
gdalImage BASECONSTVARIABLEAPI CreategdalImageDouble(QString& img_path, int height, int width, int band_num, bool overwrite = false, bool isEnvi = false);
-gdalImage BASECONSTVARIABLEAPI CreategdalImageDouble(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 BASECONSTVARIABLEAPI 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, GDALDataType datetype = GDT_Float32);
-
-gdalImage BASECONSTVARIABLEAPI 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);
+gdalImage BASECONSTVARIABLEAPI CreategdalImageFloat(QString& img_path, int height, int width, int band_num, bool overwrite = false, bool isEnvi = false);
+gdalImage BASECONSTVARIABLEAPI CreategdalImageDouble(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 BASECONSTVARIABLEAPI 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, GDALDataType datetype = GDT_Float32);
+gdalImage BASECONSTVARIABLEAPI 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 BASECONSTVARIABLEAPI 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 BASECONSTVARIABLEAPI CreategdalImageComplexNoProj(const QString& img_path, int height, int width, int band_num, bool overwrite = true);
@@ -314,8 +316,29 @@ void BASECONSTVARIABLEAPI testOutDataArr(QString filename, long* data, long row
void BASECONSTVARIABLEAPI CreateSARIntensityByLookTable(QString IntensityRasterPath, QString LookTableRasterPath, QString SARIntensityPath, long min_rid, long max_rid, long min_cid, long max_cid, std::function processBarShow = {});
+bool BASECONSTVARIABLEAPI ConvertVrtToEnvi(QString vrtPath, QString outPath);
+
+
+
+
+void BASECONSTVARIABLEAPI MultiLookRaster(QString inRasterPath, QString outRasterPath, long looklineNumrow, long looklineNumCol);
+ErrorCode BASECONSTVARIABLEAPI Complex2PhaseRaster(QString inComplexPath, QString outRasterPath);
+ErrorCode BASECONSTVARIABLEAPI Complex2dBRaster(QString inComplexPath, QString outRasterPath);
+ErrorCode BASECONSTVARIABLEAPI Complex2AmpRaster(QString inComplexPath, QString outRasterPath);
+ErrorCode BASECONSTVARIABLEAPI amp2dBRaster(QString inPath, QString outRasterPath);
+
+
+ErrorCode BASECONSTVARIABLEAPI ResampleDEM(QString indemPath, QString outdemPath, double gridx, double gridy);
+
+
+void BASECONSTVARIABLEAPI CloseAllGDALRaster();
+
+
//--------------------- 图像文件读写 ------------------------------
+
+
+
template
inline std::shared_ptr readDataArr(gdalImage& imgds, long start_row, long start_col, long& rows_count, long& cols_count, int band_ids, GDALREADARRCOPYMETHOD method)
{
@@ -513,7 +536,7 @@ inline std::shared_ptr readDataArr(gdalImage& imgds, long start_row, long sta
omp_destroy_lock(&lock); // 劫伙拷斤拷
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return result;
-}
+};
template
inline std::shared_ptr readDataArrComplex(gdalImageComplex& imgds, long start_row, long start_col, long& rows_count, long& cols_count, int band_ids, GDALREADARRCOPYMETHOD method)
@@ -586,7 +609,7 @@ inline std::shared_ptr readDataArrComplex(gdalImageComplex& imgds, long start
omp_destroy_lock(&lock); // 劫伙拷斤拷
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return result;
-}
+};
diff --git a/BaseCommonLibrary/BaseTool/MergeRasterOperator.cpp b/BaseCommonLibrary/BaseTool/MergeRasterOperator.cpp
new file mode 100644
index 0000000..f8144e7
--- /dev/null
+++ b/BaseCommonLibrary/BaseTool/MergeRasterOperator.cpp
@@ -0,0 +1,488 @@
+#include "stdafx.h"
+#include "ImageOperatorBase.h"
+#include "BaseTool.h"
+#include "GeoOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include "FileOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include // OGRSpatialReference ڿռοת
+#include // GDALWarp
+
+ErrorCode MergeRasterProcess(QVector filepaths, QString outfileptah, QString mainString, MERGEMODE mergecode, bool isENVI, ShowProessAbstract* dia)
+{
+ //
+ if (!isExists(mainString)) {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND)) << "\t" << mainString;
+ return ErrorCode::FILENOFOUND;
+ }
+ else {}
+ gdalImage mainimg(mainString);
+ QVector imgdslist(filepaths.count());
+ for (long i = 0; i < filepaths.count(); i++) {
+ if (!isExists(filepaths[i])) {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND)) << "\t" << filepaths[i];
+ return ErrorCode::FILENOFOUND;
+ }
+ else {
+ imgdslist[i] = gdalImage(filepaths[i]);
+ if (imgdslist[i].band_num != mainimg.band_num) {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTERBAND_NOTEQUAL)) << "\t" << imgdslist[i].band_num << " != " << mainimg.band_num;
+ return ErrorCode::RASTERBAND_NOTEQUAL;
+ }
+ }
+ }
+
+ // ϵǷͳһ
+ long EPSGCode = GetEPSGFromRasterFile(mainString);
+ long tempCode = 0;
+ for (long i = 0; i < filepaths.count(); i++) {
+ tempCode = GetEPSGFromRasterFile(filepaths[i]);
+ if (EPSGCode != tempCode) {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSAME)) << "\t" << EPSGCode << "!=" << tempCode;
+ return ErrorCode::EPSGCODE_NOTSAME;
+ }
+ }
+
+ // ӰǷͳһ
+ GDALDataType mainType = mainimg.getDataType();
+ for (long i = 0; i < imgdslist.count(); i++) {
+ if (mainType != imgdslist[i].getDataType()) {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTER_DATETYPE_NOTSAME)) << "\t" << mainType << "!=" << imgdslist[i].getDataType();
+ return ErrorCode::RASTER_DATETYPE_NOTSAME;
+ }
+ }
+
+ Eigen::MatrixXd maingt = mainimg.getGeoTranslation();
+ Eigen::MatrixXd rgt = Eigen::MatrixXd::Zero(2, 3);
+ RasterExtend mainExtend = mainimg.getExtend();
+ rgt(0, 1) = (mainExtend.max_x - mainExtend.min_x) / (mainimg.width - 1); //dx
+ rgt(1, 2) = -1 * std::abs(((mainExtend.max_y - mainExtend.min_y) / (mainimg.height - 1)));//dy
+ QVector extendlist(imgdslist.count());
+ for (long i = 0; i < imgdslist.count(); i++) {
+ extendlist[i] = imgdslist[i].getExtend();
+ mainExtend.min_x = mainExtend.min_x < extendlist[i].min_x ? mainExtend.min_x : extendlist[i].min_x;
+ mainExtend.max_x = mainExtend.max_x > extendlist[i].max_x ? mainExtend.max_x : extendlist[i].max_x;
+ mainExtend.min_y = mainExtend.min_y < extendlist[i].min_y ? mainExtend.min_y : extendlist[i].min_y;
+ mainExtend.max_y = mainExtend.max_y > extendlist[i].max_y ? mainExtend.max_y : extendlist[i].max_y;
+ }
+
+ rgt(0, 0) = mainExtend.min_x;
+ rgt(1, 0) = mainExtend.max_y;
+
+ //
+
+ long width = std::ceil((mainExtend.max_x - mainExtend.min_x) / rgt(0, 1) + 1);
+ long height = std::ceil(std::abs((mainExtend.min_y - mainExtend.max_y) / rgt(1, 2)) + 1);
+ OGRSpatialReference oSRS;
+ if (oSRS.importFromEPSG(EPSGCode) != OGRERR_NONE) {
+ qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSUPPORT)) << "\t" << EPSGCode;
+ return ErrorCode::EPSGCODE_NOTSUPPORT;
+ }
+
+ gdalImage resultImage = CreategdalImage(outfileptah, height, width, mainimg.band_num, rgt, EPSGCode, mainType, true, true, isENVI);
+
+ QString resultMaskString = addMaskToFileName(outfileptah, QString("_MASK"));
+ gdalImage maskImage = CreategdalImage(resultMaskString, height, width, 1, rgt, EPSGCode, GDT_Int32, true, true, isENVI);
+
+ // ʼ
+ long resultline = Memory1MB * 500 / 8 / resultImage.width;
+ resultline = resultline < 10000 ? resultline : 10000; // 100
+ resultline = resultline > 0 ? resultline : 2;
+ long bandnum = resultImage.band_num + 1;
+ long starti = 0;
+ long rasterCount = imgdslist.count();
+
+
+
+ QProgressDialog progressDialog(u8"ʼӰ", u8"ֹ", 0, resultImage.height);
+ progressDialog.setWindowTitle(u8"ʼӰ");
+ progressDialog.setWindowModality(Qt::WindowModal);
+ progressDialog.setAutoClose(true);
+ progressDialog.setValue(0);
+ progressDialog.setMaximum(resultImage.height);
+ progressDialog.setMinimum(0);
+ progressDialog.show();
+
+
+ for (starti = 0; starti < resultImage.height; starti = starti + resultline) {
+ long blocklines = resultline;
+ blocklines = starti + blocklines < resultImage.height ? blocklines : resultImage.height - starti;
+ for (long b = 1; b < bandnum; b++) {
+ Eigen::MatrixXd data = resultImage.getData(starti, 0, blocklines, resultImage.width, b);
+ Eigen::MatrixXi maskdata = maskImage.getDatai(starti, 0, blocklines, resultImage.width, b);
+ data = data.array() * 0;
+ maskdata = maskdata.array() * 0;
+ resultImage.saveImage(data, starti, 0, b);
+ maskImage.saveImage(maskdata, starti, 0, b);
+ }
+ if (nullptr != dia) {
+ dia->showProcess(starti * 1.0 / resultImage.height, u8"ʼӰ");
+ }
+ progressDialog.setValue(starti + blocklines);
+ }
+ progressDialog.close();
+
+
+
+ switch (mergecode)
+ {
+ case MERGE_GEOCODING:
+ return MergeRasterInGeoCoding(imgdslist, resultImage, maskImage, dia);
+ default:
+ break;
+ }
+
+
+ return ErrorCode::SUCCESS;
+}
+
+ErrorCode MergeRasterInGeoCoding(QVector imgdslist, gdalImage resultimg, gdalImage maskimg, ShowProessAbstract* dia)
+{
+ omp_set_num_threads(Paral_num_thread);
+ // ϲ
+ QVector extendlist(imgdslist.count());
+ for (long i = 0; i < imgdslist.count(); i++) {
+ extendlist[i] = imgdslist[i].getExtend();
+ imgdslist[i].InitInv_gt();
+ }
+
+ // ֿ
+ long resultline = Memory1MB * 1000 / 8 / resultimg.width;
+ resultline = resultline < 300 ? resultline : 300; // 100
+
+ long bandnum = resultimg.band_num + 1;
+ long starti = 0;
+ long rasterCount = imgdslist.count();
+
+ long processNumber = 0;
+ QProgressDialog progressDialog(u8"ϲӰ", u8"ֹ", 0, resultimg.height);
+ progressDialog.setWindowTitle(u8"ϲӰ");
+ progressDialog.setWindowModality(Qt::WindowModal);
+ progressDialog.setAutoClose(true);
+ progressDialog.setValue(0);
+ progressDialog.setMaximum(resultimg.height);
+ progressDialog.setMinimum(0);
+ progressDialog.show();
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+
+#pragma omp parallel for
+ for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
+ long blocklines = resultline;
+ blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
+
+ long rid = starti;
+ long cid = 0;
+
+ Landpoint pp = { 0,0,0 };
+ Landpoint lp = { 0,0,0 };
+
+ for (long ir = 0; ir < rasterCount; ir++) {// Ӱ
+ long minRid = imgdslist[ir].height;
+ long maxRid = 0;
+
+ Eigen::MatrixXd ridlist = resultimg.getData(starti, 0, blocklines, resultimg.width, 1);
+ ridlist = ridlist.array() * 0;
+ Eigen::MatrixXd cidlist = ridlist.array() * 0;
+
+ for (long i = 0; i < blocklines; i++) {// к
+ rid = starti + i;
+ for (long j = 0; j < resultimg.width; j++) {// к
+ cid = j;
+ resultimg.getLandPoint(rid, cid, 0, pp);
+ lp = imgdslist[ir].getRow_Col(pp.lon, pp.lat); // ȡ
+ ridlist(i, j) = lp.lat;
+ cidlist(i, j) = lp.lon;
+ }
+ }
+
+ //ImageShowDialogClass* dialog = new ImageShowDialogClass;
+ //dialog->show();
+ //dialog->load_double_MatrixX_data(cidlist, u8"");
+
+ //dialog->exec();
+
+
+ if (ridlist.maxCoeff() < 0 || ridlist.minCoeff() >= imgdslist[ir].height) {
+ continue;
+ }
+
+ if (cidlist.maxCoeff() < 0 || cidlist.minCoeff() >= imgdslist[ir].width) {
+ continue;
+ }
+
+ minRid = std::floor(ridlist.minCoeff());
+ maxRid = std::ceil(ridlist.maxCoeff());
+ minRid = minRid < 0 ? 0 : minRid;
+ maxRid = maxRid < imgdslist[ir].height ? maxRid : imgdslist[ir].height - 1;
+
+ long rowlen = maxRid - minRid + 1;
+ if (rowlen <= 0) {
+ continue;
+ }
+ // ȡ
+ Landpoint p0{ 0,0,0 }, p11{ 0,0,0 }, p21{ 0,0,0 }, p12{ 0,0,0 }, p22{ 0,0,0 }, p{ 0,0,0 };
+
+
+ long rowcount = 0;
+ long colcount = 0;
+ double ridtemp = 0, cidtemp = 0;
+
+ long lastr = 0, nextr = 0;
+ long lastc = 0, nextc = 0;
+
+ double r0 = 0, c0 = 0;
+
+ for (long b = 1; b < bandnum; b++) {
+ Eigen::MatrixXd resultdata = resultimg.getData(starti, 0, blocklines, resultimg.width, b);
+ Eigen::MatrixXi resultmask = maskimg.getDatai(starti, 0, blocklines, resultimg.width, b);
+ Eigen::MatrixXd data = imgdslist[ir].getData(minRid, 0, rowlen, imgdslist[ir].width, b);
+
+ double nodata = imgdslist[ir].getNoDataValue(b);
+ for (long ii = 0; ii < data.rows(); ii++) {
+ for (long jj = 0; jj < data.cols(); jj++) {
+ if (std::abs(data(ii, jj) - nodata) < 1e-6) {
+ data(ii, jj) = 0;
+ }
+ }
+ }
+ rowcount = ridlist.rows();
+ colcount = ridlist.cols();
+ double Bileanervalue = 0;
+ for (long i = 0; i < rowcount; i++) {
+ for (long j = 0; j < colcount; j++) {
+ ridtemp = ridlist(i, j);
+ cidtemp = cidlist(i, j);
+
+ lastr = std::floor(ridtemp);
+ nextr = std::ceil(ridtemp);
+ lastc = std::floor(cidtemp);
+ nextc = std::ceil(cidtemp);
+
+ if (lastr < 0 || lastr >= imgdslist[ir].height
+ || nextr < 0 || nextr >= imgdslist[ir].height
+ || lastc < 0 || lastc >= imgdslist[ir].width
+ || nextc < 0 || nextc >= imgdslist[ir].width) {
+ continue;
+ }
+ else {}
+
+ r0 = ridtemp - std::floor(ridtemp);
+ c0 = cidtemp - std::floor(cidtemp);
+
+ lastr = lastr - minRid;
+ nextr = nextr - minRid;
+
+ p0 = Landpoint{ c0,r0,0 };
+ p11 = Landpoint{ 0,0,data(lastr,lastc) };
+ p21 = Landpoint{ 0,1,data(nextr,lastc) };
+ p12 = Landpoint{ 1,0,data(lastr,nextc) };
+ p22 = Landpoint{ 1,1,data(nextr,nextc) };
+ Bileanervalue = Bilinear_interpolation(p0, p11, p21, p12, p22);
+ if (std::abs(Bileanervalue) < 1e-6 || resultmask(i, j) > 0) {
+ continue;
+ }
+ resultdata(i, j) = resultdata(i, j) + Bileanervalue;
+ resultmask(i, j) = resultmask(i, j) + 1;
+ }
+ }
+ resultimg.saveImage(resultdata, starti, 0, b);
+ maskimg.saveImage(resultmask, starti, 0, b);
+ }
+ }
+
+ omp_set_lock(&lock);
+ processNumber = processNumber + blocklines;
+ qDebug() << "\rprocess bar:\t" << processNumber * 100.0 / resultimg.height << " % " << "\t\t\t";
+ if (nullptr != dia) {
+ dia->showProcess(processNumber * 1.0 / resultimg.height, u8"ϲͼ");
+ }
+ if (progressDialog.maximum() <= processNumber) {
+ processNumber = progressDialog.maximum() - 1;
+ }
+ progressDialog.setValue(processNumber);
+ omp_unset_lock(&lock);
+ }
+ omp_destroy_lock(&lock);
+
+ progressDialog.setWindowTitle(u8"ӰĤ");
+ progressDialog.setLabelText(u8"ӰĤ");
+ for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
+ long blocklines = resultline;
+ blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
+ for (long b = 1; b < bandnum; b++) {
+ Eigen::MatrixXd data = resultimg.getData(starti, 0, blocklines, resultimg.width, b);
+ Eigen::MatrixXd maskdata = maskimg.getData(starti, 0, blocklines, maskimg.width, b);
+
+ for (long i = 0; i < data.rows(); i++) {
+ for (long j = 0; j < data.cols(); j++) {
+ if (maskdata(i, j) == 0) {
+ data(i, j) = -9999;
+ continue;
+ }
+ data(i, j) = data(i, j) / maskdata(i, j);
+ }
+ }
+
+ resultimg.saveImage(data, starti, 0, b);
+ maskimg.saveImage(maskdata, starti, 0, b);
+ }
+ if (nullptr != dia) {
+ dia->showProcess((starti + blocklines) * 1.0 / resultimg.height, u8"ӰĤ");
+ }
+ progressDialog.setValue(starti + blocklines);
+ }
+
+
+ progressDialog.close();
+ return ErrorCode::SUCCESS;
+}
+
+
+void MergeTiffs(QList inputFiles, QString outputFile) {
+ GDALAllRegister();
+
+ if (inputFiles.isEmpty()) {
+ fprintf(stderr, "No input files provided.\n");
+ return;
+ }
+
+ // Open the first file to determine the data type and coordinate system
+ GDALDataset* poFirstDS = (GDALDataset*)GDALOpen(inputFiles.first().toUtf8().constData(), GA_ReadOnly);
+ if (poFirstDS == nullptr) {
+ fprintf(stderr, "Failed to open the first file %s\n", inputFiles.first().toUtf8().constData());
+ return;
+ }
+
+ double adfGeoTransform[6];
+ CPLErr eErr = poFirstDS->GetGeoTransform(adfGeoTransform);
+ if (eErr != CE_None) {
+ fprintf(stderr, "Failed to get GeoTransform for the first file %s\n", inputFiles.first().toUtf8().constData());
+ GDALClose(poFirstDS);
+ return;
+ }
+
+ int nXSize = 0;
+ int nYSize = 0;
+ double minX = std::numeric_limits::max();
+ double minY = std::numeric_limits::max();
+ double maxX = std::numeric_limits::lowest();
+ double maxY = std::numeric_limits::lowest();
+ double pixelWidth = adfGeoTransform[1];
+ double pixelHeight = adfGeoTransform[5];
+
+ // Determine the bounding box and size of the output raster
+ for (const QString& inputFile : inputFiles) {
+ GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile.toUtf8().constData(), GA_ReadOnly);
+ if (poSrcDS == nullptr) {
+ fprintf(stderr, "Failed to open %s\n", inputFile.toUtf8().constData());
+ continue;
+ }
+
+ double adfThisTransform[6];
+ eErr = poSrcDS->GetGeoTransform(adfThisTransform);
+ if (eErr != CE_None) {
+ fprintf(stderr, "Failed to get GeoTransform for %s\n", inputFile.toUtf8().constData());
+ GDALClose(poSrcDS);
+ continue;
+ }
+
+ minX = std::min(minX, adfThisTransform[0]);
+ minY = std::min(minY, adfThisTransform[3] + adfThisTransform[5] * poSrcDS->GetRasterYSize());
+ maxX = std::max(maxX, adfThisTransform[0] + adfThisTransform[1] * poSrcDS->GetRasterXSize());
+ maxY = std::max(maxY, adfThisTransform[3]);
+
+ GDALClose(poSrcDS);
+ }
+
+ nXSize = static_cast(std::ceil((maxX - minX) / pixelWidth));
+ nYSize = static_cast(std::ceil((maxY - minY) / (-pixelHeight)));
+
+ GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ if (poDriver == nullptr) {
+ fprintf(stderr, "GTiff driver not available.\n");
+ GDALClose(poFirstDS);
+ return;
+ }
+
+ char** papszOptions = nullptr;
+ GDALDataset* poDstDS = poDriver->Create(outputFile.toUtf8().constData(), nXSize, nYSize, 1, poFirstDS->GetRasterBand(1)->GetRasterDataType(), papszOptions);
+ if (poDstDS == nullptr) {
+ fprintf(stderr, "Creation of output file failed.\n");
+ GDALClose(poFirstDS);
+ return;
+ }
+
+ poDstDS->SetGeoTransform(adfGeoTransform);
+
+ const OGRSpatialReference* oSRS = poFirstDS->GetSpatialRef();
+ poDstDS->SetSpatialRef(oSRS);
+
+ float fillValue = std::numeric_limits::quiet_NaN();
+ void* pafScanline = CPLMalloc(GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * nXSize);
+ memset(pafScanline, 0, GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * nXSize);
+
+ // Initialize all pixels to NaN
+ for (int iY = 0; iY < nYSize; ++iY) {
+ GDALRasterBand* poBand = poDstDS->GetRasterBand(1);
+ poBand->RasterIO(GF_Write, 0, iY, nXSize, 1, pafScanline, nXSize, 1, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
+ }
+
+ CPLFree(pafScanline);
+
+ // Read each source image and write into the destination image
+ for (const QString& inputFile : inputFiles) {
+ GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile.toUtf8().constData(), GA_ReadOnly);
+ if (poSrcDS == nullptr) {
+ fprintf(stderr, "Failed to open %s\n", inputFile.toUtf8().constData());
+ continue;
+ }
+
+ double adfThisTransform[6];
+ poSrcDS->GetGeoTransform(adfThisTransform);
+
+ int srcXSize = poSrcDS->GetRasterXSize();
+ int srcYSize = poSrcDS->GetRasterYSize();
+
+ int dstXOffset = static_cast(std::round((adfThisTransform[0] - minX) / pixelWidth));
+ int dstYOffset = static_cast(std::round((maxY - adfThisTransform[3]) / (-pixelHeight)));
+
+ GDALRasterBand* poSrcBand = poSrcDS->GetRasterBand(1);
+ GDALRasterBand* poDstBand = poDstDS->GetRasterBand(1);
+
+ void* pafBuffer = CPLMalloc(GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * srcXSize * srcYSize);
+ poSrcBand->RasterIO(GF_Read, 0, 0, srcXSize, srcYSize, pafBuffer, srcXSize, srcYSize, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
+
+ poDstBand->RasterIO(GF_Write, dstXOffset, dstYOffset, srcXSize, srcYSize, pafBuffer, srcXSize, srcYSize, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
+
+ CPLFree(pafBuffer);
+ GDALClose(poSrcDS);
+ }
+
+ GDALClose(poDstDS);
+ GDALClose(poFirstDS);
+
+}
+
+
+
diff --git a/BaseCommonLibrary/BaseTool/RasterToolBase.cpp b/BaseCommonLibrary/BaseTool/RasterToolBase.cpp
index 82a50be..5debe6f 100644
--- a/BaseCommonLibrary/BaseTool/RasterToolBase.cpp
+++ b/BaseCommonLibrary/BaseTool/RasterToolBase.cpp
@@ -15,6 +15,8 @@
#include
#include
#include
+#include "SARSimulationImageL1.h"
+#include
namespace RasterToolBase {
long getProjectEPSGCodeByLon_Lat(double lon, double lat, ProjectStripDelta stripState)
@@ -261,8 +263,14 @@ namespace RasterToolBase {
} else if(oSRS.IsProjected()) {
return CoordinateSystemType::ProjectCoordinateSystem;
}
+ else {
+ return CoordinateSystemType::UNKNOW;
+ }
} else {
return CoordinateSystemType::UNKNOW;
}
}
-} // namespace RasterToolBase
\ No newline at end of file
+}; // namespace RasterToolBase
+
+
+
diff --git a/BaseCommonLibrary/BaseTool/RasterToolBase.h b/BaseCommonLibrary/BaseTool/RasterToolBase.h
index bc6842f..254d4bb 100644
--- a/BaseCommonLibrary/BaseTool/RasterToolBase.h
+++ b/BaseCommonLibrary/BaseTool/RasterToolBase.h
@@ -12,40 +12,42 @@
#include "BaseConstVariable.h"
#include "gdal_priv.h"
#include
+#include "LogInfoCls.h"
+
namespace RasterToolBase {
- static bool GDALAllRegisterEnable=false;
+ static bool GDALAllRegisterEnable = false;
- enum ProjectStripDelta{
+ enum ProjectStripDelta {
Strip_6, // 6度带
Strip_3
};
- enum CoordinateSystemType{ // 坐标系类型
+ enum CoordinateSystemType { // 坐标系类型
GeoCoordinateSystem,
ProjectCoordinateSystem,
UNKNOW
};
- struct PointRaster{ // 影像坐标点
+ struct PointRaster { // 影像坐标点
double x;
double y;
double z;
};
- struct PointXYZ{
- double x,y,z;
+ struct PointXYZ {
+ double x, y, z;
};
- struct PointGeo{
- double lon,lat,ati;
+ struct PointGeo {
+ double lon, lat, ati;
};
- struct PointImage{
- double pixel_x,pixel_y;
+ struct PointImage {
+ double pixel_x, pixel_y;
};
/// 根据经纬度获取
@@ -56,14 +58,14 @@ namespace RasterToolBase {
/// \param lat 纬度
/// \return 对应投影坐标系统的 EPSG编码,-1 表示计算错误
long BASECONSTVARIABLEAPI getProjectEPSGCodeByLon_Lat(double long, double lat,
- ProjectStripDelta stripState = ProjectStripDelta::Strip_3);
+ ProjectStripDelta stripState = ProjectStripDelta::Strip_3);
long BASECONSTVARIABLEAPI getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat);
long BASECONSTVARIABLEAPI getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat);
- QString BASECONSTVARIABLEAPI GetProjectionNameFromEPSG(long epsgCode) ;
+ QString BASECONSTVARIABLEAPI GetProjectionNameFromEPSG(long epsgCode);
long BASECONSTVARIABLEAPI GetEPSGFromRasterFile(QString filepath);
@@ -72,9 +74,23 @@ namespace RasterToolBase {
CoordinateSystemType BASECONSTVARIABLEAPI getCoordinateSystemTypeByEPSGCode(long EPSGCODE);
+};// namespace RasterProjectConvertor
+
+
+
+
+// 遥感类常用数据
+
+
+
+
+
+
+
+
+
-} // namespace RasterProjectConvertor
#endif // LAMPCAE_RASTERTOOLBASE_H
diff --git a/BaseCommonLibrary/BaseTool/SARSimulationImageL1.cpp b/BaseCommonLibrary/BaseTool/SARSimulationImageL1.cpp
index 16134d2..2fbfbb8 100644
--- a/BaseCommonLibrary/BaseTool/SARSimulationImageL1.cpp
+++ b/BaseCommonLibrary/BaseTool/SARSimulationImageL1.cpp
@@ -106,7 +106,7 @@ ErrorCode SARSimulationImageL1Dataset::OpenOrNew(QString folder, QString filenam
this->saveToXml();
}
- if (this->Rasterlevel!=RasterL2||QFile(this->GPSPointFilePath).exists() == false) {
+ if (this->Rasterlevel==RasterL2||QFile(this->GPSPointFilePath).exists() == false) {
// ļ
omp_lock_t lock;
omp_init_lock(&lock);
@@ -122,7 +122,7 @@ ErrorCode SARSimulationImageL1Dataset::OpenOrNew(QString folder, QString filenam
omp_destroy_lock(&lock);
}
- if (this->Rasterlevel == RasterLevel::RasterSLC || QFile(this->ImageRasterPath).exists() == false) {
+ else if (this->Rasterlevel == RasterLevel::RasterSLC || QFile(this->ImageRasterPath).exists() == false) {
// ļ
omp_lock_t lock;
@@ -140,7 +140,7 @@ ErrorCode SARSimulationImageL1Dataset::OpenOrNew(QString folder, QString filenam
}
- if (this->Rasterlevel != RasterLevel::RasterSLC || QFile(this->ImageRasterPath).exists() == false) {
+ else if (this->Rasterlevel == RasterLevel::RasterL1B || QFile(this->ImageRasterPath).exists() == false) {
// ļ
omp_lock_t lock;
@@ -150,11 +150,11 @@ ErrorCode SARSimulationImageL1Dataset::OpenOrNew(QString folder, QString filenam
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
- std::shared_ptr poDstDS(poDriver->Create(this->ImageRasterPath.toUtf8().constData(), colCount, rowCount, 1, GDT_CFloat32, NULL));
+ std::shared_ptr poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 19, rowCount, 1, GDT_Float64, NULL));
GDALFlushCache((GDALDatasetH)poDstDS.get());
poDstDS.reset();
- omp_unset_lock(&lock); //
- omp_destroy_lock(&lock); //
+ omp_unset_lock(&lock);
+ omp_destroy_lock(&lock);
}
@@ -203,7 +203,7 @@ ErrorCode SARSimulationImageL1Dataset::Open(QString xmlPath)
QFileInfo fileInfo(xmlPath);
QString fileName = fileInfo.fileName(); // ȡļ
QString fileSuffix = fileInfo.suffix(); // ȡ
- QString fileNameWithoutSuffix = fileInfo.baseName(); // ȡļ
+ QString fileNameWithoutSuffix = fileInfo.completeBaseName(); // ȡļ
QString directoryPath = fileInfo.path(); // ȡļ·
if (fileSuffix.toLower() == "xml" || fileSuffix.toLower() == ".xml") {
return this->Open(directoryPath, fileNameWithoutSuffix);
@@ -246,12 +246,13 @@ void SARSimulationImageL1Dataset::saveToXml()
xmlWriter.writeTextElement("RowCount", QString::number(this->rowCount));
xmlWriter.writeTextElement("ColCount", QString::number(this->colCount));
- xmlWriter.writeTextElement("Rnear", QString::number(this->Rnear));
- xmlWriter.writeTextElement("Rfar", QString::number(this->Rfar));
- xmlWriter.writeTextElement("Rref", QString::number(this->Rref));
- xmlWriter.writeTextElement("CenterFreq", QString::number(this->centerFreq));
- xmlWriter.writeTextElement("Fs", QString::number(this->Fs));
- xmlWriter.writeTextElement("CenterAngle", QString::number(this->CenterAngle));
+ xmlWriter.writeTextElement("Rnear", QString::number(this->Rnear, 'e', 18));
+ xmlWriter.writeTextElement("Rfar", QString::number(this->Rfar, 'e', 18));
+ xmlWriter.writeTextElement("Rref", QString::number(this->Rref, 'e', 18));
+ xmlWriter.writeTextElement("CenterFreq", QString::number(this->centerFreq, 'e', 18));
+ xmlWriter.writeTextElement("Fs", QString::number(this->Fs, 'e', 18));
+ xmlWriter.writeTextElement("PRF", QString::number(this->prf, 'e', 18));
+ xmlWriter.writeTextElement("CenterAngle", QString::number(this->CenterAngle, 'e', 18));
xmlWriter.writeTextElement("LookSide", this->LookSide);
// sateantpos
@@ -259,60 +260,60 @@ void SARSimulationImageL1Dataset::saveToXml()
for (long i = 0; i < this->sateposes.count(); i++) {
xmlWriter.writeStartElement("AntPosParam");
- xmlWriter.writeTextElement("time", QString::number(this->sateposes[i].time)); // time
- xmlWriter.writeTextElement("Px", QString::number(this->sateposes[i].Px)); // Px
- xmlWriter.writeTextElement("Py", QString::number(this->sateposes[i].Py)); // Py
- xmlWriter.writeTextElement("Pz", QString::number(this->sateposes[i].Pz)); // Pz
- xmlWriter.writeTextElement("Vx", QString::number(this->sateposes[i].Vx)); // Vx
- xmlWriter.writeTextElement("Vy", QString::number(this->sateposes[i].Vy)); // Vy
- xmlWriter.writeTextElement("Vz", QString::number(this->sateposes[i].Vz)); // Vz
- xmlWriter.writeTextElement("AntDirectX", QString::number(this->sateposes[i].AntDirectX)); // AntDirectX
- xmlWriter.writeTextElement("AntDirectY", QString::number(this->sateposes[i].AntDirectY)); // AntDirectY
- xmlWriter.writeTextElement("AntDirectZ", QString::number(this->sateposes[i].AntDirectZ)); // AntDirectZ
- xmlWriter.writeTextElement("AVx", QString::number(this->sateposes[i].AVx)); // AVx
- xmlWriter.writeTextElement("AVy", QString::number(this->sateposes[i].AVy)); // AVy
- xmlWriter.writeTextElement("AVz", QString::number(this->sateposes[i].AVz)); // AVz
- xmlWriter.writeTextElement("lon", QString::number(this->sateposes[i].lon)); // lon
- xmlWriter.writeTextElement("lat", QString::number(this->sateposes[i].lat)); // lat
- xmlWriter.writeTextElement("ati", QString::number(this->sateposes[i].ati)); // ati
+ xmlWriter.writeTextElement("time", QString::number(this->sateposes[i].time, 'e', 18)); // time
+ xmlWriter.writeTextElement("Px", QString::number(this->sateposes[i].Px, 'e', 18)); // Px
+ xmlWriter.writeTextElement("Py", QString::number(this->sateposes[i].Py, 'e', 18)); // Py
+ xmlWriter.writeTextElement("Pz", QString::number(this->sateposes[i].Pz, 'e', 18)); // Pz
+ xmlWriter.writeTextElement("Vx", QString::number(this->sateposes[i].Vx, 'e', 18)); // Vx
+ xmlWriter.writeTextElement("Vy", QString::number(this->sateposes[i].Vy, 'e', 18)); // Vy
+ xmlWriter.writeTextElement("Vz", QString::number(this->sateposes[i].Vz, 'e', 18)); // Vz
+ xmlWriter.writeTextElement("AntDirectX", QString::number(this->sateposes[i].AntDirectX, 'e', 18)); // AntDirectX
+ xmlWriter.writeTextElement("AntDirectY", QString::number(this->sateposes[i].AntDirectY, 'e', 18)); // AntDirectY
+ xmlWriter.writeTextElement("AntDirectZ", QString::number(this->sateposes[i].AntDirectZ, 'e', 18)); // AntDirectZ
+ xmlWriter.writeTextElement("AVx", QString::number(this->sateposes[i].AVx, 'e', 18)); // AVx
+ xmlWriter.writeTextElement("AVy", QString::number(this->sateposes[i].AVy, 'e', 18)); // AVy
+ xmlWriter.writeTextElement("AVz", QString::number(this->sateposes[i].AVz, 'e', 18)); // AVz
+ xmlWriter.writeTextElement("lon", QString::number(this->sateposes[i].lon, 'e', 18)); // lon
+ xmlWriter.writeTextElement("lat", QString::number(this->sateposes[i].lat, 'e', 18)); // lat
+ xmlWriter.writeTextElement("ati", QString::number(this->sateposes[i].ati, 'e', 18)); // ati
xmlWriter.writeEndElement(); //
}
- xmlWriter.writeTextElement("ImageStartTime", QString::number(this->startImageTime));
- xmlWriter.writeTextElement("ImageEndTime", QString::number(this->EndImageTime));
+ xmlWriter.writeTextElement("ImageStartTime", QString::number(this->startImageTime, 'e', 18));
+ xmlWriter.writeTextElement("ImageEndTime", QString::number(this->EndImageTime, 'e', 18));
- xmlWriter.writeTextElement("incidenceAngleNearRange", QString::number(this->incidenceAngleNearRange));
- xmlWriter.writeTextElement("incidenceAngleFarRange", QString::number(this->incidenceAngleFarRange));
- xmlWriter.writeTextElement("TotalProcessedAzimuthBandWidth", QString::number(this->TotalProcessedAzimuthBandWidth));
- xmlWriter.writeTextElement("DopplerParametersReferenceTime", QString::number(this->DopplerParametersReferenceTime));
+ xmlWriter.writeTextElement("incidenceAngleNearRange", QString::number(this->incidenceAngleNearRange, 'e', 18));
+ xmlWriter.writeTextElement("incidenceAngleFarRange", QString::number(this->incidenceAngleFarRange, 'e', 18));
+ xmlWriter.writeTextElement("TotalProcessedAzimuthBandWidth", QString::number(this->TotalProcessedAzimuthBandWidth, 'e', 18));
+ xmlWriter.writeTextElement("DopplerParametersReferenceTime", QString::number(this->DopplerParametersReferenceTime, 'e', 18));
xmlWriter.writeStartElement("DopplerCentroidCoefficients");
- xmlWriter.writeTextElement("d0", QString::number(this->d0));
- xmlWriter.writeTextElement("d1", QString::number(this->d1));
- xmlWriter.writeTextElement("d2", QString::number(this->d2));
- xmlWriter.writeTextElement("d3", QString::number(this->d3));
- xmlWriter.writeTextElement("d4", QString::number(this->d4));
+ xmlWriter.writeTextElement("d0", QString::number(this->d0, 'e', 18));
+ xmlWriter.writeTextElement("d1", QString::number(this->d1, 'e', 18));
+ xmlWriter.writeTextElement("d2", QString::number(this->d2, 'e', 18));
+ xmlWriter.writeTextElement("d3", QString::number(this->d3, 'e', 18));
+ xmlWriter.writeTextElement("d4", QString::number(this->d4, 'e', 18));
xmlWriter.writeEndElement(); // DopplerCentroidCoefficients
xmlWriter.writeStartElement("DopplerRateValuesCoefficients");
- xmlWriter.writeTextElement("r0", QString::number(this->r0));
- xmlWriter.writeTextElement("r1", QString::number(this->r1));
- xmlWriter.writeTextElement("r2", QString::number(this->r2));
- xmlWriter.writeTextElement("r3", QString::number(this->r3));
- xmlWriter.writeTextElement("r4", QString::number(this->r4));
+ xmlWriter.writeTextElement("r0", QString::number(this->r0, 'e', 18));
+ xmlWriter.writeTextElement("r1", QString::number(this->r1, 'e', 18));
+ xmlWriter.writeTextElement("r2", QString::number(this->r2, 'e', 18));
+ xmlWriter.writeTextElement("r3", QString::number(this->r3, 'e', 18));
+ xmlWriter.writeTextElement("r4", QString::number(this->r4, 'e', 18));
xmlWriter.writeEndElement(); // DopplerRateValuesCoefficients
- xmlWriter.writeTextElement("latitude_center", QString::number(this->latitude_center));
- xmlWriter.writeTextElement("longitude_center", QString::number(this->longitude_center));
- xmlWriter.writeTextElement("latitude_topLeft", QString::number(this->latitude_topLeft));
- xmlWriter.writeTextElement("longitude_topLeft", QString::number(this->longitude_topLeft));
- xmlWriter.writeTextElement("latitude_topRight", QString::number(this->latitude_topRight));
- xmlWriter.writeTextElement("longitude_topRight", QString::number(this->longitude_topRight));
- xmlWriter.writeTextElement("latitude_bottomLeft", QString::number(this->latitude_bottomLeft));
- xmlWriter.writeTextElement("longitude_bottomLeft", QString::number(this->longitude_bottomLeft));
- xmlWriter.writeTextElement("latitude_bottomRight", QString::number(this->latitude_bottomRight));
- xmlWriter.writeTextElement("longitude_bottomRight", QString::number(this->longitude_bottomRight));
+ xmlWriter.writeTextElement("latitude_center", QString::number(this->latitude_center, 'e', 18));
+ xmlWriter.writeTextElement("longitude_center", QString::number(this->longitude_center, 'e', 18));
+ xmlWriter.writeTextElement("latitude_topLeft", QString::number(this->latitude_topLeft, 'e', 18));
+ xmlWriter.writeTextElement("longitude_topLeft", QString::number(this->longitude_topLeft, 'e', 18));
+ xmlWriter.writeTextElement("latitude_topRight", QString::number(this->latitude_topRight, 'e', 18));
+ xmlWriter.writeTextElement("longitude_topRight", QString::number(this->longitude_topRight, 'e', 18));
+ xmlWriter.writeTextElement("latitude_bottomLeft", QString::number(this->latitude_bottomLeft, 'e', 18));
+ xmlWriter.writeTextElement("longitude_bottomLeft", QString::number(this->longitude_bottomLeft, 'e', 18));
+ xmlWriter.writeTextElement("latitude_bottomRight", QString::number(this->latitude_bottomRight, 'e', 18));
+ xmlWriter.writeTextElement("longitude_bottomRight", QString::number(this->longitude_bottomRight, 'e', 18));
xmlWriter.writeEndElement(); // Parameters
xmlWriter.writeEndDocument();
@@ -371,6 +372,9 @@ ErrorCode SARSimulationImageL1Dataset::loadFromXml()
else if (xmlReader.name() == "Fs") {
this->Fs = xmlReader.readElementText().toDouble();
}
+ else if (xmlReader.name() == "PRF") {
+ this->prf = xmlReader.readElementText().toDouble();
+ }
else if(xmlReader.name() == "ImageStartTime"){
this->startImageTime = xmlReader.readElementText().toDouble();
}
@@ -873,6 +877,7 @@ QVector SARSimulationImageL1Dataset::getDopplerParams()
result[2] = d2;
result[3] = d3;
result[4] = d4;
+
return result;
}
@@ -888,11 +893,13 @@ void SARSimulationImageL1Dataset::setDopplerParams(double id0, double id1, doubl
QVector SARSimulationImageL1Dataset::getDopplerCenterCoff()
{
QVector result(5);
- result[0]=r0;
- result[1]=r1;
- result[2]=r2;
- result[3]=r3;
- result[4]=r4;
+ result[0] = r0;
+ result[1] = r1;
+ result[2] = r2;
+ result[3] = r3;
+ result[4] = r4;
+
+
return result;
}
diff --git a/BaseCommonLibrary/BaseTool/SARSimulationImageL1.h b/BaseCommonLibrary/BaseTool/SARSimulationImageL1.h
index c2abecf..6d8325a 100644
--- a/BaseCommonLibrary/BaseTool/SARSimulationImageL1.h
+++ b/BaseCommonLibrary/BaseTool/SARSimulationImageL1.h
@@ -143,9 +143,11 @@ public:
double getDopplerParametersReferenceTime();
void setDopplerParametersReferenceTime(double v);
+ // ղ
QVector getDopplerParams();
void setDopplerParams(double d0, double d1, double d2, double d3, double d4);
+ // ϵ
QVector getDopplerCenterCoff();
void setDopplerCenterCoff(double r0, double r1, double r2, double r3, double r4);
@@ -208,3 +210,9 @@ private:
};
+
+
+
+
+
+
diff --git a/BaseCommonLibrary/BaseTool/ShowProessAbstract.cpp b/BaseCommonLibrary/BaseTool/ShowProessAbstract.cpp
new file mode 100644
index 0000000..e862722
--- /dev/null
+++ b/BaseCommonLibrary/BaseTool/ShowProessAbstract.cpp
@@ -0,0 +1,18 @@
+#include "stdafx.h"
+#include "ShowProessAbstract.h"
+#include "BaseTool.h"
+#include "GeoOperator.h"
+#include "FileOperator.h"
+#include
+
+
+void ShowProessAbstract::showProcess(double precent, QString tip)
+{
+
+}
+
+void ShowProessAbstract::showToolInfo(QString tip)
+{
+}
+
+
diff --git a/BaseCommonLibrary/BaseTool/TestImageOperator.cpp b/BaseCommonLibrary/BaseTool/TestImageOperator.cpp
new file mode 100644
index 0000000..d2c9819
--- /dev/null
+++ b/BaseCommonLibrary/BaseTool/TestImageOperator.cpp
@@ -0,0 +1,146 @@
+#include "stdafx.h"
+#include "ImageOperatorBase.h"
+#include "BaseTool.h"
+#include "GeoOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include "FileOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include // OGRSpatialReference ڿռοת
+#include // GDALWarp
+
+
+
+void testOutAmpArr(QString filename, float* amp, long rowcount, long colcount)
+{
+
+
+ Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
+
+ for (long hii = 0; hii < rowcount; hii++) {
+ for (long hjj = 0; hjj < colcount; hjj++) {
+ h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
+ }
+ }
+ QString ampPath = getDebugDataPath(filename);
+ saveEigenMatrixXd2Bin(h_amp_img, ampPath);
+ qDebug() << filename.toLocal8Bit().constData();
+ qDebug() << "max:\t" << h_amp_img.maxCoeff();
+ qDebug() << "min:\t" << h_amp_img.minCoeff();
+}
+
+void testOutAmpArr(QString filename, double* amp, long rowcount, long colcount)
+{
+
+
+ Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
+
+ for (long hii = 0; hii < rowcount; hii++) {
+ for (long hjj = 0; hjj < colcount; hjj++) {
+ h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
+ }
+ }
+ QString ampPath = getDebugDataPath(filename);
+ saveEigenMatrixXd2Bin(h_amp_img, ampPath);
+ qDebug() << filename.toLocal8Bit().constData();
+ qDebug() << "max:\t" << h_amp_img.maxCoeff();
+ qDebug() << "min:\t" << h_amp_img.minCoeff();
+}
+
+
+void testOutClsArr(QString filename, long* amp, long rowcount, long colcount) {
+
+ Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
+
+ for (long hii = 0; hii < rowcount; hii++) {
+ for (long hjj = 0; hjj < colcount; hjj++) {
+ h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
+ }
+ }
+ QString ampPath = getDebugDataPath(filename);
+ saveEigenMatrixXd2Bin(h_amp_img, ampPath);
+ qDebug() << filename.toLocal8Bit().constData();
+ qDebug() << "max:\t" << h_amp_img.maxCoeff();
+ qDebug() << "min:\t" << h_amp_img.minCoeff();
+
+}
+
+void testOutComplexDoubleArr(QString filename, std::complex* data, long rowcount, long colcount)
+{
+ QString ampPath = getDebugDataPath(filename);
+ gdalImageComplex compleximg = CreateEchoComplex(ampPath, rowcount, colcount, 1);
+ compleximg.saveImage(data, 0, 0, rowcount, colcount, 1);
+
+ return;
+}
+
+void testOutDataArr(QString filename, double* data, long rowcount, long colcount)
+{
+ return testOutAmpArr(filename, data, rowcount, colcount);
+}
+
+void testOutDataArr(QString filename, float* data, long rowcount, long colcount)
+{
+ return testOutAmpArr(filename, data, rowcount, colcount);
+}
+
+void testOutDataArr(QString filename, long* data, long rowcount, long colcount)
+{
+ return testOutClsArr(filename, data, rowcount, colcount);
+}
+
+void testOutAntPatternTrans(QString antpatternfilename, double* antPatternArr,
+ double starttheta, double deltetheta,
+ double startphi, double deltaphi,
+ long thetanum, long phinum)
+{
+
+
+ Eigen::MatrixXd antPatternMatrix(thetanum, phinum);
+ for (long t = 0; t < thetanum; ++t) {
+ for (long p = 0; p < phinum; ++p) {
+ long index = t * phinum + p;
+ if (index < thetanum * phinum) {
+ antPatternMatrix(t, p) = static_cast(antPatternArr[index]); // Copy to Eigen matrix
+ }
+ }
+ }
+
+ Eigen::MatrixXd gt(2, 3);
+ gt(0, 0) = startphi;//x
+ gt(0, 1) = deltaphi;
+ gt(0, 2) = 0;
+
+ gt(1, 0) = starttheta;
+ gt(1, 1) = 0;
+ gt(1, 2) = deltetheta;
+
+ QString antpatternfilepath = getDebugDataPath(antpatternfilename);
+ gdalImage ds = CreategdalImageDouble(antpatternfilepath, thetanum, phinum, 1, gt, "", true, true, true);
+ ds.saveImage(antPatternMatrix, 0, 0, 1);
+}
+
+
+
+
+
+
diff --git a/BaseCommonLibrary/BaseTool/gdalImageComplexOperator.cpp b/BaseCommonLibrary/BaseTool/gdalImageComplexOperator.cpp
new file mode 100644
index 0000000..7d0df86
--- /dev/null
+++ b/BaseCommonLibrary/BaseTool/gdalImageComplexOperator.cpp
@@ -0,0 +1,547 @@
+#include "stdafx.h"
+#include "ImageOperatorBase.h"
+#include "BaseTool.h"
+#include "GeoOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include "FileOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include // OGRSpatialReference ڿռοת
+#include // GDALWarp
+
+
+
+
+
+
+
+
+
+gdalImageComplex::gdalImageComplex(const QString& raster_path)
+{
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ this->img_path = raster_path;
+
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
+ raster_path.toUtf8().constData(), GA_ReadOnly));
+ this->width = rasterDataset->GetRasterXSize();
+ this->height = rasterDataset->GetRasterYSize();
+ this->band_num = rasterDataset->GetRasterCount();
+
+ double* gt = new double[6];
+ rasterDataset->GetGeoTransform(gt);
+ this->gt = Eigen::MatrixXd(2, 3);
+ this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5];
+
+ double a = this->gt(0, 0);
+ double b = this->gt(0, 1);
+ double c = this->gt(0, 2);
+ double d = this->gt(1, 0);
+ double e = this->gt(1, 1);
+ double f = this->gt(1, 2);
+
+ this->projection = rasterDataset->GetProjectionRef();
+
+ // ͷͶӰ
+ GDALFlushCache((GDALDatasetH)rasterDataset);
+ GDALClose((GDALDatasetH)rasterDataset);
+ rasterDataset = NULL; // ָÿ�
+ this->InitInv_gt();
+ delete[] gt;
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+}
+
+gdalImageComplex::~gdalImageComplex() {}
+
+void gdalImageComplex::setData(Eigen::MatrixXcd data)
+{
+ this->data = data;
+}
+
+void gdalImageComplex::saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids)
+{
+
+
+
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
+ QString tip = u8"file path: " + this->img_path;
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, GDT_CFloat64, NULL); // ���
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ //delete[] gt_ptr;
+ }
+
+ int datarows = data.rows();
+ int datacols = data.cols();
+
+ if (this->getDataType() == GDT_CFloat64)
+ {
+
+ double* databuffer = new double[data.size() * 2];
+ for (int i = 0; i < data.rows(); i++) {
+ for (int j = 0; j < data.cols(); j++) {
+ databuffer[i * data.cols() * 2 + j * 2] = data(i, j).real();
+ databuffer[i * data.cols() * 2 + j * 2 + 1] = data(i, j).imag();
+ }
+ }
+
+ // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
+ // datarows, GDT_Float32,band_ids, num,0,0,0);
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, GDT_CFloat64, 0, 0);
+ GDALFlushCache(poDstDS);
+ delete databuffer;
+
+ }
+ else if (this->getDataType() == GDT_CFloat32) {
+
+ float* databuffer = new float[data.size() * 2];
+ for (int i = 0; i < data.rows(); i++) {
+ for (int j = 0; j < data.cols(); j++) {
+ databuffer[i * data.cols() * 2 + j * 2] = float(data(i, j).real());
+ databuffer[i * data.cols() * 2 + j * 2 + 1] =float( data(i, j).imag());
+ }
+ }
+
+ // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
+ // datarows, GDT_Float32,band_ids, num,0,0,0);
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, GDT_CFloat32, 0, 0);
+ GDALFlushCache(poDstDS);
+ delete databuffer;
+ }
+ else {
+ throw std::exception("gdalImageComplex::saveImage: data type error");
+ }
+
+
+
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ omp_unset_lock(&lock); //
+ omp_destroy_lock(&lock); //
+}
+
+void gdalImageComplex::saveImage(std::shared_ptr> data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
+{
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + rowCount > this->height || start_col + colCount > this->width) {
+ QString tip = u8"file path: " + this->img_path;
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ return;
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, GDT_CFloat64, NULL); // ���
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ //delete[] gt_ptr;
+ }
+
+ double* databuffer = new double[rowCount * colCount * 2];
+ for (long i = 0; i < rowCount; i++) {
+ for (long j = 0; j < colCount; j++) {
+ databuffer[i * colCount * 2 + j * 2] = data.get()[i * colCount + j].real();
+ databuffer[i * colCount * 2 + j * 2 + 1] = data.get()[i * colCount + j].imag();
+ }
+ }
+
+ // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
+ // datarows, GDT_Float32,band_ids, num,0,0,0);
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, colCount, rowCount,
+ databuffer, colCount, rowCount, GDT_CFloat64, 0, 0);
+ GDALFlushCache(poDstDS);
+ delete databuffer;
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ omp_unset_lock(&lock); //
+ omp_destroy_lock(&lock); //
+
+}
+
+void gdalImageComplex::saveImage(std::complex* data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
+{
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + rowCount > this->height || start_col + colCount > this->width) {
+ QString tip = u8"file path: " + this->img_path;
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ return;
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, GDT_CFloat64, NULL); // ���
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ //delete[] gt_ptr;
+ }
+
+ double* databuffer = new double[rowCount * colCount * 2];
+ for (long i = 0; i < rowCount; i++) {
+ for (long j = 0; j < colCount; j++) {
+ databuffer[i * colCount * 2 + j * 2] = data[i * colCount + j].real();
+ databuffer[i * colCount * 2 + j * 2 + 1] = data[i * colCount + j].imag();
+ }
+ }
+
+ // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
+ // datarows, GDT_Float32,band_ids, num,0,0,0);
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, colCount, rowCount,
+ databuffer, colCount, rowCount, GDT_CFloat64, 0, 0);
+ GDALFlushCache(poDstDS);
+ delete databuffer;
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ omp_unset_lock(&lock); //
+ omp_destroy_lock(&lock); //
+}
+
+Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col, int rows_count,
+ int cols_count, int band_ids)
+{
+ GDALDataset* poDataset;
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+
+ // TIFFļ
+ poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
+ if (poDataset == nullptr) {
+ QMessageBox::warning(nullptr, u8"", u8"ļ" + this->img_path);
+ qDebug() << u8"ļ" + this->img_path;
+ }
+
+ // ȡݼĵһ
+ GDALRasterBand* poBand;
+ poBand = poDataset->GetRasterBand(1);
+ rows_count = start_row + rows_count <= this->height ? rows_count : this->height - start_row;
+ cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
+ // ȡϢǸ
+ int nXSize = cols_count; poBand->GetXSize();
+ int nYSize = rows_count; poBand->GetYSize();
+ Eigen::MatrixXcd rasterData(nYSize, nXSize); // ʹEigenMatrixXcd
+ if (this->getDataType() == GDT_CFloat64)
+ {
+ long long pixelCount =long long( nXSize) *long long( nYSize);
+ double* databuffer = new double[pixelCount * 2];
+ poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count, rows_count, GDT_CFloat64, 0, 0);
+ GDALClose((GDALDatasetH)poDataset);
+
+ for (long long i = 0; i < nYSize; i++) {
+ for (long long j = 0; j < nXSize; j++) {
+ rasterData(i, j) = std::complex(databuffer[i * nXSize * 2 + j * 2],
+ databuffer[i * nXSize * 2 + j * 2 + 1]);
+ }
+ }
+
+ delete[] databuffer;
+ }
+ else if(this->getDataType()==GDT_CFloat32)
+ {
+ long long pixelCount = long long(nXSize) * long long(nYSize);
+ float* databuffer = new float[pixelCount * 2];
+ poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count, rows_count, GDT_CFloat32, 0, 0);
+ GDALClose((GDALDatasetH)poDataset);
+
+ for (long long i = 0; i < nYSize; i++) {
+ for (long long j = 0; j < nXSize; j++) {
+ rasterData(i, j) = std::complex(databuffer[i * nXSize * 2 + j * 2],
+ databuffer[i * nXSize * 2 + j * 2 + 1]);
+ }
+ }
+
+ delete[] databuffer;
+ }
+
+
+ return rasterData;
+}
+
+std::shared_ptr> gdalImageComplex::getDataComplexSharePtr(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
+{
+ GDALDataset* poDataset;
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+
+ // TIFFļ
+ poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
+ if (poDataset == nullptr) {
+ QMessageBox::warning(nullptr, u8"", u8"ļ" + this->img_path);
+ qDebug() << u8"ļ" + this->img_path;
+ }
+
+ // ȡݼĵһ
+ GDALRasterBand* poBand;
+ poBand = poDataset->GetRasterBand(1);
+
+ // ȡϢǸ
+
+
+ double* databuffer = new double[rows_count * cols_count * 2];
+ poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count,
+ rows_count, GDT_CFloat64, 0, 0);
+ GDALClose((GDALDatasetH)poDataset);
+
+ std::shared_ptr> rasterData(new std::complex[rows_count * cols_count], delArrPtr);
+
+ for (size_t i = 0; i < rows_count; i++) {
+ for (size_t j = 0; j < cols_count; j++) {
+ rasterData.get()[i * cols_count + j] = std::complex(databuffer[i * cols_count * 2 + j * 2],
+ databuffer[i * cols_count * 2 + j * 2 + 1]);
+ }
+ }
+
+ delete[] databuffer;
+ return rasterData;
+}
+
+void gdalImageComplex::saveComplexImage()
+{
+ this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
+}
+void gdalImageComplex::savePreViewImage()
+{
+ qDebug() << "void gdalImageComplex::savePreViewImage()";
+ Eigen::MatrixXd data_abs = Eigen::MatrixXd::Zero(this->height, this->width);
+ data_abs = (this->data.array().real().pow(2) + this->data.array().imag().pow(2))
+ .array()
+ .log10() * 10.0; //
+
+ double min_abs = data_abs.minCoeff(); // ֵ
+ double max_abs = data_abs.maxCoeff(); // Сֵ
+ double delta = (max_abs - min_abs) / 1000; // 1000λ
+ Eigen::MatrixX data_idx =
+ ((data_abs.array() - min_abs).array() / delta).array().floor().cast();
+
+ // ʼ
+ double hist[1001];
+ for (size_t i = 0; i < 1001; i++) {
+ hist[i] = 0; // ʼ
+ }
+ for (size_t i = 0; i < this->height; i++) {
+ for (size_t j = 0; j < this->width; j++) {
+ hist[data_idx(i, j)]++;
+ }
+ }
+
+ // ͳ
+ size_t count = this->height * this->width;
+ double precent = 0;
+ size_t curCount = 0;
+ double pre2 = 0;
+ bool findprec_2 = true;
+ double pre98 = 0;
+ bool findprec_98 = true;
+ for (size_t i = 0; i < 1001; i++) {
+ precent = precent + hist[i];
+ if (findprec_2 & precent / count > 0.02) {
+ pre2 = i * delta + min_abs;
+ findprec_2 = false;
+ }
+ if (findprec_98 & precent / count > 0.98) {
+ pre98 = (i - 1) * delta + min_abs;
+ findprec_98 = false;
+ }
+ }
+ //
+ delta = (pre98 - pre2) / 200;
+ data_idx =
+ ((data_abs.array() - pre2).array() / delta).array().floor().cast();
+
+ for (size_t i = 0; i < this->height; i++) {
+ for (size_t j = 0; j < this->width; j++) {
+ if (data_idx(i, j) < 0) {
+ data_idx(i, j) = 0;
+ }
+ else if (data_idx(i, j) > 255) {
+ data_idx(i, j) = 255;
+ }
+ else {
+
+ }
+ }
+ }
+
+ // ֵ
+ QString filePath = this->img_path;
+ QFile file(filePath);
+ QFileInfo fileInfo(file);
+
+ QString directory = fileInfo.absolutePath();
+ qDebug() << "ļĿ¼" << directory;
+ QString baseName = fileInfo.completeBaseName();
+ qDebug() << "ļ" << baseName;
+
+ // ļ·
+ QString previewImagePath = JoinPath(directory, baseName + "_preview.png");
+ cv::Mat img(this->height, this->width, CV_8U, cv::Scalar(0));
+
+ for (size_t i = 0; i < this->height; i++) {
+ for (size_t j = 0; j < this->width; j++) {
+ img.at(i, j) = (uchar)(data_idx(i, j));
+ }
+ }
+ //std::string outimgpath=previewImagePath.toUtf8().data();
+ cv::imwrite(previewImagePath.toUtf8().data(), img);
+}
+
+
+
+gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width,
+ int band_num, Eigen::MatrixXd gt, QString projection,
+ bool need_gt, bool overwrite)
+{
+ if (exists_test(img_path.toUtf8().constData())) {
+ if (overwrite) {
+ gdalImageComplex result_img(img_path);
+ return result_img;
+ }
+ else {
+ throw "file has exist!!!";
+ exit(1);
+ }
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
+ GDT_CFloat64, NULL);
+ if (need_gt) {
+ poDstDS->SetProjection(projection.toUtf8().constData());
+
+ // ����ת������
+ double gt_ptr[6] = { 0 };
+ for (int i = 0; i < gt.rows(); i++) {
+ for (int j = 0; j < gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ }
+ //for(int i = 1; i <= band_num; i++) {
+ // poDstDS->GetRasterBand(i)->SetNoDataValue(0); // ز
+ //}
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ gdalImageComplex result_img(img_path);
+ return result_img;
+}
+
+gdalImageComplex CreategdalImageComplexNoProj(const QString& img_path, int height, int width, int band_num, bool overwrite)
+{
+
+ // ļ
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+
+ GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = (poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_CFloat64, NULL));
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ //poDstDS.reset();
+ omp_unset_lock(&lock); //
+ omp_destroy_lock(&lock); //
+
+ gdalImageComplex result_img(img_path);
+ return result_img;
+}
+
+gdalImageComplex CreateEchoComplex(const QString& img_path, int height, int width, int band_num)
+{
+ // ͼ
+ Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
+ //Xgeo = GeoTransform[0] + Xpixel * GeoTransform[1] + Ypixel * GeoTransform[2]
+ //Ygeo = GeoTransform[3] + Xpixel * GeoTransform[4] + Ypixel * GeoTransform[5]
+ // X
+ gt(0, 0) = 0; gt(0, 2) = 1; gt(0, 2) = 0;
+ gt(1, 0) = 0; gt(1, 1) = 0; gt(1, 2) = 1;
+ // Y
+ QString projection = "";
+ gdalImageComplex echodata = CreategdalImageComplex(img_path, height, width, 1, gt, projection, false, true);
+ return echodata;
+
+}
+
diff --git a/BaseCommonLibrary/BaseTool/gdalImageOperator.cpp b/BaseCommonLibrary/BaseTool/gdalImageOperator.cpp
new file mode 100644
index 0000000..9fe8bfd
--- /dev/null
+++ b/BaseCommonLibrary/BaseTool/gdalImageOperator.cpp
@@ -0,0 +1,1530 @@
+#include "stdafx.h"
+#include "ImageOperatorBase.h"
+#include "BaseTool.h"
+#include "GeoOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include "FileOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include // OGRSpatialReference ڿռοת
+#include // GDALWarp
+
+
+
+
+
+gdalImage::gdalImage()
+{
+ this->height = 0;
+ this->width = 0;
+ this->data_band_ids = 1;
+ this->start_row = 0;
+ this->start_col = 0;
+}
+
+///
+/// �ͼȡӰ�?1?7
+///
+///
+gdalImage::gdalImage(const QString& raster_path)
+{
+ omp_lock_t lock;
+ omp_init_lock(&lock); // ��ʼ��
+ omp_set_lock(&lock); // ��û�?1?7
+ this->img_path = raster_path;
+
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // עʽ��?1?7
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ // ��DEMӰ��
+ GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.toUtf8().constData(), GA_ReadOnly));
+ if (nullptr == rasterDataset || NULL == rasterDataset) {
+ QMessageBox::warning(nullptr, u8"", QString(u8"ļ") + QString(raster_path));
+ exit(1);
+ }
+ this->width = rasterDataset->GetRasterXSize();
+ this->height = rasterDataset->GetRasterYSize();
+ this->band_num = rasterDataset->GetRasterCount();
+
+ double* gt = new double[6];
+ // ��÷���
+ rasterDataset->GetGeoTransform(gt);
+ this->gt = Eigen::MatrixXd(2, 3);
+ this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5];
+
+ this->projection = rasterDataset->GetProjectionRef();
+ // ��
+ // double* inv_gt = new double[6];;
+ // GDALInvGeoTransform(gt, inv_gt); // ��
+ // �ͶӰ
+ GDALFlushCache((GDALDatasetH)rasterDataset);
+ GDALClose((GDALDatasetH)rasterDataset);
+ rasterDataset = NULL; // ָÿ�
+ this->InitInv_gt();
+ delete[] gt;
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+}
+
+gdalImage::~gdalImage() {}
+
+void gdalImage::setHeight(int height)
+{
+ this->height = height;
+}
+
+void gdalImage::setWidth(int width)
+{
+ this->width = width;
+}
+
+void gdalImage::setTranslationMatrix(Eigen::MatrixXd gt)
+{
+ this->gt = gt;
+}
+
+void gdalImage::setData(Eigen::MatrixXd, int data_band_ids)
+{
+ this->data = data;
+ this->data_band_ids = data_band_ids;
+}
+
+Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, int cols_count,
+ int band_ids = 1)
+{
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
+ this->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 <= this->height ? rows_count : this->height - start_row;
+ cols_count = start_col + cols_count <= this->width ? cols_count : this->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);
+ for (int i = 0; i < rows_count; i++) {
+ for (int j = 0; j < cols_count; j++) {
+ datamatrix(i, j) = long(temp[i * cols_count + j]) * 1.0;
+ }
+ }
+ 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);
+ for (int i = 0; i < rows_count; i++) {
+ for (int j = 0; j < cols_count; j++) {
+ datamatrix(i, j) = long(temp[i * cols_count + j]) * 1.0;
+ }
+ }
+ 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);
+ for (int i = 0; i < rows_count; i++) {
+ for (int j = 0; j < cols_count; j++) {
+ datamatrix(i, j) = long(temp[i * cols_count + j]) * 1.0;
+ }
+ }
+ 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);
+ 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_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);
+ for (int i = 0; i < rows_count; i++) {
+ for (int j = 0; j < cols_count; j++) {
+ datamatrix(i, j) = long(temp[i * cols_count + j]) * 1.0;
+ }
+ }
+ 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);
+
+ 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_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);
+
+ 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 {
+ }
+ GDALClose((GDALDatasetH)rasterDataset);
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ return datamatrix;
+}
+
+Eigen::MatrixXf gdalImage::getDataf(int start_row, int start_col, int rows_count, int cols_count,
+ int band_ids = 1)
+{
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
+ this->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 <= this->height ? rows_count : this->height - start_row;
+ cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
+
+ Eigen::MatrixXf 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);
+ 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_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);
+ 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_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);
+ 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_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);
+ 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_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);
+ 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_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);
+
+ 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_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);
+
+ 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 {
+ }
+ GDALClose((GDALDatasetH)rasterDataset);
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ return datamatrix;
+}
+
+
+
+Eigen::MatrixXi gdalImage::getDatai(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
+{
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
+ this->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 <= this->height ? rows_count : this->height - start_row;
+ cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
+
+ Eigen::MatrixXi 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);
+ 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_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);
+ 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_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);
+ 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_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);
+ 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_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);
+ 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_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);
+
+ 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_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);
+
+ 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 {
+ }
+ GDALClose((GDALDatasetH)rasterDataset);
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ return datamatrix;
+}
+
+ErrorCode gdalImage::getData(double* data, int start_row, int start_col, int rows_count, int cols_count, int band_ids)
+{
+ ErrorCode state = ErrorCode::SUCCESS;
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->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 <= this->height ? rows_count : this->height - start_row;
+ cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
+
+
+
+ if (gdal_datatype == GDT_Float64) {
+ demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, data, cols_count, rows_count, gdal_datatype, 0, 0);
+ }
+ else {
+ state = ErrorCode::FAIL;
+ }
+ GDALClose((GDALDatasetH)rasterDataset);
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+ return state;
+}
+
+ErrorCode gdalImage::getData(long* data, int start_row, int start_col, int rows_count, int cols_count, int band_ids)
+{
+ ErrorCode state = ErrorCode::SUCCESS;
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->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 <= this->height ? rows_count : this->height - start_row;
+ cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
+
+ if (gdal_datatype == GDT_Int32) {
+ demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, data, cols_count, rows_count, gdal_datatype, 0, 0);
+ }
+ else {
+ state = ErrorCode::FAIL;
+ }
+ GDALClose((GDALDatasetH)rasterDataset);
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+ return state;
+}
+
+Eigen::MatrixXd gdalImage::getGeoTranslation()
+{
+ return this->gt;
+}
+
+GDALDataType gdalImage::getDataType()
+{
+ GDALDataset* rasterDataset =
+ (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));
+ GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
+ return gdal_datatype;
+}
+
+///
+///
+///
+///
+///
+///
+///
+void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0,
+ int band_ids = 1)
+{
+ GDALDataType datetype = this->getDataType();
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
+ QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+
+ QString filesuffer = getFileExtension(this->img_path).toLower();
+ bool isTiff = filesuffer.contains("tif");
+ GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, datetype, NULL); // ���
+
+ if (nullptr == poDstDS) {
+ QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ return;
+ }
+
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ //delete gt_ptr;
+ }
+
+
+ int datarows = data.rows();
+ int datacols = data.cols();
+ void* databuffer = nullptr;
+ if (datetype == GDT_Float32) {
+ databuffer = new float[datarows * datacols];
+ for (int i = 0; i < datarows; i++) {
+ for (int j = 0; j < datacols; j++) {
+ ((float*)databuffer)[i * datacols + j] = float(data(i, j));
+ }
+ }
+
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, datetype, 0, 0);
+ }
+ else if (datetype == GDT_Float64) {
+ databuffer = new double[datarows * datacols];
+ for (int i = 0; i < datarows; i++) {
+ for (int j = 0; j < datacols; j++) {
+ ((double*)databuffer)[i * datacols + j] = double(data(i, j));
+ }
+ }
+
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, datetype, 0, 0);
+ }
+ else {
+
+ }
+ GDALFlushCache(poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ delete[] databuffer;
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+}
+
+void gdalImage::saveImage(Eigen::MatrixXf data, int start_row = 0, int start_col = 0,
+ int band_ids = 1)
+{
+ GDALDataType datetype = this->getDataType();
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
+ QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ QString filesuffer = getFileExtension(this->img_path).toLower();
+ bool isTiff = filesuffer.contains("tif");
+ GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, datetype, NULL); // ���
+
+ if (nullptr == poDstDS) {
+ QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ return;
+ }
+
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ //delete gt_ptr;
+ }
+
+
+ int datarows = data.rows();
+ int datacols = data.cols();
+ void* databuffer = nullptr;
+ if (datetype == GDT_Float32) {
+ databuffer = new float[datarows * datacols];
+ for (int i = 0; i < datarows; i++) {
+ for (int j = 0; j < datacols; j++) {
+ ((float*)databuffer)[i * datacols + j] = float(data(i, j));
+ }
+ }
+
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, datetype, 0, 0);
+ }
+ else if (datetype == GDT_Float64) {
+ databuffer = new double[datarows * datacols];
+ for (int i = 0; i < datarows; i++) {
+ for (int j = 0; j < datacols; j++) {
+ ((double*)databuffer)[i * datacols + j] = double(data(i, j));
+ }
+ }
+
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, datetype, 0, 0);
+ }
+ else {
+
+ }
+ GDALFlushCache(poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ delete[] databuffer;
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+}
+
+
+void gdalImage::saveImage(Eigen::MatrixXi data, int start_row, int start_col, int band_ids)
+{
+ GDALDataType datetype = this->getDataType();
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
+ QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ QString filesuffer = getFileExtension(this->img_path).toLower();
+ bool isTiff = filesuffer.contains("tif");
+ GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, GDT_Float32, NULL); // ���
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ //delete gt_ptr;
+ }
+
+ long datarows = data.rows();
+ long datacols = data.cols();
+
+ long* databuffer = new long[datarows * datacols]; // (float*)malloc(datarows * datacols * sizeof(float));
+
+ for (long i = 0; i < datarows; i++) {
+ for (long j = 0; j < datacols; j++) {
+ databuffer[i * datacols + j] = data(i, j);
+ }
+ }
+ // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
+ // datarows, GDT_Float32,band_ids, num,0,0,0);
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, datetype, 0, 0);
+
+ GDALFlushCache(poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ delete[] databuffer;
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+}
+
+void gdalImage::saveImage(std::shared_ptr data, int start_row, int start_col, int rowcount, int colcount, int band_ids)
+{
+ GDALDataType datetype = this->getDataType();
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + rowcount > this->height || start_col + colcount > this->width) {
+ QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + rowcount) + ", " + QString::number(start_col + colcount) + ") ";
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ QString filesuffer = getFileExtension(this->img_path).toLower();
+ bool isTiff = filesuffer.contains("tif");
+ GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, GDT_Float64, NULL); // ���
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ }
+
+ long datarows = rowcount;
+ long datacols = colcount;
+ double* databuffer = new double[datarows * datacols];
+ if (datetype == GDT_Float64) {
+ memcpy(databuffer, data.get(), sizeof(double) * datarows * datacols);
+ }
+ else {
+ for (long i = 0; i < datarows; i++) {
+ for (long j = 0; j < datacols; j++) {
+ databuffer[i * datacols + j] = data.get()[i * datacols + j];
+ }
+ }
+ }
+
+ // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
+ // datarows, GDT_Float32,band_ids, num,0,0,0);
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, datetype, 0, 0);
+
+ GDALFlushCache(poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ delete[] databuffer;
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+}
+
+void gdalImage::saveImage(std::shared_ptr data, int start_row, int start_col, int rowcount, int colcount, int band_ids)
+{
+ GDALDataType datetype = this->getDataType();
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + rowcount > this->height || start_col + colcount > this->width) {
+ QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + rowcount) + ", " + QString::number(start_col + colcount) + ") ";
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ QString filesuffer = getFileExtension(this->img_path).toLower();
+ bool isTiff = filesuffer.contains("tif");
+ GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, GDT_Float32, NULL); // ���
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ }
+
+ long datarows = rowcount;
+ long datacols = colcount;
+ float* databuffer = new float[datarows * datacols];
+ if (datetype == GDT_Float32) {
+ memcpy(databuffer, data.get(), sizeof(float) * datarows * datacols);
+ }
+ else {
+ for (long i = 0; i < datarows; i++) {
+ for (long j = 0; j < datacols; j++) {
+ databuffer[i * datacols + j] = data.get()[i * datacols + j];
+ }
+ }
+ }
+
+ // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
+ // datarows, GDT_Float32,band_ids, num,0,0,0);
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, datetype, 0, 0);
+
+ GDALFlushCache(poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ //delete poDstDS;
+ //poDstDS = nullptr;
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ delete[] databuffer;
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+}
+
+
+void gdalImage::saveImage(std::shared_ptr data, int start_row, int start_col, int rowcount, int colcount, int band_ids)
+{
+ GDALDataType datetype = this->getDataType();
+ omp_lock_t lock;
+ omp_init_lock(&lock);
+ omp_set_lock(&lock);
+ if (start_row + rowcount > this->height || start_col + colcount > this->width) {
+ QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + rowcount) + ", " + QString::number(start_col + colcount) + ") ";
+ qDebug() << tip;
+ throw std::exception(tip.toUtf8().constData());
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
+ QString filesuffer = getFileExtension(this->img_path).toLower();
+ bool isTiff = filesuffer.contains("tif");
+ GDALDriver* poDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
+ GDALDataset* poDstDS = nullptr;
+ if (exists_test(this->img_path)) {
+ poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
+ }
+ else {
+ poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
+ this->band_num, GDT_Float32, NULL); // ���
+ poDstDS->SetProjection(this->projection.toUtf8().constData());
+
+ double gt_ptr[6];
+ for (int i = 0; i < this->gt.rows(); i++) {
+ for (int j = 0; j < this->gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = this->gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ }
+
+ long datarows = rowcount;
+ long datacols = colcount;
+ int* databuffer = new int[datarows * datacols];
+ if (datetype == GDT_Int32) {
+ memcpy(databuffer, data.get(), sizeof(int) * datarows * datacols);
+ }
+ else {
+ for (long i = 0; i < datarows; i++) {
+ for (long j = 0; j < datacols; j++) {
+ databuffer[i * datacols + j] = data.get()[i * datacols + j];
+ }
+ }
+ }
+
+ // poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
+ // datarows, GDT_Float32,band_ids, num,0,0,0);
+ poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
+ databuffer, datacols, datarows, datetype, 0, 0);
+
+ GDALFlushCache(poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ delete[] databuffer;
+ omp_unset_lock(&lock); // �ͷŻ��
+ omp_destroy_lock(&lock); // ٻ��
+}
+
+
+void gdalImage::saveImage()
+{
+ this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
+}
+
+
+
+
+void gdalImage::setNoDataValue(double nodatavalue = -9999, int band_ids = 1)
+{
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // עʽ��?1?7
+ // GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
+ poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue);
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+}
+
+void gdalImage::setNoDataValuei(int nodatavalue, int band_ids)
+{
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // עʽ��?1?7
+ // GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
+ poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue);
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+}
+
+double gdalImage::getNoDataValue(int band_ids)
+{
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // עʽ��?1?7
+ // GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
+ double v = poDstDS->GetRasterBand(band_ids)->GetNoDataValue();
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ return v;
+}
+
+int gdalImage::getNoDataValuei(int band_ids)
+{
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // עʽ��?1?7
+ // GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
+ int v = poDstDS->GetRasterBand(band_ids)->GetNoDataValue();
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ return v;
+}
+
+
+int gdalImage::InitInv_gt()
+{
+ // 1 lon lat = x
+ // 1 lon lat = y
+ Eigen::MatrixXd temp = Eigen::MatrixXd::Zero(2, 3);
+ this->inv_gt = temp;
+ double a = this->gt(0, 0);
+ double b = this->gt(0, 1);
+ double c = this->gt(0, 2);
+ double d = this->gt(1, 0);
+ double e = this->gt(1, 1);
+ double f = this->gt(1, 2);
+ double g = 1;
+ double det_gt = b * f - c * e;
+ if (det_gt == 0) {
+ return 0;
+ }
+ this->inv_gt(0, 0) = (c * d - a * f) / det_gt; // 2
+ this->inv_gt(0, 1) = f / det_gt; // lon
+ this->inv_gt(0, 2) = -c / det_gt; // lat
+ this->inv_gt(1, 0) = (a * e - b * d) / det_gt; // 1
+ this->inv_gt(1, 1) = -e / det_gt; // lon
+ this->inv_gt(1, 2) = b / det_gt; // lat
+ return 1;
+}
+
+Landpoint gdalImage::getRow_Col(double lon, double lat)
+{
+ Landpoint p{ 0, 0, 0 };
+ p.lon = this->inv_gt(0, 0) + lon * this->inv_gt(0, 1) + lat * this->inv_gt(0, 2); // x
+ p.lat = this->inv_gt(1, 0) + lon * this->inv_gt(1, 1) + lat * this->inv_gt(1, 2); // y
+ return p;
+}
+
+Landpoint gdalImage::getLandPoint(double row, double col, double ati = 0)
+{
+ Landpoint p{ 0, 0, 0 };
+ p.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); // x
+ p.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); // y
+ p.ati = ati;
+ return p;
+}
+
+void gdalImage::getLandPoint(double row, double col, double ati, Landpoint& Lp)
+{
+ Lp.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); // x
+ Lp.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); // y
+ Lp.ati = ati;
+
+}
+
+
+
+
+
+
+double gdalImage::mean(int bandids)
+{
+ double mean_value = 0;
+ double count = this->height * this->width;
+ int line_invert = 100;
+ int start_ids = 0;
+ do {
+ Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
+ mean_value = mean_value + sar_a.sum() / count;
+ start_ids = start_ids + line_invert;
+ } while (start_ids < this->height);
+ return mean_value;
+}
+
+double gdalImage::BandmaxValue(int bandids)
+{
+ double max_value = 0;
+ bool state_max = true;
+ int line_invert = 100;
+ int start_ids = 0;
+ double temp_max = 0;
+ do {
+ Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
+ if (state_max) {
+ state_max = false;
+ max_value = sar_a.maxCoeff();
+ }
+ else {
+ temp_max = sar_a.maxCoeff();
+ if (max_value < temp_max) {
+ max_value = temp_max;
+ }
+ }
+ start_ids = start_ids + line_invert;
+ } while (start_ids < this->height);
+ return max_value;
+}
+
+double gdalImage::BandminValue(int bandids)
+{
+ double min_value = 0;
+ bool state_min = true;
+ int line_invert = 100;
+ int start_ids = 0;
+ double temp_min = 0;
+ do {
+ Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
+ if (state_min) {
+ state_min = false;
+ min_value = sar_a.minCoeff();
+ }
+ else {
+ temp_min = sar_a.minCoeff();
+ if (min_value < temp_min) {
+ min_value = temp_min;
+ }
+ }
+ start_ids = start_ids + line_invert;
+ } while (start_ids < this->height);
+ return min_value;
+}
+
+GDALRPCInfo gdalImage::getRPC()
+{
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
+ CPLSetConfigOption("GDAL_DATA", "./data");
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // ע���
+ // ���
+ GDALDataset* pDS = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
+ // ��Ԫ��л�ȡRPC��Ϣ
+ char** papszRPC = pDS->GetMetadata("RPC");
+
+ // �ȡ��RPC��Ϣ�ɽṹ�?1?7
+ GDALRPCInfo oInfo;
+ GDALExtractRPCInfo(papszRPC, &oInfo);
+
+ GDALClose((GDALDatasetH)pDS);
+
+ return oInfo;
+}
+
+Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points)
+{
+ if (points.cols() != 3) {
+ throw new std::exception("the size of points is equit 3!!!");
+ }
+
+ Eigen::MatrixXd result(points.rows(), 3);
+ result.col(2) = points.col(2); // �߳�
+ points.col(2) = points.col(2).array() * 0 + 1;
+ points.col(0).swap(points.col(2)); // �
+ Eigen::MatrixXd gts(3, 2);
+ gts.col(0) = this->gt.row(0);
+ gts.col(1) = this->gt.row(1);
+
+ result.block(0, 0, points.rows(), 2) = points * gts;
+ return result;
+}
+
+Eigen::MatrixXd gdalImage::getHist(int bandids)
+{
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // עʽ��?1?7
+ // ��DEMӰ��
+ GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
+ this->img_path.toUtf8().constData(), GA_ReadOnly)); // ��ֻ�ʽ��ȡ�Ӱ��
+
+ GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
+ GDALRasterBand* xBand = rasterDataset->GetRasterBand(bandids);
+
+ double dfMin = this->BandminValue(bandids);
+ double dfMax = this->BandmaxValue(bandids);
+ int count = int((dfMax - dfMin) / 0.01);
+ count = count > 255 ? count : 255;
+ GUIntBig* panHistogram = new GUIntBig[count];
+ xBand->GetHistogram(dfMin, dfMax, count, panHistogram, TRUE, FALSE, NULL, NULL);
+ Eigen::MatrixXd result(count, 2);
+ double delta = (dfMax - dfMin) / count;
+ for (int i = 0; i < count; i++) {
+ result(i, 0) = dfMin + i * delta;
+ result(i, 1) = double(panHistogram[i]);
+ }
+ delete[] panHistogram;
+ GDALClose((GDALDatasetH)rasterDataset);
+ return result;
+}
+
+RasterExtend gdalImage::getExtend()
+{
+ RasterExtend extend{ 0,0,0,0 };
+ double x1 = this->gt(0, 0);
+ double y1 = this->gt(1, 0);
+
+ double x2 = this->gt(0, 0) + (this->width - 1) * gt(0, 1) + (0) * gt(0, 2); //
+ double y2 = this->gt(1, 0) + (this->width - 1) * gt(1, 1) + (0) * gt(1, 2); // γ
+
+ double x3 = this->gt(0, 0) + (0) * gt(0, 1) + (this->height - 1) * gt(0, 2);
+ double y3 = this->gt(1, 0) + (0) * gt(1, 1) + (this->height - 1) * gt(1, 2);
+
+ double x4 = this->gt(0, 0) + (this->width - 1) * gt(0, 1) + (this->height - 1) * gt(0, 2);
+ double y4 = this->gt(1, 0) + (this->width - 1) * gt(1, 1) + (this->height - 1) * gt(1, 2);
+
+
+ extend.min_x = x1 < x2 ? x1 : x2;
+ extend.max_x = x1 < x2 ? x2 : x1;
+ extend.min_y = y1 < y2 ? y1 : y2;
+ extend.max_y = y1 < y2 ? y2 : y1;
+
+
+ extend.min_x = extend.min_x < x3 ? extend.min_x : x3;
+ extend.max_x = extend.max_x > x3 ? extend.max_x : x3;
+ extend.min_y = extend.min_y < y3 ? extend.min_y : y3;
+ extend.max_y = extend.max_y > y3 ? extend.max_y : y3;
+
+
+ extend.min_x = extend.min_x < x4 ? extend.min_x : x4;
+ extend.max_x = extend.max_x > x4 ? extend.max_x : x4;
+ extend.min_y = extend.min_y < y4 ? extend.min_y : y4;
+ extend.max_y = extend.max_y > y4 ? extend.max_y : y4;
+
+ return extend;
+}
+
+
+
+
+
+
+gdalImage CreategdalImageDouble(QString& img_path, int height, int width, int band_num, bool overwrite, bool isEnvi)
+{
+
+ if (exists_test(img_path.toUtf8().constData())) {
+ if (overwrite) {
+ gdalImage result_img(img_path);
+ return result_img;
+ }
+ else {
+ throw "file has exist!!!";
+ exit(1);
+ }
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // ע���ʽ�����?1?7
+ GDALDriver* poDriver = nullptr;
+ if (isEnvi) {
+ poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ }
+ else {
+ poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ }
+
+
+ GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_Float64, NULL); // ������
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ gdalImage result_img(img_path);
+ return result_img;
+
+}
+
+gdalImage CreategdalImageFloat(QString& img_path, int height, int width, int band_num, bool overwrite, bool isEnvi)
+{
+
+ if (exists_test(img_path.toUtf8().constData())) {
+ if (overwrite) {
+ gdalImage result_img(img_path);
+ return result_img;
+ }
+ else {
+ throw "file has exist!!!";
+ exit(1);
+ }
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // ע���ʽ�����?1?7
+ GDALDriver* poDriver = nullptr;
+ if (isEnvi) {
+ poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ }
+ else {
+ poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ }
+
+
+ GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_Float32, NULL); // ������
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ gdalImage result_img(img_path);
+ return result_img;
+}
+
+gdalImage CreategdalImageDouble(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite, bool isEnvi)
+{
+ if (exists_test(img_path.toUtf8().constData())) {
+ if (overwrite) {
+ gdalImage result_img(img_path);
+ return result_img;
+ }
+ else {
+ throw "file has exist!!!";
+ exit(1);
+ }
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // ע���ʽ�����?1?7
+ GDALDriver* poDriver = nullptr;
+ if (isEnvi) {
+ poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ }
+ else {
+ poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ }
+
+
+ GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
+ GDT_Float64, NULL); // ������
+ if (need_gt) {
+ if (!projection.isEmpty()) {
+ poDstDS->SetProjection(projection.toUtf8().constData());
+ }
+ double gt_ptr[6] = { 0 };
+ for (int i = 0; i < gt.rows(); i++) {
+ for (int j = 0; j < gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ }
+ for (int i = 1; i <= band_num; i++) {
+ poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
+ }
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ gdalImage result_img(img_path);
+ return result_img;
+}
+
+gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num,
+ Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite, bool isEnvi, GDALDataType datetype)
+{
+ if (exists_test(img_path.toUtf8().constData())) {
+ if (overwrite) {
+ gdalImage result_img(img_path);
+ return result_img;
+ }
+ else {
+ throw "file has exist!!!";
+ exit(1);
+ }
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // ע���ʽ�����?1?7
+ GDALDriver* poDriver = nullptr;
+ if (isEnvi) {
+ poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
+ }
+ else {
+ poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+ }
+
+
+ GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
+ datetype, NULL); // ������
+
+ if (NULL == poDstDS)
+ {
+ qDebug() << "Create image failed!";
+ throw "Create image failed!";
+ exit(1);
+ }
+
+
+ if (need_gt) {
+ if (!projection.isEmpty()) {
+ poDstDS->SetProjection(projection.toUtf8().constData());
+ }
+ double gt_ptr[6] = { 0 };
+ for (int i = 0; i < gt.rows(); i++) {
+ for (int j = 0; j < gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ }
+ for (int i = 1; i <= band_num; i++) {
+ poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
+ }
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ gdalImage result_img(img_path);
+ return result_img;
+}
+
+gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long epsgCode, GDALDataType eType, bool need_gt, bool overwrite, bool isENVI)
+{
+ if (exists_test(img_path.toUtf8().constData())) {
+ if (overwrite) {
+ gdalImage result_img(img_path);
+ return result_img;
+ }
+ else {
+ throw "file has exist!!!";
+ exit(1);
+ }
+ }
+ GDALAllRegister();
+ CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // ע���ʽ�����?1?7
+ GDALDriver* poDriver = isENVI ? GetGDALDriverManager()->GetDriverByName("ENVI") : GetGDALDriverManager()->GetDriverByName("GTiff");
+ GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, eType, NULL); // ������
+ if (need_gt) {
+ OGRSpatialReference oSRS;
+
+ if (oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) {
+ qDebug() << "Failed to import EPSG code " << epsgCode;
+ throw "Failed to import EPSG code ";
+ exit(1);
+ }
+ char* pszWKT = NULL;
+ oSRS.exportToWkt(&pszWKT);
+ qDebug() << "WKT of EPSG:" << epsgCode << " :\n" << pszWKT;
+ poDstDS->SetProjection(pszWKT);
+ double gt_ptr[6] = { 0 };
+ for (int i = 0; i < gt.rows(); i++) {
+ for (int j = 0; j < gt.cols(); j++) {
+ gt_ptr[i * 3 + j] = gt(i, j);
+ }
+ }
+ poDstDS->SetGeoTransform(gt_ptr);
+ }
+ for (int i = 1; i <= band_num; i++) {
+ poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
+ }
+ GDALFlushCache((GDALDatasetH)poDstDS);
+ GDALClose((GDALDatasetH)poDstDS);
+ GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
+ gdalImage result_img(img_path);
+ return result_img;
+}
+
+
+bool CopyProjectTransformMatrixFromRasterAToRasterB(QString RasterAPath, QString RasterBPath) {
+ // עGDAL
+ GDALAllRegister();
+
+ // ӰAֻģʽ
+ GDALDataset* ds_a = (GDALDataset*)GDALOpen(RasterAPath.toUtf8().constData(), GA_ReadOnly);
+ if (ds_a == nullptr) {
+ std::cerr << "ӰA" << std::endl;
+ return false;
+ }
+
+ // ȡAķͶӰϢ
+ double geotransform[6];
+ ds_a->GetGeoTransform(geotransform); // 任
+ const char* projection = ds_a->GetProjectionRef(); // WKTʽͶӰ
+
+ // ӰBģʽ
+ GDALDataset* ds_b = (GDALDataset*)GDALOpen(RasterBPath.toUtf8().constData(), GA_Update);
+ if (ds_b == nullptr) {
+ std::cerr << "ӰB" << std::endl;
+ GDALClose(ds_a);
+ return false;
+ }
+
+ // ÷
+ if (ds_b->SetGeoTransform(geotransform) != CE_None) {
+ std::cerr << "÷ʧ" << std::endl;
+ }
+
+ // ͶӰϵ
+ if (ds_b->SetProjection(projection) != CE_None) {
+ std::cerr << "ͶӰʧ" << std::endl;
+ }
+
+ // ͷԴ
+ GDALClose(ds_a);
+ GDALClose(ds_b);
+
+ return true;
+
+
+
+
+
+}
+
+
+
+
diff --git a/BaseCommonLibrary/ImageOperatorFuntion.cpp b/BaseCommonLibrary/ImageOperatorFuntion.cpp
new file mode 100644
index 0000000..ee94390
--- /dev/null
+++ b/BaseCommonLibrary/ImageOperatorFuntion.cpp
@@ -0,0 +1,2138 @@
+#include "stdafx.h"
+#include "ImageOperatorBase.h"
+#include "BaseTool.h"
+#include "GeoOperator.h"
+#include
+#include
+#include
+#include
+#include
+#include