From 78d58e678e7fc465a0e8a7c440334af3c1a3f9db Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=E5=89=91=E5=8F=A4=E6=95=9B=E9=94=8B?= <3045316072@qq.com>
Date: Fri, 1 Sep 2023 12:45:38 +0800
Subject: [PATCH] init V2.2

---
 SIMOrthoProgram.cpp |  6 ++---
 baseTool.cpp        | 26 +-------------------
 baseTool.h          |  2 ++
 simptsn.cpp         | 60 +++++++++++++++++++++++++++++++++++++++------
 simptsn.h           | 16 ++++++------
 5 files changed, 67 insertions(+), 43 deletions(-)
diff --git a/SIMOrthoProgram.cpp b/SIMOrthoProgram.cpp
index 43c78a1..df00f2d 100644
--- a/SIMOrthoProgram.cpp
+++ b/SIMOrthoProgram.cpp
@@ -204,9 +204,9 @@ void createRPC_lon_lat(int argc, char* argv[]) {
 void Scatter2Grid_lon_lat(int argc, char* argv[]) {
 	std::cout << "mode 10";
 	std::cout << "SIMOrthoProgram.exe 10 lon_lat_path  data_tiff  grid_path  space";
-	std::string lon_lat_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\RD_ori_sim.tif";
-	std::string data_tiff = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1B_HH_L10003923848_db.tif";
-	std::string grid_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1B_HH_L10003923848__db.tif";
+	std::string lon_lat_path = "F:\\orthtest\\ori_sim_preprocessed.tif";
+	std::string data_tiff = "F:\\orthtest\\SoilMoistureProduct_geo.tif";
+	std::string grid_path = "F:\\orthtest\\SoilMoistureProduct_geo_test.tif";
 	double space = 5;
 
 	lon_lat_path = argv[2]; 
diff --git a/baseTool.cpp b/baseTool.cpp
index aea74d8..6b89b76 100644
--- a/baseTool.cpp
+++ b/baseTool.cpp
@@ -762,20 +762,6 @@ void clipGdalImage(string in_path, string out_path, DemBox box, double pixelinte
 }
 
 
-
-/***
-* 遥锟斤拷影锟斤拷锟截诧拷锟斤拷   (要锟斤拷影锟斤拷锟斤拷锟斤拷锟酵队帮拷锟斤拷锟斤拷锟斤拷卟锟酵?1锟?7)
-* @param pszSrcFile        锟斤拷锟斤拷锟侥硷拷锟斤拷路锟斤拷
-* @param pszOutFile        写锟斤拷慕锟斤拷图锟斤拷锟铰凤拷锟?1锟?7
-* @param eResample         锟斤拷锟斤拷模式锟斤拷锟斤拷锟斤拷锟街ｏ拷锟斤拷锟斤拷渭锟紾DALResampleAlg锟斤拷锟藉，默锟斤拷为双锟斤拷锟斤拷锟节诧拷
-* @param gt             X转锟斤拷锟斤拷锟斤拷锟饺ｏ拷默锟较达拷小为1.0锟斤拷锟斤拷锟斤拷1图锟斤拷锟斤拷小锟斤拷1锟斤拷示图锟斤拷锟斤拷小锟斤拷锟斤拷值锟较碉拷锟节诧拷锟斤拷锟斤拷图锟斤拷目锟斤拷群筒锟斤拷锟角巴硷拷锟斤拷锟饺的憋拷
-* @param new_width
-* @param new_hegiht
-* @retrieve     0   锟缴癸拷
-* @retrieve     -1  锟斤拷源锟侥硷拷失锟斤拷
-* @retrieve     -2  锟斤拷锟斤拷锟斤拷锟侥硷拷失锟斤拷
-* @retrieve     -3  锟斤拷锟斤拷锟斤拷锟斤拷锟叫筹拷锟斤拷
-*/
 int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample)
 {
 	GDALAllRegister();
@@ -941,11 +927,6 @@ int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResa
 }
 
 
-/// <summary>
-/// 锟斤拷锟斤拷转锟街凤拷锟斤拷
-/// </summary>
-/// <param name="Num"></param>
-/// <returns></returns>
 string Convert(float Num)
 {
 	ostringstream oss;
@@ -955,12 +936,7 @@ string Convert(float Num)
 }
 
 
-/// <summary>
-/// 路锟斤拷拼锟斤拷
-/// </summary>
-/// <param name="path"></param>
-/// <param name="filename"></param>
-/// <returns></returns>
+
 std::string JoinPath(const std::string& path, const std::string& filename)
 {
 	int dir_size = path.size();
diff --git a/baseTool.h b/baseTool.h
index a0639a9..64e53cf 100644
--- a/baseTool.h
+++ b/baseTool.h
@@ -29,6 +29,8 @@ using namespace Eigen;
 
 #define PI_180 180/3.141592653589793238462643383279
 #define T180_PI 3.141592653589793238462643383279/180
+#define LIGHTSPEED 299792458
+
 
 #define Radians2Degrees(Radians) Radians*PI_180
 #define Degrees2Radians(Degrees) Degrees*T180_PI
diff --git a/simptsn.cpp b/simptsn.cpp
index a6cb6c5..4e49d39 100644
--- a/simptsn.cpp
+++ b/simptsn.cpp
@@ -43,6 +43,28 @@ PSTNAlgorithm::PSTNAlgorithm(string infile_path)
 	std::cout << "=========================================================================" << endl;
 	std::cout << "read  parameters :\t" << infile_path << endl;
 	std::cout << "parameters files :\t" << infile_path << endl;
+	std::cout << "params file Format:" << endl;
+	std::cout << "height --- image height" << endl;
+	std::cout << "width --- image width" << endl;
+	std::cout << "R0  --- image nearRange" << endl;
+	std::cout << "near_in_angle " << endl;
+	std::cout << "far incident angle " << endl;
+	std::cout << "LightSpeed " << endl;
+	std::cout << "lamda" << endl;
+	std::cout << "imgStartTime" << endl;
+	std::cout << "PRF" << endl;
+	std::cout << "refrange" << endl;
+	std::cout << "fs" << endl;
+	std::cout << "doppler_paramenter_number " << endl;
+	std::cout << "doppler_paramenter 1 " << endl;
+	std::cout << "doppler_paramenter 2" << endl;
+	std::cout << "doppler_paramenter 3 " << endl;
+	std::cout << "doppler_paramenter 4 " << endl;
+	std::cout << "doppler_paramenter 5 " << endl;
+	std::cout << "ispolySatelliteModel" << endl;
+	std::cout << "sate polynum " << endl;
+	std::cout << "SatelliteModelStartTime" << endl;
+	std::cout << "orbit params list " << endl;
 	std::cout << "============================================================================" << endl;
 
 	// 解析文件
@@ -64,9 +86,9 @@ PSTNAlgorithm::PSTNAlgorithm(string infile_path)
 		getline(infile, buf); this->LightSpeed = stod(buf);		std::cout << "LightSpeed\t" << this->LightSpeed << endl;
 		getline(infile, buf); this->lamda = stod(buf);			std::cout << "lamda\t" << this->lamda << endl;
 		getline(infile, buf); this->imgStartTime = stod(buf);	std::cout << "imgStartTime\t" << this->imgStartTime << endl;
-		getline(infile, buf); this->PRF = stod(buf);				std::cout << "PRF\t" << this->PRF << endl;
+		getline(infile, buf); this->PRF = stod(buf);			std::cout << "PRF\t" << this->PRF << endl;
 		getline(infile, buf); this->refrange = stod(buf);		std::cout << "refrange\t" << this->refrange << endl;
-		getline(infile, buf); this->widthspace = stod(buf);		std::cout << "widthspace\t" << this->widthspace << endl;
+		getline(infile, buf); this->fs = stod(buf);				std::cout << "fs\t" << this->fs << endl;
 
 		getline(infile, buf); this->doppler_paramenter_number = stoi(buf); std::cout << "doppler_paramenter_number\t" << this->doppler_paramenter_number << endl;
 
@@ -96,6 +118,7 @@ PSTNAlgorithm::PSTNAlgorithm(string infile_path)
 				for (int j = 0; j < polynum; j++) {
 					getline(infile, buf);
 					polySatellitePara(j, i) = stod(buf);
+					std::cout << "orbit params " << i << " , " << j << buf << endl;
 				}
 			}
 		}
@@ -165,7 +188,7 @@ int PSTNAlgorithm::dem2SAR_RC(std::string dem_path, std::string sim_rc_path)
 
 Eigen::MatrixX<double> PSTNAlgorithm::calNumericalDopplerValue(Eigen::MatrixX<double> R)
 {
-	Eigen::MatrixX<double> t = (R.array() - this->refrange) * (1000000 / this->LightSpeed);
+	Eigen::MatrixX<double> t = (R.array() - this->refrange) * (1000000 / LIGHTSPEED);
 	Eigen::MatrixX<double> tlist(t.rows(), this->doppler_paramenter_number);// nx5
 	tlist.col(0) = Eigen::MatrixX<double>::Ones(t.rows(), 1).col(0);
 	for (int i = 1; i < this->doppler_paramenter_number; i++)
@@ -180,7 +203,7 @@ Eigen::MatrixX<double> PSTNAlgorithm::calNumericalDopplerValue(Eigen::MatrixX<do
 
 double PSTNAlgorithm::calNumericalDopplerValue(double R)
 {
-	double t = (R - this->refrange) * (1000000 / this->LightSpeed);
+	double t = (R - this->refrange) * (1e6 / LIGHTSPEED);
 	return 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;
 }
 
@@ -198,7 +221,6 @@ Eigen::MatrixX<double> PSTNAlgorithm::PSTN(Eigen::MatrixX<double> landpoints, do
 	double S_T = 0;
 	double R_T = 0;
 	double cos_theta = 0;
-	double delta_col_ref = ((this->refrange - this->R0) * 100000000 / this->LightSpeed) / this->widthspace;
 
 	Eigen::MatrixX<double> landpoint = Eigen::MatrixX<double>::Ones(1, 6);// 卫星状态
 	Eigen::MatrixX<double> satestate = Eigen::MatrixX<double>::Ones(1, 6);// 卫星状态
@@ -256,7 +278,7 @@ Eigen::MatrixX<double> PSTNAlgorithm::PSTN(Eigen::MatrixX<double> landpoints, do
 		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) = (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;
 	return SARPoint;
@@ -719,7 +741,7 @@ int simProcess::dem2SAR()
 		for (int max_rows_ids = 0; max_rows_ids < dem_r.height; max_rows_ids = max_rows_ids + line_invert) {
 			Eigen::MatrixXd dem_r_block = dem_r.getData(max_rows_ids, 0, line_invert, dem_r.width, 1);  // 行号
 			Eigen::MatrixXd slant_range = plant_slant_range.getData(max_rows_ids, 0, line_invert, dem_r.width, 1); // 地形
-			Eigen::MatrixXd dem_c = (slant_range.array() - this->pstn.R0) / this->pstn.widthspace;
+			Eigen::MatrixXd dem_c = getRangeColumn(slant_range, this->pstn.R0, this->pstn.fs, this->pstn.lamda);//  (slant_range.array() - this->pstn.R0) / this->pstn.widthspace;
 			omp_set_lock(&lock); //获得互斥器
 			dem_rc.saveImage(dem_r_block, max_rows_ids, 0, 1);
 			dem_rc.saveImage(dem_c, max_rows_ids, 0, 2);
@@ -2113,7 +2135,7 @@ int simProcess::ASF(std::string in_dem_path, std::string out_latlon_path, ASFOrt
 				double beta0 = 50 * d2r;// (180-getAngle(sata, sata-initp))* d2r;
 
 				for (int j = 0; j < this->pstn.width; j++) {
-					double R = this->pstn.R0 + j * this->pstn.widthspace;
+					double R = getRangebyColumn(j, this->pstn.R0, this->pstn.fs, this->pstn.lamda);
 					double ati_H = sar_dem(i - row_min, j);
 					if (!isnan(ati_H) && !isinf(ati_H)) {
 						ati_H = aveageAti;
@@ -3140,6 +3162,28 @@ void simProcess::CorrectionFromSLC2Geo(string in_lon_lat_path, string in_radar_p
 /// <returns></returns>
 
 
+double getRangeColumn(long double R, long double NearRange, long double Fs, long double lamda)
+{
+	double result = (double)(2 * (R - NearRange) / lamda * Fs);
+	return result;
+}
+
+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>();
+}
+
+double getRangebyColumn(double j, long double NearRange, long double Fs, long double lamda)
+{
+	return (j *  lamda / 2) /  Fs + NearRange;
+ 
+}
+
+Eigen::MatrixXd getRangebyColumn(Eigen::MatrixXd j, long double NearRange, long double Fs, long double lamda)
+{
+	return ((j.array() * lamda / 2) / Fs + NearRange).array().cast<double>();
+}
+
 double getlocalIncAngle(Landpoint& satepoint, Landpoint& landpoint, Landpoint& slopeVector) {
 	Landpoint Rsc = satepoint - landpoint; // AB=B-A
 	//double R = getlength(Rsc);
diff --git a/simptsn.h b/simptsn.h
index d2a639e..fd69252 100644
--- a/simptsn.h
+++ b/simptsn.h
@@ -18,6 +18,13 @@
 using namespace std;
 using namespace Eigen;
 
+////////////// 常用函数 ////////////////////////
+
+double getRangeColumn(long double 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);
+
+double getRangebyColumn(double j, long double NearRange, long double Fs, long double lamda);
+Eigen::MatrixXd getRangebyColumn(Eigen::MatrixXd j, long double NearRange, long double Fs, long double lamda);
 /////////////  ptsn 算法 /////////////////////
 
 class PSTNAlgorithm {
@@ -29,7 +36,7 @@ public:
 
 	int dem2lat_lon(std::string dem_path, std::string lon_path, std::string lat_path);
 	int dem2SAR_RC(std::string dem_path, std::string sim_rc_path);
-
+	
 	// 模拟算法函数
 	Eigen::MatrixX<double> calNumericalDopplerValue(Eigen::MatrixX<double> R);// 数值模拟法计算多普勒频移值
 	double calNumericalDopplerValue(double R);// 数值模拟法计算多普勒频移值
@@ -45,11 +52,6 @@ public: // WGS84 
 	double Dut1 = 0.0366742;
 	double Dat = 37;
 public:
-//	std::string dem_path; // 模拟用的DEM影像
-//	std::string lon_path;// 经度影像
-//	std::string lat_path;// 纬度影像
-
-//	std::string simsar_path;// 仿真影像坐标
 	int height; // 影像的高
 	int width;
 
@@ -58,7 +60,7 @@ public:
 	double far_in_angle;// 远入射角
 
 	// SAR的成像参数
-	double widthspace;// 距离向分辨率
+	double fs;// 距离向采样率
 	double R0;//近斜距
 	double LightSpeed; // 光速
 	double lamda;// 波长
