KJ135
/
SIMOrthoProgram-Orth_LT1AB-Strip
mirror of http://172.16.0.12:5000/WBJY/SIMOrthoProgram-Orth_LT1AB-Strip.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

LT参数适配

LT1ABOrth
剑古敛锋 2023-11-17 20:44:00 +08:00
parent 6290ac10dd
commit fef48a75b1
2 changed files with 248 additions and 210 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

59
simorthoprogram-orth_L_sar-strip/SIMOrthoProgram.cpp
View File

@ -34,17 +34,17 @@ void PreProcess(int argc, char* argv[])
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
std::cout << "预处理计算结果，可以计算出DEM 范围 " << endl; std::cout << "预处理计算结果，可以计算出DEM 范围 " << endl;
std::cout << "SIMOrthoProgram.exe 1 in_parameter_path in_dem_path in_ori_sar_path in_work_path in_taget_path "; std::cout << "SIMOrthoProgram.exe 1 in_parameter_path in_dem_path in_ori_sar_path in_work_path in_taget_path ";
std::string parameter_path = argv[2]; // 参数文件 //std::string parameter_path = argv[2]; // 参数文件
std::string dem_path = argv[3]; // dem 文件 //std::string dem_path = argv[3]; // dem 文件
std::string in_sar_path = argv[4]; // 输入SAR文件 //std::string in_sar_path = argv[4]; // 输入SAR文件
std::string work_path = argv[5]; // 目标空间文件 //std::string work_path = argv[5]; // 目标空间文件
std::string taget_path = argv[6]; // 输出坐标映射文件 //std::string taget_path = argv[6]; // 输出坐标映射文件
//std::string parameter_path = "D:\\micro\\LWork\\Ortho\\Temporary\\package\\orth_para.txt"; // 参数文件 std::string parameter_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\package\\orth_para.txt"; // 参数文件
//std::string dem_path = "D:\\micro\\LWork\\Ortho\\Temporary\\TestDEM\\mergedDEM.tif"; // dem 文件 std::string dem_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\TestDEM\\mergedDEM.tif"; // dem 文件
//std::string in_sar_path = "D:\\micro\\LWork\\Ortho\\Temporary\\unpack\\LT1B_MONO_MYC_STRIP4_005860_E130.9_N47.7_20230327_SLC_AHV_L1A_0000086966\\LT1B_MONO_MYC_STRIP4_005860_E130.9_N47.7_20230327_SLC_HH_L1A_0000086966.tiff"; // 输入SAR文件 std::string in_sar_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\unpack\\LT1B_MONO_KSC_STRIP2_007065_E100.8_N26.6_20230616_SLC_HH_L1A_0000135014\\LT1B_MONO_KSC_STRIP2_007065_E100.8_N26.6_20230616_SLC_HH_L1A_0000135014.tiff"; // 输入SAR文件
//std::string work_path = "D:\\micro\\LWork\\Ortho\\Temporary"; // 目标空间文件 std::string work_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary"; // 目标空间文件
//std::string taget_path = "D:\\micro\\LWork\\Ortho\\Temporary\\package"; // 输出坐标映射文件 std::string taget_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\package"; // 输出坐标映射文件
//std::string Incident_path = argv[7];// 输出入射角文件 //std::string Incident_path = argv[7];// 输出入射角文件
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
@ -58,6 +58,23 @@ void PreProcess(int argc, char* argv[])
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
process.InitSimulationSAR(parameter_path, work_path, taget_path, dem_path, in_sar_path); process.InitSimulationSAR(parameter_path, work_path, taget_path, dem_path, in_sar_path);
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
// ----- 测试最终插值
//std::string parameter_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\package\\orth_para.txt";
std::string in_rc_wgs84_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\package\\RD_sim_ori.tif";
std::string in_ori_sar_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\in_sar_power.tiff";
std::string out_orth_sar_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\in_sar_power_GTC.tiff";
//simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
process.pstn = pstn;
std::cout << "==========================================================================" << endl;
process.interpolation_GTC_sar_sigma(in_rc_wgs84_path, in_ori_sar_path, out_orth_sar_path, pstn);
} }
//mode 2 //mode 2
@ -178,14 +195,14 @@ void calInterpolation_cubic_Wgs84_rc_sar_sigma(int argc, char* argv[]) {
std::cout << "mode 9: interpolation(cubic convolution) orth sar value by rc_wgs84 and ori_sar image and model:\n "; std::cout << "mode 9: interpolation(cubic convolution) orth sar value by rc_wgs84 and ori_sar image and model:\n ";
std::cout << "SIMOrthoProgram.exe 9 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path"; std::cout << "SIMOrthoProgram.exe 9 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path";
std::string parameter_path = "D:\\micro\\LWork\\Ortho\\Temporary\\package\\orth_para.txt"; std::string parameter_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\package\\orth_para.txt";
std::string in_rc_wgs84_path = "D:\\micro\\LWork\\Ortho\\Temporary\\package\\RD_sim_ori.tif"; std::string in_rc_wgs84_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\package\\RD_sim_ori.tif";
std::string in_ori_sar_path = "D:\\micro\\LWork\\Ortho\\Temporary\\in_sar_power.tiff"; std::string in_ori_sar_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\in_sar_power.tiff";
std::string out_orth_sar_path = "D:\\micro\\LWork\\Ortho\\Temporary\\in_sar_power_GTC.tiff"; std::string out_orth_sar_path = "D:\\micro\\WorkSpace\\Ortho\\Temporary\\in_sar_power_GTC.tiff";
parameter_path = argv[2]; //parameter_path = argv[2];
in_rc_wgs84_path = argv[3]; //in_rc_wgs84_path = argv[3];
in_ori_sar_path = argv[4]; //in_ori_sar_path = argv[4];
out_orth_sar_path = argv[5]; //out_orth_sar_path = argv[5];
simProcess process; simProcess process;
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
@ -339,7 +356,7 @@ int main(int argc, char* argv[])
Landpoint p2 = { -3421843,5089485,3630606 }; Landpoint p1 = { -2609414,4763328,3332879 }; Landpoint p2 = { -3421843,5089485,3630606 }; Landpoint p1 = { -2609414,4763328,3332879 };
cout << getIncAngle(p2, p1) << endl;; cout << getIncAngle(p2, p1) << endl;;
std::cout << "program start:\t" << getCurrentTimeString() << endl;; std::cout << "program start:\t" << getCurrentTimeString() << endl;;
int mode = 9; int mode = 1;
GDALAllRegister(); GDALAllRegister();
if (argc == 0) { // 测试参数 if (argc == 0) { // 测试参数
// 算法说明 // 算法说明
@ -370,12 +387,12 @@ int main(int argc, char* argv[])
} }
//calInterpolation_cubic_Wgs84_rc_sar(argc, argv); //calInterpolation_cubic_Wgs84_rc_sar(argc, argv);
} }
else if (argc > 1) { // 预处理模块 else if (argc >= 1) { // 预处理模块
std::cout << "=============================description V2.0 =============================================" << endl; std::cout << "=============================description V2.0 =============================================" << endl;
std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl; std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl;
std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl; std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl;
mode = stoi(argv[1]); //mode = stoi(argv[1]);
if (mode == 0) { if (mode == 0) {
test_main(mode, argv[2]); test_main(mode, argv[2]);
} }

47
simorthoprogram-orth_L_sar-strip/simptsn.cpp
View File

@ -203,9 +203,13 @@ Eigen::MatrixX<double> PSTNAlgorithm::calNumericalDopplerValue(Eigen::MatrixX<do
double PSTNAlgorithm::calNumericalDopplerValue(double R) double PSTNAlgorithm::calNumericalDopplerValue(double R)
{ {
//double t = (R - this->refrange) * (1e6 / LIGHTSPEED); GF3 //double t = (R - this->refrange) * (1e6 / LIGHTSPEED); GF3
double t = R / LIGHTSPEED - this->doppler_reference_time; //S_SAR¹«Ê½ //double t = R / LIGHTSPEED - this->doppler_reference_time; //S_SAR公式
double t = R * 2 / LIGHTSPEED;// -this->doppler_reference_time; //LT1 公式1
//double t = R; //LT1 公式3
//double t = c/RSR; //LT1 公式4
//double temp = this->doppler_paras(0) + this->doppler_paras(1) * t + this->doppler_paras(2) * t * t + this->doppler_paras(3) * t * t * t + this->doppler_paras(4) * t * t * t * t; // GF3 HJ //double temp = this->doppler_paras(0) + this->doppler_paras(1) * t + this->doppler_paras(2) * t * t + this->doppler_paras(3) * t * t * t + this->doppler_paras(4) * t * t * t * t; // GF3 HJ
double temp = this->doppler_paras(0) + this->doppler_paras(1) * t + this->doppler_paras(2) * t * t + this->doppler_paras(3) * t * t * t; // LT double temp = 0; this->doppler_paras(0) + this->doppler_paras(1) * t + this->doppler_paras(2) * t * t + this->doppler_paras(3) * t * t * t; // LT
return temp; return temp;
} }
@ -270,7 +274,8 @@ Eigen::MatrixX<double> PSTNAlgorithm::PSTN(Eigen::MatrixX<double> landpoints, do
break; break;
} }
} }
SARPoint(j, 0) = ti;
SARPoint(j, 0) = FdTheory2;
SARPoint(j, 1) = (ti - this->imgStartTime) * this->PRF; SARPoint(j, 1) = (ti - this->imgStartTime) * this->PRF;
// 尝试使用入射角切分列号 // 尝试使用入射角切分列号
@ -278,9 +283,6 @@ Eigen::MatrixX<double> PSTNAlgorithm::PSTN(Eigen::MatrixX<double> landpoints, do
satepoint = satestate(Eigen::all, { 0,1,2 }); satepoint = satestate(Eigen::all, { 0,1,2 });
sate_land = satepoint - landpoint; sate_land = satepoint - landpoint;
R1 = (sate_land.array().pow(2).array().rowwise().sum().array().sqrt())[0];; R1 = (sate_land.array().pow(2).array().rowwise().sum().array().sqrt())[0];;
//SARPoint(j, 2) = getIncAngle(satePoints, landPoint);
//SARPoint(j, 2) = (R1 - this->R0 - (this->refrange - this->R0) / 3) / this->widthspace;
SARPoint(j, 2) = getRangeColumn(R1, this->R0, this->fs, this->lamda);// (R1 - this->R0) / this->widthspace; // SARPoint(j, 2) = getRangeColumn(R1, this->R0, this->fs, this->lamda);// (R1 - this->R0) / this->widthspace; //
} }
//std::cout <<"SARPoints(200, 0):\t" << SARPoint(200, 0) << "\t" << SARPoint(200, 1) << "\t" << SARPoint(200, 2) << "\t" << endl; //std::cout <<"SARPoints(200, 0):\t" << SARPoint(200, 0) << "\t" << SARPoint(200, 1) << "\t" << SARPoint(200, 2) << "\t" << endl;
@ -402,7 +404,7 @@ int simProcess::InitSimulationSAR(std::string paras_path, std::string workspace_
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
this->CreateSARDEM(); this->CreateSARDEM();
this->dem2SAR_row(); // 获取行号 this->dem2SAR_row(); // 获取行号
this->SARIncidentAngle(); //this->SARIncidentAngle();
//this->SARSimulationWGS(); //this->SARSimulationWGS();
//this->SARSimulation(); //this->SARSimulation();
this->in_sar_power(); this->in_sar_power();
@ -420,6 +422,11 @@ int simProcess::InitSimulationSAR(std::string paras_path, std::string workspace_
////插值计算行列号 ////插值计算行列号
//this->ReflectTable_WGS2Range(); //this->ReflectTable_WGS2Range();
//this->SARIncidentAngle_RPC(); //this->SARIncidentAngle_RPC();
return 0; return 0;
} }
int simProcess::InitASFSAR(std::string paras_path, std::string workspace_path, std::string out_dir_path, std::string in_dem_path) int simProcess::InitASFSAR(std::string paras_path, std::string workspace_path, std::string out_dir_path, std::string in_dem_path)
@ -589,12 +596,24 @@ int simProcess::CreateSARDEM()
std::cout << "min_lon-max_lon:\t" << box.min_lon << "\t" << box.max_lon << endl; std::cout << "min_lon-max_lon:\t" << box.min_lon << "\t" << box.max_lon << endl;
std::cout << "cal box of sar over" << endl; std::cout << "cal box of sar over" << endl;
// 计算采样之后的影像大小 // 计算采样之后的影像大小
//int width = 1.8 * this->pstn.width; // 原采样参数
//int height = 1.8 * this->pstn.height;
int width = 1.8 * this->pstn.width; int width = 1.8 * this->pstn.width;
int height = 1.8 * this->pstn.height; int height = 1.8 * this->pstn.height;
width = this->pstn.width; // 为了测试
height = this->pstn.height;
double heightspace = (box.max_lat - box.min_lat) / height; double heightspace = (box.max_lat - box.min_lat) / height;
double widthspace = (box.max_lon - box.min_lon) / width; double widthspace = (box.max_lon - box.min_lon) / width;
double space = heightspace > widthspace ? heightspace : widthspace;
std::cout << "resample dem:\n" << getCurrentTimeString() << std::endl; std::cout << "resample dem:\n" << getCurrentTimeString() << std::endl;
std::cout << "in_dem:\t" << this->in_dem_path << endl; std::cout << "in_dem:\t" << this->in_dem_path << endl;
std::cout << "out_dem:\t" << this->dem_path << endl; std::cout << "out_dem:\t" << this->dem_path << endl;
@ -2686,7 +2705,7 @@ void simProcess::interpolation_GTC_sar_sigma(std::string in_rc_wgs84_path, std::
// 正式计算插值 // 正式计算插值
{ {
int line_invert = 600; int line_invert = 2000;
line_invert = line_invert > 10 ? line_invert : 10; line_invert = line_invert > 10 ? line_invert : 10;
int start_ids = 0; int start_ids = 0;
do { do {
@ -2694,6 +2713,7 @@ void simProcess::interpolation_GTC_sar_sigma(std::string in_rc_wgs84_path, std::
std::cout << "rows:\t" << start_ids << "/" << sar_rc.height << "\t computing.....\t" << getCurrentTimeString() << endl; std::cout << "rows:\t" << start_ids << "/" << sar_rc.height << "\t computing.....\t" << getCurrentTimeString() << endl;
Eigen::MatrixXd sar_r = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 1); Eigen::MatrixXd sar_r = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 1);
Eigen::MatrixXd sar_c = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 2); Eigen::MatrixXd sar_c = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 2);
sar_r = sar_r.array();
// 处理行列号获取范围 // 处理行列号获取范围
int min_r = floor(sar_r.array().minCoeff()); int min_r = floor(sar_r.array().minCoeff());
@ -2701,7 +2721,7 @@ void simProcess::interpolation_GTC_sar_sigma(std::string in_rc_wgs84_path, std::
int min_c = floor(sar_c.array().minCoeff()); int min_c = floor(sar_c.array().minCoeff());
int max_c = ceil(sar_c.array().maxCoeff()); int max_c = ceil(sar_c.array().maxCoeff());
if (max_r < 0 || max_c < 0 || min_r >= this->pstn.height || min_c >= this->pstn.width) { if (max_r < 0 || max_c < 0||min_r>max_r||min_c>max_c || min_r >= this->pstn.height || min_c >= this->pstn.width) {
start_ids = start_ids + line_invert; start_ids = start_ids + line_invert;
if (start_ids < sar_rc.height) { if (start_ids < sar_rc.height) {
continue; continue;
@ -2735,9 +2755,9 @@ void simProcess::interpolation_GTC_sar_sigma(std::string in_rc_wgs84_path, std::
//Eigen::MatrixXd sar_b = sar_img.getData(start_ids, 0, line_invert, sar_rc.width, 2); //Eigen::MatrixXd sar_b = sar_img.getData(start_ids, 0, line_invert, sar_rc.width, 2);
sar_r = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 1); sar_r = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 1);
sar_c = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 2); sar_c = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 2);
row_count = sar_r.rows()-100; row_count = sar_r.rows();
col_count = sar_r.cols(); col_count = sar_r.cols();
#pragma omp parallel for num_threads(8) // NEW ADD //#pragma omp parallel for num_threads(8) // NEW ADD
for (int i = 0; i < row_count; i++) { for (int i = 0; i < row_count; i++) {
int r0, r1, r2, r3, c0, c1, c2, c3; int r0, r1, r2, r3, c0, c1, c2, c3;
double r, c; double r, c;
@ -2759,6 +2779,7 @@ void simProcess::interpolation_GTC_sar_sigma(std::string in_rc_wgs84_path, std::
if (r0<-1 || c0<-1 || c3>this->pstn.width || r3>this->pstn.height) { if (r0<-1 || c0<-1 || c3>this->pstn.width || r3>this->pstn.height) {
continue; continue;
} }
// 边界插值计算，插值模块权重 // 边界插值计算，插值模块权重
Eigen::MatrixX<complex<double>> img_block(4, 4); Eigen::MatrixX<complex<double>> img_block(4, 4);
if (r0 == -1 || c0 == -1) { if (r0 == -1 || c0 == -1) {
@ -2767,7 +2788,7 @@ void simProcess::interpolation_GTC_sar_sigma(std::string in_rc_wgs84_path, std::
img_block(3, 3) = ori(r3 - min_r, c3); img_block(3, 3) = ori(r3 - min_r, c3);
img_block(3, 2) = ori(r3 - min_r, c2); img_block(3, 2) = ori(r3 - min_r, c2);
img_block(3, 1) = ori(r3 - min_r, c1); img_block(3, 1) = ori(r3 - min_r, c1);
img_block(3, 0) = ori(r3 - min_r, c1) * complex<double>{3, 0} - complex<double>{3,0} *ori(r3 - min_r, c2) + ori(r3 - min_r, c3); // k,-1 img_block(3, 0) = ori(r3 - min_r, c1) * complex<double>{3, 0} - complex<double>{3, 0} *ori(r3 - min_r, c2) + ori(r3 - min_r, c3); // k,-1
img_block(2, 3) = ori(r2 - min_r, c3); img_block(2, 3) = ori(r2 - min_r, c3);
img_block(2, 2) = ori(r2 - min_r, c2); img_block(2, 2) = ori(r2 - min_r, c2);
@ -3414,7 +3435,7 @@ double getRangeColumn(long double R, long double NearRange, long double Fs, long
Eigen::MatrixXd getRangeColumn(Eigen::MatrixXd R, long double NearRange, long double Fs, long double lamda) Eigen::MatrixXd getRangeColumn(Eigen::MatrixXd R, long double NearRange, long double Fs, long double lamda)
{ {
return ((R.array() - NearRange).array() / lamda * Fs).array().cast<double>(); return (2*(R.array() - NearRange).array() / lamda * Fs).array().cast<double>();
} }
double getRangebyColumn(double j, long double NearRange, long double Fs, long double lamda) double getRangebyColumn(double j, long double NearRange, long double Fs, long double lamda)