KJ135
/
RasterProcessTool
mirror of http://172.16.0.12:5000/LAMPSARToolSoftware/RasterProcessTool.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

修改子模块

pull/1/head
陈增辉 2024-12-24 10:50:25 +08:00
parent 6ffbb622f6
commit c1b7c11560
24 changed files with 9337 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

163
BaseTool/BaseConstVariable.h Normal file
View File

@ -0,0 +1,163 @@
#pragma once
#ifndef BASECONSTVARIABLE_H
#define BASECONSTVARIABLE_H
#define EIGEN_USE_MKL_ALL
//#define EIGEN_NO_DEBUG
#define EIGEN_USE_BLAS
#define EIGEN_USE_LAPACK
#define EIGEN_VECTORIZE_SSE2
//#define DEBUGSHOWDIALOG
#define __CUDANVCC___ // 定义CUDA函数
//#define __PRFDEBUG__
//#include <mkl.h>
#include <complex>
#include <math.h>
#include <complex>
#define MATPLOTDRAWIMAGE
#define r2d 180/3.141592653589793238462643383279
#define d2r 3.141592653589793238462643383279/180
#define LIGHTSPEED 299792458
#define PRECISIONTOLERANCE 1e-6
#define Radians2Degrees(Radians) Radians*PI_180
#define Degrees2Radians(Degrees) Degrees*T180_PI
#define EARTHWE 0.000072292115
#define PI 3.141592653589793238462643383279
#define earthRoute 0.000072292115
#define Memory1GB 1073741824
#define Memory1MB 1048576
#define Memory1KB 1024
const std::complex<double> imagI(0, 1);
const double epsilon = 0.000000000000001;
const double pi = 3.14159265358979323846;
const double a = 6378137.0; //椭球长半轴
const double ae = 6378137.0; //椭球长半轴
const double ee = 0.0818191910428;// 第一偏心率
const double f_inverse = 298.257223563; //扁率倒数
const double b = a - a / f_inverse;
const double eSquare = (a * a - b * b) / (a * a);
const double e = sqrt(eSquare);
const double earth_Re = 6378136.49;
const double earth_Rp = (1 - 1 / f_inverse) * earth_Re;
const double earth_We = 0.000072292115;
/*********************************************** openMap参数 ********************************************************************/
static long Paral_num_thread = 14;
/*********************************************** 基础枚举 ********************************************************************/
enum POLARTYPEENUM { // 极化类型
POLARHH,
POLARHV,
POLARVH,
POLARVV,
POLARUNKOWN
};
/*********************************************** 基础结构体区域 ********************************************************************/
/// <summary>
/// 地理坐标点
/// </summary>
struct Landpoint // 点 SAR影像的像素坐标；
{
double lon; // 经度x lon pixel_col
double lat; // 纬度y lat pixel_row
double ati; // 高程z ati pixel_time
};
struct Point2 {
double x = 0;
double y = 0;
};
struct Point3 {
double x = 0;
double y = 0;
double z = 0;
};
struct DemBox {
double min_lon;
double max_lon;
double min_lat;
double max_lat;
};
/*********************************************** FEKO仿真参数 ********************************************************************/
struct SatellitePos {
double time;
double Px;
double Py;
double Pz;
double Vx;
double Vy;
double Vz;
};
struct SatelliteAntPos {
double time; // 0
double Px;
double Py;
double Pz;
double Vx;
double Vy;
double Vz; //7
double AntDirectX;
double AntDirectY;
double AntDirectZ;
double AVx;
double AVy;
double AVz;//13
double ZeroAntDiectX;
double ZeroAntDiectY;
double ZeroAntDiectZ;
double lon;
double lat;
double ati; // 19
};
/*********************************************** 指针回收区域 ********************************************************************/
inline void delArrPtr(void* p)
{
delete[] p;
p = nullptr;
}
inline void delPointer(void* p)
{
delete p;
p = nullptr;
}
#endif

629
BaseTool/BaseTool.cpp Normal file
View File

@ -0,0 +1,629 @@
#include "stdafx.h"
#define EIGEN_USE_BLAS
#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
////#include <mkl.h>
#include <omp.h>
#include <QDir>
#include < io.h >
#include < stdio.h >
#include < stdlib.h >
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
//#include <ogr_geos.h>
#include <ogrsf_frmts.h> //#include "ogrsf_frmts.h"
#include <fstream>
#include <proj.h>
#include "GeoOperator.h"
#include <iostream>
#include <boost/date_time/posix_time/posix_time.hpp>
#include "baseTool.h"
#include <gsl/gsl_fit.h> // 多项式拟合
#include <gsl/gsl_cdf.h> /* 提供了 gammaq 函数 */
#include <gsl/gsl_vector.h> /* 提供了向量结构*/
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include <qcoreapplication.h>
QString getCurrentTimeString() {
struct tm ConversionTime;
std::time_t t = std::time(NULL);
char mbstr[100];
_localtime64_s(&ConversionTime, &t);
std::strftime(mbstr, sizeof(mbstr), "%Y-%m-%d %H:%M:%S", &ConversionTime);
QString strTime = mbstr;
return strTime;
}
QString getCurrentShortTimeString() {
struct tm ConversionTime;
std::time_t t = std::time(NULL);
char mbstr[100];
_localtime64_s(&ConversionTime, &t);
std::strftime(mbstr, sizeof(mbstr), "%Y-%m-%d %H:%M:%S", &ConversionTime);
QString strTime = mbstr;
return strTime;
}
std::vector<QString> splitString(const QString& str, char delimiter)
{
QStringList tokens = str.split(delimiter);
return convertQStringListToStdVector(tokens);
}
std::complex<double> Cubic_Convolution_interpolation(double u, double v, Eigen::MatrixX<std::complex<double>> img)
{
if (img.rows() != 4 || img.cols() != 4) {
throw std::exception("the size of img's block is not right");
}
// 斤拷锟斤拷模锟斤拷
Eigen::MatrixX<std::complex<double>> wrc(1, 4);// 使锟斤拷 complex<double> 斤拷锟斤拷要原斤拷为锟剿伙拷取值
Eigen::MatrixX<std::complex<double>> wcr(4, 1);//
for (int i = 0; i < 4; i++) {
wrc(0, i) = Cubic_kernel_weight(u + 1 - i); // u+1,u,u-1,u-2
wcr(i, 0) = Cubic_kernel_weight(v + 1 - i);
}
Eigen::MatrixX<std::complex<double>> interValue = wrc * img * wcr;
return interValue(0, 0);
}
std::complex<double> Cubic_kernel_weight(double s)
{
s = abs(s);
if (s <= 1) {
return std::complex<double>(1.5 * s * s * s - 2.5 * s * s + 1, 0);
}
else if (s <= 2) {
return std::complex<double>(-0.5 * s * s * s + 2.5 * s * s - 4 * s + 2, 0);
}
else {
return std::complex<double>(0, 0);
}
}
/// <summary>
/// p11 p12 -- x
/// p0(u,v)
/// p21 p22
/// |
/// y
/// p11(0,0)
/// p21(0,1)
/// P12(1,0)
/// p22(1,1)
/// </summary>
/// <param name="p0">x,y,z</param>
/// <param name="p1">x,y,z</param>
/// <param name="p2">x,y,z</param>
/// <param name="p3">x,y,z</param>
/// <param name="p4">x,y,z</param>
/// <returns></returns>
double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12, Landpoint p22)
{
return p11.ati * (1 - p0.lon) * (1 - p0.lat) +
p12.ati * p0.lon * (1 - p0.lat) +
p21.ati * (1 - p0.lon) * p0.lat +
p22.ati * p0.lon * p0.lat;
}
bool onSegment(Point3 Pi, Point3 Pj, Point3 Q)
{
if ((Q.x - Pi.x) * (Pj.y - Pi.y) == (Pj.x - Pi.x) * (Q.y - Pi.y) //叉乘
//保证Q点坐标在pi,pj之间
&& std::min(Pi.x, Pj.x) <= Q.x && Q.x <= std::max(Pi.x, Pj.x)
&& std::min(Pi.y, Pj.y) <= Q.y && Q.y <= std::max(Pi.y, Pj.y))
return true;
else
return false;
}
Point3 invBilinear(Point3 p, Point3 a, Point3 b, Point3 c, Point3 d)
{
Point3 res;
Point3 e = b - a;
Point3 f = d - a;
Point3 g = a - b + c - d;
Point3 h = p - a;
double k2 = cross2d(g, f);
double k1 = cross2d(e, f) + cross2d(h, g);
double k0 = cross2d(h, e);
double u, v;
// if edges are parallel, this is a linear equation
if (abs(k2) < 0.001)
{
v = -k0 / k1;
u = (h.x - f.x * v) / (e.x + g.x * v);
p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v;
return p;
}
// otherwise, it's a quadratic
else
{
float w = k1 * k1 - 4.0 * k0 * k2;
if (w < 0.0) {
// 可能在边界上
if (onSegment(a, b, p)) {
Point3 tt = b - a;
Point3 ttpa = p - a;
double scater = ttpa / tt;
if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; }
p.z = a.z + scater * tt.z;
return p;
}
else if (onSegment(b, c, p)) {
Point3 tt = c - b;
Point3 ttpa = p - b;
double scater = ttpa / tt;
if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; }
p.z = b.z + scater * tt.z;
return p;
}
else if (onSegment(c, d, p)) {
Point3 tt = d - c;
Point3 ttpa = p - c;
double scater = ttpa / tt;
if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; }
p.z = c.z + scater * tt.z;
return p;
}
else if (onSegment(d, a, p)) {
Point3 tt = a - d;
Point3 ttpa = p - d;
double scater = ttpa / tt;
if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; }
p.z = d.z + scater * tt.z;
return p;
}
return { -9999,-9999,-9999 };
}
else {
w = sqrt(w);
float ik2 = 0.5 / k2;
float v = (-k1 - w) * ik2;
float u = (h.x - f.x * v) / (e.x + g.x * v);
if (u < 0.0 || u>1.0 || v < 0.0 || v>1.0)
{
v = (-k1 + w) * ik2;
u = (h.x - f.x * v) / (e.x + g.x * v);
}
p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v;
return p;
}
}
p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v;
return p;
}
double sind(double degree)
{
return sin(degree * d2r);
}
double cosd(double d)
{
return cos(d * d2r);
}
std::string Convert(float Num)
{
std::ostringstream oss;
oss << Num;
std::string str(oss.str());
return str;
}
QString JoinPath(const QString& path, const QString& filename)
{
QDir dir(path);
// Ensure the path ends with the appropriate separator
if (!QDir::isAbsolutePath(path))
dir.makeAbsolute();
return dir.filePath(filename);
}
std::vector<QString> convertQStringListToStdVector(const QStringList& qStringList)
{
std::vector<QString> stdVector;
for (const QString& str : qStringList) {
stdVector.push_back(str);
}
return stdVector;
};
ErrorCode polyfit(const double* x, const double* y, int xyLength, int poly_n, std::vector<double>& out_factor, double& out_chisq) {
/*
* x：自变量，视差
* y：因变量，距离
* xyLength: x、y长度
* poly_n：拟合的阶次
* out_factor：拟合的系数结果，从0阶到poly_n阶的系数
* out_chisq：拟合曲线与数据点的优值函数最小值 ,χ2 检验
*/
gsl_matrix* XX = gsl_matrix_alloc(xyLength, poly_n + 1);
gsl_vector* c = gsl_vector_alloc(poly_n + 1);
gsl_matrix* cov = gsl_matrix_alloc(poly_n + 1, poly_n + 1);
gsl_vector* vY = gsl_vector_alloc(xyLength);
for (size_t i = 0; i < xyLength; i++)
{
gsl_matrix_set(XX, i, 0, 1.0);
gsl_vector_set(vY, i, y[i]);
for (int j = 1; j <= poly_n; j++)
{
gsl_matrix_set(XX, i, j, pow(x[i], j));
}
}
gsl_multifit_linear_workspace* workspace = gsl_multifit_linear_alloc(xyLength, poly_n + 1);
int r = gsl_multifit_linear(XX, vY, c, cov, &out_chisq, workspace);
gsl_multifit_linear_free(workspace);
out_factor.resize(c->size, 0);
for (size_t i = 0; i < c->size; ++i)
{
out_factor[i] = gsl_vector_get(c, i);
}
gsl_vector_free(vY);
gsl_matrix_free(XX);
gsl_matrix_free(cov);
gsl_vector_free(c);
return GSLState2ErrorCode(r);
}
Point3 crossProduct(const Point3& a, const Point3& b)
{
return Point3{
a.y * b.z - a.z * b.y, // C_x
a.z * b.x - a.x * b.z, // C_y
a.x * b.y - a.y * b.x // C_z
};
}
// 计算绕任意轴旋转的旋转矩阵
Eigen::Matrix3d rotationMatrix(const Eigen::Vector3d& axis, double angle) {
// 确保旋转轴是单位向量
Eigen::Vector3d u = axis.normalized();
double cos_theta = cos(angle);
double sin_theta = sin(angle);
// 计算旋转矩阵
Eigen::Matrix3d R;
R <<
cos_theta + u.x() * u.x() * (1 - cos_theta), u.x()* u.y()* (1 - cos_theta) - u.z() * sin_theta, u.x()* u.z()* (1 - cos_theta) + u.y() * sin_theta,
u.y()* u.x()* (1 - cos_theta) + u.z() * sin_theta, cos_theta + u.y() * u.y() * (1 - cos_theta), u.y()* u.z()* (1 - cos_theta) - u.x() * sin_theta,
u.z()* u.x()* (1 - cos_theta) - u.y() * sin_theta, u.z()* u.y()* (1 - cos_theta) + u.x() * sin_theta, cos_theta + u.z() * u.z() * (1 - cos_theta);
return R;
}
long double convertToMilliseconds(const std::string& dateTimeStr) {
// 定义日期时间字符串
std::string dateTimeString = dateTimeStr;
// 使用 Boost.Date_Time 解析日期时间字符串
boost::posix_time::ptime dateTime = boost::posix_time::time_from_string(dateTimeString);
// 将 ptime 转换为自 epoch (1970年1月1日) 以来的毫秒数
boost::posix_time::ptime epoch(boost::gregorian::date(1970, 1, 1));
boost::posix_time::time_duration duration = dateTime - epoch;
long double milliseconds = duration.total_milliseconds() / 1000.0;
return milliseconds;
}
long FindValueInStdVector(std::vector<double>& list, double& findv)
{
if (list.size() == 0) {
return -1;
}
else if (list.size() == 1) {
if (abs(list[0] - findv) < 1e-9) {
return 0;
}
else {
return -1;
}
}
else {}
if (abs(list[0] - findv) < 1e-9) {
return 0;
}
else if (abs(list[list.size() - 1] - findv) < 1e-9) {
return list.size() - 1;
}
else if (list[0] > findv) {
return -1;
}
else if (list[list.size() - 1] < findv) {
return -1;
}
else {}
long low = 0;
long high = list.size() - 1;
while (low <= high) {
long mid = (low + high) / 2;
if (abs(list[mid] - findv) < 1e-9) {
return mid;
}
else if (list[mid] < findv) {
low = mid + 1;
}
else if (list[mid] > findv) {
high = mid - 1;
}
else {}
}
return -1;
}
long InsertValueInStdVector(std::vector<double>& list, double insertValue, bool repeatValueInsert)
{
if (list.size() == 0) {
list.insert(list.begin(), insertValue);
return 0;
}
else if (list.size() == 1) {
if (std::abs(list[0] - insertValue) < PRECISIONTOLERANCE) {
if (repeatValueInsert) {
list.push_back(insertValue);
return 1;
}
else {
return -1;
}
}
else if (insertValue > list[0]) {
list.push_back(insertValue);
return 1;
}
else if (insertValue < list[0]) {
list.push_back(insertValue);
return 0;
}
}
else {
long low = 0;
long high = list.size() - 1;
long mid = -1;
// 处理边界
if (list[high] <= insertValue)
{
if (std::abs(list[high] - insertValue) < PRECISIONTOLERANCE) {
if (repeatValueInsert) {
list.push_back(insertValue);
}
else {}
}
else {
list.push_back(insertValue);
}
return list.size() - 1;
}
if (list[low] >= insertValue) {
if (std::abs(list[low] - insertValue) < PRECISIONTOLERANCE) {
if (repeatValueInsert) {
list.insert(list.begin(), insertValue);
}
else {}
}
else {
list.insert(list.begin(), insertValue);
}
return 0;
}
// low<insertValue < high
while ((high - low) > 1) {
mid = (low + high) / 2;
if (std::abs(list[mid] - insertValue) < PRECISIONTOLERANCE) {
if (repeatValueInsert) {
list.insert(list.begin() + mid, insertValue);
}
return mid;
}
else if (insertValue < list[mid]) {
high = mid;
}
else if (list[mid] < insertValue) {
low = mid;
}
}
if (list[low] <= insertValue && insertValue <= list[high])
{
if (std::abs(list[high] - insertValue) < PRECISIONTOLERANCE) {
if (repeatValueInsert) {
list.insert(list.begin() + high, insertValue);
}
else {}
}
else {
list.insert(list.begin() + high, insertValue);
}
return low;
}
else {
return -1;
}
}
return -1;
}
long FindValueInStdVectorLast(std::vector<double>& list, double& insertValue)
{
if (list.size() == 0 || list.size() == 1) {
return -1;
}
else {
long low = 0;
long high = list.size() - 1;
long mid = -1;
// 处理边界
if (list[high]+ PRECISIONTOLERANCE < insertValue)
{
return -1;
}
else if (std::abs(list[high] - insertValue) < PRECISIONTOLERANCE)
{
return high;
}
else if (list[low] > insertValue) {
return -1;
}
else if (std::abs(list[low] - insertValue) < PRECISIONTOLERANCE) {
return low;
}
else {}
// low<insertValue < high
while ((high - low) > 1) {
mid = (low + high) / 2;
if (insertValue < list[mid]) {
high = mid;
}
else if (list[mid] < insertValue) {
low = mid;
}
else {// insertValue==list[mid] , list[mid]===insertValue
return mid;//
}
}
if (list[low] <= insertValue && insertValue <= list[high])
{
return low;
}
else {
return -1;
}
}
return -1;
}
ErrorCode polynomial_fit(const std::vector<double>& x, const std::vector<double>& y, int degree, std::vector<double>& out_factor, double& out_chisq) {
int xyLength = x.size();
double* xdata = new double[xyLength];
double* ydata = new double[xyLength];
for (long i = 0; i < xyLength; i++) {
xdata[i] = x[i];
ydata[i] = y[i];
}
ErrorCode state = polyfit(xdata, ydata, xyLength, degree, out_factor, out_chisq);
delete[] xdata;
delete[] ydata; // 释放内存
return state;
}
QVector<SatelliteAntPos> SatellitePos2SatelliteAntPos(QVector<SatellitePos> poses)
{
QVector<SatelliteAntPos> antposes(poses.count());
for (long i = 0; i < poses.count(); i++) {
antposes[i].time = poses[i].time;
antposes[i].Px = poses[i].Px;
antposes[i].Py = poses[i].Py;
antposes[i].Pz = poses[i].Pz;
antposes[i].Vx = poses[i].Vx;
antposes[i].Vy = poses[i].Vy;
antposes[i].Vz = poses[i].Vz;
}
return antposes;
}
QVector<SatellitePos> SatelliteAntPos2SatellitePos(QVector<SatelliteAntPos> poses)
{
QVector<SatellitePos> antposes(poses.count());
for (long i = 0; i < poses.count(); i++) {
antposes[i].time = poses[i].time;
antposes[i].Px = poses[i].Px;
antposes[i].Py = poses[i].Py;
antposes[i].Pz = poses[i].Pz;
antposes[i].Vx = poses[i].Vx;
antposes[i].Vy = poses[i].Vy;
antposes[i].Vz = poses[i].Vz;
}
return antposes;
}
QString getDebugDataPath(QString filename)
{
QString folderName = "debugdata";
QString appDir = QCoreApplication::applicationDirPath();
QString folderpath = JoinPath(appDir, folderName);
if (!QDir(folderpath).exists()) {
QDir(appDir).mkdir(folderName);
}
QString datapath = JoinPath(folderpath, filename);
QFile datafile(datapath);
if (datafile.exists()) {
datafile.remove();
}
return datapath;
}
std::vector<std::string> split(const std::string& str, char delimiter) {
std::vector<std::string> tokens;
std::string token;
std::istringstream tokenStream(str);
while (std::getline(tokenStream, token, delimiter)) {
tokens.push_back(token);
}
return tokens;
}
Eigen::VectorXd linspace(double start, double stop, int num) {
Eigen::VectorXd result(num);
double step = (stop - start) / (num - 1); // 计算步长
for (int i = 0; i < num; ++i) {
result[i] = start + i * step; // 生成等间距数值
}
return result;
}

117
BaseTool/BaseTool.h Normal file
View File

@ -0,0 +1,117 @@
#pragma once
#ifndef BASETOOL_H
#define BASETOOL_H
#include "BaseConstVariable.h"
///
/// 基本类、基本函数
///
// //#include <mkl.h>
#include <complex>
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#include <string>
#include <omp.h>
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <ogrsf_frmts.h>
#include <fstream>
#include "GeoOperator.h"
#include <vector>
#include <string>
#include <QString>
#include <QStringList>
#include "LogInfoCls.h"
///////////////////////////////////// 运行时间打印
/////////////////////////////////////////////////////////////
QString getCurrentTimeString();
QString getCurrentShortTimeString();
std::vector<QString> splitString(const QString& str, char delimiter);
std::vector<QString> convertQStringListToStdVector(const QStringList& qStringList);
/////////////////////////////// 基本图像类 结束
/////////////////////////////////////////////////////////////
std::string Convert(float Num);
QString JoinPath(const QString& path, const QString& filename);
////////////////////////////// 坐标部分基本方法 //////////////////////////////////////////
////////////////////////////// 坐标部分基本方法 //////////////////////////////////////////
////////////////////////////// 插值 ////////////////////////////////////////////
std::complex<double> Cubic_Convolution_interpolation(double u, double v,
Eigen::MatrixX<std::complex<double>> img);
std::complex<double> Cubic_kernel_weight(double s);
double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12,Landpoint p22);
bool onSegment(Point3 Pi, Point3 Pj, Point3 Q);
Point3 invBilinear(Point3 p, Point3 a, Point3 b, Point3 c, Point3 d);
//
// WGS84 到J2000 坐标系的变换
// 参考网址：https://blog.csdn.net/hit5067/article/details/116894616
// 资料网址：http://celestrak.org/spacedata/
// 参数文件：
// a. Earth Orientation Parameter 文件： http://celestrak.org/spacedata/EOP-Last5Years.csv
// b. Space Weather Data 文件： http://celestrak.org/spacedata/SW-Last5Years.csv
// 备注：上述文件是自2017年-五年内
/**
在wgs84 坐标系转到J2000 坐标系 主要 涉及到坐标的相互转换。一般给定的WGS坐标为 给定时刻的 t ,BLH
转换步骤：
step 1: WGS 84 转换到协议地球坐标系
step 2: 协议地球坐标系 转换为瞬时地球坐标系
step 3: 瞬时地球坐标系 转换为 瞬时真天球坐标系
step 4: 瞬时真天球坐标系 转到瞬时平天球 坐标系
step 5: 瞬时平天球坐标系转换为协议天球坐标系（J2000）
**/
double sind(double degree);
double cosd(double d);
// 插值
ErrorCode polyfit(const double* x, const double* y, int xyLength, int poly_n, std::vector<double>& out_factor, double& out_chisq);
// 叉乘
Point3 crossProduct(const Point3& a, const Point3& b);
Eigen::Matrix3d rotationMatrix(const Eigen::Vector3d& axis, double angle);
long double convertToMilliseconds(const std::string& dateTimeStr);
/// <summary>
/// list 应该是按照从小到大的顺序排好
/// </summary>
/// <param name="list"></param>
/// <param name="findv"></param>
/// <returns></returns>
long FindValueInStdVector(std::vector<double>& list,double& findv);
long InsertValueInStdVector(std::vector<double>& list, double insertValue, bool repeatValueInsert = false);
long FindValueInStdVectorLast(std::vector<double>& list, double& findv);
ErrorCode polynomial_fit(const std::vector<double>& x, const std::vector<double>& y, int degree, std::vector<double>& out_factor, double& out_chisq);
QVector<SatelliteAntPos> SatellitePos2SatelliteAntPos(QVector<SatellitePos> poses);
QVector<SatellitePos> SatelliteAntPos2SatellitePos(QVector<SatelliteAntPos> poses);
QString getDebugDataPath(QString filename);
std::vector<std::string> split(const std::string& str, char delimiter);
Eigen::VectorXd linspace(double start, double stop, int num);
#endif

634
BaseTool/EchoDataFormat.cpp Normal file
View File

@ -0,0 +1,634 @@
#include "stdafx.h"
#include "EchoDataFormat.h"
#include <QFile>
#include <QFileInfo>
#include <QDir>
#include <QString>
#include <QFile>
#include <QXmlStreamWriter>
#include <QXmlStreamReader>
#include <QDebug>
#include "ImageOperatorBase.h"
std::shared_ptr<PluseAntPos> CreatePluseAntPosArr(long pluseCount)
{
return std::shared_ptr<PluseAntPos>(new PluseAntPos[pluseCount],delArrPtr);
}
long getPluseDataSize(PluseData& pluseData)
{
long datasize = sizeof(long) + sizeof(double) + sizeof(double) * 6 + sizeof(long);// PRFId,time,px-vz,pluseCount
datasize = datasize + pluseData.plusePoints * sizeof(double) * 2; // 数据块
return datasize;
}
ErrorCode getPluseDataFromBuffer(char* buffer, PluseData & data)
{
long seekid = 0;
std::memcpy(&data.id, buffer + seekid, sizeof(data.id)); seekid += sizeof(data.id);
std::memcpy(&data.time, buffer + seekid, sizeof(data.time)); seekid += sizeof(data.time);
std::memcpy(&data.Px, buffer + seekid, sizeof(data.Px)); seekid += sizeof(data.Px);
std::memcpy(&data.Py, buffer + seekid, sizeof(data.Py)); seekid += sizeof(data.Py);
std::memcpy(&data.Pz, buffer + seekid, sizeof(data.Pz)); seekid += sizeof(data.Pz);
std::memcpy(&data.Vx, buffer + seekid, sizeof(data.Vx)); seekid += sizeof(data.Vx);
std::memcpy(&data.Vy, buffer + seekid, sizeof(data.Vy)); seekid += sizeof(data.Vy);
std::memcpy(&data.Vz, buffer + seekid, sizeof(data.Vz)); seekid += sizeof(data.Vz);
std::memcpy(&data.plusePoints, buffer + seekid, sizeof(data.plusePoints)); seekid += sizeof(data.plusePoints);
Eigen::MatrixXd echoData_real = Eigen::MatrixXd::Zero(1, data.plusePoints);
Eigen::MatrixXd echoData_imag = Eigen::MatrixXd::Zero(1, data.plusePoints);
std::memcpy(echoData_real.data(), buffer + seekid, sizeof(data.plusePoints)); seekid += data.plusePoints * sizeof(double);
std::memcpy(echoData_imag.data(), buffer + seekid, sizeof(data.plusePoints));
std::shared_ptr<Eigen::MatrixXcd> echoData = std::make_shared<Eigen::MatrixXcd>(1, data.plusePoints);
echoData->real() = echoData_real.array();
echoData->imag() = echoData_imag.array();
return ErrorCode::SUCCESS;
}
std::shared_ptr<PluseData> CreatePluseDataArr(long pluseCount)
{
return std::shared_ptr<PluseData>(new PluseData[pluseCount],delArrPtr);
}
std::shared_ptr<std::complex<double>> CreateEchoData(long plusePoints)
{
return std::shared_ptr<std::complex<double>>(new std::complex<double>[plusePoints],delArrPtr);
}
EchoL0Dataset::EchoL0Dataset()
{
this->folder="";
this->filename="";
this->xmlname="";
this->GPSPointFilename="";
this->echoDataFilename="";
this->xmlFilePath="";
this->GPSPointFilePath="";
this->echoDataFilePath="";
this->simulationTaskName = "";
this->PluseCount = 0;
this->PlusePoints = 0;
this->NearRange = 0;
this->FarRange = 0;
this->centerFreq = 0;
this->Fs = 0;
this->folder.clear();
this->filename.clear();
this->xmlname.clear();
this->GPSPointFilename.clear();
this->echoDataFilename.clear();
this->xmlFilePath.clear();
this->GPSPointFilePath.clear();
this->echoDataFilePath.clear();
this->simulationTaskName.clear();
}
EchoL0Dataset::~EchoL0Dataset()
{
}
ErrorCode EchoL0Dataset::OpenOrNew(QString folder, QString filename, long PluseCount, long PlusePoints)
{
qDebug() << "--------------Echo Data OpenOrNew ---------------------------------------";
qDebug() << "Folder: " << folder;
qDebug() << "Filename: " << filename;
QDir dir(folder);
if (dir.exists() == false)
{
dir.mkpath(folder);
}
else {}
if (dir.exists() == false) {
return ErrorCode::FOLDER_NOT_EXIST;
}
else {}
QString filePath = dir.filePath(filename); // 生成完整的文件路径
this->folder = folder;
this->filename = filename;
this->simulationTaskName = filename;
this->xmlname=filename+".xml";
this->GPSPointFilename=filename+".gpspos.data";
this->echoDataFilename = filename + ".L0echo.data";
this->xmlFilePath = dir.filePath(this->xmlname);
this->GPSPointFilePath = dir.filePath(this->GPSPointFilename);
this->echoDataFilePath = dir.filePath(this->echoDataFilename);
this->PluseCount = PluseCount;
this->PlusePoints = PlusePoints;
//
if (QFile(this->xmlFilePath).exists())
{
this->loadFromXml();
}
else
{
this->saveToXml();
}
if (QFile(this->GPSPointFilePath).exists() == false) {
// 创建新文件
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 19, PluseCount, 1, GDT_Float64, NULL));
GDALFlushCache((GDALDatasetH)poDstDS.get());
poDstDS.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
if (QFile(this->echoDataFilePath).exists() == false) {
// 创建新文件
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(this->echoDataFilePath.toUtf8().constData(), PlusePoints, PluseCount, 1, GDT_CFloat64, NULL));
GDALFlushCache((GDALDatasetH)poDstDS.get());
poDstDS.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
return ErrorCode::SUCCESS;
}
ErrorCode EchoL0Dataset::Open(QString xmlfilepath)
{
QFileInfo fileInfo(xmlfilepath);
QString fileName = fileInfo.fileName(); // 获取文件名
QString fileSuffix = fileInfo.suffix(); // 获取后缀名
QString fileNameWithoutSuffix = fileInfo.baseName(); // 获取不带后缀的文件名
QString directoryPath = fileInfo.path(); // 获取文件夹路径
if (fileSuffix.toLower() == "xml" || fileSuffix.toLower() == ".xml") {
return this->Open(directoryPath,fileNameWithoutSuffix);
}
else {
return ErrorCode::ECHO_L0DATA_XMLNAMEERROR;
}
return ErrorCode::SUCCESS;
}
ErrorCode EchoL0Dataset::Open(QString folder, QString filename)
{
QDir dir(folder);
if (dir.exists() == false)
{
return ErrorCode::FOLDER_NOT_EXIST;
}
else {}
if (dir.exists() == false) {
return ErrorCode::FOLDER_NOT_EXIST;
}
else {}
QString filePath = dir.filePath(filename); // 生成完整的文件路径
this->folder = folder;
this->filename = filename;
this->simulationTaskName = filename;
this->xmlname = filename + ".xml";
this->GPSPointFilename = filename + ".gpspos.data";
this->echoDataFilename = filename + ".L0echo.data";
this->xmlFilePath = dir.filePath(this->xmlname);
this->GPSPointFilePath = dir.filePath(this->GPSPointFilename);
this->echoDataFilePath = dir.filePath(this->echoDataFilename);
this->PluseCount = PluseCount;
this->PlusePoints = PlusePoints;
//
if (QFile(this->xmlFilePath).exists())
{
this->loadFromXml();
}
else
{
return ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR;
}
}
QString EchoL0Dataset::getxmlName()
{
return xmlname;
}
QString EchoL0Dataset::getGPSPointFilename()
{
return GPSPointFilename;
}
QString EchoL0Dataset::getEchoDataFilename()
{
return GPSPointFilePath;
}
void EchoL0Dataset::initEchoArr(std::complex<double> init0)
{
long blockline = Memory1MB * 2000 / 8 / 2 / this->PlusePoints;
long start = 0;
for (start = 0; start < this->PluseCount; start = start + blockline) {
long templine = start + blockline < this->PluseCount ? blockline : this->PluseCount - start;
std::shared_ptr<std::complex<double>> echotemp = this->getEchoArr(start, templine);
for (long i = 0; i < templine; i++) {
for (long j = 0; j < this->PlusePoints; j++) {
echotemp.get()[i * this->PlusePoints + j] = init0;
}
}
this->saveEchoArr(echotemp, start, templine);
qDebug() << "echo init col : " << start << "\t-\t" << start + templine << "\t/\t" << this->PluseCount;
}
}
// Getter 和 Setter 方法实现
long EchoL0Dataset::getPluseCount() { return this->PluseCount; }
void EchoL0Dataset::setPluseCount(long pulseCount) { this->PluseCount = pulseCount; }
long EchoL0Dataset::getPlusePoints() { return this->PlusePoints; }
void EchoL0Dataset::setPlusePoints(long pulsePoints) { this->PlusePoints = pulsePoints; }
double EchoL0Dataset::getNearRange() { return this->NearRange; }
void EchoL0Dataset::setNearRange(double nearRange) { this->NearRange = nearRange; }
double EchoL0Dataset::getFarRange() { return this->FarRange; }
void EchoL0Dataset::setFarRange(double farRange) { this->FarRange = farRange; }
double EchoL0Dataset::getCenterFreq() { return this->centerFreq; }
void EchoL0Dataset::setCenterFreq(double freq) { this->centerFreq = freq; }
double EchoL0Dataset::getFs() { return this->Fs; }
void EchoL0Dataset::setFs(double samplingFreq) { this->Fs = samplingFreq; }
QString EchoL0Dataset::getSimulationTaskName() { return this->simulationTaskName; }
void EchoL0Dataset::setSimulationTaskName(const QString& taskName) { this->simulationTaskName = taskName; }
double EchoL0Dataset::getCenterAngle()
{
return this->CenterAngle;
}
void EchoL0Dataset::setCenterAngle(double angle)
{
this->CenterAngle = angle;
}
QString EchoL0Dataset::getLookSide()
{
return this->LookSide;
}
void EchoL0Dataset::setLookSide(QString lookside)
{
this->LookSide = lookside;
}
SatelliteAntPos EchoL0Dataset::getSatelliteAntPos(long prf_id)
{
std::shared_ptr<double> antpos = this->getAntPos();
SatelliteAntPos prfpos{};
prfpos.time = antpos.get()[prf_id *19 + 0];
prfpos.Px = antpos.get()[prf_id *19 + 1];
prfpos.Py = antpos.get()[prf_id *19 + 2];
prfpos.Pz = antpos.get()[prf_id *19 + 3];
prfpos.Vx = antpos.get()[prf_id *19 + 4];
prfpos.Vy = antpos.get()[prf_id *19 + 5];
prfpos.Vz = antpos.get()[prf_id *19 + 6];
prfpos.AntDirectX = antpos.get()[prf_id *19 + 7];
prfpos.AntDirectY = antpos.get()[prf_id *19 + 8];
prfpos.AntDirectZ = antpos.get()[prf_id *19 + 9];
prfpos.AVx = antpos.get()[prf_id *19 + 10];
prfpos.AVy = antpos.get()[prf_id *19 + 11];
prfpos.AVz =antpos.get()[prf_id *19 + 12];
prfpos.ZeroAntDiectX = antpos.get()[prf_id *19 + 13];
prfpos.ZeroAntDiectY = antpos.get()[prf_id *19 + 14];
prfpos.ZeroAntDiectZ = antpos.get()[prf_id *19 + 15];
prfpos.lon =antpos.get()[prf_id *19 + 16];
prfpos.lat =antpos.get()[prf_id *19 + 17];
prfpos.ati =antpos.get()[prf_id *19 + 18];
return prfpos;
}
// 打印信息的实现
void EchoL0Dataset::printInfo() {
std::cout << "Simulation Task Name: " << this->simulationTaskName.toStdString() << std::endl;
std::cout << "Pulse Count: " << this->PluseCount << std::endl;
std::cout << "Pulse Points: " << this->PlusePoints << std::endl;
std::cout << "Near Range: " << this->NearRange << std::endl;
std::cout << "Far Range: " << this->FarRange << std::endl;
std::cout << "Center Frequency: " << this->centerFreq << std::endl;
std::cout << "Sampling Frequency: " << this->Fs << std::endl;
}
// xml文件读写
void EchoL0Dataset::saveToXml() {
QString filePath = this->xmlFilePath;
QFile file(filePath);
if (!file.open(QIODevice::WriteOnly)) {
qWarning() << "Cannot open file for writing:" << filePath;
return;
}
QXmlStreamWriter xmlWriter(&file);
xmlWriter.setAutoFormatting(true);
xmlWriter.writeStartDocument();
xmlWriter.writeStartElement("SimulationConfig");
xmlWriter.writeTextElement("PluseCount", QString::number(this->PluseCount));
xmlWriter.writeTextElement("PlusePoints", QString::number(this->PlusePoints));
xmlWriter.writeTextElement("NearRange", QString::number(this->NearRange));
xmlWriter.writeTextElement("FarRange", QString::number(this->FarRange));
xmlWriter.writeTextElement("CenterFreq", QString::number(this->centerFreq));
xmlWriter.writeTextElement("Fs", QString::number(this->Fs));
xmlWriter.writeTextElement("CenterAngle", QString::number(this->CenterAngle));
xmlWriter.writeTextElement("LookSide", this->LookSide);
xmlWriter.writeTextElement("SimulationTaskName", this->simulationTaskName);
xmlWriter.writeTextElement("Filename", this->filename);
xmlWriter.writeTextElement("Xmlname", this->xmlname);
xmlWriter.writeTextElement("GPSPointFilename", this->GPSPointFilename);
xmlWriter.writeTextElement("EchoDataFilename", this->echoDataFilename);
xmlWriter.writeEndElement(); // SimulationConfig
xmlWriter.writeEndDocument();
file.close();
}
ErrorCode EchoL0Dataset::loadFromXml() {
QString filePath = this->xmlFilePath;
QFile file(filePath);
if (!file.open(QIODevice::ReadOnly)) {
qWarning() << "Cannot open file for reading:" << filePath;
return ErrorCode::FILEOPENFAIL;
}
bool PluseCountflag = false;
bool PlusePointsflag = false;
bool NearRangeflag = false;
bool FarRangeflag = false;
bool CenterFreqflag = false;
bool Fsflag = false;
QXmlStreamReader xmlReader(&file);
while (!xmlReader.atEnd() && !xmlReader.hasError()) {
xmlReader.readNext();
if (xmlReader.isStartElement()) {
QString elementName = xmlReader.name().toString();
if (elementName == "PluseCount") {
this->PluseCount = xmlReader.readElementText().toLong();
PluseCountflag = true;
}
else if (elementName == "PlusePoints") {
this->PlusePoints = xmlReader.readElementText().toLong();
PlusePointsflag = true;
}
else if (elementName == "NearRange") {
this->NearRange = xmlReader.readElementText().toDouble();
NearRangeflag = true;
}
else if (elementName == "FarRange") {
this->FarRange = xmlReader.readElementText().toDouble();
FarRangeflag = true;
}
else if (elementName == "CenterFreq") {
this->centerFreq = xmlReader.readElementText().toDouble();
CenterFreqflag = true;
}
else if (elementName == "Fs") {
this->Fs = xmlReader.readElementText().toDouble();
Fsflag = true;
}
else if (elementName == "SimulationTaskName") {
this->simulationTaskName = xmlReader.readElementText();
}
else if (elementName == "Filename") {
this->filename = xmlReader.readElementText();
}
else if (elementName == "Xmlname") {
this->xmlname = xmlReader.readElementText();
}
else if (elementName == "GPSPointFilename") {
this->GPSPointFilename = xmlReader.readElementText();
}
else if (elementName == "EchoDataFilename") {
this->echoDataFilename = xmlReader.readElementText();
}
else if (elementName == "CenterAngle") {
this->CenterAngle = xmlReader.readElementText().toDouble();
}
else if (elementName == "LookSide") {
this->LookSide = xmlReader.readElementText();
}
}
}
if (xmlReader.hasError()) {
qWarning() << "XML error:" << xmlReader.errorString();
file.close();
return ErrorCode::ECHO_L0DATA_XMLFILENOTOPEN;
}
if (!(PlusePointsflag && PlusePointsflag && FarRangeflag && NearRangeflag && CenterFreqflag && Fsflag)) {
file.close();
return ErrorCode::ECHO_L0DATA_XMLFILENOTOPEN;
}
file.close();
return ErrorCode::SUCCESS;
}
std::shared_ptr<double> EchoL0Dataset::getAntPos()
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
// 读取文件
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(this->GPSPointFilePath, GDALAccess::GA_ReadOnly);
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(1);
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
std::shared_ptr<double> temp = nullptr;
if (gdal_datatype == GDT_Float64) {
temp=std::shared_ptr<double>(new double[this->PluseCount * 19],delArrPtr);
demBand->RasterIO(GF_Read, 0, 0, 19, this->PluseCount, temp.get(), 19, this->PluseCount, gdal_datatype, 0, 0);
}
else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)) ;
}
rasterDataset.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return temp;
}
std::shared_ptr<std::complex<double>> EchoL0Dataset::getEchoArr(long startPRF, long PRFLen)
{
if (!(startPRF < this->PluseCount)) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE))<<startPRF<<" "<<this->PluseCount;
return nullptr;
}
else {}
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
std::shared_ptr<GDALDataset> 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 nullptr;
}
else {}
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
std::shared_ptr<std::complex<double>> temp = nullptr;
if (height != this->PluseCount || width != this->PlusePoints) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
}
else {
if (gdal_datatype == GDT_CFloat64) {
temp= std::shared_ptr<std::complex<double>>(new std::complex<double>[PRFLen * width ], delArrPtr);
poBand->RasterIO(GF_Read, 0, startPRF, width, PRFLen, temp.get(), width, PRFLen, GDT_CFloat64, 0, 0);
GDALFlushCache((GDALDatasetH)rasterDataset.get());
}
else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
}
}
rasterDataset.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return temp;
}
std::shared_ptr<std::complex<double>> EchoL0Dataset::getEchoArr()
{
return this->getEchoArr(0,this->PluseCount);
}
ErrorCode EchoL0Dataset::saveAntPos(std::shared_ptr<double> ptr)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
// 读取文件
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(this->GPSPointFilePath, GDALAccess::GA_Update);
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(1);
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
if (gdal_datatype == GDT_Float64) {
demBand->RasterIO(GF_Write, 0, 0, 19, this->PluseCount, ptr.get(), 19, this->PluseCount, gdal_datatype, 0, 0);
}
else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR));
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR;
}
rasterDataset.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return ErrorCode::SUCCESS;
}
ErrorCode EchoL0Dataset::saveEchoArr(std::shared_ptr<std::complex<double>> echoPtr, long startPRF, long PRFLen)
{
if (!(startPRF < this->PluseCount)) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE));
return ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE;
}
else {}
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(this->echoDataFilePath, GDALAccess::GA_Update);
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 ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR;
}
else {}
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
if (height != this->PluseCount || width != this->PlusePoints) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR;
}
else {
if (gdal_datatype == GDT_CFloat64) {
poBand->RasterIO(GF_Write, 0, startPRF, width, PRFLen, echoPtr.get(), width, PRFLen, GDT_CFloat64, 0, 0);
GDALFlushCache((GDALDatasetH)rasterDataset.get());
}
else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
}
}
rasterDataset.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return ErrorCode::SUCCESS;
}

193
BaseTool/EchoDataFormat.h Normal file
View File

@ -0,0 +1,193 @@
#pragma once
/*****************************************************************//**
* \file EchoDataFormat.h
* \brief 存储仿真过程中的各级产品数据格式，主要包含回波数据格式、仿真参数数据库
*
* \author 30453
* \date October 2024
*********************************************************************/
#include "BaseConstVariable.h"
#include <complex>
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <omp.h>
#include <memory>
#include <vector>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <QString>
#include "BaseTool.h"
#include <iostream>
#include <fstream>
#include <complex>
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <string>
#include <memory>
#include <QString>
#include <cpl_conv.h> // for CPLMalloc()
//========================================================================
// 成像回波格式
// file type:
// PRFCOUNT
// PRFPOINT
// nearRange
// farRange
// PRF1,time,Px,Py,Pz,Vx,Vy,Vz
// PRF2,time,Px,Py,Pz,Vx,Vy,Vz
// ------------------ 脉冲文件部分 --------------------------------------
// PRF1,time,Px,Py,Pz,Vx,Vy,Vz,PRFPOINT,RealData,imagData
// PRF2,time,Px,Py,Pz,Vx,Vy,Vz,PRFPOINT,RealData,imagData
// --------------------------------------------------------------
// 注意Bp并不关心脉冲的顺序，只是关注脉冲的坐标位置，默认按照脉冲的解析顺序进行组织
//========================================================================
// 成像用的回波脉冲文件格式
/// <summary>
/// 回波数据--废弃
/// </summary>
struct PluseData {
long id; // PRF id
double time;// 时间
double Px; // 坐标
double Py;
double Pz;
double Vx; // 速度
double Vy;
double Vz;
long plusePoints; // 脉冲点数
std::shared_ptr<std::complex<double>> echoData; // 回波数据
};
long getPluseDataSize(PluseData& pluseData);
ErrorCode getPluseDataFromBuffer(char* buffer, PluseData& data);
std::shared_ptr<PluseData> CreatePluseDataArr(long pluseCount);
std::shared_ptr<std::complex<double>> CreateEchoData(long plusePoints);
/// <summary>
/// 姿态数据
/// </summary>
struct PluseAntPos {
long id; // PRF id
double time;// 时间
double Px; // 坐标
double Py;
double Pz;
double Vx; // 速度
double Vy;
double Vz;
};
std::shared_ptr<PluseAntPos> CreatePluseAntPosArr(long pluseCount);
// 定义L0级数据
class EchoL0Dataset {
public:
EchoL0Dataset();
~EchoL0Dataset();
public:
ErrorCode OpenOrNew(QString folder, QString filename,long PluseCount,long PlusePoints);
ErrorCode Open(QString xmlfilepath);
ErrorCode Open(QString folder, QString filename);
QString getxmlName();
QString getGPSPointFilename();
QString getEchoDataFilename();
void initEchoArr(std::complex<double> init0);
private: // 产品名称设置
QString folder;
QString filename;
QString xmlname;
QString GPSPointFilename;
QString echoDataFilename;
//
QString xmlFilePath;
QString GPSPointFilePath;
QString echoDataFilePath;
public: // 文件处理部分
// Getter 和 Setter 方法
long getPluseCount() ;
void setPluseCount(long pulseCount);
long getPlusePoints() ;
void setPlusePoints(long pulsePoints);
double getNearRange() ;
void setNearRange(double nearRange);
double getFarRange() ;
void setFarRange(double farRange);
double getCenterFreq() ;
void setCenterFreq(double freq);
double getFs() ;
void setFs(double samplingFreq);
QString getSimulationTaskName() ;
void setSimulationTaskName(const QString& taskName);
double getCenterAngle();
void setCenterAngle(double angle);
QString getLookSide();
void setLookSide(QString lookside);
SatelliteAntPos getSatelliteAntPos(long plusePRFID);
// 打印信息的成员函数
void printInfo() ;
private: // 回波变量
long PluseCount;
long PlusePoints;
double NearRange;
double FarRange;
double centerFreq; // 入射中心频率
double Fs; // 等效采样频率
QString simulationTaskName;
double CenterAngle;
QString LookSide;
public: // 读写 XML 的函数
void saveToXml();
ErrorCode loadFromXml();
public:
// 读取文件
std::shared_ptr<double> getAntPos();
std::shared_ptr<std::complex<double>> getEchoArr(long startPRF, long PRFLen);
std::shared_ptr<std::complex<double>> getEchoArr();
//保存文件
ErrorCode saveAntPos(std::shared_ptr<double> ptr); // 注意这个方法很危险，请写入前检查数据是否正确
ErrorCode saveEchoArr(std::shared_ptr<std::complex<double>> echoPtr, long startPRF, long PRFLen);
};

253
BaseTool/FileOperator.cpp Normal file
View File

@ -0,0 +1,253 @@
#include "stdafx.h"
#include "FileOperator.h"
#include <boost/filesystem.hpp>
#include <string>
#include <memory.h>
#include <memory>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <QDebug>
#include <QString>
#include <QMessageBox>
QString addMaskToFileName(const QString& filePath,QString _zzname) {
// 获取文件路径和文件名
QFileInfo fileInfo(filePath);
// 获取文件名和扩展名
QString baseName = fileInfo.baseName(); // 不包含扩展名的文件名
QString extension = fileInfo.suffix(); // 文件扩展名（例如 ".txt", ".jpg"）
// 拼接新的文件名
QString newFileName = baseName + _zzname; // 在文件名中增加 "_mask"
if (!extension.isEmpty()) {
newFileName += "." + extension; // 如果有扩展名，添加扩展名
}
// 返回新的文件路径
QString newFilePath = fileInfo.path() + "/" + newFileName;
return newFilePath;
}
std::vector<QString> getFilelist(const QString& folderpath, const QString& filenameExtension, int (*logfun)(QString logtext, int value))
{
QString filenameExtensionStr = filenameExtension;
filenameExtensionStr = filenameExtensionStr.remove(0, 1);
std::vector<QString> filenames(0);
if (isExists(folderpath) && isDirectory(folderpath)) {
QDir directory(folderpath);
if (directory.exists() && directory.isReadable()) {
QFileInfoList fileList = directory.entryInfoList(QDir::Files | QDir::NoDotAndDotDot);
for (const QFileInfo& fileInfo : fileList) {
// qDebug() << fileInfo.filePath() + "\nExtension: (" + filenameExtensionStr + ", " + fileInfo.suffix() + ")";
if (filenameExtensionStr == u8"*" || filenameExtensionStr == fileInfo.suffix()) {
filenames.push_back(fileInfo.filePath());
}
if (logfun) {
logfun(fileInfo.filePath() + "\nExtension: (" + filenameExtensionStr + ", " + fileInfo.suffix() + ")", 4);
}
}
}
else {
qWarning() << "Folder does not exist or is not readable: " << folderpath;
}
return filenames;
}
else {
return std::vector<QString>(0);
}
}
QString getParantFolderNameFromPath(const QString& path)
{
QDir directory(path);
directory.cdUp();
QString parentPath = directory.absolutePath();
return directory.dirName();
}
QString getParantFromPath(const QString& path)
{
//qDebug() << path;
QDir directory(path);
directory.cdUp();
QString parentPath = directory.absolutePath();
//qDebug() << parentPath;
return parentPath;
}
QString getFileNameFromPath(const QString& path)
{
QFileInfo fileInfo(path);
return fileInfo.fileName();
}
QString getFileNameWidthoutExtend(QString path)
{
QFileInfo fileInfo(path);
QString fileNameWithoutExtension = fileInfo.baseName(); // 获取无后缀文件名
return fileNameWithoutExtension;
}
bool isDirectory(const QString& path)
{
QFileInfo fileinfo(path);
return fileinfo.isDir();
}
bool isExists(const QString& path)
{
QFileInfo fileinfo(path);
return fileinfo.exists();
}
bool isFile(const QString& path)
{
QFileInfo fileinfo(path);
return fileinfo.isFile();
}
int write_binfile(char* filepath, char* data, size_t data_len)
{
FILE* pd = fopen(filepath, "w");
if (NULL == pd) {
return 2;
}
//数据块首地址: "&a"，元素大小: "sizeof(unsigned __int8)"， 元素个数: "10"， 文件指针："pd"
fwrite(data, sizeof(char), data_len, pd);
fclose(pd);
return -1;
}
char* read_textfile(char* text_path, int* length)
{
char* data = NULL;
FILE* fp1 = fopen(text_path, "r");
if (fp1 == NULL) {
return NULL;
}
else {}
// 读取文件
fseek(fp1, 0, SEEK_END);
int data_length = ftell(fp1);
data = (char*)malloc((data_length + 1) * sizeof(char));
rewind(fp1);
if (data_length == fread(data, sizeof(char), data_length, fp1)) {
data[data_length] = '\0'; // 文件尾
}
else {
free(data);
fclose(fp1);
return NULL;
}
fclose(fp1);
*length = data_length + 1;
return data;
}
bool exists_test(const QString& name)
{
return isExists(name);
}
size_t fsize(FILE* fp)
{
size_t n;
fpos_t fpos; // 当前位置
fgetpos(fp, &fpos); // 获取当前位置
fseek(fp, 0, SEEK_END);
n = ftell(fp);
fsetpos(fp, &fpos); // 恢复之前的位置
return n;
}
void removeFile(const QString& filePath)
{
QFile file(filePath);
if (file.exists()) {
if (file.remove()) {
qDebug() << "File removed successfully: " << filePath;
}
else {
qWarning() << "Failed to remove file: " << filePath;
}
}
else {
qDebug() << "File does not exist: " << filePath;
}
}
unsigned long convertToULong(const QString& input) {
bool ok; // Used to check if the conversion was successful
unsigned long result = input.toULong(&ok);
if (!ok) {
qWarning() << "Conversion to unsigned long failed for input: " << input;
}
return result;
}
void copyFile(const QString& sourcePath, const QString& destinationPath) {
QFile sourceFile(sourcePath);
QFile destinationFile(destinationPath);
qDebug() << QString("copy file ready !! from ") + sourcePath+" to "+destinationPath ;
if (sourceFile.exists()) {
if (sourceFile.copy(destinationPath)) {
qDebug() << QString("copy file sucessfully !! from ") + sourcePath+" to "+destinationPath ;
// 复制成功
//QMessageBox::information(nullptr, u8"成功", u8"文件复制成功");
}
else {
// 复制失败
if(sourceFile.exists()){
QMessageBox::critical(nullptr, QObject::tr("error"), QObject::tr("The file already exists !!"));
}
else{
QMessageBox::critical(nullptr, QObject::tr("error"), QObject::tr("file copy error"));
}
}
}
else {
// 源文件不存在
QMessageBox::warning(nullptr, QObject::tr("warning"), QObject::tr("Source file not found"));
}
}
bool copyAndReplaceFile(const QString& sourceFilePath, const QString& destinationFilePath) {
// 检查源文件是否存在
if (!QFile::exists(sourceFilePath)) {
qDebug() << "Source file does not exist:" << sourceFilePath;
return false;
}
// 如果目标文件存在，则删除它
if (QFile::exists(destinationFilePath)) {
if (!QFile::remove(destinationFilePath)) {
qDebug() << "Failed to remove existing destination file:" << destinationFilePath;
return false;
}
}
// 复制文件
if (!QFile::copy(sourceFilePath, destinationFilePath)) {
qDebug() << "Failed to copy file from" << sourceFilePath << "to" << destinationFilePath;
return false;
}
qDebug() << "File copied successfully from" << sourceFilePath << "to" << destinationFilePath;
return true;
}

56
BaseTool/FileOperator.h Normal file
View File

@ -0,0 +1,56 @@
#pragma once
#ifndef FILEOPERATOR_H
#define FILEOPERATOR_H
#include "BaseConstVariable.h"
#include <string.h>
#include <memory.h>
#include <memory>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <string>
#include <vector>
#include <QString>
#include <QDir>
#include <QFile>
#include <QDebug>
bool isDirectory(const QString& path);
bool isExists(const QString& path);
bool isFile(const QString& path);
void removeFile(const QString& filePath);
unsigned long convertToULong(const QString& input);
/// <summary>
/// 获取文件(绝对路径）
/// </summary>
/// <param name="folderpath"></param>
/// <param name="FilenameExtension"></param>
/// <returns></returns>
std::vector<QString> getFilelist(const QString& folderpath, const QString& FilenameExtension = ".*",int (*logfun)(QString logtext,int value)=nullptr);
QString getParantFolderNameFromPath(const QString& path);
QString getFileNameFromPath(const QString& path);
QString getFileNameWidthoutExtend(QString path);
int write_binfile(char* filepath, char* data, size_t data_len);
char* read_textfile(char* text_path, int* length);
bool exists_test(const QString& name);
size_t fsize(FILE* fp);
QString getParantFromPath(const QString& path);
void copyFile(const QString& sourcePath, const QString& destinationPath);
QString addMaskToFileName(const QString& filePath, QString _zzname);
// QT FileOperator
bool copyAndReplaceFile(const QString& sourceFilePath, const QString& destinationFilePath);
#endif

953
BaseTool/GPUTool.cu Normal file
View File

@ -0,0 +1,953 @@
#include <iostream>
#include <memory>
#include <cmath>
#include <complex>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#ifdef __CUDANVCC___
#define CUDAMEMORY Memory1MB*100
#define LAMP_CUDA_PI 3.141592653589793238462643383279
// 定义参数
__device__ cuComplex cuCexpf(cuComplex x)
{
float factor = exp(x.x);
return make_cuComplex(factor * cos(x.y), factor * sin(x.y));
}
// 定义仿真所需参数
__device__ float GPU_getSigma0dB(CUDASigmaParam param, float theta) {
float sigma= param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
return sigma;
}
__device__ CUDAVector GPU_VectorAB(CUDAVector A, CUDAVector B) {
CUDAVector C;
C.x = B.x - A.x;
C.y = B.y - A.y;
C.z = B.z - A.z;
return C;
}
__device__ float GPU_VectorNorm2(CUDAVector A) {
return sqrtf(A.x * A.x + A.y * A.y + A.z * A.z);
}
__device__ float GPU_dotVector(CUDAVector A, CUDAVector B) {
return A.x * B.x + A.y * B.y + A.z * B.z;
}
__device__ float GPU_CosAngle_VectorA_VectorB(CUDAVector A, CUDAVector B) {
return GPU_dotVector(A, B) / (GPU_VectorNorm2(A) * GPU_VectorNorm2(B));
}
__device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(float RstX, float RstY, float RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ
) {
CUDAVectorEllipsoidal result{ 0,0,-1 };
float Xst = -1 * RstX; // 卫星 --> 地面
float Yst = -1 * RstY;
float Zst = -1 * RstZ;
float AntXaxisX = antXaxisX;
float AntXaxisY = antXaxisY;
float AntXaxisZ = antXaxisZ;
float AntYaxisX = antYaxisX;
float AntYaxisY = antYaxisY;
float AntYaxisZ = antYaxisZ;
float AntZaxisX = antZaxisX;
float AntZaxisY = antZaxisY;
float AntZaxisZ = antZaxisZ;
// 天线指向在天线坐标系下的值
float Xant = (Xst * (AntYaxisY * AntZaxisZ - AntYaxisZ * AntZaxisY) + Xst * (AntXaxisZ * AntZaxisY - AntXaxisY * AntZaxisZ) + Xst * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
float Yant = (Yst * (AntYaxisZ * AntZaxisX - AntYaxisX * AntZaxisZ) + Yst * (AntXaxisX * AntZaxisZ - AntXaxisZ * AntZaxisX) + Yst * (AntYaxisX * AntXaxisZ - AntXaxisX * AntYaxisZ)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
float Zant = (Zst * (AntYaxisX * AntZaxisY - AntYaxisY * AntZaxisX) + Zst * (AntXaxisY * AntZaxisX - AntXaxisX * AntZaxisY) + Zst * (AntXaxisX * AntYaxisY - AntYaxisX * AntXaxisY)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
// 计算theta 与 phi
float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
float ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
float YsinTheta = Yant / sinf(ThetaAnt);
float PhiAnt = (YsinTheta / abs(YsinTheta)) * acosf(Xant / (Norm * sinf(ThetaAnt)));
result.theta = ThetaAnt;
result.phi = PhiAnt;
result.pho = Norm;
return result;
}
/**
天线方向图插值方法，以双线性插值算法为基础，由theta与phi组合得到的矩阵图为基础数据，通过插值计算的方法获取目标点的数据。
其中行是theta、列是phi
*/
__device__ float GPU_BillerInterpAntPattern(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
long thetapoints, long phipoints,
float searththeta, float searchphi) {
float stheta = searththeta;
float sphi = searchphi;
if (stheta > 90) {
return 0;
}
else {}
float pthetaid = (stheta - starttheta) / dtheta;//
float pphiid = (sphi - startphi) / dphi;
long lasttheta = floorf(pthetaid);
long nextTheta = lasttheta + 1;
long lastphi = floorf(pphiid);
long nextPhi = lastphi + 1;
if (lasttheta < 0 || nextTheta < 0 || lastphi < 0 || nextPhi < 0 ||
lasttheta >= thetapoints || nextTheta >= thetapoints || lastphi >= phipoints || nextPhi >= phipoints)
{
return 0;
}
else {
float x = stheta;
float y = sphi;
float x1 = lasttheta * dtheta + starttheta;
float x2 = nextTheta * dtheta + starttheta;
float y1 = lastphi * dphi + startphi;
float y2 = nextPhi * dphi + startphi;
float z11 = antpattern[lasttheta * phipoints + lastphi];
float z12 = antpattern[lasttheta * phipoints + nextPhi];
float z21 = antpattern[nextTheta * phipoints + lastphi];
float z22 = antpattern[nextTheta * phipoints + nextPhi];
//z11 = powf(10, z11 / 10); // dB-> 线性
//z12 = powf(10, z12 / 10);
//z21 = powf(10, z21 / 10);
//z22 = powf(10, z22 / 10);
float GainValue = (z11 * (x2 - x) * (y2 - y)
+ z21 * (x - x1) * (y2 - y)
+ z12 * (x2 - x) * (y - y1)
+ z22 * (x - x1) * (y - y1));
GainValue = GainValue / ((x2 - x1) * (y2 - y1));
return GainValue;
}
}
__device__ cuComplex GPU_calculationEcho(float sigma0, float TransAnt, float ReciveAnt,
float localangle, float R, float slopeangle, float Pt, float lamda) {
float r = R;
float amp = Pt * TransAnt * ReciveAnt;
amp = amp * sigma0;
amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // 反射强度
float phi = (-4 * LAMP_CUDA_PI / lamda) * r;
cuComplex echophi = make_cuComplex(0, phi);
cuComplex echophiexp = cuCexpf(echophi);
cuComplex echo;
echo.x = echophiexp.x * amp;
echo.y = echophiexp.y * amp;
return echo;
}
__global__ void CUDA_DistanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* R, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
R[idx] = sqrtf(powf(Ax[idx] - Bx[idx], 2) + powf(Ay[idx] - By[idx], 2) + powf(Az[idx] - Bz[idx], 2));
}
}
__global__ void CUDA_B_DistanceA(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
R[idx] = sqrtf(powf(Ax[idx] - Bx, 2) + powf(Ay[idx] - By, 2) + powf(Az[idx] - Bz, 2));
}
}
__global__ void CUDA_make_VectorA_B(float sX, float sY, float sZ, float* tX, float* tY, float* tZ, float* RstX, float* RstY, float* RstZ, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
RstX[idx] = sX - tX[idx]; // 地面->天
RstY[idx] = sY - tY[idx];
RstZ[idx] = sZ - tZ[idx];
}
}
__global__ void CUDA_Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
R[idx] = sqrtf(powf(Vx[idx], 2) + powf(Vy[idx], 2) + powf(Vz[idx], 2));
}
}
__global__ void CUDA_cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float tAx = Ax[idx];
float tAy = Ay[idx];
float tAz = Az[idx];
float tBx = Bx[idx];
float tBy = By[idx];
float tBz = Bz[idx];
float AR = sqrtf(powf(tAx, 2) + powf(tAy, 2) + powf(tAz, 2));
float BR = sqrtf(powf(tBx, 2) + powf(tBy, 2) + powf(tBz, 2));
float dotAB = tAx * tBx + tAy * tBy + tAz * tBz;
float result = acosf(dotAB / (AR * BR));
anglecos[idx] = result;
}
}
__global__ void CUDA_SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* thetaAnt, float* phiAnt
, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float Xst = -1 * RstX[idx]; // 卫星 --> 地面
float Yst = -1 * RstY[idx];
float Zst = -1 * RstZ[idx];
float AntXaxisX = antXaxisX;
float AntXaxisY = antXaxisY;
float AntXaxisZ = antXaxisZ;
float AntYaxisX = antYaxisX;
float AntYaxisY = antYaxisY;
float AntYaxisZ = antYaxisZ;
float AntZaxisX = antZaxisX;
float AntZaxisY = antZaxisY;
float AntZaxisZ = antZaxisZ;
// 归一化
float RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
float AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
float AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
float AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
float Rx = Xst / RstNorm;
float Ry = Yst / RstNorm;
float Rz = Zst / RstNorm;
float Xx = AntXaxisX / AntXaxisNorm;
float Xy = AntXaxisY / AntXaxisNorm;
float Xz = AntXaxisZ / AntXaxisNorm;
float Yx = AntYaxisX / AntYaxisNorm;
float Yy = AntYaxisY / AntYaxisNorm;
float Yz = AntYaxisZ / AntYaxisNorm;
float Zx = AntZaxisX / AntZaxisNorm;
float Zy = AntZaxisY / AntZaxisNorm;
float Zz = AntZaxisZ / AntZaxisNorm;
float Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
// 计算theta 与 phi
float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
float ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
float PhiAnt = atanf(Yant / Xant); // -pi/2 ~pi/2
if (abs(Yant) < PRECISIONTOLERANCE) { // X轴上
PhiAnt = 0;
}
else if (abs(Xant) < PRECISIONTOLERANCE) { // Y轴上，原点
if (Yant > 0) {
PhiAnt = PI / 2;
}
else {
PhiAnt = -PI / 2;
}
}
else if (Xant < 0) {
if (Yant > 0) {
PhiAnt = PI + PhiAnt;
}
else {
PhiAnt = -PI+PhiAnt ;
}
}
else { // Xant>0 X 正轴
}
if (isnan(PhiAnt)) {
printf("V=[%f,%f,%f];norm=%f;thetaAnt=%f;phiAnt=%f;\n", Xant, Yant, Zant,Norm, ThetaAnt, PhiAnt);
}
//if (abs(ThetaAnt - 0) < PRECISIONTOLERANCE) {
// PhiAnt = 0;
//}
//else {}
thetaAnt[idx] = ThetaAnt*r2d;
phiAnt[idx] = PhiAnt*r2d;
//printf("Rst=[%f,%f,%f];AntXaxis = [%f, %f, %f];AntYaxis=[%f,%f,%f];AntZaxis=[%f,%f,%f];phiAnt=%f;thetaAnt=%f;\n", Xst, Yst, Zst
// , AntXaxisX, AntXaxisY, AntXaxisZ
// , AntYaxisX, AntYaxisY, AntYaxisZ
// , AntZaxisX, AntZaxisY, AntZaxisZ
// , phiAnt[idx]
// , thetaAnt[idx]
//);
}
}
__global__ void CUDA_BillerInterpAntPattern(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
long thetapoints, long phipoints,
float* searththeta, float* searchphi, float* searchantpattern,
long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float stheta = searththeta[idx];
float sphi = searchphi[idx];
float pthetaid = (stheta - starttheta) / dtheta;//
float pphiid = (sphi - startphi) / dphi;
long lasttheta = floorf(pthetaid);
long nextTheta = lasttheta + 1;
long lastphi = floorf(pphiid);
long nextPhi = lastphi + 1;
if (lasttheta < 0 || nextTheta < 0 || lastphi < 0 || nextPhi < 0 ||
lasttheta >= thetapoints || nextTheta >= thetapoints || lastphi >= phipoints || nextPhi >= phipoints)
{
searchantpattern[idx] = 0;
}
else {
float x = stheta;
float y = sphi;
float x1 = lasttheta * dtheta + starttheta;
float x2 = nextTheta * dtheta + starttheta;
float y1 = lastphi * dphi + startphi;
float y2 = nextPhi * dphi + startphi;
float z11 = antpattern[lasttheta * phipoints + lastphi];
float z12 = antpattern[lasttheta * phipoints + nextPhi];
float z21 = antpattern[nextTheta * phipoints + lastphi];
float z22 = antpattern[nextTheta * phipoints + nextPhi];
z11 = powf(10, z11 / 10);
z12 = powf(10, z12 / 10);
z21 = powf(10, z21 / 10);
z22 = powf(10, z22 / 10);
float GainValue = (z11 * (x2 - x) * (y2 - y)
+ z21 * (x - x1) * (y2 - y)
+ z12 * (x2 - x) * (y - y1)
+ z22 * (x - x1) * (y - y1));
GainValue = GainValue / ((x2 - x1) * (y2 - y1));
searchantpattern[idx] = GainValue;
}
}
}
__global__ void CUDA_Test_HelloWorld(float a, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
printf("\nidx:\t %d %d \n", idx, len);
}
__global__ void CUDA_RTPC(
float antPx, float antPy, float antPz,// 天线坐标
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* demx, float* demy, float* demz, long* demcls,
float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,
float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, long Tthetapoints, long Tphipoints,
float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, long Rthetapoints, long Rphipoints,
float lamda, float fs, float nearrange, float Pt, long Freqnumbers, // 参数
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby
cuComplex* outecho, int* d_echoAmpFID,
int linecount,int plusepoint) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
//printf("\nidx:\t %d %d %d\n", idx, linecount, plusepoint);
if (idx < linecount* plusepoint) {
long clsid = demcls[idx];
CUDAVector Rs{ antPx,antPy,antPz };
CUDAVector Rt{ demx[idx],demy[idx],demz[idx] };
CUDAVector Rst{ Rs.x - Rt.x,Rs.y - Rt.y,Rs.z - Rt.z };
CUDAVector Vslope{ demslopex[idx],demslopey[idx],demslopez[idx] };
float R = GPU_VectorNorm2(Rst); // 斜距
float slopeangle = demslopeangle[idx];
CUDAVectorEllipsoidal Rtanttheta = GPU_SatelliteAntDirectNormal( // 地面目标在天线的位置
Rst.x, Rst.y, Rst.z,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ);
float localangle = GPU_CosAngle_VectorA_VectorB(Rst, Vslope); // 距地入射角
float sigma = GPU_getSigma0dB(sigma0Paramslist[clsid], localangle * r2d);
sigma = powf(10.0, sigma / 10.0);// 后向散射系数
//printf("\ntheta: %f\t,%f ,%f ,%f ,%f ,%f ,%f \n", localangle * r2d, sigma0Paramslist[clsid].p1, sigma0Paramslist[clsid].p2, sigma0Paramslist[clsid].p3,
// sigma0Paramslist[clsid].p4, sigma0Paramslist[clsid].p5, sigma0Paramslist[clsid].p6);
// 发射方向图
float transPattern = GPU_BillerInterpAntPattern(Tantpattern,
Tstarttheta, Tstartphi, Tdtheta, Tdphi, Tthetapoints, Tphipoints,
Rtanttheta.theta, Rtanttheta.phi) * r2d;
// 接收方向图
float receivePattern = GPU_BillerInterpAntPattern(Rantpattern,
Rstarttheta, Rstartphi, Rdtheta, Rdphi, Rthetapoints, Rphipoints,
Rtanttheta.theta, Rtanttheta.phi) * r2d;
// 计算振幅、相位
float amp = Pt * transPattern * receivePattern * sigma * (1 / cos(slopeangle) * sin(localangle));
amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4));
float phi = (-4 * LAMP_CUDA_PI / lamda) * R;
// 构建回波
cuComplex echophi = make_cuComplex(0, phi);
cuComplex echophiexp = cuCexpf(echophi);
float timeR = 2 * (R - nearrange) / LIGHTSPEED * fs;
long timeID = floorf(timeR);
if (timeID < 0 || timeID >= Freqnumbers) {
timeID = 0;
amp = 0;
}
else {}
cuComplex echo;
echo.x = echophiexp.x * amp;
echo.y = echophiexp.y * amp;
outecho[idx] = echo;
d_echoAmpFID[idx] = timeID;
}
}
__global__ void CUDA_TBPImage(
float* antPx, float* antPy, float* antPz,
float* imgx, float* imgy, float* imgz,
cuComplex* echoArr, cuComplex* imgArr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount,
long prfid, long freqcount
) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
//printf("\nidx:\t %d %d %d\n", idx, linecount, plusepoint);
if (idx < rowcount * colcount) {
float R = sqrtf(powf(antPx[prfid] - imgx[idx], 2) + powf(antPy[prfid] - imgy[idx], 2) + powf(antPz[prfid] - imgz[idx], 2));
float Ridf = ((R - Rnear) * 2 / LIGHTSPEED) * fs;
long Rid = floorf(Ridf);
if(Rid <0|| Rid >= freqcount){}
else {
float factorj = freq * 4 * PI / LIGHTSPEED;
cuComplex Rphi =cuCexpf(make_cuComplex(0, factorj * R));// 校正项
imgArr[idx] = cuCaddf(imgArr[idx], cuCmulf(echoArr[Rid] , Rphi));// 矫正
}
}
}
__global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R, float* slopeangle,
float nearRange, float Fs, float Pt, float lamda, long FreqIDmax,
cuComplex* echoArr, long* FreqID,
long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float r = R[idx];
float amp = Pt * TransAnt[idx] * ReciveAnt[idx];
amp = amp * sigma0[idx];
amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // 反射强度
// 处理相位
float phi = (-4 * LAMP_CUDA_PI / lamda) * r;
cuComplex echophi = make_cuComplex(0, phi);
cuComplex echophiexp = cuCexpf(echophi);
float timeR = 2 * (r - nearRange) / LIGHTSPEED * Fs;
long timeID = floorf(timeR);
if (timeID < 0 || timeID >= FreqIDmax) {
timeID = 0;
amp = 0;
}
cuComplex echo;
echo.x = echophiexp.x * amp;
echo.y = echophiexp.y * amp;
echoArr[idx] = echo;
FreqID[idx] = timeID;
}
}
__global__ void CUDA_AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float temptheta = anttheta[idx];
float tempphi = antphi[idx];
float antPatternGain = GPU_BillerInterpAntPattern(antpattern,
starttheta, startphi, dtheta, dphi, thetapoints, phipoints,
temptheta, tempphi) ;
gain[idx] = antPatternGain;
}
}
//__global__ void Sigma0InterpPixel(long* demcls, float* demslopeangle, CUDASigmaParam* sigma0Paramslist, float* localangle, float* sigma0list, long sigmaparamslistlen, long len)
//{
// long idx = blockIdx.x * blockDim.x + threadIdx.x;
// if (idx < len) {
// long clsid = demcls[idx];
// if(clsid<=)
// sigma0list[idx] = 0;
// }
//}
__global__ void CUDA_InterpSigma(
long* demcls, float* sigmaAmp, float* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
long clsid = demcls[idx];
float localangle = localanglearr[idx] * r2d;
CUDASigmaParam tempsigma = sigma0Paramslist[clsid];
//printf("cls:%d;localangle=%f;\n",clsid, localangle);
if (localangle < 0 || localangle >= 90) {
sigmaAmp[idx] = 0;
}
else {}
if (abs(tempsigma.p1)< PRECISIONTOLERANCE&&
abs(tempsigma.p2) < PRECISIONTOLERANCE &&
abs(tempsigma.p3) < PRECISIONTOLERANCE &&
abs(tempsigma.p4) < PRECISIONTOLERANCE&&
abs(tempsigma.p5) < PRECISIONTOLERANCE&&
abs(tempsigma.p6) < PRECISIONTOLERANCE
) {
sigmaAmp[idx] = 0;
}
else {
float sigma = GPU_getSigma0dB(tempsigma, localangle);
sigma = powf(10.0, sigma / 10.0);// 后向散射系数
//printf("cls:%d;localangle=%f;sigma0=%f;\n", clsid, localangle, sigma);
sigmaAmp[idx] = sigma;
}
}
}
//错误提示
void checkCudaError(cudaError_t err, const char* msg) {
if (err != cudaSuccess) {
std::cerr << "CUDA error: " << msg << " (" << cudaGetErrorString(err) << ")" << std::endl;
exit(EXIT_FAILURE);
}
}
// 主机参数内存声明
extern "C" void* mallocCUDAHost(long memsize) {
void* ptr;
cudaMallocHost(&ptr, memsize);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("mallocCUDAHost CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
return ptr;
}
// 主机参数内存释放
extern "C" void FreeCUDAHost(void* ptr) {
cudaFreeHost(ptr);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("FreeCUDAHost CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
// GPU参数内存声明
extern "C" void* mallocCUDADevice(long memsize) {
void* ptr;
cudaMalloc(&ptr, memsize);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("mallocCUDADevice CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
return ptr;
}
// GPU参数内存释放
extern "C" void FreeCUDADevice(void* ptr) {
cudaFree(ptr);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("FreeCUDADevice CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
// GPU 内存数据转移
extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize) {
cudaMemcpy(deviceptr, hostptr, memsize, cudaMemcpyHostToDevice);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("HostToDevice CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize) {
cudaMemcpy(hostptr, deviceptr, memsize, cudaMemcpyDeviceToHost);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("DeviceToHost CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void CUDATestHelloWorld(float a,long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_Test_HelloWorld << <numBlocks, blockSize >> > (a, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("FreeCUDADevice CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
void CUDATBPImage(float* antPx, float* antPy, float* antPz,
float* imgx, float* imgy, float* imgz,
cuComplex* echoArr, cuComplex* imgArr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount,
long prfid, long freqcount)
{
int blockSize = 256; // 每个块的线程数
int numBlocks = (rowcount * colcount + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
//printf("\nCUDA_RTPC_SiglePRF blockSize:%d ,numBlock:%d\n",blockSize,numBlocks);
// 调用 CUDA 核函数 CUDA_RTPC_Kernel
CUDA_TBPImage << <numBlocks, blockSize >> > (
antPx, antPy, antPz,
imgx, imgy, imgz,
echoArr, imgArr,
freq, fs, Rnear, Rfar,
rowcount, colcount,
prfid, freqcount
);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDATBPImage CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void distanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* R, long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_DistanceAB << <numBlocks, blockSize >> > (Ax, Ay, Az, Bx, By, Bz, R, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void BdistanceAs(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_B_DistanceA << <numBlocks, blockSize >> > (Ax, Ay, Az, Bx, By, Bz, R, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void make_VectorA_B(float sX, float sY, float sZ, float* tX, float* tY, float* tZ, float* RstX, float* RstY, float* RstZ, long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_make_VectorA_B << <numBlocks, blockSize >> > (sX, sY, sZ, tX, tY, tZ, RstX, RstY, RstZ, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_Norm_Vector << <numBlocks, blockSize >> > (Vx, Vy, Vz, R, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos, long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_cosAngle_VA_AB << <numBlocks, blockSize >> > (Ax, Ay, Az, Bx, By, Bz, anglecos, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* thetaAnt, float* phiAnt
, long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_SatelliteAntDirectNormal << <numBlocks, blockSize >> > (RstX, RstY, RstZ,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ,
thetaAnt, phiAnt
, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
//printf("\nCUDA_RTPC_SiglePRF blockSize:%d ,numBlock:%d\n", blockSize, numBlocks);
CUDA_AntPatternInterpGain << <numBlocks, blockSize >> > ( anttheta,antphi, gain,
antpattern,
starttheta, startphi, dtheta, dphi, thetapoints, phipoints,
len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void CUDARTPCPRF(float antPx, long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
printf("\nCUDA_RTPC_SiglePRF blockSize:%d ,numBlock:%d\n", blockSize, numBlocks);
CUDA_Test_HelloWorld << <numBlocks, blockSize >> > (antPx, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R, float* slopeangle,
float nearRange, float Fs, float pt, float lamda, long FreqIDmax,
cuComplex* echoAmp, long* FreqID,
long len)
{
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_calculationEcho << <numBlocks, blockSize >> > (sigma0, TransAnt, ReciveAnt,
localangle, R, slopeangle,
nearRange, Fs, pt, lamda, FreqIDmax,
echoAmp, FreqID,
len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void CUDA_RTPC_SiglePRF(
float antPx, float antPy, float antPZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* demx, float* demy, float* demz, long* demcls,
float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,
float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, int Tthetapoints, int Tphipoints,
float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, int Rthetapoints, int Rphipoints,
float lamda, float fs, float nearrange, float Pt, int Freqnumbers,
CUDASigmaParam* sigma0Paramslist, int sigmaparamslistlen,
cuComplex* outecho, int* d_echoAmpFID,
int linecount,int colcount) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (linecount* colcount + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
//printf("\nCUDA_RTPC_SiglePRF blockSize:%d ,numBlock:%d\n",blockSize,numBlocks);
// 调用 CUDA 核函数 CUDA_RTPC_Kernel
CUDA_RTPC << <numBlocks, blockSize >> > (
antPx, antPy, antPZ,// 天线坐标
antXaxisX, antXaxisY, antXaxisZ, // 天线坐标系
antYaxisX, antYaxisY, antYaxisZ, //
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ,// 天线指向
demx, demy, demz,
demcls, // 地面坐标
demslopex, demslopey, demslopez, demslopeangle,// 地面坡度
Tantpattern, Tstarttheta, Tstartphi, Tdtheta, Tdphi, Tthetapoints, Tphipoints,// 天线方向图相关
Rantpattern, Rstarttheta, Rstartphi, Rdtheta, Rdphi, Rthetapoints, Rphipoints,// 天线方向图相关
lamda, fs, nearrange, Pt, Freqnumbers, // 参数
sigma0Paramslist, sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby
outecho, d_echoAmpFID,
linecount, colcount
);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void CUDAInterpSigma(
long* demcls,float* sigmaAmp, float* localanglearr,long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) {// 地表覆盖类型-sigma插值对应函数-ulaby
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_InterpSigma << <numBlocks, blockSize >> > (
demcls, sigmaAmp, localanglearr, len,
sigma0Paramslist, sigmaparamslistlen
);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RTPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
#endif

452
BaseTool/GPUTool.cuh Normal file
View File

@ -0,0 +1,452 @@
#ifndef GPUTOOL_H
#define GPUTOOL_H
#ifdef __CUDANVCC___
#include "BaseConstVariable.h"
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#define __CUDADEBUG__
// 默认显存分布
enum LAMPGPUDATETYPE {
LAMP_LONG,
LAMP_FLOAT,
LAMP_COMPLEXFLOAT
};
extern "C" struct CUDASigmaParam {
float p1;
float p2;
float p3;
float p4;
float p5;
float p6;
};
extern "C" struct CUDAVector {
float x;
float y;
float z;
};
extern "C" struct CUDAVectorEllipsoidal {
float theta;
float phi;
float pho;
};
// GPU 内存函数
extern "C" void* mallocCUDAHost( long memsize); // 主机内存声明
extern "C" void FreeCUDAHost(void* ptr);
extern "C" void* mallocCUDADevice( long memsize); // GPU内存声明
extern "C" void FreeCUDADevice(void* ptr);
extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移 设备 -> GPU
extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移 GPU -> 设备
// 仿真所需的常用函数
extern "C" void distanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* R, long member);
extern "C" void BdistanceAs(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long member);
extern "C" void make_VectorA_B(float sX, float sY, float sZ, float* tX, float* tY, float* tZ, float* RstX, float* RstY, float* RstZ, long member);
extern "C" void Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long member);
extern "C" void cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos, long len);
extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ, float antXaxisX, float antXaxisY, float antXaxisZ, float antYaxisX, float antYaxisY, float antYaxisZ, float antZaxisX, float antZaxisY, float antZaxisZ, float antDirectX, float antDirectY, float antDirectZ, float* thetaAnt, float* phiAnt, long len);
extern "C" void AntPatternInterpGain(float* anttheta, float* antphi, float* gain, float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints,long len);
extern "C" void CUDA_RTPC_SiglePRF(
float antPx, float antPy, float antPZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* demx, float* demy, float* demz, long* demcls,
float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,
float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, int Tthetapoints, int Tphipoints,
float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, int Rthetapoints, int Rphipoints,
float lamda, float fs, float nearrange, float Pt, int Freqnumbers,
CUDASigmaParam* sigma0Paramslist, int sigmaparamslistlen,
cuComplex* outecho, int* d_echoAmpFID,
int linecount, int colcount
);
extern "C" void CUDARTPCPRF(float antPx, long len);
extern "C" void CUDATestHelloWorld(float a, long len);
extern "C" void CUDATBPImage(
float* antPx,
float* antPy,
float* antPz,
float* imgx,
float* imgy,
float* imgz,
cuComplex* echoArr,
cuComplex* imgArr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount,
long prfid, long freqcount
);
extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R, float* slopeangle,
float nearRange, float Fs, float pt, float lamda, long FreqIDmax,
cuComplex* echoAmp, long* FreqID,
long len);
extern "C" void CUDAInterpSigma(
long* demcls, float* sigmaAmp, float* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen);
#endif
#endif
/**
*
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<double>(databuffer[i * nXSize * 2 + j * 2],
databuffer[i * nXSize * 2 + j * 2 + 1]);
}
}
delete[] databuffer;
gdalImage demxyz(demxyzPath);// 地面点坐标
gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型
gdalImage demsloperxyz(this->demsloperPath);// 地面坡向
omp_lock_t lock; // 定义锁
omp_init_lock(&lock); // 初始化锁
long start_ids = 1250;
for (start_ids = 1; start_ids < demxyz.height; start_ids = start_ids + line_invert) { // 8+ 17 + 0.3 MB
QDateTime current = QDateTime::currentDateTime();
long pluseStep = Memory1MB * 100 / 3 / PlusePoint;
if (pluseStep * num_thread * 3 > this->PluseCount) {
pluseStep = this->PluseCount / num_thread / 3;
}
pluseStep = pluseStep > 50 ? pluseStep : 50;
qDebug() << current.toString("yyyy-MM-dd HH:mm:ss.zzz") << " \tstart \t " << start_ids << " - " << start_ids + line_invert << "\t" << demxyz.height << "\t pluseCount:\t" << pluseStep;
// 文件读取
Eigen::MatrixXd dem_x = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 1); //
Eigen::MatrixXd dem_y = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 2); //
Eigen::MatrixXd dem_z = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 3); //
// 地表覆盖
std::shared_ptr<long> dem_landcls = readDataArr<long>(demlandcls, start_ids - 1, 0, line_invert + 1, demxyz.width, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD); // 地表覆盖类型
long* dem_landcls_ptr = dem_landcls.get();
double localAngle = 30.0;
bool sigmaNoZeroFlag = true;
for (long ii = 0; ii < dem_x.rows(); ii++) {
for (long jj = 0; jj < dem_y.cols(); jj++) {
if (0 != this->SigmaDatabasePtr->getAmp(dem_landcls_ptr[dem_x.cols() * ii + jj], localAngle, polartype)) {
sigmaNoZeroFlag = false;
break;
}
}
if (!sigmaNoZeroFlag) {
break;
}
}
if (sigmaNoZeroFlag) {
continue;
}
//#ifdef DEBUGSHOWDIALOG
// dialog->load_double_MatrixX_data(dem_z, "dem_z");
//#endif
Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 1); //
Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 2); //
Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 3); //
Eigen::MatrixXd sloperAngle = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 4); //
sloperAngle = sloperAngle.array() * T180_PI;
long dem_rows = dem_x.rows();
long dem_cols = dem_x.cols();
long freqidx = 0;//
#ifdef DEBUGSHOWDIALOG
ImageShowDialogClass* dialog = new ImageShowDialogClass(nullptr);
dialog->show();
Eigen::MatrixXd landaArr = Eigen::MatrixXd::Zero(dem_rows, dem_cols);
for (long i = 0; i < dem_rows; i++) {
for (long j = 0; j < dem_cols; j++) {
landaArr(i, j) = dem_landcls.get()[i * dem_cols + j];
}
}
dialog->load_double_MatrixX_data(landaArr, "landCover");
#endif
//qDebug() << " pluse bolck size :\t " << pluseStep << " all size:\t" << this->PluseCount;
long processNumber = 0;
#pragma omp parallel for
for (long startprfidx = 0; startprfidx < this->PluseCount; startprfidx = startprfidx + pluseStep) { // 17 + 0.3 MB
long prfcount_step = startprfidx + pluseStep < this->PluseCount ? pluseStep : this->PluseCount - startprfidx;
Eigen::MatrixXcd echoPluse = Eigen::MatrixXcd::Zero(prfcount_step, PlusePoint); // 当前脉冲的回波积分情况
// 内存预分配
Eigen::MatrixXd Rst_x = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Rst_y = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Rst_z = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd R = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd localangle = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Vst_x = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Vst_y = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Vst_z = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd fde = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd fr = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Rx = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd sigam = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd echoAmp = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols()).array() + Pt;
Eigen::MatrixXd Rphi = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd TimeRange = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd TransAnt = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd ReciveAnt = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd AntTheta = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd AntPhi = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
double minR = 0, maxR = 0;
double minLocalAngle = 0, maxLocalAngle = 0;
Vector3D Rt = { 0,0,0 };
SatelliteOribtNode oRs = SatelliteOribtNode{ 0 };;
Vector3D p0 = {}, slopeVector = {}, sateAntDirect = {};
Vector3D Rs = {}, Vs = {}, Ast = {};
SatelliteAntDirect antdirectNode = {};
std::complex<double> echofreq;
std::complex<double> Imag1(0, 1);
double TAntPattern = 1; // 发射天线方向图
double RAntPanttern = 1;// 接收天线方向图
double maxechoAmp = 1;
double tempAmp = 1;
for (long prfidx = 0; prfidx < prfcount_step; prfidx++)
{
oRs = sateOirbtNodes[prfidx + startprfidx];
// 计算天线方向图
for (long jj = 1; jj < dem_cols - 1; jj++) {
for (long ii = 1; ii < dem_rows - 1; ii++) {
p0.x = dem_x(ii, jj);
p0.y = dem_y(ii, jj);
p0.z = dem_z(ii, jj);
this->TaskSetting->getSatelliteAntDirectNormal(oRs, p0, antdirectNode);
//antdirectNode.ThetaAnt = antdirectNode.ThetaAnt * r2d;
//antdirectNode.PhiAnt = antdirectNode.PhiAnt * r2d;
AntTheta(ii, jj) = antdirectNode.ThetaAnt * r2d;
AntPhi(ii, jj) = antdirectNode.PhiAnt * r2d;
}
}
// 计算发射天线方向图
for (long jj = 1; jj < dem_cols - 1; jj++) {
for (long ii = 1; ii < dem_rows - 1; ii++) {
TransformPattern->getGainLinear(AntTheta(ii, jj), AntPhi(ii, jj), TransAnt(ii, jj));
//TransAnt(ii, jj) = TAntPattern;
}
}
// 计算接收天线方向图
for (long jj = 1; jj < dem_cols - 1; jj++) {
for (long ii = 1; ii < dem_rows - 1; ii++) {
TransformPattern->getGainLinear(AntTheta(ii, jj), AntPhi(ii, jj), ReciveAnt(ii, jj));
//ReciveAnt(ii, jj) = RAntPanttern;
}
}
// 计算经过增益的能量
echoAmp = Pt * TransAnt.array() * ReciveAnt.array();
maxechoAmp = echoAmp.maxCoeff();
if (std::abs(maxechoAmp) < PRECISIONTOLERANCE) { // 这种情况下，不在合成孔径范围中
continue;
}
Rs.x = sateOirbtNodes[prfidx + startprfidx].Px; // 卫星位置
Rs.y = sateOirbtNodes[prfidx + startprfidx].Py;
Rs.z = sateOirbtNodes[prfidx + startprfidx].Pz;
Vs.x = sateOirbtNodes[prfidx + startprfidx].Vx; // 卫星速度
Vs.y = sateOirbtNodes[prfidx + startprfidx].Vy;
Vs.z = sateOirbtNodes[prfidx + startprfidx].Vz;
Ast.x = sateOirbtNodes[prfidx + startprfidx].AVx;// 卫星加速度
Ast.y = sateOirbtNodes[prfidx + startprfidx].AVy;
Ast.z = sateOirbtNodes[prfidx + startprfidx].AVz;
Rst_x = Rs.x - dem_x.array(); // Rst = Rs - Rt;
Rst_y = Rs.y - dem_y.array();
Rst_z = Rs.z - dem_z.array();
R = (Rst_x.array().pow(2) + Rst_y.array().pow(2) + Rst_z.array().pow(2)).array().sqrt(); // R
minR = R.minCoeff();
maxR = R.maxCoeff();
//qDebug() << "minR:\t" << minR << " maxR:\t" << maxR;
if (maxR<NearRange || minR>FarRange) {
continue;
}
else {}
// getCosAngle
// double c = dot(a, b) / (getlength(a) * getlength(b));
// return acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d;
// localangle = getCosAngle(Rst, slopeVector) * T180_PI; // 注意这个只能实时计算，因为非实时计算代价太大
localangle = (Rst_x.array() * demsloper_x.array() + Rst_y.array() * demsloper_y.array() + Rst_z.array() * demsloper_z.array()).array(); // dot(a, b)
localangle = localangle.array() / R.array();
localangle = localangle.array() / (demsloper_x.array().pow(2) + demsloper_y.array().pow(2) + demsloper_z.array().pow(2)).array().sqrt().array();
localangle = localangle.array().acos(); // 弧度值
minLocalAngle = localangle.minCoeff();
maxLocalAngle = localangle.maxCoeff();
if (maxLocalAngle<0 || minLocalAngle>PI / 2) {
continue;
}
else {}
//Vst_x = Vs.x + 1 * earthRoute * dem_y.array(); // Vst = Vs - Vt;
//Vst_y = Vs.y - 1 * earthRoute * dem_x.array();
//Vst_z = Vs.z - Eigen::MatrixXd::Zero(dem_x.rows(), dem_y.cols()).array();
//// 计算多普勒中心频率 Rst, Vst : ( - 2 / lamda) * dot(Rs - Rt, Vs - Vt) / R; // 星载合成孔径雷达原始回波数据模拟研究 3.18
//fde = (-2 / lamda) * (Rst_x.array() * Vst_x.array() + Rst_y.array() * Vst_y.array() + Rst_z.array() * Vst_z.array()).array() / (R.array());
//// 计算多普勒频率斜率 // 星载合成孔径雷达原始回波数据模拟研究 3.19
//// -(2/lamda)*( dot(Vs - Vt, Vs - Vt)/R + dot(Ast, Rs - Rt)/R - std::pow(dot(Vs - Vt, Rs - Rt),2 )/std::pow(R,3));
//fr = (-2 / lamda) *
// (Vst_x.array() * Vst_x.array() + Vst_y.array() * Vst_y.array() + Vst_z.array() * Vst_z.array()).array() / (R.array()) +
// (-2 / lamda) *
// (Ast.x * Rst_x.array() + Ast.y * Rst_y.array() + Ast.z * Rst_z.array()).array() / (R.array()) -
// (-2 / lamda) *
// (Vst_x.array() * Rst_x.array() + Vst_y.array() * Rst_y.array() + Vst_z.array() * Rst_z.array()).array().pow(2) / (R.array().pow(3));
// 计算回波
Rx = R;// -(lamda / 2) * (fde * TRx + 0.5 * fr * TRx * TRx); // 斜距历程值
// 逐点计算 this->SigmaDatabasePtr->getAmp(covercls, localangle, polartype); // 后向散射系数 HH
for (long ii = 0; ii < dem_x.rows(); ii++) {
for (long jj = 0; jj < dem_y.cols(); jj++) {
sigam(ii, jj) = this->SigmaDatabasePtr->getAmp(dem_landcls_ptr[dem_x.cols() * ii + jj], localangle(ii, jj) * r2d, polartype);
}
}
if (sigam.maxCoeff() > 0) {}
else {
continue;
}
// projArea = 1 / std::cos(sloperAngle) * std::sin(localangle); // 投影面积系数，单位投影面积 1m x 1m --注意这里是假设，后期再补充
// echoAmp = projArea*TAntPattern * RAntPanttern * sigam / (4 * PI * R * R);
echoAmp = echoAmp.array() * sigam.array() * (1 / sloperAngle.array().cos() * localangle.array().sin()); // 反射强度
echoAmp = echoAmp.array() / (4 * PI * R.array().pow(2)); // 距离衰减
Rphi = -4 * PI / lamda * Rx.array();// 距离徙动相位
// 积分
TimeRange = ((2 * R.array() / LIGHTSPEED - TimgNearRange).array() * Fs).array();
double localAnglepoint = -1;
long prf_freq_id = 0;
for (long jj = 1; jj < dem_cols - 1; jj++) {
for (long ii = 1; ii < dem_rows - 1; ii++) {
prf_freq_id = std::floor(TimeRange(ii, jj));
if (prf_freq_id < 0 || prf_freq_id >= PlusePoint || localangle(ii, jj) < 0 || localangle(ii, jj) > PI / 2 || echoAmp(ii, jj) == 0) {
continue;
}
echofreq = echoAmp(ii, jj) * std::exp(Rphi(ii, jj) * Imag1);
echoPluse(prfidx, prf_freq_id) = echoPluse(prfidx, prf_freq_id) + echofreq;
}
}
#ifdef DEBUGSHOWDIALOG
ImageShowDialogClass* localangledialog = new ImageShowDialogClass(dialog);
localangledialog->show();
localangledialog->load_double_MatrixX_data(localangle.array() * r2d, "localangle");
ImageShowDialogClass* sigamdialog = new ImageShowDialogClass(dialog);
sigamdialog->show();
sigamdialog->load_double_MatrixX_data(TimeRange, "TimeRange");
ImageShowDialogClass* ampdialog = new ImageShowDialogClass(dialog);
ampdialog->show();
ampdialog->load_double_MatrixX_data(echoAmp, "echoAmp");
Eigen::MatrixXd echoPluseamp = echoPluse.array().abs().cast<double>().array();
ImageShowDialogClass* echoampdialog = new ImageShowDialogClass(dialog);
echoampdialog->show();
echoampdialog->load_double_MatrixX_data(echoPluseamp, "echoPluseamp");
dialog->exec();
#endif
//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << " end " << prfidx;
}
//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz")<<" step "<< prfcount_step;
omp_set_lock(&lock); // 回波整体赋值处理
for (long prfidx = 0; prfidx < prfcount_step; prfidx++) {
for (long freqidx = 0; freqidx < PlusePoint; freqidx++)
{
//qDebug() << prfidx << " " << freqidx << " " << echoPluse(prfidx, freqidx).real() << " + " << echoPluse(prfidx, freqidx).imag() << " j";
echo.get()[(prfidx + startprfidx) * PlusePoint + freqidx] = echo.get()[(prfidx + startprfidx) * PlusePoint + freqidx] + echoPluse(prfidx, freqidx);
}
}
//this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
omp_unset_lock(&lock); // 解锁
//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << " step 2" << prfcount_step;
}
omp_set_lock(&lock); // 保存文件
processNumber = processNumber + pluseStep;
this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
omp_unset_lock(&lock); // 解锁
qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << " \t " << start_ids << "\t--\t " << start_ids + line_invert << "\t/\t" << demxyz.height;
}
omp_destroy_lock(&lock); // 销毁锁
*/

410
BaseTool/GeoOperator.cpp Normal file
View File

@ -0,0 +1,410 @@
#include "stdafx.h"
#include "GeoOperator.h"
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
////#include <mkl.h>
#include <string>
#include <omp.h>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
//#include <ogr_geos.h>
#include <ogrsf_frmts.h> //#include "ogrsf_frmts.h"
#include <fstream>
#include <proj.h>
#include "GeoOperator.h"
Landpoint operator +(const Landpoint& p1, const Landpoint& p2)
{
return Landpoint{ p1.lon + p2.lon,p1.lat + p2.lat,p1.ati + p2.ati };
}
Landpoint operator -(const Landpoint& p1, const Landpoint& p2)
{
return Landpoint{ p1.lon - p2.lon,p1.lat - p2.lat,p1.ati - p2.ati };
}
bool operator ==(const Landpoint& p1, const Landpoint& p2)
{
return p1.lat == p2.lat && p1.lon == p2.lon && p1.ati == p2.ati;
}
Landpoint operator *(const Landpoint& p, double scale)
{
return Landpoint{
p.lon * scale,
p.lat * scale,
p.ati * scale
};
}
Landpoint LLA2XYZ(const Landpoint& LLA) {
double L = LLA.lon * d2r;
double B = LLA.lat * d2r;
double H = LLA.ati;
double sinB = sin(B);
double cosB = cos(B);
//double N = a / sqrt(1 - e * e * sin(B) * sin(B));
double N = a / sqrt(1 - eSquare * sinB * sinB);
Landpoint result = { 0,0,0 };
result.lon = (N + H) * cosB * cos(L);
result.lat = (N + H) * cosB * sin(L);
//result.z = (N * (1 - e * e) + H) * sin(B);
result.ati = (N * (1 - 1 / f_inverse) * (1 - 1 / f_inverse) + H) * sinB;
return result;
}
void LLA2XYZ(const Landpoint& LLA, Point3& XYZ)
{
double L = LLA.lon * d2r;
double B = LLA.lat * d2r;
double H = LLA.ati;
double sinB = sin(B);
double cosB = cos(B);
//double N = a / sqrt(1 - e * e * sin(B) * sin(B));
double N = a / sqrt(1 - eSquare * sinB * sinB);
Landpoint result = { 0,0,0 };
XYZ.x = (N + H) * cosB * cos(L);
XYZ.y = (N + H) * cosB * sin(L);
//result.z = (N * (1 - e * e) + H) * sin(B);
XYZ.z = (N * (1 - 1 / f_inverse) * (1 - 1 / f_inverse) + H) * sinB;
}
Eigen::MatrixXd LLA2XYZ(Eigen::MatrixXd landpoint)
{
landpoint.col(0) = landpoint.col(0).array() * d2r; // lon
landpoint.col(1) = landpoint.col(1).array() * d2r; // lat
Eigen::MatrixXd sinB = (landpoint.col(1).array().sin());//lat
Eigen::MatrixXd cosB = (landpoint.col(1).array().cos());//lat
Eigen::MatrixXd N = a / ((1 - sinB.array().pow(2) * eSquare).array().sqrt());
Eigen::MatrixXd result(landpoint.rows(), 3);
result.col(0) = (N.array() + landpoint.col(2).array()) * cosB.array() * Eigen::cos(landpoint.col(0).array()).array(); //x
result.col(1) = (N.array() + landpoint.col(2).array()) * cosB.array() * Eigen::sin(landpoint.col(0).array()).array(); //y
result.col(2) = (N.array() * (1 - 1 / f_inverse) * (1 - 1 / f_inverse) + landpoint.col(2).array()) * sinB.array(); //z
return result;
}
Landpoint XYZ2LLA(const Landpoint& XYZ) {
double tmpX = XYZ.lon;//
double temY = XYZ.lat;//
double temZ = XYZ.ati;
double curB = 0;
double N = 0;
double sqrtTempXY = sqrt(tmpX * tmpX + temY * temY);
double calB = atan2(temZ, sqrtTempXY);
int counter = 0;
double sinCurB = 0;
while (abs(curB - calB) * r2d > epsilon && counter < 25)
{
curB = calB;
sinCurB = sin(curB);
N = a / sqrt(1 - eSquare * sinCurB * sinCurB);
calB = atan2(temZ + N * eSquare * sinCurB, sqrtTempXY);
counter++;
}
Landpoint result = { 0,0,0 };
result.lon = atan2(temY, tmpX) * r2d;
result.lat = curB * r2d;
result.ati = temZ / sinCurB - N * (1 - eSquare);
return result;
}
double getAngle(const Landpoint& a, const Landpoint& b)
{
double c = dot(a, b) / (getlength(a) * getlength(b));
if (a.lon * b.lat - a.lat * b.lon >= 0) {
return acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d;
}
else {
return 360 - acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d;
}
}
double dot(const Landpoint& p1, const Landpoint& p2)
{
return p1.lat * p2.lat + p1.lon * p2.lon + p1.ati * p2.ati;
}
double getlength(const Landpoint& p1) {
return sqrt(dot(p1, p1));
}
Landpoint crossProduct(const Landpoint& a, const Landpoint& b) {
return Landpoint{
a.lat * b.ati - a.ati * b.lat,//x
a.ati * b.lon - a.lon * b.ati,//y
a.lon * b.lat - a.lat * b.lon//z
};
}
float cross2d(Point3 a, Point3 b)
{
return a.x * b.y - a.y * b.x;
}
Point3 operator -(Point3 a, Point3 b)
{
return Point3{ a.x - b.x, a.y - b.y, a.z - b.z };
}
Point3 operator +(Point3 a, Point3 b)
{
return Point3{ a.x + b.x, a.y + b.y, a.z + b.z };
}
double operator /(Point3 a, Point3 b)
{
return sqrt(pow(a.x, 2) + pow(a.y, 2)) / sqrt(pow(b.x, 2) + pow(b.y, 2));
}
Landpoint getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landpoint& p2, const Landpoint& p3, const Landpoint& p4,bool inLBH) {
Landpoint n0, n1, n2, n3, n4;
if (inLBH) {
n0 = LLA2XYZ(p0);
n1 = LLA2XYZ(p1);
n2 = LLA2XYZ(p2);
n3 = LLA2XYZ(p3);
n4 = LLA2XYZ(p4);
}
else {
n0 = p0;
n1 = p1;
n2 = p2;
n3 = p3;
n4 = p4;
}
Landpoint n01 = n1 - n0, n02 = n2 - n0, n03 = n3 - n0, n04 = n4 - n0;
// ��n01��������
Landpoint np01 = p1 - p0, np02 = p2 - p0, np03 = p3 - p0, np04 = p4 - p0;
double a2 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np02.lon,np02.lat,0 });// 01->02
double a3 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np03.lon,np03.lat,0 });// 01->03
double a4 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np04.lon,np04.lat,0 });// 01->04
//qDebug() << a2 << "\t" << a3 << "\t" << a4 << endl;
a2 = 360 - a2;
a3 = 360 - a3;
a4 = 360 - a4;
Landpoint N, W, S, E;
N = n01;
if (a2 >= a3 && a2 >= a4) {
W = n02;
if (a3 >= a4) {
S = n03;
E = n04;
}
else {
S = n04;
E = n03;
}
}
else if (a3 >= a2 && a3 >= a4) {
W = n03;
if (a2 >= a4) {
S = n02;
E = n04;
}
else {
S = n04;
E = n02;
}
}
else if (a4 >= a2 && a4 >= a3)
{
W = n04;
if (a2 >= a3) {
S = n02;
E = n03;
}
else {
S = n03;
E = n02;
}
}
return (crossProduct(N, W) + crossProduct(W, S) + crossProduct(S, E) + crossProduct(E, N)) * 0.25;
}
double distance(const Vector3D& p1, const Vector3D& p2)
{
double dx = p1.x - p2.x;
double dy = p1.y - p2.y;
double dz = p1.z - p2.z;
return std::sqrt(dx * dx + dy * dy + dz * dz);
}
double pointToLineDistance(const Vector3D& point, const Vector3D& linePoint, const Vector3D& lineDirection)
{
Vector3D pointToLine = { point.x - linePoint.x, point.y - linePoint.y, point.z - linePoint.z };
// 计算点到直线的投影点的ä½<C3A4>ç½?
double t = (pointToLine.x * lineDirection.x + pointToLine.y * lineDirection.y + pointToLine.z * lineDirection.z) /
(lineDirection.x * lineDirection.x + lineDirection.y * lineDirection.y + lineDirection.z * lineDirection.z);
// 计算投影�
Vector3D projection = { linePoint.x + t * lineDirection.x, linePoint.y + t * lineDirection.y, linePoint.z + t * lineDirection.z };
// 计算点到直线的è·<C3A8>ç¦?
return distance(point, projection);
}
Vector3D operator+(const Vector3D& p1, const Vector3D& p2)
{
return Vector3D{ p1.x + p2.x,p1.y + p2.y,p1.z + p2.z };
}
Vector3D operator-(const Vector3D& p1, const Vector3D& p2)
{
return Vector3D{ p1.x - p2.x,p1.y - p2.y,p1.z - p2.z };
}
bool operator==(const Vector3D& p1, const Vector3D& p2)
{
return p1.x == p2.x && p1.y == p2.y && p1.z == p2.z;
}
Vector3D operator*(const Vector3D& p, double scale)
{
return Vector3D{
p.x * scale,
p.y * scale,
p.z * scale
};
}
Vector3D operator*(double scale, const Vector3D& p)
{
return Vector3D{
p.x * scale,
p.y * scale,
p.z * scale
};
}
double getAngle(const Vector3D& a, const Vector3D& b)
{
double c = dot(a, b) / (getlength(a) * getlength(b));
if (a.x * b.y - a.y * b.x >= 0) {
return acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d;
}
else {
return 360 - acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d;
}
}
double getCosAngle(const Vector3D& a, const Vector3D& b)
{
double c = dot(a, b) / (getlength(a) * getlength(b));
return acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d;
}
double dot(const Vector3D& p1, const Vector3D& p2)
{
return p1.y * p2.y + p1.x * p2.x + p1.z * p2.z;
}
double getlength(const Vector3D& p1)
{
return std::sqrt(std::pow(p1.x, 2) + std::pow(p1.y, 2) + std::pow(p1.z, 2));
}
Vector3D crossProduct(const Vector3D& a, const Vector3D& b)
{
return Vector3D{
a.y * b.z - a.z * b.y,//x
a.z * b.x - a.x * b.z,//y
a.x * b.y - a.y * b.x//z
};
}
/// <summary>
/// n1
/// n4 n0 n2
/// n3
/// </summary>
Vector3D getSlopeVector(const Vector3D& n0, const Vector3D& n1, const Vector3D& n2, const Vector3D& n3, const Vector3D& n4)
{
Vector3D n01 = n1 - n0, n02 = n2 - n0, n03 = n3 - n0, n04 = n4 - n0;
return (crossProduct(n01, n04) + crossProduct(n04, n03) + crossProduct(n03, n02) + crossProduct(n02, n01)) * 0.25;
}
SphericalCoordinates cartesianToSpherical(const CartesianCoordinates& cartesian)
{
SphericalCoordinates spherical;
spherical.r = std::sqrt(cartesian.x * cartesian.x + cartesian.y * cartesian.y + cartesian.z * cartesian.z);
spherical.theta = std::acos(cartesian.z / spherical.r);
spherical.phi = std::atan2(cartesian.y, cartesian.x);
return spherical;
}
CartesianCoordinates sphericalToCartesian(const SphericalCoordinates& spherical)
{
CartesianCoordinates cartesian;
cartesian.x = spherical.r * std::sin(spherical.theta) * std::cos(spherical.phi);
cartesian.y = spherical.r * std::sin(spherical.theta) * std::sin(spherical.phi);
cartesian.z = spherical.r * std::cos(spherical.theta);
return cartesian;
}
double getlocalIncAngle(Landpoint& satepoint, Landpoint& landpoint, Landpoint& slopeVector) {
Landpoint Rsc = satepoint - landpoint; // AB=B-A
//double R = getlength(Rsc);
//std::cout << R << endl;
double angle = getAngle(Rsc, slopeVector);
if (angle >= 180) {
return 360 - angle;
}
else {
return angle;
}
}
double getlocalIncAngle(Vector3D& satepoint, Vector3D& landpoint, Vector3D& slopeVector){
Vector3D Rsc = satepoint - landpoint; // AB=B-A
//double R = getlength(Rsc);
//std::cout << R << endl;
double angle = getAngle(Rsc, slopeVector);
if (angle >= 180) {
return 360 - angle;
}
else {
return angle;
}
}

126
BaseTool/GeoOperator.h Normal file
View File

@ -0,0 +1,126 @@
#pragma once
#ifndef _GEOOPERATOR_H
#define _GEOOPERATOR_H
#include "BaseConstVariable.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <string.h>
#include <iostream>
/// <summary>
/// 将经纬度转换为地固参心坐标系
/// </summary>
/// <param name="XYZP">经纬度点--degree</param>
/// <returns>投影坐标系点</returns>
Landpoint LLA2XYZ(const Landpoint& LLA);
void LLA2XYZ(const Landpoint& LLA,Point3& XYZ);
Eigen::MatrixXd LLA2XYZ(Eigen::MatrixXd landpoint);
/// <summary>
/// 将地固参心坐标系转换为经纬度
/// </summary>
/// <param name="XYZ">固参心坐标系</param>
/// <returns>经纬度--degree</returns>
Landpoint XYZ2LLA(const Landpoint& XYZ);
Landpoint operator +(const Landpoint& p1, const Landpoint& p2);
Landpoint operator -(const Landpoint& p1, const Landpoint& p2);
bool operator ==(const Landpoint& p1, const Landpoint& p2);
Landpoint operator *(const Landpoint& p, double scale);
double getAngle(const Landpoint& a, const Landpoint& b);
double dot(const Landpoint& p1, const Landpoint& p2);
double getlength(const Landpoint& p1);
Landpoint crossProduct(const Landpoint& a, const Landpoint& b);
Landpoint getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landpoint& p2, const Landpoint& p3, const Landpoint& p4, bool inLBH=true);
double getlocalIncAngle(Landpoint& satepoint, Landpoint& landpoint, Landpoint& slopeVector);
float cross2d(Point3 a, Point3 b);
Point3 operator -(Point3 a, Point3 b);
Point3 operator +(Point3 a, Point3 b);
double operator /(Point3 a, Point3 b);
// 矢量计算
struct Vector3D {
double x, y, z;
};
// 计算两点之间的距离
double distance(const Vector3D& p1, const Vector3D& p2);
// 计算点到直线的最短距离
double pointToLineDistance(const Vector3D& point, const Vector3D& linePoint, const Vector3D& lineDirection);
Vector3D operator +(const Vector3D& p1, const Vector3D& p2);
Vector3D operator -(const Vector3D& p1, const Vector3D& p2);
bool operator ==(const Vector3D& p1, const Vector3D& p2);
Vector3D operator *(const Vector3D& p, double scale);
Vector3D operator *(double scale,const Vector3D& p );
double getAngle(const Vector3D& a, const Vector3D& b);
double getCosAngle(const Vector3D& a, const Vector3D& b);
double dot(const Vector3D& p1, const Vector3D& p2);
double getlength(const Vector3D& p1);
Vector3D crossProduct(const Vector3D& a, const Vector3D& b);
/// <summary>
/// n1
/// n4 n0 n2
/// n3
/// </summary>
/// <param name="n0"></param>
/// <param name="n1"></param>
/// <param name="n2"></param>
/// <param name="n3"></param>
/// <param name="n4"></param>
/// <returns></returns>
Vector3D getSlopeVector(const Vector3D& n0, const Vector3D& n1, const Vector3D& n2, const Vector3D& n3, const Vector3D& n4);
struct CartesianCoordinates {
double x, y, z;
};
struct SphericalCoordinates {
double r, theta, phi;
};
SphericalCoordinates cartesianToSpherical(const CartesianCoordinates& cartesian);
CartesianCoordinates sphericalToCartesian(const SphericalCoordinates& spherical);
double getlocalIncAngle(Vector3D& satepoint, Vector3D& landpoint, Vector3D& slopeVector);
#endif

2810
BaseTool/ImageOperatorBase.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

558
BaseTool/ImageOperatorBase.h Normal file
View File

@ -0,0 +1,558 @@
#pragma once
/**
* 影像数据操作项
* 所有数据相关的操作，都是基于ENVI的数据格式
* 影像读写操作时基于 GDAL 库
* **/
#ifndef IMAGEOPERATORBASE_H
#define IMAGEOPERATORBASE_H
#include "BaseConstVariable.h"
#include "GeoOperator.h"
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <complex>
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <string>
#include <memory>
#include <QString>
#include <cpl_conv.h> // for CPLMalloc()
#include "LogInfoCls.h"
#include <QObject>
enum ProjectStripDelta {
Strip_6, // 6度带
Strip_3
};
enum CoordinateSystemType { // 坐标系类型
GeoCoordinateSystem,
ProjectCoordinateSystem,
UNKNOW
};
struct PointRaster { // 影像坐标点
double x;
double y;
double z;
};
struct PointXYZ {
double x, y, z;
};
struct PointGeo {
double lon, lat, ati;
};
struct PointImage {
double pixel_x, pixel_y;
};
struct ImageGEOINFO {
int width;
int height;
int bandnum;
};
enum GDALREADARRCOPYMETHOD {
MEMCPYMETHOD, // 直接拷贝
VARIABLEMETHOD // 变量赋值
};
class ShowProessAbstract{
public:
virtual void showProcess(double precent,QString tip);
virtual void showToolInfo( QString tip) ;
};
/// 根据经纬度获取
/// EPSG代码，根据经纬度返回对应投影坐标系统，其中如果在中华人民共和国境内，默认使用
/// CGCS2000坐标系统(EPSG 4502 ~ 4512 6度带,EPSG 4534 ~ 4554 3度带)，其余地方使用WGS坐标系统，
/// 投影方法 高斯克吕格（国内）， 高斯克吕格
/// \param long 经度
/// \param lat 纬度
/// \return 对应投影坐标系统的 EPSG编码,-1 表示计算错误
long getProjectEPSGCodeByLon_Lat(double long, double lat,ProjectStripDelta stripState );
long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat);
long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat);
QString GetProjectionNameFromEPSG(long epsgCode);
long GetEPSGFromRasterFile(QString filepath);
std::shared_ptr<PointRaster> GetCenterPointInRaster(QString filepath);
CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long EPSGCODE);
// 文件打开 // 当指令销毁时，调用GDALClose 销毁类型
std::shared_ptr<GDALDataset> OpenDataset(const QString& in_path, GDALAccess rwmode= GA_ReadOnly);
void CloseDataset(GDALDataset* ptr);
// 数据格式转换
int TIFF2ENVI(QString in_tiff_path,QString out_envi_path);
int ENVI2TIFF(QString in_envi_path,QString out_tiff_path);
// 保存影像数据 --直接保存 ENVI 文件
int CreateDataset(QString new_file_path, int height, int width, int band_num, double* gt, QString projection, GDALDataType gdal_dtype, bool need_gt); // 创建文件
int saveDataset(QString new_file_path, int start_line, int start_cols, int band_ids, int datacols, int datarows, void* databuffer);
// 根据限制条件估算分块大小
int block_num_pre_memory(int width, int height, GDALDataType gdal_dtype,double memey_size);
// 将结果转换为复数 或者 实数
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> ReadComplexMatrixData(int start_line,int width, int line_num, std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype);
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> ReadMatrixDoubleData(int start_line, int width, int line_num, std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype,int band_idx);
Eigen::MatrixXd getGeoTranslationArray(QString in_path);
ImageGEOINFO getImageINFO(QString in_path);
GDALDataType getGDALDataType(QString fileptah);
struct RasterExtend {
double min_x; //纬度
double min_y;//经度
double max_x;//纬度
double max_y;//经度
};
/// <summary>
/// gdalImage图像操作类
/// </summary>
class gdalImage
{
public: // 方法
gdalImage();
gdalImage(const QString& raster_path);
~gdalImage();
virtual void setHeight(int);
virtual void setWidth(int);
virtual void setTranslationMatrix(Eigen::MatrixXd gt);
virtual void setData(Eigen::MatrixXd,int data_band_ids=1);
virtual Eigen::MatrixXd getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
virtual Eigen::MatrixXi getDatai(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
virtual Eigen::MatrixXd getGeoTranslation();
virtual GDALDataType getDataType();
virtual void saveImage(Eigen::MatrixXd, int start_row, int start_col, int band_ids);
virtual void saveImage(Eigen::MatrixXi, int start_row, int start_col, int band_ids);
virtual void saveImage();
virtual void setNoDataValue(double nodatavalue, int band_ids);
virtual void setNoDataValuei(int nodatavalue, int band_ids);
virtual double getNoDataValue(int band_ids);
virtual int getNoDataValuei(int band_ids);
virtual int InitInv_gt();
virtual Landpoint getRow_Col(double lon, double lat);
virtual Landpoint getLandPoint(double i, double j, double ati);
virtual void getLandPoint(double i, double j, double ati, Landpoint& Lp);
virtual double mean(int bandids = 1);
double BandmaxValue(int bandids = 1);
double BandminValue(int bandids = 1);
virtual GDALRPCInfo getRPC();
virtual Eigen::MatrixXd getLandPoint(Eigen::MatrixXd points);
virtual Eigen::MatrixXd getHist(int bandids);
virtual RasterExtend getExtend();
public:
QString img_path; // 图像文件
int height; // 高
int width; // 宽
int band_num;// 波段数
int start_row;//
int start_col;//
int data_band_ids;
Eigen::MatrixXd gt; // 变换矩阵
Eigen::MatrixXd inv_gt; // 逆变换矩阵
Eigen::MatrixXd data;
QString projection;
};
/// <summary>
/// gdalImage图像操作类
/// </summary>
class gdalImageComplex:public gdalImage
{
public: // 方法
gdalImageComplex(const QString& raster_path);
~gdalImageComplex();
void setData(Eigen::MatrixXcd);
void saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids);
Eigen::MatrixXcd getDataComplex(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
void saveImage() override;
void savePreViewImage();
public:
Eigen::MatrixXcd data;
};
// 创建影像
gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false, bool isEnvi = false);
gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long espgcode, GDALDataType eType=GDT_Float32, bool need_gt = true, bool overwrite = false,bool isENVI=false);
gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false);
gdalImageComplex CreateEchoComplex(const QString& img_path, int height, int width, int band_num);
ErrorCode DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QString demsloperPath);
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample);
void resampleRaster(const char* inputRaster, const char* outputRaster, double targetPixelSizeX, double targetPixelSizeY);
void cropRasterByLatLon(const char* inputFile, const char* outputFile, double minLon, double maxLon, double minLat, double maxLat);
int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample = GRIORA_Bilinear);
void transformRaster(const char* inputFile, const char* outputFile, int sourceEPSG, int targetEPSG);
ErrorCode transformCoordinate(double x, double y, int sourceEPSG, int targetEPSG, Point2& p);
int alignRaster(QString inputPath, QString referencePath, QString outputPath, GDALResampleAlg eResample);
//--------------------- 保存文博 -------------------------------
int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path);
//----------------------------------------------------
//--------------------- 图像合并流程 ------------------------------
enum MERGEMODE
{
MERGE_GEOCODING,
};
ErrorCode MergeRasterProcess(QVector<QString> filepath, QString outfileptah, QString mainString, MERGEMODE mergecode = MERGEMODE::MERGE_GEOCODING, bool isENVI = false, ShowProessAbstract* dia=nullptr);
ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> inimgs, gdalImage resultimg,gdalImage maskimg, ShowProessAbstract* dia = nullptr);
// 保存矩阵转换为envi文件，默认数据格式为double
bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath);
template<typename T>
std::shared_ptr<T> readDataArr(gdalImage& imgds, int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method)
{
std::shared_ptr<T> result = nullptr;
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(imgds.img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids);
rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row;
cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col;
//Eigen::MatrixXd datamatrix(rows_count, cols_count);
if (gdal_datatype == GDT_Byte) {
char* temp = new char[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_UInt16) {
unsigned short* temp = new unsigned short[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_Int16) {
short* temp = new short[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_UInt32) {
unsigned int* temp = new unsigned int[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_Int32) {
int* temp = new int[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
}
delete[] temp;
}
// else if (gdal_datatype == GDT_UInt64) {
// unsigned long* temp = new unsigned long[rows_count * cols_count];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
//rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { for (int j = 0; j <
//cols_count; j++) { datamatrix(i, j) = temp[i * cols_count + j];
// }
// }
// delete[] temp;
// }
// else if (gdal_datatype == GDT_Int64) {
// long* temp = new long[rows_count * cols_count];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
//rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) { for (int j = 0; j <
//cols_count; j++) { datamatrix(i, j) = temp[i * cols_count + j];
// }
// }
// delete[] temp;
// }
else if (gdal_datatype == GDT_Float32) {
float* temp = new float[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_Float64) {
double* temp = new double[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
}
delete[] temp;
}
//else if (gdal_datatype == GDT_CFloat32) {
// if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
// float* temp = new float[rows_count * cols_count * 2];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
// result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
// if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
// std::memcpy(result.get(), temp, rows_count * cols_count);
// }
// else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
// long count = rows_count * cols_count;
// for (long i = 0; i < count; i++) {
// result.get()[i] = T(temp[i * 2],
// temp[i * 2 + 1]);
// }
// }
// delete[] temp;
// }
// else {
// result = nullptr;
// }
//}
//else if (gdal_datatype == GDT_CFloat64 ) {
// if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
// double* temp = new double[rows_count * cols_count * 2];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
// result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
// if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
// std::memcpy(result.get(), temp, rows_count * cols_count);
// }
// else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
// long count = rows_count * cols_count;
// for (long i = 0; i < count; i++) {
// result.get()[i] = T(temp[i * 2],
// temp[i * 2 + 1]);
// }
// }
// delete[] temp;
// }
// else {
// result = nullptr;
// }
//}
else {
}
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return result;
}
template<typename T>
std::shared_ptr<T> readDataArrComplex(gdalImageComplex& imgds, int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method)
{
std::shared_ptr<T> result = nullptr;
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(imgds.img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids);
rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row;
cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col;
//Eigen::MatrixXd datamatrix(rows_count, cols_count);
if (gdal_datatype == GDT_CFloat32) {
if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
float* temp = new float[rows_count * cols_count * 2];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i * 2],
temp[i * 2 + 1]);
}
}
delete[] temp;
}
else {
result = nullptr;
}
}
else if (gdal_datatype == GDT_CFloat64) {
if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
double* temp = new double[rows_count * cols_count * 2];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i * 2],
temp[i * 2 + 1]);
}
}
delete[] temp;
}
else {
result = nullptr;
}
}
else {
}
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return result;
}
//--------------------- 图像分块 ------------------------------
#endif

233
BaseTool/LogInfoCls.cpp Normal file
View File

@ -0,0 +1,233 @@
#include "stdafx.h"
#include "LogInfoCls.h"
std::string errorCode2errInfo(ErrorCode e)
{
switch (e)
{
_CASE_STR(SUCCESS);
_CASE_STR(VIRTUALABSTRACT);
_CASE_STR(FAIL);
_CASE_STR(FILENOFOUND);
_CASE_STR(OrbitNodeNotEnough);
_CASE_STR(XYDataPointNotEqual);
_CASE_STR(FILEOPENFAIL );
_CASE_STR(XMLPARSEFAIL );
_CASE_STR(XMLNOTFOUNDElEMENT );
_CASE_STR(FILEPATHISEMPTY);
_CASE_STR(FOLDER_NOT_EXIST);
_CASE_STR(FILE_NOT_EXIST);
_CASE_STR(FIND_ID_ERROR);
_CASE_STR(INSERT_ID_ERROR);
_CASE_STR(EPSGCODE_NOTSAME);
_CASE_STR(EPSGCODE_NOTSUPPORT);
_CASE_STR(RASTERBAND_NOTEQUAL);
_CASE_STR(RASTER_DATETYPE_NOTSAME);
//GSL 1xx
_CASE_STR(Error_GSL_FAILURE );
_CASE_STR(Error_GSL_CONTINUE );
_CASE_STR(Error_GSL_EDOM );
_CASE_STR(Error_GSL_ERANGE );
_CASE_STR(Error_GSL_EFAULT );
_CASE_STR(Error_GSL_EINVAL );
_CASE_STR(Error_GSL_EFAILED );
_CASE_STR(Error_GSL_EFACTOR );
_CASE_STR(Error_GSL_ESANITY );
_CASE_STR(Error_GSL_ENOMEM );
_CASE_STR(Error_GSL_EBADFUNC );
_CASE_STR(Error_GSL_ERUNAWAY );
_CASE_STR(Error_GSL_EMAXITER );
_CASE_STR(Error_GSL_EZERODIV );
_CASE_STR(Error_GSL_EBADTOL );
_CASE_STR(Error_GSL_ETOL );
_CASE_STR(Error_GSL_EUNDRFLW );
_CASE_STR(Error_GSL_EOVRFLW );
_CASE_STR(Error_GSL_ELOSS );
_CASE_STR(Error_GSL_EROUND );
_CASE_STR(Error_GSL_EBADLEN );
_CASE_STR(Error_GSL_ENOTSQR );
_CASE_STR(Error_GSL_ESING );
_CASE_STR(Error_GSL_EDIVERGE );
_CASE_STR(Error_GSL_EUNSUP );
_CASE_STR(Error_GSL_EUNIMPL );
_CASE_STR(Error_GSL_ECACHE );
_CASE_STR(Error_GSL_ETABLE );
_CASE_STR(Error_GSL_ENOPROG );
_CASE_STR(Error_GSL_ENOPROGJ );
_CASE_STR(Error_GSL_ETOLF );
_CASE_STR(Error_GSL_ETOLX );
_CASE_STR(Error_GSL_ETOLG );
_CASE_STR(Error_GSL_EOF );
// RTPC
_CASE_STR(RTPC_PARAMSISEMPTY);
_CASE_STR(ECHO_L0DATA_NOTOPEN);
_CASE_STR(ECHO_L0DATA_ROW_COL_NOEQUAL);
_CASE_STR(ECHO_L0DATA_PRFIDXOUTRANGE);
_CASE_STR(ECHO_L0DATA_GPSFILEFORMATERROR);
_CASE_STR(ECHO_L0DATA_ECHOFILEFORMATERROR);
_CASE_STR(ECHO_L0DATA_ECHOFILENOTOPEN);
_CASE_STR(ECHO_L0DATA_GPSFILEFNOTOPEN);
_CASE_STR(ECHO_L0DATA_XMLFILENOTOPEN);
_CASE_STR(OUTOFRANGE);
_CASE_STR(ECHO_L0DATA_XMLNAMEERROR);
//
_CASE_STR(TBP_L0OPENFAIL);
//
_CASE_STR(IMAGE_L1DATA_XMLNAMEERROR);
_CASE_STR(IMAGE_L1DATA_XMLNAMEOPENERROR);
default:
break;
}
return "UNKNOW_EVENT!";
}
ErrorCode GSLState2ErrorCode(int gslState)
{
switch (gslState)
{
case 0:
return ErrorCode::SUCCESS;
break;
case -1:
return ErrorCode::Error_GSL_FAILURE;
break;
case -2:
return ErrorCode::Error_GSL_CONTINUE;
break;
case 1:
return ErrorCode::Error_GSL_EDOM; // 输入域错误，例如 sqrt(-1)
break;
case 2:
return ErrorCode::Error_GSL_ERANGE; // 输出范围错误，例如 exp(1e100)
break;
case 3:
return ErrorCode::Error_GSL_EFAULT; // 无效指针
break;
case 4:
return ErrorCode::Error_GSL_EINVAL; // 用户提供的无效参数
break;
case 5:
return ErrorCode::Error_GSL_EFAILED; // 通用失败
break;
case 6:
return ErrorCode::Error_GSL_EFACTOR; // 因子分解失败
break;
case 7:
return ErrorCode::Error_GSL_ESANITY; // 理智检查失败，通常不应该发生
break;
case 8:
return ErrorCode::Error_GSL_ENOMEM; // 内存分配失败
break;
case 9:
return ErrorCode::Error_GSL_EBADFUNC; // 用户提供的函数存在问题
break;
case 10:
return ErrorCode::Error_GSL_ERUNAWAY; // 迭代过程失控
break;
case 11:
return ErrorCode::Error_GSL_EMAXITER; // 超过最大迭代次数
break;
case 12:
return ErrorCode::Error_GSL_EZERODIV; // 尝试进行零除
break;
case 13:
return ErrorCode::Error_GSL_EBADTOL; // 用户指定的容忍度无效
break;
case 14:
return ErrorCode::Error_GSL_ETOL; // 未能达到指定的容忍度
break;
case 15:
return ErrorCode::Error_GSL_EUNDRFLW; // 下溢
break;
case 16:
return ErrorCode::Error_GSL_EOVRFLW; // 上溢
break;
case 17:
return ErrorCode::Error_GSL_ELOSS; // 精度损失
break;
case 18:
return ErrorCode::Error_GSL_EROUND; // 因舍入错误而失败
break;
case 19:
return ErrorCode::Error_GSL_EBADLEN; // 矩阵或向量长度不匹配
break;
case 20:
return ErrorCode::Error_GSL_ENOTSQR; // 矩阵不是方阵
break;
case 21:
return ErrorCode::Error_GSL_ESING; // 检测到明显的奇异性
break;
case 22:
return ErrorCode::Error_GSL_EDIVERGE; // 积分或级数发散
break;
case 23:
return ErrorCode::Error_GSL_EUNSUP; // 请求的特性不受硬件支持
break;
case 24:
return ErrorCode::Error_GSL_EUNIMPL; // 请求的特性尚未实现
break;
case 25:
return ErrorCode::Error_GSL_ECACHE; // 超过缓存限制
break;
case 26:
return ErrorCode::Error_GSL_ETABLE; // 超过表限制
break;
case 27:
return ErrorCode::Error_GSL_ENOPROG; // 迭代未能朝向解决方案取得进展
break;
case 28:
return ErrorCode::Error_GSL_ENOPROGJ; // 雅可比评估未能改善解决方案
break;
case 29:
return ErrorCode::Error_GSL_ETOLF; // 无法达到 F 的指定容忍度
break;
case 30:
return ErrorCode::Error_GSL_ETOLX; // 无法达到 X 的指定容忍度
break;
case 31:
return ErrorCode::Error_GSL_ETOLG; // 无法达到梯度的指定容忍度
break;
case 32:
return ErrorCode::Error_GSL_EOF; // 文件结束
break;
default:
return ErrorCode::Error_GSL_FAILURE; // 未知错误，返回一般失败
break;
}
}

102
BaseTool/LogInfoCls.h Normal file
View File

@ -0,0 +1,102 @@
#pragma once
/*****************************************************************//**
* \file LogInfoCls.h
* \brief 错误枚举类，以及返回错误信息
*
* \author 陈增辉
* \date October 2024
*********************************************************************/
#include <string>
// 定义变换
#define _CASE_STR(x) case x : return #x; break;
enum ErrorCode {
SUCCESS = -1,// 执行成功
VIRTUALABSTRACT = -2,// virtual abstract function ，not implement
FAIL=0,
FILENOFOUND = 1,
OrbitNodeNotEnough = 2,
XYDataPointNotEqual = 3,
FILEOPENFAIL = 4,
XMLPARSEFAIL = 5,
XMLNOTFOUNDElEMENT = 6,
FILEPATHISEMPTY = 7,
FOLDER_NOT_EXIST = 8,
FILE_NOT_EXIST = 9,
FIND_ID_ERROR = 10,
INSERT_ID_ERROR = 11,
EPSGCODE_NOTSAME = 12,
EPSGCODE_NOTSUPPORT = 13,
RASTERBAND_NOTEQUAL = 14,
RASTER_DATETYPE_NOTSAME = 15,
//GSL 1xx
Error_GSL_FAILURE = -101,
Error_GSL_CONTINUE = -102, /* iteration has not converged */
Error_GSL_EDOM = 101, /* input domain error, e.g sqrt(-1) */
Error_GSL_ERANGE = 102, /* output range error, e.g. exp(1e100) */
Error_GSL_EFAULT = 103, /* invalid pointer */
Error_GSL_EINVAL = 104, /* invalid argument supplied by user */
Error_GSL_EFAILED = 105, /* generic failure */
Error_GSL_EFACTOR = 106, /* factorization failed */
Error_GSL_ESANITY = 107, /* sanity check failed - shouldn't happen */
Error_GSL_ENOMEM = 108, /* malloc failed */
Error_GSL_EBADFUNC = 109, /* problem with user-supplied function */
Error_GSL_ERUNAWAY = 110, /* iterative process is out of control */
Error_GSL_EMAXITER = 111, /* exceeded max number of iterations */
Error_GSL_EZERODIV = 112, /* tried to divide by zero */
Error_GSL_EBADTOL = 113, /* user specified an invalid tolerance */
Error_GSL_ETOL = 114, /* failed to reach the specified tolerance */
Error_GSL_EUNDRFLW = 115, /* underflow */
Error_GSL_EOVRFLW = 116, /* overflow */
Error_GSL_ELOSS = 117, /* loss of accuracy */
Error_GSL_EROUND = 118, /* failed because of roundoff error */
Error_GSL_EBADLEN = 119, /* matrix, vector lengths are not conformant */
Error_GSL_ENOTSQR = 120, /* matrix not square */
Error_GSL_ESING = 121, /* apparent singularity detected */
Error_GSL_EDIVERGE = 122, /* integral or series is divergent */
Error_GSL_EUNSUP = 123, /* requested feature is not supported by the hardware */
Error_GSL_EUNIMPL = 124, /* requested feature not (yet) implemented */
Error_GSL_ECACHE = 125, /* cache limit exceeded */
Error_GSL_ETABLE = 126, /* table limit exceeded */
Error_GSL_ENOPROG = 127, /* iteration is not making progress towards solution */
Error_GSL_ENOPROGJ = 128, /* jacobian evaluations are not improving the solution */
Error_GSL_ETOLF = 129, /* cannot reach the specified tolerance in F */
Error_GSL_ETOLX = 130, /* cannot reach the specified tolerance in X */
Error_GSL_ETOLG = 131, /* cannot reach the specified tolerance in gradient */
Error_GSL_EOF = 132, /* end of file */
// RTPC
RTPC_PARAMSISEMPTY = 201,
// L0 数据
ECHO_L0DATA_NOTOPEN = 202,
ECHO_L0DATA_ROW_COL_NOEQUAL = 203, // 行或者列数量错误
ECHO_L0DATA_PRFIDXOUTRANGE = 204, // PRF 索引超出范围
ECHO_L0DATA_GPSFILEFORMATERROR = 205, // GPS文件错误
ECHO_L0DATA_ECHOFILEFORMATERROR = 206, // 回波文件格式错误
ECHO_L0DATA_ECHOFILENOTOPEN = 207, // 回波文件打开错误
ECHO_L0DATA_GPSFILEFNOTOPEN = 208, // GPS文件打开错误
ECHO_L0DATA_XMLFILENOTOPEN = 209, // xml文件打开错误
OUTOFRANGE = 210, // 超出范围
ECHO_L0DATA_XMLNAMEERROR = 211, // 超出范围
// BP成像
TBP_L0OPENFAIL = 301, // 0级文件打开错误
// L1图像
IMAGE_L1DATA_XMLNAMEERROR = 401,
IMAGE_L1DATA_XMLNAMEOPENERROR = 402,
IMAGE_L1DATA_XMLNAMEPARASEERROR = 403,
};
std::string errorCode2errInfo(ErrorCode code);
ErrorCode GSLState2ErrorCode(int gslState);

24
BaseTool/QToolProcessBarDialog.cpp Normal file
View File

@ -0,0 +1,24 @@
#include "QToolProcessBarDialog.h"
QToolProcessBarDialog::QToolProcessBarDialog(QWidget *parent)
: QDialog(parent)
{
ui.setupUi(this);
ui.progressBar->setRange(0, 100);
}
QToolProcessBarDialog::~QToolProcessBarDialog()
{}
void QToolProcessBarDialog::showProcess(double precent, QString tip)
{
ui.progressBar->setValue(std::ceil(precent * 100));
ui.labelTip->setText(tip);
this->update();
}
void QToolProcessBarDialog::showToolInfo(QString tip)
{
ui.textEditTip->append("\n"+tip);
}

19
BaseTool/QToolProcessBarDialog.h Normal file
View File

@ -0,0 +1,19 @@
#pragma once
#include <QDialog>
#include "ui_QToolProcessBarDialog.h"
#include "ImageOperatorBase.h"
class QToolProcessBarDialog : public QDialog, public ShowProessAbstract
{
Q_OBJECT
public:
QToolProcessBarDialog(QWidget *parent = nullptr);
~QToolProcessBarDialog();
private:
Ui::QToolProcessBarDialogClass ui;
public:
virtual void showProcess(double precent, QString tip) override;
virtual void showToolInfo(QString tip) override;
};

73
BaseTool/QToolProcessBarDialog.ui Normal file
View File

@ -0,0 +1,73 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>QToolProcessBarDialogClass</class>
<widget class="QDialog" name="QToolProcessBarDialogClass">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>600</width>
<height>400</height>
</rect>
</property>
<property name="windowTitle">
<string>QToolProcessBarDialog</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout">
<item>
<widget class="QWidget" name="widget" native="true">
<layout class="QGridLayout" name="gridLayout">
<item row="1" column="2">
<widget class="QPushButton" name="OKpushButton">
<property name="minimumSize">
<size>
<width>120</width>
<height>26</height>
</size>
</property>
<property name="maximumSize">
<size>
<width>120</width>
<height>26</height>
</size>
</property>
<property name="text">
<string>退出</string>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="labelTip">
<property name="text">
<string>提示</string>
</property>
</widget>
</item>
<item row="2" column="0" colspan="2">
<widget class="QProgressBar" name="progressBar">
<property name="value">
<number>24</number>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item>
<widget class="QGroupBox" name="groupBox">
<property name="title">
<string/>
</property>
<layout class="QVBoxLayout" name="verticalLayout_2">
<item>
<widget class="QTextEdit" name="textEditTip"/>
</item>
</layout>
</widget>
</item>
</layout>
</widget>
<layoutdefault spacing="6" margin="11"/>
<resources/>
<connections/>
</ui>

268
BaseTool/RasterToolBase.cpp Normal file
View File

@ -0,0 +1,268 @@
/**
* @file RasterProjectBase.cpp
* @brief None
* @author 陈增辉 (3045316072@qq.com)
* @version 2.5.0
* @date 24-6-4
* @copyright Copyright (c) Since 2024 中科卫星应用研究院 All rights reserved.
*/
#include <QDebug>
#include "RasterToolBase.h"
#include "gdal_priv.h"
#include "ogr_spatialref.h"
#include "cpl_conv.h" // for CPLMalloc()
#include <QTextCodec>
#include <iostream>
#include <QFile>
namespace RasterToolBase {
long getProjectEPSGCodeByLon_Lat(double lon, double lat, ProjectStripDelta stripState)
{
long EPSGCode = 0;
switch(stripState) {
case ProjectStripDelta::Strip_3: {
break;
};
case ProjectStripDelta::Strip_6: {
break;
}
default: {
EPSGCode = -1;
break;
}
}
qDebug() << QString(" EPSG code : %1").arg(EPSGCode);
return EPSGCode;
}
long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat)
{
// EPSG 4534 ~ 4554 3 度带
// 首先判断是否是在 中国带宽范围
// 中心经度范围 ：75E ~ 135E 实际范围 73.5E ~ 136.5E,
// 纬度范围 3N ~ 54N，放宽到 0N~ 60N
if(73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内
long code = trunc((lon - 73.5) / 3) + 4534;
return code;
} else { // 非中国境内 使用 高斯克吕格
bool isSouth = lat < 0; // 简单判断南北半球，这里仅为示例，实际应用可能需要更细致的逻辑
long prefix = isSouth ? 327000 : 326000;
// std::string prefix = isSouth ? "327" : "326";
lon = fmod(lon + 360.0, 360.0);
long zone = std::floor((lon + 180.0) / 3.0);
prefix = prefix + zone;
return prefix;
}
return 0;
}
long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat)
{
// EPSG 4502 ~ 4512 6度带
// 首先判断是否是在 中国带宽范围
// 中心经度范围 ：75E ~ 135E 实际范围 72.0E ~ 138E,
// 纬度范围 3N ~ 54N，放宽到 0N~ 60N
if(73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内
long code = trunc((lon - 72.0) / 6) + 4502;
return code;
} else { // 非中国境内 使用 UTM// 确定带号，6度带从1开始到60，每6度一个带
int zone = static_cast<int>((lon + 180.0) / 6.0) + 1;
bool isSouth = lon < 0; // 判断是否在南半球
long epsgCodeBase = isSouth ? 32700 : 32600; // 计算EPSG代码
long prefix = static_cast<int>(epsgCodeBase + zone);
return prefix;
}
return 0;
}
QString GetProjectionNameFromEPSG(long epsgCode)
{
qDebug() << "============= GetProjectionNameFromEPSG ======================";
OGRSpatialReference oSRS;
// 设置EPSG代码
if(oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) {
qDebug() << "epsgcode not recognition";
return "";
}
// 获取并输出坐标系的描述（名称）
const char* pszName = oSRS.GetAttrValue("GEOGCS");
if(pszName) {
qDebug() << "Coordinate system name for EPSG " + QString::number(epsgCode)
<< " is: " + QString::fromStdString(pszName);
return QString::fromStdString(pszName);
} else {
qDebug() << "Unable to retrieve the name for EPSG " + QString::number(epsgCode);
return "";
}
// char* wkt = NULL;
// // 转换为WKT格式
// oSRS.exportToWkt(&wkt);
//
// qDebug() << wkt;
//
// // 从WKT中解析投影名称，这里简化处理，实际可能需要更复杂的逻辑来准确提取名称
// std::string wktStr(wkt);
// long start = wktStr.find("PROJCS[\"") + 8; // 找到"PROJCS["后的第一个双引号位置
// // 从start位置开始找下一个双引号，这之间的内容即为投影名称
// int end = wktStr.find('\"', start);
// QString projName = QString::fromStdString(wktStr.substr(start, end - start));
//
// // 释放WKT字符串内存
// CPLFree(wkt);
// return projName;
}
long GetEPSGFromRasterFile(QString filepath)
{
qDebug() << "============= GetEPSGFromRasterFile ======================";
// QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器
// QByteArray byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray
//,这个应该会自动释放 const char* charArray = byteArray.constData(); //
// 获取QByteArray的const char*指针
{
if(QFile(filepath).exists()){
qDebug() << "info: the image found.\n";
}else{
return -1;
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注册GDAL驱动
// std::cout<<filepath.toLocal8Bit().constData()<<std::endl;
// 打开影像文件
GDALDataset* poDataset;
poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
if(NULL==poDataset) {
qDebug() << "Error: Unable to open the image.\n";
return -1;
}
// 获取影像的投影信息
const char* pszProjection = poDataset->GetProjectionRef();
qDebug() << QString::fromUtf8(pszProjection);
// 创建SpatialReference对象
OGRSpatialReference oSRS;
if(oSRS.importFromWkt((char**)&pszProjection) != OGRERR_NONE) {
qDebug() << ("Error: Unable to import projection information.\n");
GDALClose(poDataset);
return -1;
}
long epsgCode = atoi(oSRS.GetAuthorityCode(nullptr)); // 获取EPSG代码
if(epsgCode != 0) {
GDALClose(poDataset);
qDebug() << QString("file %1 :epsg Code %2").arg(filepath).arg(epsgCode);
return epsgCode;
} else {
qDebug() << "EPSG code could not be determined from the spatial reference.";
GDALClose(poDataset);
return -1;
}
}
}
std::shared_ptr<PointRaster> GetCenterPointInRaster(QString filepath)
{
qDebug() << "============= GetCenterPointInRaster ======================";
// QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器
// QByteArray byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray
//,这个应该会自动释放 const char* charArray = byteArray.constData(); //
// 获取QByteArray的const char*指针
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// std::cout<<filepath.toLocal8Bit().constData()<<std::endl;
GDALDataset* poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
if(nullptr==poDataset||NULL==poDataset) {
qDebug() << "Could not open dataset";
return nullptr; // 表示读取失败
}
double x, y, z;
bool flag = false;
{
double adfGeoTransform[6];
if(poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
qDebug() << "Failed to get GeoTransform";
return nullptr;
}
double dfWidth = poDataset->GetRasterXSize();
double dfHeight = poDataset->GetRasterYSize();
// 计算中心点坐标（像素坐标）
double dfCenterX = adfGeoTransform[0] + dfWidth * adfGeoTransform[1] / 2.0
+ dfHeight * adfGeoTransform[2] / 2.0;
double dfCenterY = adfGeoTransform[3] + dfWidth * adfGeoTransform[4] / 2.0
+ dfHeight * adfGeoTransform[5] / 2.0;
OGRSpatialReference oSRS;
oSRS.importFromWkt(poDataset->GetProjectionRef());
if(oSRS.IsGeographic()) {
qDebug() << "Center coords (already in geographic): (" + QString::number(dfCenterX)
+ ", " + QString::number(dfCenterY) + ")";
flag = true;
x = dfCenterX;
y = dfCenterY;
} else {
// 如果不是地理坐标系，转换到WGS84
OGRSpatialReference oSRS_WGS84;
oSRS_WGS84.SetWellKnownGeogCS("WGS84");
OGRCoordinateTransformation* poCT =
OGRCreateCoordinateTransformation(&oSRS, &oSRS_WGS84);
if(poCT == nullptr) {
qDebug() << "Failed to create coordinate transformation";
return nullptr;
}
// double dfLon, dfLat;
if(poCT->Transform(1, &dfCenterX, &dfCenterY)) {
qDebug() << "Center coords (transformed to WGS84): ("
+ QString::number(dfCenterX) + ", " + QString::number(dfCenterY)
<< ")";
flag = true;
x = dfCenterX;
y = dfCenterY;
} else {
qDebug() << "Transformation failed.";
}
OGRCoordinateTransformation::DestroyCT(poCT);
}
}
if(nullptr==poDataset||NULL==poDataset){}else{
GDALClose(poDataset);
}
if(flag) {
std::shared_ptr<PointRaster> RasterCenterPoint = std::make_shared<PointRaster>();
RasterCenterPoint->x = x;
RasterCenterPoint->y = y;
RasterCenterPoint->z = 0;
return RasterCenterPoint;
} else {
return nullptr;
}
}
CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long epsg_code)
{
OGRSpatialReference oSRS;
if(oSRS.importFromEPSG(epsg_code) == OGRERR_NONE) {
if(oSRS.IsGeographic()) {
return CoordinateSystemType::GeoCoordinateSystem;
} else if(oSRS.IsProjected()) {
return CoordinateSystemType::ProjectCoordinateSystem;
}
} else {
return CoordinateSystemType::UNKNOW;
}
}
} // namespace RasterToolBase

79
BaseTool/RasterToolBase.h Normal file
View File

@ -0,0 +1,79 @@
/**
* @file RasterProjectBase.h
* @brief None
* @author 陈增辉 (3045316072@qq.com)
* @version 2.5.0
* @date 24-6-4
* @copyright Copyright (c) Since 2024 中科卫星应用研究院 All rights reserved.
*/
#ifndef LAMPCAE_RASTERTOOLBASE_H
#define LAMPCAE_RASTERTOOLBASE_H
#include "gdal_priv.h"
#include <memory>
namespace RasterToolBase {
static bool GDALAllRegisterEnable=false;
enum ProjectStripDelta{
Strip_6, // 6度带
Strip_3
};
enum CoordinateSystemType{ // 坐标系类型
GeoCoordinateSystem,
ProjectCoordinateSystem,
UNKNOW
};
struct PointRaster{ // 影像坐标点
double x;
double y;
double z;
};
struct PointXYZ{
double x,y,z;
};
struct PointGeo{
double lon,lat,ati;
};
struct PointImage{
double pixel_x,pixel_y;
};
/// 根据经纬度获取
/// EPSG代码，根据经纬度返回对应投影坐标系统，其中如果在中华人民共和国境内，默认使用
/// CGCS2000坐标系统(EPSG 4502 ~ 4512 6度带,EPSG 4534 ~ 4554 3度带)，其余地方使用WGS坐标系统，
/// 投影方法 高斯克吕格（国内）， 高斯克吕格
/// \param long 经度
/// \param lat 纬度
/// \return 对应投影坐标系统的 EPSG编码,-1 表示计算错误
long getProjectEPSGCodeByLon_Lat(double long, double lat,
ProjectStripDelta stripState = ProjectStripDelta::Strip_3);
long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat);
long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat);
QString GetProjectionNameFromEPSG(long epsgCode) ;
long GetEPSGFromRasterFile(QString filepath);
std::shared_ptr<PointRaster> GetCenterPointInRaster(QString filepath);
CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long EPSGCODE);
} // namespace RasterProjectConvertor
#endif // LAMPCAE_RASTERTOOLBASE_H

985
BaseTool/SARSimulationImageL1.cpp Normal file
View File

@ -0,0 +1,985 @@
#include "stdafx.h"
#include "SARSimulationImageL1.h"
#include "LogInfoCls.h"
#include <memory>
#include <QFile>
#include <QDomDocument>
#include <QXmlStreamReader>
SARSimulationImageL1Dataset::SARSimulationImageL1Dataset(RasterLevel level)
{
this->Rasterlevel = level;
}
SARSimulationImageL1Dataset::~SARSimulationImageL1Dataset()
{
}
RasterLevel SARSimulationImageL1Dataset::getRasterLevel()
{
return this->Rasterlevel;
}
QString SARSimulationImageL1Dataset::getoutFolderPath()
{
return outFolderPath;
}
QString SARSimulationImageL1Dataset::getDatesetName()
{
return L1DatasetName;
}
QString SARSimulationImageL1Dataset::getxmlfileName()
{
return xmlfileName;
}
QString SARSimulationImageL1Dataset::getxmlFilePath()
{
return xmlFilePath;
}
QString SARSimulationImageL1Dataset::getImageRasterName()
{
return ImageRasterName;
}
QString SARSimulationImageL1Dataset::getImageRasterPath()
{
return ImageRasterPath;
}
QString SARSimulationImageL1Dataset::getGPSPointFilename()
{
if (this->Rasterlevel == RasterLevel::RasterL2) {
return "";
}
return GPSPointFilename;
}
QString SARSimulationImageL1Dataset::getGPSPointFilePath()
{
if (this->Rasterlevel == RasterLevel::RasterL2) {
return "";
}
return GPSPointFilePath;
}
ErrorCode SARSimulationImageL1Dataset::OpenOrNew(QString folder, QString filename, long inrowCount, long incolCount)
{
qDebug() << "--------------Image Raster L1 Data OpenOrNew ---------------------------------------";
qDebug() << "Folder: " << folder;
qDebug() << "Filename: " << filename;
QDir dir(folder);
if (dir.exists() == false)
{
dir.mkpath(folder);
}
else {}
if (dir.exists() == false) {
return ErrorCode::FOLDER_NOT_EXIST;
}
else {}
QString filePath = dir.filePath(filename); // 生成完整的文件路径
this->rowCount = inrowCount;
this->colCount = incolCount;
this->outFolderPath = folder;
this->L1DatasetName = filename;
this->xmlfileName = filename + ".xml";
this->GPSPointFilename = filename + ".gpspos.data";
this->ImageRasterName = filename + ".bin";
this->xmlFilePath = dir.filePath(this->xmlfileName);
this->GPSPointFilePath = dir.filePath(this->GPSPointFilename);
this->ImageRasterPath = dir.filePath(this->ImageRasterName);
if (QFile(this->xmlFilePath).exists())
{
this->loadFromXml();
}
else
{
this->saveToXml();
}
if (this->Rasterlevel!=RasterL2||QFile(this->GPSPointFilePath).exists() == false) {
// 创建新文件
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
std::shared_ptr<GDALDataset> 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);
}
if (this->Rasterlevel == RasterLevel::RasterSLC || QFile(this->ImageRasterPath).exists() == false) {
// 创建新文件
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(this->ImageRasterPath.toUtf8().constData(), colCount, rowCount, 1, GDT_CFloat64, NULL));
GDALFlushCache((GDALDatasetH)poDstDS.get());
poDstDS.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
if (this->Rasterlevel != RasterLevel::RasterSLC || QFile(this->ImageRasterPath).exists() == false) {
// 创建新文件
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(this->ImageRasterPath.toUtf8().constData(), colCount, rowCount, 1, GDT_CFloat32, NULL));
GDALFlushCache((GDALDatasetH)poDstDS.get());
poDstDS.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
return ErrorCode::SUCCESS;
}
ErrorCode SARSimulationImageL1Dataset::Open(QString folder, QString filename)
{
QDir dir(folder);
if (dir.exists() == false)
{
return ErrorCode::FOLDER_NOT_EXIST;
}
else {}
if (dir.exists() == false) {
return ErrorCode::FOLDER_NOT_EXIST;
}
else {}
QString filePath = dir.filePath(filename); // 生成完整的文件路径
this->outFolderPath = folder;
this->L1DatasetName = filename;
this->xmlfileName = filename + ".xml";
this->GPSPointFilename = filename + ".gpspos.data";
this->ImageRasterName = filename + ".bin";
this->xmlFilePath = dir.filePath(this->xmlfileName);
this->GPSPointFilePath = dir.filePath(this->GPSPointFilename);
this->ImageRasterPath = dir.filePath(this->ImageRasterName);
if (QFile(this->xmlFilePath).exists())
{
this->loadFromXml();
}
else
{
return ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR;
}
}
ErrorCode SARSimulationImageL1Dataset::Open(QString xmlPath)
{
QFileInfo fileInfo(xmlPath);
QString fileName = fileInfo.fileName(); // 获取文件名
QString fileSuffix = fileInfo.suffix(); // 获取后缀名
QString fileNameWithoutSuffix = fileInfo.baseName(); // 获取不带后缀的文件名
QString directoryPath = fileInfo.path(); // 获取文件夹路径
if (fileSuffix.toLower() == "xml" || fileSuffix.toLower() == ".xml") {
return this->Open(directoryPath, fileNameWithoutSuffix);
}
else {
return ErrorCode::IMAGE_L1DATA_XMLNAMEERROR;
}
return ErrorCode::SUCCESS;
}
void SARSimulationImageL1Dataset::saveToXml()
{
QFile file(this->xmlFilePath);
if (!file.open(QIODevice::WriteOnly | QIODevice::Text)) {
qDebug() << "Cannot open file for writing:" << file.errorString();
return;
}
QXmlStreamWriter xmlWriter(&file);
xmlWriter.setAutoFormatting(true);
xmlWriter.writeStartDocument();
xmlWriter.writeStartElement("Parameters");
switch (this->Rasterlevel)
{
case(RasterLevel::RasterSLC):
xmlWriter.writeTextElement("rasterlevel", "SLC");
break;
case(RasterLevel::RasterL1B):
xmlWriter.writeTextElement("rasterlevel", "L1B");
break;
case(RasterLevel::RasterL2):
xmlWriter.writeTextElement("rasterlevel", "L2");
break;
default:
break;
}
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("LookSide", this->LookSide);
// 保存sateantpos
xmlWriter.writeStartElement("AntPos");
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.writeEndElement(); // 结束 <AntPosParam>
}
xmlWriter.writeTextElement("ImageStartTime", QString::number(this->startImageTime));
xmlWriter.writeTextElement("ImageEndTime", QString::number(this->EndImageTime));
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.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.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.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.writeEndElement(); // Parameters
xmlWriter.writeEndDocument();
file.close();
}
ErrorCode SARSimulationImageL1Dataset::loadFromXml()
{
QFile file(this->xmlFilePath);
if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) {
qDebug() << "Cannot open file for reading:" << file.errorString();
return ErrorCode::IMAGE_L1DATA_XMLNAMEOPENERROR;
}
QXmlStreamReader xmlReader(&file);
while (!xmlReader.atEnd() && !xmlReader.hasError()) {
QXmlStreamReader::TokenType token = xmlReader.readNext();
if (token == QXmlStreamReader::StartDocument) {
continue;
}
if (token == QXmlStreamReader::StartElement) {
if (xmlReader.name() == "Parameters") {
continue;
}
else if (xmlReader.name() == "rasterlevel") {
QString rasterlevelstr = xmlReader.readElementText();
if (rasterlevelstr == "SLC") {
this->Rasterlevel = RasterLevel::RasterSLC;
}
else if (rasterlevelstr == "L1B") {
this->Rasterlevel = RasterLevel::RasterL1B;
}
else if (rasterlevelstr == "L2") {
this->Rasterlevel = RasterLevel::RasterL2;
}
}
else if (xmlReader.name() == "RowCount") {
this->rowCount = xmlReader.readElementText().toLong();
}
else if (xmlReader.name() == "ColCount") {
this->colCount = xmlReader.readElementText().toLong();
}
else if (xmlReader.name() == "Rnear") {
this->Rnear = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Rfar") {
this->Rfar = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Rref") {
this->Rref = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "CenterFreq") {
this->centerFreq = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Fs") {
this->Fs = xmlReader.readElementText().toDouble();
}
else if(xmlReader.name() == "ImageStartTime"){
this->startImageTime = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "ImageEndTime") {
this->EndImageTime = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "CenterAngle") {
this->CenterAngle = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "LookSide") {
this->LookSide = xmlReader.readElementText();
}
else if (xmlReader.name() == "AntPosParam") {
SatelliteAntPos antPosParam;
while (!(xmlReader.isEndElement() && xmlReader.name() == "AntPosParam")) {
if (xmlReader.isStartElement()) {
if (xmlReader.name() == "time") {
antPosParam.time = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Px") {
antPosParam.Px = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Py") {
antPosParam.Py = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Pz") {
antPosParam.Pz = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Vx") {
antPosParam.Vx = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Vy") {
antPosParam.Vy = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "Vz") {
antPosParam.Vz = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "AntDirectX") {
antPosParam.AntDirectX = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "AntDirectY") {
antPosParam.AntDirectY = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "AntDirectZ") {
antPosParam.AntDirectZ = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "AVx") {
antPosParam.AVx = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "AVy") {
antPosParam.AVy = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "AVz") {
antPosParam.AVz = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "lon") {
antPosParam.lon = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "lat") {
antPosParam.lat = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "ati") {
antPosParam.ati = xmlReader.readElementText().toDouble();
}
}
xmlReader.readNext();
}
sateposes.append(antPosParam); // 将解析的数据添加到列表中
}
else if (xmlReader.name() == "incidenceAngleNearRange") {
incidenceAngleNearRange = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "incidenceAngleFarRange") {
incidenceAngleFarRange = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "TotalProcessedAzimuthBandWidth") {
this->TotalProcessedAzimuthBandWidth = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "DopplerParametersReferenceTime") {
this->DopplerParametersReferenceTime = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "DopplerCentroidCoefficients") {
while (!(xmlReader.tokenType() == QXmlStreamReader::EndElement && xmlReader.name() == "DopplerCentroidCoefficients")) {
if (xmlReader.tokenType() == QXmlStreamReader::StartElement) {
if (xmlReader.name() == "d0") {
this->d0 = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "d1") {
this->d1 = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "d2") {
this->d2 = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "d3") {
this->d3 = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "d4") {
this->d4 = xmlReader.readElementText().toDouble();
}
}
xmlReader.readNext();
}
}
else if (xmlReader.name() == "DopplerRateValuesCoefficients") {
while (!(xmlReader.tokenType() == QXmlStreamReader::EndElement && xmlReader.name() == "DopplerRateValuesCoefficients")) {
if (xmlReader.tokenType() == QXmlStreamReader::StartElement) {
if (xmlReader.name() == "r0") {
this->r0 = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "r1") {
this->r1 = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "r2") {
this->r2 = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "r3") {
this->r3 = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "r4") {
this->r4 = xmlReader.readElementText().toDouble();
}
}
xmlReader.readNext();
}
}
else if (xmlReader.name() == "latitude_center") {
this->latitude_center = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "longitude_center") {
this->longitude_center = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "latitude_topLeft") {
this->latitude_topLeft = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "longitude_topLeft") {
this->longitude_topLeft = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "latitude_topRight") {
this->latitude_topRight = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "longitude_topRight") {
this->longitude_topRight = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "latitude_bottomLeft") {
this->latitude_bottomLeft = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "longitude_bottomLeft") {
this->longitude_bottomLeft = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "latitude_bottomRight") {
this->latitude_bottomRight = xmlReader.readElementText().toDouble();
}
else if (xmlReader.name() == "longitude_bottomRight") {
this->longitude_bottomRight = xmlReader.readElementText().toDouble();
}
}
}
if (xmlReader.hasError()) {
qDebug() << "XML error:" << xmlReader.errorString();
return ErrorCode::IMAGE_L1DATA_XMLNAMEPARASEERROR;
}
file.close();
return ErrorCode::SUCCESS;
}
std::shared_ptr<double> SARSimulationImageL1Dataset::getAntPos()
{
if (this->Rasterlevel == RasterLevel::RasterL2) {
return nullptr;
}
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
// 读取文件
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(this->GPSPointFilePath, GDALAccess::GA_ReadOnly);
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(1);
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
std::shared_ptr<double> temp = nullptr;
if (gdal_datatype == GDT_Float64) {
temp = std::shared_ptr<double>(new double[this->rowCount * 19], delArrPtr);
demBand->RasterIO(GF_Read, 0, 0, 19, this->rowCount, temp.get(), 19, this->rowCount, gdal_datatype, 0, 0);
}
else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR));
}
rasterDataset.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return temp;
}
ErrorCode SARSimulationImageL1Dataset::saveAntPos(std::shared_ptr<double> ptr)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
// 读取文件
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(this->GPSPointFilePath, GDALAccess::GA_Update);
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(1);
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
if (gdal_datatype == GDT_Float64) {
demBand->RasterIO(GF_Write, 0, 0, 19, this->rowCount, ptr.get(), 19, this->rowCount, gdal_datatype, 0, 0);
}
else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR));
omp_unset_lock(&lock);
omp_destroy_lock(&lock);
return ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR;
}
rasterDataset.reset();
omp_unset_lock(&lock);
omp_destroy_lock(&lock);
return ErrorCode::SUCCESS;
}
std::shared_ptr<std::complex<double>> SARSimulationImageL1Dataset::getImageRaster()
{
if (this->Rasterlevel != RasterLevel::RasterSLC) {
return nullptr;
}
return this->getImageRaster(0, this->rowCount);
}
ErrorCode SARSimulationImageL1Dataset::saveImageRaster(std::shared_ptr<std::complex<double>> echoPtr, long startPRF, long PRFLen)
{
if (!(startPRF < this->rowCount)) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE));
return ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE;
}
else {}
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(this->ImageRasterPath, GDALAccess::GA_Update);
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 ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR;
}
else {}
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
if (height != this->rowCount || width != this->colCount) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR;
}
else {
if (gdal_datatype == GDT_CFloat64) {
poBand->RasterIO(GF_Write, 0, startPRF, width, PRFLen, echoPtr.get(), width, PRFLen, GDT_CFloat64, 0, 0);
GDALFlushCache((GDALDatasetH)rasterDataset.get());
}
else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
}
}
rasterDataset.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return ErrorCode::SUCCESS;
}
std::shared_ptr<std::complex<double>> SARSimulationImageL1Dataset::getImageRaster(long startPRF, long PRFLen)
{
if (this->Rasterlevel != RasterLevel::RasterSLC) {
return nullptr;
}
if (!(startPRF < this->rowCount)) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE)) << startPRF << " " << this->rowCount;
return nullptr;
}
else {}
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(this->ImageRasterPath, 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 nullptr;
}
else {}
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
std::shared_ptr<std::complex<double>> temp = nullptr;
if (height != this->rowCount || width != this->colCount) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
}
else {
if (gdal_datatype == GDT_CFloat64) {
temp = std::shared_ptr<std::complex<double>>(new std::complex<double>[PRFLen * width], delArrPtr);
poBand->RasterIO(GF_Read, 0, startPRF, width, PRFLen, temp.get(), width, PRFLen, GDT_CFloat64, 0, 0);
GDALFlushCache((GDALDatasetH)rasterDataset.get());
}
else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
}
}
rasterDataset.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
return temp;
}
Eigen::MatrixXcd SARSimulationImageL1Dataset::getImageRasterMatrix()
{
if (this->Rasterlevel != RasterLevel::RasterSLC) {
return Eigen::MatrixXcd::Zero(0,0);
}
std::shared_ptr<std::complex<double>> data = this->getImageRaster();
Eigen::MatrixXcd dataimg = Eigen::MatrixXcd::Zero(this->rowCount, this->colCount);
for (long i = 0; i < this->rowCount; i++) {
for (long j = 0; j < this->colCount; j++) {
dataimg(i, j) = data.get()[i*colCount+j];
}
}
return Eigen::MatrixXcd();
}
ErrorCode SARSimulationImageL1Dataset::saveImageRaster(Eigen::MatrixXcd& dataimg, long startPRF)
{
std::shared_ptr<std::complex<double>> data(new std::complex<double>[dataimg.rows()* dataimg.cols()]);
long dataimgrows = dataimg.rows();
long dataimgcols = dataimg.cols();
for (long i = 0; i < dataimgrows; i++) {
for (long j = 0; j < dataimgcols; j++) {
data.get()[i * dataimgcols + j] = dataimg(i, j);
}
}
this->saveImageRaster(data, startPRF,dataimgrows);
return ErrorCode();
}
gdalImage SARSimulationImageL1Dataset::getImageRasterGdalImage()
{
return gdalImage(this->ImageRasterPath);
}
void SARSimulationImageL1Dataset::setStartImageTime(double imageTime)
{
this->startImageTime = imageTime;
}
double SARSimulationImageL1Dataset::getStartImageTime()
{
return this->startImageTime;
}
void SARSimulationImageL1Dataset::setEndImageTime(double imageTime)
{
this->EndImageTime = imageTime;
}
double SARSimulationImageL1Dataset::getEndImageTime()
{
return this->EndImageTime;
}
QVector<SatelliteAntPos> SARSimulationImageL1Dataset::getXmlSateAntPos()
{
if (this->Rasterlevel == RasterLevel::RasterL2) {
return QVector<SatelliteAntPos>(0);
}
return this->sateposes;
}
void SARSimulationImageL1Dataset::setSateAntPos(QVector<SatelliteAntPos> pos)
{
this->sateposes = pos;
}
// Getter 和 Setter 方法实现
long SARSimulationImageL1Dataset::getrowCount() { return this->rowCount; }
void SARSimulationImageL1Dataset::setrowCount(long rowCountin) { this->rowCount = rowCountin; }
long SARSimulationImageL1Dataset::getcolCount() { return this->colCount; }
void SARSimulationImageL1Dataset::setcolCount(long pulsePointsin) { this->colCount = pulsePointsin; }
double SARSimulationImageL1Dataset::getNearRange() { return this->Rnear; }
void SARSimulationImageL1Dataset::setNearRange(double nearRange) { this->Rnear = nearRange; }
double SARSimulationImageL1Dataset::getFarRange() { return this->Rfar; }
void SARSimulationImageL1Dataset::setFarRange(double farRange) { this->Rfar = farRange; }
double SARSimulationImageL1Dataset::getRefRange()
{
return this->Rref;
}
void SARSimulationImageL1Dataset::setRefRange(double refRange)
{
this->Rref = refRange;
}
double SARSimulationImageL1Dataset::getCenterFreq() { return this->centerFreq; }
void SARSimulationImageL1Dataset::setCenterFreq(double freq) { this->centerFreq = freq; }
double SARSimulationImageL1Dataset::getFs() { return this->Fs; }
void SARSimulationImageL1Dataset::setFs(double samplingFreq) { this->Fs = samplingFreq; }
double SARSimulationImageL1Dataset::getPRF()
{
return this->prf;
}
void SARSimulationImageL1Dataset::setPRF(double PRF)
{
this->prf = PRF;
}
double SARSimulationImageL1Dataset::getCenterAngle()
{
return this->CenterAngle;
}
void SARSimulationImageL1Dataset::setCenterAngle(double angle)
{
this->CenterAngle = angle;
}
QString SARSimulationImageL1Dataset::getLookSide()
{
return this->LookSide;
}
void SARSimulationImageL1Dataset::setLookSide(QString looksides)
{
this->LookSide = looksides;
}
double SARSimulationImageL1Dataset::getTotalProcessedAzimuthBandWidth()
{
return TotalProcessedAzimuthBandWidth;
}
void SARSimulationImageL1Dataset::setTotalProcessedAzimuthBandWidth(double v)
{
TotalProcessedAzimuthBandWidth = v;
}
double SARSimulationImageL1Dataset::getDopplerParametersReferenceTime()
{
return DopplerParametersReferenceTime;
}
void SARSimulationImageL1Dataset::setDopplerParametersReferenceTime(double v)
{
DopplerParametersReferenceTime = v;
}
QVector<double> SARSimulationImageL1Dataset::getDopplerParams()
{
QVector<double> result(5);
result[0] = d0;
result[1] = d1;
result[2] = d2;
result[3] = d3;
result[4] = d4;
return result;
}
void SARSimulationImageL1Dataset::setDopplerParams(double id0, double id1, double id2, double id3, double id4)
{
this->d0 = id0;
this->d1 = id1;
this->d2 = id2;
this->d3 = id3;
this->d4 = id4;
}
QVector<double> SARSimulationImageL1Dataset::getDopplerCenterCoff()
{
QVector<double> result(5);
result[0]=r0;
result[1]=r1;
result[2]=r2;
result[3]=r3;
result[4]=r4;
return result;
}
void SARSimulationImageL1Dataset::setDopplerCenterCoff(double ir0, double ir1, double ir2, double ir3, double ir4)
{
this->r0 = ir0;
this->r1 = ir1;
this->r2 = ir2;
this->r3 = ir3;
this->r4 = ir4;
}
double SARSimulationImageL1Dataset::getIncidenceAngleNearRange()
{
return incidenceAngleNearRange;
}
void SARSimulationImageL1Dataset::setIncidenceAngleNearRangeet(double angle)
{
this->incidenceAngleNearRange = angle;
}
double SARSimulationImageL1Dataset::getIncidenceAngleFarRange()
{
return incidenceAngleFarRange;
}
void SARSimulationImageL1Dataset::setIncidenceAngleFarRange(double angle)
{
this->incidenceAngleFarRange = angle;
}
double SARSimulationImageL1Dataset::getLatitudeCenter() { return latitude_center; }
void SARSimulationImageL1Dataset::setLatitudeCenter(double value) { latitude_center = value; }
double SARSimulationImageL1Dataset::getLongitudeCenter() { return longitude_center; }
void SARSimulationImageL1Dataset::setLongitudeCenter(double value) { longitude_center = value; }
double SARSimulationImageL1Dataset::getLatitudeTopLeft() { return latitude_topLeft; }
void SARSimulationImageL1Dataset::setLatitudeTopLeft(double value) { latitude_topLeft = value; }
double SARSimulationImageL1Dataset::getLongitudeTopLeft() { return longitude_topLeft; }
void SARSimulationImageL1Dataset::setLongitudeTopLeft(double value) { longitude_topLeft = value; }
double SARSimulationImageL1Dataset::getLatitudeTopRight() { return latitude_topRight; }
void SARSimulationImageL1Dataset::setLatitudeTopRight(double value) { latitude_topRight = value; }
double SARSimulationImageL1Dataset::getLongitudeTopRight() { return longitude_topRight; }
void SARSimulationImageL1Dataset::setLongitudeTopRight(double value) { longitude_topRight = value; }
double SARSimulationImageL1Dataset::getLatitudeBottomLeft() { return latitude_bottomLeft; }
void SARSimulationImageL1Dataset::setLatitudeBottomLeft(double value) { latitude_bottomLeft = value; }
double SARSimulationImageL1Dataset::getLongitudeBottomLeft() { return longitude_bottomLeft; }
void SARSimulationImageL1Dataset::setLongitudeBottomLeft(double value) { longitude_bottomLeft = value; }
double SARSimulationImageL1Dataset::getLatitudeBottomRight() { return latitude_bottomRight; }
void SARSimulationImageL1Dataset::setLatitudeBottomRight(double value) { latitude_bottomRight = value; }
double SARSimulationImageL1Dataset::getLongitudeBottomRight() { return longitude_bottomRight; }
void SARSimulationImageL1Dataset::setLongitudeBottomRight(double value) { longitude_bottomRight = value; }
DemBox SARSimulationImageL1Dataset::getExtend()
{
double minlon = 0, maxlon = 0;
double minlat = 0, maxlat = 0;
minlon = this->longitude_bottomLeft < this->longitude_bottomRight ? this->longitude_bottomLeft : this->longitude_bottomRight;
minlon = minlon < this->longitude_topLeft ? minlon : this->longitude_topLeft;
minlon = minlon < this->longitude_topRight ? minlon : this->longitude_topRight;
minlat = this->latitude_bottomLeft < this->latitude_bottomRight ? this->latitude_bottomLeft : this->latitude_bottomRight;
minlat = minlat < this->latitude_topLeft ? minlat : this->latitude_topLeft;
minlat = minlat < this->latitude_bottomRight ? minlat : this->latitude_bottomRight;
maxlon = this->longitude_bottomLeft > this->longitude_bottomRight ? this->longitude_bottomLeft : this->longitude_bottomRight;
maxlon = maxlon > this->longitude_topLeft ? maxlon : this->longitude_topLeft;
maxlon = maxlon > this->longitude_topRight ? maxlon : this->longitude_topRight;
maxlat = this->latitude_bottomLeft > this->latitude_bottomRight ? this->latitude_bottomLeft : this->latitude_bottomRight;
maxlat = maxlat > this->latitude_topLeft ? maxlat : this->latitude_topLeft;
maxlat = maxlat > this->latitude_bottomRight ? maxlat : this->latitude_bottomRight;
DemBox box;
box.min_lat = minlat;
box.min_lon = minlon;
box.max_lat = maxlat;
box.max_lon = maxlon;
return box;
}

199
BaseTool/SARSimulationImageL1.h Normal file
View File

@ -0,0 +1,199 @@
#pragma once
#include "BaseConstVariable.h"
#include "BaseTool.h"
#include "ImageOperatorBase.h"
#include "GeoOperator.h"
#include "FileOperator.h"
#include "LogInfoCls.h"
enum RasterLevel {
RasterSLC,
RasterL1B,
RasterL2
};
class SARSimulationImageL1Dataset
{
public:
SARSimulationImageL1Dataset(RasterLevel Rasterlevel= RasterLevel::RasterSLC);
~SARSimulationImageL1Dataset();
private:
RasterLevel Rasterlevel;
public:
RasterLevel getRasterLevel();
private :
QString outFolderPath;
QString L1DatasetName;
QString xmlfileName;
QString xmlFilePath;
QString ImageRasterName;
QString ImageRasterPath;
QString GPSPointFilename;
QString GPSPointFilePath;
public:
QString getoutFolderPath();
QString getDatesetName();
QString getxmlfileName();
QString getxmlFilePath();
QString getImageRasterName();
QString getImageRasterPath();
QString getGPSPointFilename(); // GPS点
QString getGPSPointFilePath();
public:
ErrorCode OpenOrNew(QString folder, QString filename, long rowCount, long colCount);
ErrorCode Open(QString folderPath, QString Name);
ErrorCode Open(QString xmlPath);
public:
void saveToXml();
ErrorCode loadFromXml();
std::shared_ptr<double> getAntPos();
ErrorCode saveAntPos(std::shared_ptr<double> ptr); // 注意这个方法很危险，请写入前检查数据是否正确
std::shared_ptr<std::complex<double>> getImageRaster();
ErrorCode saveImageRaster(std::shared_ptr<std::complex<double>> echoPtr, long startPRF, long PRFLen);
std::shared_ptr<std::complex<double>> getImageRaster(long startPRF, long PRFLen);
Eigen::MatrixXcd getImageRasterMatrix();
ErrorCode saveImageRaster(Eigen::MatrixXcd& data, long startPRF);
gdalImage getImageRasterGdalImage();
private: // xml中参数
long rowCount;
long colCount;
double Rnear;
double Rfar;
double Rref;
double centerFreq;
double Fs;
double prf;
double CenterAngle;
QString LookSide;
QVector<SatelliteAntPos> sateposes;
double startImageTime;
double EndImageTime;
public:
void setStartImageTime(double imageTime);
double getStartImageTime();
void setEndImageTime(double imageTime);
double getEndImageTime();
QVector<SatelliteAntPos> getXmlSateAntPos();
void setSateAntPos(QVector<SatelliteAntPos> pos);
long getrowCount();
void setrowCount(long rowCount);
long getcolCount();
void setcolCount(long pulsePoints);
double getNearRange();
void setNearRange(double nearRange);
double getFarRange();
void setFarRange(double farRange);
double getRefRange();
void setRefRange(double refRange);
double getCenterFreq();
void setCenterFreq(double freq);
double getFs();
void setFs(double samplingFreq);
double getPRF();
void setPRF(double PRF);
double getCenterAngle();
void setCenterAngle(double angle);
QString getLookSide();
void setLookSide(QString lookside);
public:// 多普勒参数
double TotalProcessedAzimuthBandWidth, DopplerParametersReferenceTime;
double d0, d1, d2, d3, d4;
double r0, r1, r2, r3, r4;
double DEM;
double incidenceAngleNearRange, incidenceAngleFarRange;
public:
double getTotalProcessedAzimuthBandWidth();
void setTotalProcessedAzimuthBandWidth(double v);
double getDopplerParametersReferenceTime();
void setDopplerParametersReferenceTime(double v);
QVector<double> getDopplerParams();
void setDopplerParams(double d0, double d1, double d2, double d3, double d4);
QVector<double> getDopplerCenterCoff();
void setDopplerCenterCoff(double r0, double r1, double r2, double r3, double r4);
double getIncidenceAngleNearRange();
void setIncidenceAngleNearRangeet(double angle);
double getIncidenceAngleFarRange();
void setIncidenceAngleFarRange(double angle);
private:
double latitude_center, longitude_center,
latitude_topLeft, longitude_topLeft,
latitude_topRight, longitude_topRight,
latitude_bottomLeft, longitude_bottomLeft,
latitude_bottomRight, longitude_bottomRight;
public:
// Getter and Setter functions
double getLatitudeCenter();
void setLatitudeCenter(double value);
double getLongitudeCenter();
void setLongitudeCenter(double value);
double getLatitudeTopLeft();
void setLatitudeTopLeft(double value);
double getLongitudeTopLeft();
void setLongitudeTopLeft(double value);
double getLatitudeTopRight();
void setLatitudeTopRight(double value);
double getLongitudeTopRight();
void setLongitudeTopRight(double value);
double getLatitudeBottomLeft();
void setLatitudeBottomLeft(double value);
double getLongitudeBottomLeft();
void setLongitudeBottomLeft(double value);
double getLatitudeBottomRight();
void setLatitudeBottomRight(double value);
double getLongitudeBottomRight();
void setLongitudeBottomRight(double value);
public:
DemBox getExtend();
};

1
BaseTool/stdafx.cpp Normal file
View File

@ -0,0 +1 @@
#include "stdafx.h"

0
BaseTool/stdafx.h Normal file
View File