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

编码格式由utf-8 修改为 gb 2312

master
剑古敛锋 2023-05-05 15:27:20 +08:00
parent cb3db54958
commit 027c1484ee
16 changed files with 681 additions and 697 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
ImageMatch.cpp
View File

@ -2,7 +2,7 @@
//#define EIGEN_USE_MKL_ALL //#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2 //#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h> //#include <mkl.h>
// 本地方法 // 本地方法
#include <iostream> #include <iostream>
#include <Eigen/Core> #include <Eigen/Core>
@ -43,7 +43,7 @@ int ImageMatch::gdal2JPG(std::string gdal_path, std::string jpg_path, int band_i
Eigen::MatrixXd temp(line_invert, gdalimg.width); Eigen::MatrixXd temp(line_invert, gdalimg.width);
double min_value = 0, temp_min = 0; double min_value = 0, temp_min = 0;
double max_value = 0, temp_max = 0; double max_value = 0, temp_max = 0;
// 线性拉伸2% // 线性拉伸2%
Eigen::MatrixXd hist = gdalimg.getHist(band_ids); Eigen::MatrixXd hist = gdalimg.getHist(band_ids);
int count = 0; int count = 0;
@ -66,9 +66,9 @@ int ImageMatch::gdal2JPG(std::string gdal_path, std::string jpg_path, int band_i
count = count + hist(i, 1); count = count + hist(i, 1);
} }
// 重新缩放最大值，最小值 // 重新缩放最大值，最小值
std::cout << "min value ：\t" << min_value << "\n"; std::cout << "min value ：\t" << min_value << "\n";
std::cout << "max value ：\t" << max_value << "\n"; std::cout << "max value ：\t" << max_value << "\n";
start_ids = 0; start_ids = 0;
do { do {
line_invert = line_invert + start_ids < gdalimg.height ? line_invert : gdalimg.height - start_ids; line_invert = line_invert + start_ids < gdalimg.height ? line_invert : gdalimg.height - start_ids;
@ -89,8 +89,8 @@ int ImageMatch::gdal2JPG(std::string gdal_path, std::string jpg_path, int band_i
result1 = img; result1 = img;
//bilateralFilter(img, result1, 10, 80, 50); //bilateralFilter(img, result1, 10, 80, 50);
vector<int> compression_params; vector<int> compression_params;
compression_params.push_back(cv::IMWRITE_JPEG_QUALITY); //选择jpeg compression_params.push_back(cv::IMWRITE_JPEG_QUALITY); //选择jpeg
compression_params.push_back(100); //在这个填入你要的图片质量 compression_params.push_back(100); //在这个填入你要的图片质量
cv::Mat out = result1; cv::Mat out = result1;
//cv::resize(img, out, cv::Size(img.cols, img.rows), cv::INTER_AREA); //cv::resize(img, out, cv::Size(img.cols, img.rows), cv::INTER_AREA);
@ -100,9 +100,9 @@ int ImageMatch::gdal2JPG(std::string gdal_path, std::string jpg_path, int band_i
std::cout << "convert gdal to jpg , overing : \t" << getCurrentTimeString() << endl; std::cout << "convert gdal to jpg , overing : \t" << getCurrentTimeString() << endl;
std::cout << "=========================================\n"; std::cout << "=========================================\n";
std::cout << " convert gdal to jpg :\n"; std::cout << " convert gdal to jpg :\n";
std::cout << "input gdal img file path ：\t" << gdal_path << "\n"; std::cout << "input gdal img file path ：\t" << gdal_path << "\n";
std::cout << "input gdal img band ：\t" << band_ids << "\n"; std::cout << "input gdal img band ：\t" << band_ids << "\n";
std::cout << "out jpg file path ：\t" << jpg_path << "\n"; std::cout << "out jpg file path ：\t" << jpg_path << "\n";
std::cout << "=========================================\n"; std::cout << "=========================================\n";
return 0; return 0;
} }
@ -110,7 +110,7 @@ int ImageMatch::gdal2JPG(std::string gdal_path, std::string jpg_path, int band_i
/// <summary> /// <summary>
/// 获取模型匹配点 /// 获取模型匹配点
/// </summary> /// </summary>
/// <param name="ori_power_path"></param> /// <param name="ori_power_path"></param>
/// <param name="sim_sum_path"></param> /// <param name="sim_sum_path"></param>
@ -127,7 +127,7 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
std::cout << "input ori image path : \t" << ori_power_path << endl; std::cout << "input ori image path : \t" << ori_power_path << endl;
std::cout << "input sim image path : \t" << sim_sum_path << endl; std::cout << "input sim image path : \t" << sim_sum_path << endl;
//读取影像 //读取影像
cv::Mat ori = openJPG(ori_power_path); cv::Mat ori = openJPG(ori_power_path);
@ -169,13 +169,13 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
cv::meanStdDev(templet_mat, mean1, stddevMat); cv::meanStdDev(templet_mat, mean1, stddevMat);
double minvalue = 0; double minvalue = 0;
double maxvalue = 0; double maxvalue = 0;
cv::minMaxLoc(templet_mat, &minvalue, &maxvalue, NULL, NULL);//用于检测矩阵中最大值和最小值的位置 cv::minMaxLoc(templet_mat, &minvalue, &maxvalue, NULL, NULL);//用于检测矩阵中最大值和最小值的位置
double sig = (stddevMat.at<double>(0, 0)) / (maxvalue - minvalue); double sig = (stddevMat.at<double>(0, 0)) / (maxvalue - minvalue);
if (sig >1) { if (sig >1) {
//continue; //continue;
} }
// 构建搜索域 // 构建搜索域
int search_i = i - 2 * searchSize >= 0 ? i - 2 * searchSize >= 0 : 0; int search_i = i - 2 * searchSize >= 0 ? i - 2 * searchSize >= 0 : 0;
int search_j = j - 2 * searchSize >= 0 ? j - 2 * searchSize >= 0 : 0; int search_j = j - 2 * searchSize >= 0 ? j - 2 * searchSize >= 0 : 0;
int len_i = search_i + 4 * searchSize; int len_i = search_i + 4 * searchSize;
@ -183,14 +183,14 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
len_i = search_i + len_i < ori.rows ? len_i : ori.rows - search_i - 1; len_i = search_i + len_i < ori.rows ? len_i : ori.rows - search_i - 1;
len_j = search_j + len_j < ori.cols ? len_j : ori.rows - search_j - 1; len_j = search_j + len_j < ori.cols ? len_j : ori.rows - search_j - 1;
cv::Mat serch = ori(Rect(search_j, search_i, len_j, len_i)); cv::Mat serch = ori(Rect(search_j, search_i, len_j, len_i));
// 检索 // 检索
cv::Mat result; cv::Mat result;
matchTemplate(serch, templet_mat, result, cv::TM_CCOEFF_NORMED); matchTemplate(serch, templet_mat, result, cv::TM_CCOEFF_NORMED);
Point minLoc; Point minLoc;
Point maxLoc; Point maxLoc;
Point matchLoc; Point matchLoc;
double tempminVal = 100, tempmaxVal = 0; double tempminVal = 100, tempmaxVal = 0;
cv::minMaxLoc(result, &tempminVal, &tempmaxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置 cv::minMaxLoc(result, &tempminVal, &tempmaxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
double offset_j = maxLoc.x + search_j; double offset_j = maxLoc.x + search_j;
double offset_i = maxLoc.y + search_i; double offset_i = maxLoc.y + search_i;
@ -239,7 +239,7 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
/* /*
// step 1: 粗匹配，分块均匀匹配 // step 1: 粗匹配，分块均匀匹配
std::cout << "rough match , begining : \t" << getCurrentTimeString() << endl; std::cout << "rough match , begining : \t" << getCurrentTimeString() << endl;
double offset_x = 0, offset_y = 0; double offset_x = 0, offset_y = 0;
cv::Mat mask = sim; cv::Mat mask = sim;
@ -260,7 +260,7 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
double minVal = 100, maxVal = 0.7; double minVal = 100, maxVal = 0.7;
omp_lock_t lock; omp_lock_t lock;
std::cout << "ori_x\tori_y\tsim_x\tsim_y\tmaxval \t" << endl; std::cout << "ori_x\tori_y\tsim_x\tsim_y\tmaxval \t" << endl;
omp_init_lock(&lock); // 初始化互斥锁 omp_init_lock(&lock); // 初始化互斥锁
count = 0; count = 0;
int search_count = 0; int search_count = 0;
#pragma omp parallel for num_threads(8) #pragma omp parallel for num_threads(8)
@ -272,12 +272,12 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
templet_mat = ori(Rect(j - roughSize, i - roughSize, roughSize, roughSize)); templet_mat = ori(Rect(j - roughSize, i - roughSize, roughSize, roughSize));
matchTemplate(sim, templet_mat, result, cv::TM_CCOEFF_NORMED); matchTemplate(sim, templet_mat, result, cv::TM_CCOEFF_NORMED);
//normalize(result, result, 1, 0, cv::NORM_MINMAX); //normalize(result, result, 1, 0, cv::NORM_MINMAX);
// 通过函数 minMaxLoc 定位最匹配的位置; // 通过函数 minMaxLoc 定位最匹配的位置;
omp_set_lock(&lock); //获得互斥器 omp_set_lock(&lock); //获得互斥器
Point minLoc; Point minLoc;
Point maxLoc; Point maxLoc;
Point matchLoc; Point matchLoc;
cv::minMaxLoc(result, &tempminVal, &tempmaxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置 cv::minMaxLoc(result, &tempminVal, &tempmaxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
if (tempmaxVal >= maxVal) { if (tempmaxVal >= maxVal) {
offset_x = maxLoc.x - (j - roughSize); offset_x = maxLoc.x - (j - roughSize);
offset_y = maxLoc.y - (i - roughSize); offset_y = maxLoc.y - (i - roughSize);
@ -296,10 +296,10 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
} }
search_count = search_count + 1; search_count = search_count + 1;
std::cout << j - roughSize << "\t" << i - roughSize << "\t" << maxLoc.x << "\t" << maxLoc.y << "\t" << tempmaxVal << "\t" << search_count << "\t" << matchpoints.rows() << endl; std::cout << j - roughSize << "\t" << i - roughSize << "\t" << maxLoc.x << "\t" << maxLoc.y << "\t" << tempmaxVal << "\t" << search_count << "\t" << matchpoints.rows() << endl;
omp_unset_lock(&lock); //释放互斥器 omp_unset_lock(&lock); //释放互斥器
} }
} }
omp_destroy_lock(&lock); //销毁互斥器 omp_destroy_lock(&lock); //销毁互斥器
offset_x = offset_x*1.0 / offsetcount; offset_x = offset_x*1.0 / offsetcount;
offset_y = offset_y * 1.0 / offsetcount; offset_y = offset_y * 1.0 / offsetcount;
std::cout << "rough match point : "<< offsetcount <<"\n" << endl; std::cout << "rough match point : "<< offsetcount <<"\n" << endl;
@ -308,7 +308,7 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
std::cout << "maxVal : \t" << maxVal << endl; std::cout << "maxVal : \t" << maxVal << endl;
} }
std::cout << "rough match out : \t" << getCurrentTimeString() << endl; std::cout << "rough match out : \t" << getCurrentTimeString() << endl;
// step1.1: 粗匹配绘制结果 // step1.1: 粗匹配绘制结果
std::string rough_math_path = sim_sum_path; std::string rough_math_path = sim_sum_path;
boost::algorithm::replace_last(rough_math_path, ".", "_ori_in_sim."); boost::algorithm::replace_last(rough_math_path, ".", "_ori_in_sim.");
@ -317,13 +317,13 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
cv::rectangle(ori_in_sim, cv::Point(offset_x, offset_y), Point(offset_x + ori.cols, offset_y + ori.rows), Scalar(0, 255, 0), 2, 8, 0); cv::rectangle(ori_in_sim, cv::Point(offset_x, offset_y), Point(offset_x + ori.cols, offset_y + ori.rows), Scalar(0, 255, 0), 2, 8, 0);
vector<int> compression_params; vector<int> compression_params;
compression_params.push_back(cv::IMWRITE_JPEG_QUALITY); //选择jpeg compression_params.push_back(cv::IMWRITE_JPEG_QUALITY); //选择jpeg
compression_params.push_back(100); //在这个填入你要的图片质量 compression_params.push_back(100); //在这个填入你要的图片质量
cv::Mat out = ori_in_sim; cv::Mat out = ori_in_sim;
cv::imwrite(rough_math_path, out, compression_params); cv::imwrite(rough_math_path, out, compression_params);
if (offsetcount == 0) { // 表示全局粗匹配失败，无法进行精校准，考虑直接定位法 if (offsetcount == 0) { // 表示全局粗匹配失败，无法进行精校准，考虑直接定位法
std::cout << "there are not effective point in rought match \t" << endl; std::cout << "there are not effective point in rought match \t" << endl;
return Eigen::MatrixXd(0, 5); return Eigen::MatrixXd(0, 5);
} }
@ -341,7 +341,7 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
//std::cout << "rough match , overing : \t" << getCurrentTimeString() << endl; //std::cout << "rough match , overing : \t" << getCurrentTimeString() << endl;
//// step 2: 精匹配， //// step 2: 精匹配，
//std::cout << "Precise match , begining : \t" << getCurrentTimeString() << endl; //std::cout << "Precise match , begining : \t" << getCurrentTimeString() << endl;
//Eigen::MatrixXd matchpoints; //Eigen::MatrixXd matchpoints;
//{ //{
@ -366,8 +366,8 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
// row_count = row_count - PreciseSize; // row_count = row_count - PreciseSize;
// col_count = col_count - PreciseSize; // col_count = col_count - PreciseSize;
// omp_lock_t lock; // omp_lock_t lock;
// omp_init_lock(&lock); // 初始化互斥锁 // omp_init_lock(&lock); // 初始化互斥锁
// // 以搜索范围为核心计算 // // 以搜索范围为核心计算
// #pragma omp parallel for num_threads(8) // #pragma omp parallel for num_threads(8)
// for (int i = searchSize + offset_y; i < row_count; i = i + preciseStep) { // y // for (int i = searchSize + offset_y; i < row_count; i = i + preciseStep) { // y
// cv::Mat templet_mat; // cv::Mat templet_mat;
@ -377,13 +377,13 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
// cv::Mat result; // cv::Mat result;
// double minVal = 100, maxVal = 0; // double minVal = 100, maxVal = 0;
// for (int j = searchSize + offset_x; j < col_count; j = j + preciseStep) { //x // for (int j = searchSize + offset_x; j < col_count; j = j + preciseStep) { //x
// // 计算 起始点 // // 计算 起始点
// sim_start_x = j - searchSize; // sim_start_x = j - searchSize;
// sim_start_y = i - searchSize; // sim_start_y = i - searchSize;
// ori_start_x = (j - searchSize / 2) - PreciseSize / 2; // ori_start_x = (j - searchSize / 2) - PreciseSize / 2;
// ori_start_y = (i - searchSize / 2) - PreciseSize / 2; // ori_start_y = (i - searchSize / 2) - PreciseSize / 2;
// // 匹配模板，待匹配模板 // // 匹配模板，待匹配模板
// templet_mat = ori(Rect(ori_start_x, ori_start_y, PreciseSize, PreciseSize)); // templet_mat = ori(Rect(ori_start_x, ori_start_y, PreciseSize, PreciseSize));
// search_mat = ori(Rect(sim_start_x, sim_start_y, searchSize, searchSize)); // search_mat = ori(Rect(sim_start_x, sim_start_y, searchSize, searchSize));
@ -391,12 +391,12 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
// resample_search_mat = resampledMat(search_mat, search_size); // resample_search_mat = resampledMat(search_mat, search_size);
// matchTemplate(sim, templet_mat, result, cv::TM_CCORR_NORMED); // matchTemplate(sim, templet_mat, result, cv::TM_CCORR_NORMED);
// // 通过函数 minMaxLoc 定位最匹配的位置; // // 通过函数 minMaxLoc 定位最匹配的位置;
// cv::Point minLoc; cv::Point maxLoc; // cv::Point minLoc; cv::Point maxLoc;
// cv::Point matchLoc; // cv::Point matchLoc;
// cv::minMaxLoc(result, &minVal, &maxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置 // cv::minMaxLoc(result, &minVal, &maxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
// if (maxVal > 0.7) { // if (maxVal > 0.7) {
// omp_set_lock(&lock); //获得互斥器 // omp_set_lock(&lock); //获得互斥器
// //
// matchpoints(count, 0) = ori_start_x; //ori_x // matchpoints(count, 0) = ori_start_x; //ori_x
// matchpoints(count, 1) = ori_start_y; //ori_y // matchpoints(count, 1) = ori_start_y; //ori_y
@ -406,11 +406,11 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
// matchpoints(count, 4) = maxVal; // maxVal // matchpoints(count, 4) = maxVal; // maxVal
// count = count + 1; // count = count + 1;
// //
// omp_unset_lock(&lock); //释放互斥器 // omp_unset_lock(&lock); //释放互斥器
// } // }
// } // }
// } // }
// omp_destroy_lock(&lock); //销毁互斥器 // omp_destroy_lock(&lock); //销毁互斥器
//} //}
// //
@ -437,7 +437,7 @@ Eigen::MatrixXd ImageMatch::CreateMatchModel(Eigen::MatrixXd offsetXY_matrix)
// 0 1 2 3 4 // 0 1 2 3 4
Eigen::MatrixXd offset_x= offsetXY_matrix.col(2) - offsetXY_matrix.col(0); Eigen::MatrixXd offset_x= offsetXY_matrix.col(2) - offsetXY_matrix.col(0);
Eigen::MatrixXd offset_y = offsetXY_matrix.col(3) - offsetXY_matrix.col(1); Eigen::MatrixXd offset_y = offsetXY_matrix.col(3) - offsetXY_matrix.col(1);
// 计算最小二乘法模型 // 计算最小二乘法模型
Eigen::MatrixXd temp(offset_x.rows(), 6); Eigen::MatrixXd temp(offset_x.rows(), 6);
temp.col(0) = temp.col(0).array()*0+1; //1 temp.col(0) = temp.col(0).array()*0+1; //1
@ -488,16 +488,16 @@ int ImageMatch::outMatchModel(std::string matchmodel_path)
} }
/// <summary> /// <summary>
/// 读取jpg 文件 /// 读取jpg 文件
/// </summary> /// </summary>
/// <param name="jpg_path"></param> /// <param name="jpg_path"></param>
/// <returns></returns> /// <returns></returns>
cv::Mat openJPG(std::string jpg_path) cv::Mat openJPG(std::string jpg_path)
{ {
cv::Mat image = cv::imread(jpg_path); cv::Mat image = cv::imread(jpg_path);
if (image.data == nullptr) //nullptr是c++11新出现的空指针常量 if (image.data == nullptr) //nullptr是c++11新出现的空指针常量
{ {
throw new exception("图片文件不存在"); throw new exception("图片文件不存在");
} }
return image; return image;
} }

4
ImageMatch.h
View File

@ -2,7 +2,7 @@
//#define EIGEN_USE_MKL_ALL //#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2 //#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h> //#include <mkl.h>
// 本地方法 // 本地方法
#include <iostream> #include <iostream>
#include <Eigen/Core> #include <Eigen/Core>
@ -30,7 +30,7 @@ public:
Eigen::MatrixXd correctMatchModel(Eigen::MatrixXd r, Eigen::MatrixXd c); Eigen::MatrixXd correctMatchModel(Eigen::MatrixXd r, Eigen::MatrixXd c);
int outMatchModel(std::string matchmodel_path); int outMatchModel(std::string matchmodel_path);
//参数 //参数
Eigen::MatrixXd offsetXY_matrix; Eigen::MatrixXd offsetXY_matrix;
Eigen::MatrixXd match_model; Eigen::MatrixXd match_model;
}; };

21
LICENSE.txt
View File

@ -1,21 +0,0 @@
MIT License
Copyright (c) [year] [fullname]
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

78
OctreeNode.h
View File

@ -1,19 +1,19 @@
#pragma once #pragma once
#include <iostream> #include <iostream>
using namespace std; using namespace std;
//定义八叉树节点类 //定义八叉树节点类
template<class T> template<class T>
struct OctreeNode struct OctreeNode
{ {
T data; //节点数据 T data; //节点数据
T xmin, xmax; //节点坐标，即六面体个顶点的坐标 T xmin, xmax; //节点坐标，即六面体个顶点的坐标
T ymin, ymax; T ymin, ymax;
T zmin, zmax; T zmin, zmax;
OctreeNode <T>* top_left_front, * top_left_back; //该节点的个子结点 OctreeNode <T>* top_left_front, * top_left_back; //该节点的个子结点
OctreeNode <T>* top_right_front, * top_right_back; OctreeNode <T>* top_right_front, * top_right_back;
OctreeNode <T>* bottom_left_front, * bottom_left_back; OctreeNode <T>* bottom_left_front, * bottom_left_back;
OctreeNode <T>* bottom_right_front, * bottom_right_back; OctreeNode <T>* bottom_right_front, * bottom_right_back;
OctreeNode //节点类 OctreeNode //节点类
(T nodeValue = T(), (T nodeValue = T(),
T xminValue = T(), T xmaxValue = T(), T xminValue = T(), T xmaxValue = T(),
T yminValue = T(), T ymaxValue = T(), T yminValue = T(), T ymaxValue = T(),
@ -39,28 +39,28 @@ struct OctreeNode
bottom_right_front(bottom_right_front_Node), bottom_right_front(bottom_right_front_Node),
bottom_right_back(bottom_right_back_Node) {} bottom_right_back(bottom_right_back_Node) {}
}; };
//创建八叉树 //创建八叉树
template <class T> template <class T>
void createOctree(OctreeNode<T>*& root, int maxdepth, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) void createOctree(OctreeNode<T>*& root, int maxdepth, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
{ {
//cout<<"处理中，请稍候……"<<endl; //cout<<"处理中，请稍候……"<<endl;
maxdepth = maxdepth - 1; //每递归一次就将最大递归深度-1 maxdepth = maxdepth - 1; //每递归一次就将最大递归深度-1
if (maxdepth >= 0) if (maxdepth >= 0)
{ {
root = new OctreeNode<T>(); root = new OctreeNode<T>();
//cout << "请输入节点值:"; //cout << "请输入节点值:";
//root->data =9;//为节点赋值，可以存储节点信息，如物体可见性。由于是简单实现八叉树功能，简单赋值为9。 //root->data =9;//为节点赋值，可以存储节点信息，如物体可见性。由于是简单实现八叉树功能，简单赋值为9。
cin >> root->data; //为节点赋值 cin >> root->data; //为节点赋值
root->xmin = xmin; //为节点坐标赋值 root->xmin = xmin; //为节点坐标赋值
root->xmax = xmax; root->xmax = xmax;
root->ymin = ymin; root->ymin = ymin;
root->ymax = ymax; root->ymax = ymax;
root->zmin = zmin; root->zmin = zmin;
root->zmax = zmax; root->zmax = zmax;
double xm = (xmax - xmin) / 2;//计算节点个维度上的半边长 double xm = (xmax - xmin) / 2;//计算节点个维度上的半边长
double ym = (ymax - ymin) / 2; double ym = (ymax - ymin) / 2;
double zm = (ymax - ymin) / 2; double zm = (ymax - ymin) / 2;
//递归创建子树，根据每一个节点所处（是几号节点）的位置决定其子结点的坐标。 //递归创建子树，根据每一个节点所处（是几号节点）的位置决定其子结点的坐标。
createOctree(root->top_left_front, maxdepth, xmin, xmax - xm, ymax - ym, ymax, zmax - zm, zmax); createOctree(root->top_left_front, maxdepth, xmin, xmax - xm, ymax - ym, ymax, zmax - zm, zmax);
createOctree(root->top_left_back, maxdepth, xmin, xmax - xm, ymin, ymax - ym, zmax - zm, zmax); createOctree(root->top_left_back, maxdepth, xmin, xmax - xm, ymin, ymax - ym, zmax - zm, zmax);
createOctree(root->top_right_front, maxdepth, xmax - xm, xmax, ymax - ym, ymax, zmax - zm, zmax); createOctree(root->top_right_front, maxdepth, xmax - xm, xmax, ymax - ym, ymax, zmax - zm, zmax);
@ -72,13 +72,13 @@ void createOctree(OctreeNode<T>*& root, int maxdepth, double xmin, double xmax,
} }
} }
int i = 1; int i = 1;
//先序遍历八叉树 //先序遍历八叉树
template <class T> template <class T>
void preOrder(OctreeNode<T>*& p) void preOrder(OctreeNode<T>*& p)
{ {
if (p) if (p)
{ {
//cout << i << ".当前节点的值为：" << p->data << "\n坐标为："; //cout << i << ".当前节点的值为：" << p->data << "\n坐标为：";
//cout << "xmin: " << p->xmin << " xmax: " << p->xmax; //cout << "xmin: " << p->xmin << " xmax: " << p->xmax;
//cout << "ymin: " << p->ymin << " ymax: " << p->ymax; //cout << "ymin: " << p->ymin << " ymax: " << p->ymax;
//cout << "zmin: " << p->zmin << " zmax: " << p->zmax; //cout << "zmin: " << p->zmin << " zmax: " << p->zmax;
@ -95,7 +95,7 @@ void preOrder(OctreeNode<T>*& p)
cout << endl; cout << endl;
} }
} }
//求八叉树的深度 //求八叉树的深度
template<class T> template<class T>
int depth(OctreeNode<T>*& p) int depth(OctreeNode<T>*& p)
{ {
@ -111,7 +111,7 @@ int num(OctreeNode<T>*& p)
return 0; return 0;
return 1 + num(p->top_left_front) + num(p->top_left_back) + num(p->top_right_back) + num(p->top_right_front) + num(p->bottom_left_back) + num(p->bottom_left_front) + num(p->bottom_right_back) + num(p->bottom_right_front); return 1 + num(p->top_left_front) + num(p->top_left_back) + num(p->top_right_back) + num(p->top_right_front) + num(p->bottom_left_back) + num(p->bottom_left_front) + num(p->bottom_right_back) + num(p->bottom_right_front);
} }
//计算单位长度，为查找点做准备 //计算单位长度，为查找点做准备
int cal(int num) int cal(int num)
{ {
int result = 1; int result = 1;
@ -128,7 +128,7 @@ int cal(int num)
template<class T> template<class T>
int find(OctreeNode<T>*& p, double x, double y, double z) int find(OctreeNode<T>*& p, double x, double y, double z)
{ {
//查找点 //查找点
int maxdepth = 0; int maxdepth = 0;
int times = 0; int times = 0;
static double xmin = 0, xmax = 0, ymin = 0, ymax = 0, zmin = 0, zmax = 0; static double xmin = 0, xmax = 0, ymin = 0, ymax = 0, zmin = 0, zmax = 0;
@ -141,17 +141,17 @@ int find(OctreeNode<T>*& p, double x, double y, double z)
times++; times++;
if (x > xmax || x<xmin || y>ymax || y<ymin || z>zmax || z < zmin) if (x > xmax || x<xmin || y>ymax || y<ymin || z>zmax || z < zmin)
{ {
//cout << "该点不在场景中！" << endl; //cout << "该点不在场景中！" << endl;
return 0; return 0;
} }
if (x <= p->xmin + txm && x >= p->xmax - txm && y <= p->ymin + tym && y >= p->ymax - tym && z <= p->zmin + tzm && z >= p->zmax - tzm) if (x <= p->xmin + txm && x >= p->xmax - txm && y <= p->ymin + tym && y >= p->ymax - tym && z <= p->zmin + tzm && z >= p->zmax - tzm)
{ {
//cout << endl << "找到该点！" << "该点位于" << endl; //cout << endl << "找到该点！" << "该点位于" << endl;
//cout << "xmin: " << p->xmin << " xmax: " << p->xmax; //cout << "xmin: " << p->xmin << " xmax: " << p->xmax;
//cout << "ymin: " << p->ymin << " ymax: " << p->ymax; //cout << "ymin: " << p->ymin << " ymax: " << p->ymax;
//cout << "zmin: " << p->zmin << " zmax: " << p->zmax; //cout << "zmin: " << p->zmin << " zmax: " << p->zmax;
//cout << "节点内！" << endl; //cout << "节点内！" << endl;
//cout << "共经过" << times << "次递归！" << endl; //cout << "共经过" << times << "次递归！" << endl;
return 1; return 1;
} }
else if (x < (p->xmax - xm) && y < (p->ymax - ym) && z < (p->zmax - zm)) else if (x < (p->xmax - xm) && y < (p->ymax - ym) && z < (p->zmax - zm))
@ -187,16 +187,16 @@ int find(OctreeNode<T>*& p, double x, double y, double z)
find(p->top_right_front, x, y, z); find(p->top_right_front, x, y, z);
} }
} }
//main函数 //main函数
/* /*
int main() int main()
{ {
OctreeNode<double>* rootNode = NULL; OctreeNode<double>* rootNode = NULL;
int choiced = 0; int choiced = 0;
cout << "系统开始前请先创建八叉树" << endl; cout << "系统开始前请先创建八叉树" << endl;
cout << "请输入最大递归深度：" << endl; cout << "请输入最大递归深度：" << endl;
cin >> maxdepth; cin >> maxdepth;
cout << "请输入外包盒坐标，顺序如下：xmin,xmax,ymin,ymax,zmin,zmax" << endl; cout << "请输入外包盒坐标，顺序如下：xmin,xmax,ymin,ymax,zmin,zmax" << endl;
cin >> xmin >> xmax >> ymin >> ymax >> zmin >> zmax; cin >> xmin >> xmax >> ymin >> ymax >> zmin >> zmax;
if (maxdepth >= 0 || xmax > xmin || ymax > ymin || zmax > zmin || xmin > 0 || ymin > 0 || zmin > 0) if (maxdepth >= 0 || xmax > xmin || ymax > ymin || zmax > zmin || xmin > 0 || ymin > 0 || zmin > 0)
{ {
@ -210,12 +210,12 @@ int main()
{ {
system("cls"); system("cls");
cout << "请选择操作：\n"; cout << "请选择操作：\n";
cout << "\t1.计算空间中区域的个数\n"; cout << "\t1.计算空间中区域的个数\n";
cout << "\t2.先序遍历八叉树\n"; cout << "\t2.先序遍历八叉树\n";
cout << "\t3.查看树深度\n"; cout << "\t3.查看树深度\n";
cout << "\t4.查找节点 \n"; cout << "\t4.查找节点 \n";
cout << "\t0.退出\n"; cout << "\t0.退出\n";
cin >> choiced; cin >> choiced;
if (choiced == 0) if (choiced == 0)
@ -223,14 +223,14 @@ int main()
if (choiced == 1) if (choiced == 1)
{ {
system("cls"); system("cls");
cout << "空间区域个数" << endl; cout << "空间区域个数" << endl;
cout << num(rootNode); cout << num(rootNode);
} }
if (choiced == 2) if (choiced == 2)
{ {
system("cls"); system("cls");
cout << "先序遍历八叉树结果：/n"; cout << "先序遍历八叉树结果：/n";
i = 1; i = 1;
preOrder(rootNode); preOrder(rootNode);
cout << endl; cout << endl;
@ -240,24 +240,24 @@ int main()
{ {
system("cls"); system("cls");
int dep = depth(rootNode); int dep = depth(rootNode);
cout << "此八叉树的深度为" << dep + 1 << endl; cout << "此八叉树的深度为" << dep + 1 << endl;
system("pause"); system("pause");
} }
if (choiced == 4) if (choiced == 4)
{ {
system("cls"); system("cls");
cout << "请输入您希望查找的点的坐标，顺序如下：x,y,z\n"; cout << "请输入您希望查找的点的坐标，顺序如下：x,y,z\n";
double x, y, z; double x, y, z;
cin >> x >> y >> z; cin >> x >> y >> z;
times = 0; times = 0;
cout << endl << "开始搜寻该点……" << endl; cout << endl << "开始搜寻该点……" << endl;
find(rootNode, x, y, z); find(rootNode, x, y, z);
system("pause"); system("pause");
} }
else else
{ {
system("cls"); system("cls");
cout << "\n\n错误选择！\n"; cout << "\n\n错误选择！\n";
system("pause"); system("pause");
} }
} }

2
RPC_Correct.h
View File

@ -18,7 +18,7 @@ using namespace std;
using namespace Eigen; using namespace Eigen;
// 专门用于解析RPC // 专门用于解析RPC

60
SIMOrthoProgram.cpp
View File

@ -1,4 +1,4 @@
// SIMOrthoProgram.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。 // SIMOrthoProgram.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
// //
//#define EIGEN_USE_MKL_ALL //#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2 //#define EIGEN_VECTORIZE_SSE4_2
@ -13,7 +13,7 @@
#include <direct.h> #include <direct.h>
// gdal // gdal
#include <proj.h> #include <proj.h>
#include <string>' #include <string>
#include "gdal_priv.h" #include "gdal_priv.h"
#include "ogr_geometry.h" #include "ogr_geometry.h"
#include "gdalwarper.h" #include "gdalwarper.h"
@ -28,16 +28,16 @@ using namespace Eigen;
void PreProcess(int argc, char* argv[]) void PreProcess(int argc, char* argv[])
{ {
// .\baseTool\x64\Release\SIMOrthoProgram.exe 1 D:\MicroWorkspace\C-SAR\Ortho\Temporary\unpack\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_AHV_L10003923848\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_HH_L10003923848.tiff // .\baseTool\x64\Release\SIMOrthoProgram.exe 1 D:\MicroWorkspace\C-SAR\Ortho\Temporary\unpack\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_AHV_L10003923848\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_HH_L10003923848.tiff
//输入参数 //输入参数
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 Incident_path = argv[7];// 输出入射角文件 //std::string Incident_path = argv[7];// 输出入射角文件
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
std::cout << "in parameters:========================================================" << endl; std::cout << "in parameters:========================================================" << endl;
@ -96,7 +96,7 @@ void calInterpolation_cubic_Wgs84_rc_sar(int argc, char* argv[]) {
process.interpolation_GTC_sar(in_rc_wgs84_path, in_ori_sar_path, out_orth_sar_path, pstn); process.interpolation_GTC_sar(in_rc_wgs84_path, in_ori_sar_path, out_orth_sar_path, pstn);
} }
// mode 4 处理 RPC的入射角 // mode 4 处理 RPC的入射角
void getRPC_Incident_localIncident_angle(int argc, char* argv[]) { void getRPC_Incident_localIncident_angle(int argc, char* argv[]) {
std::cout << "mode 4: get RPC incident and local incident angle sar model:"; std::cout << "mode 4: get RPC incident and local incident angle sar model:";
std::cout << "SIMOrthoProgram.exe 4 in_parameter_path in_dem_path in_rpc_rc_path out_rpc_dem_path out_incident_angle_path out_local_incident_angle_path"; std::cout << "SIMOrthoProgram.exe 4 in_parameter_path in_dem_path in_rpc_rc_path out_rpc_dem_path out_incident_angle_path out_local_incident_angle_path";
@ -226,17 +226,17 @@ string GetExePath()
char szFilePath[MAX_PATH + 1] = { 0 }; char szFilePath[MAX_PATH + 1] = { 0 };
GetModuleFileNameA(NULL, szFilePath, MAX_PATH); GetModuleFileNameA(NULL, szFilePath, MAX_PATH);
/* /*
strrchr:函数功能：查找一个字符c在另一个字符串str中末次出现的位置（也就是从str的右侧开始查找字符c首次出现的位置）， strrchr:函数功能：查找一个字符c在另一个字符串str中末次出现的位置（也就是从str的右侧开始查找字符c首次出现的位置），
并返回这个位置的地址。如果未能找到指定字符，那么函数将返回NULL。 并返回这个位置的地址。如果未能找到指定字符，那么函数将返回NULL。
使用这个地址返回从最后一个字符c到str末尾的字符串。 使用这个地址返回从最后一个字符c到str末尾的字符串。
*/ */
(strrchr(szFilePath, '\\'))[0] = 0; // 删除文件名，只获得路径字串// (strrchr(szFilePath, '\\'))[0] = 0; // 删除文件名，只获得路径字串//
string path = szFilePath; string path = szFilePath;
return path; return path;
} }
/// <summary> /// <summary>
/// 初始化 /// 初始化
/// </summary> /// </summary>
void initProjEnv() { void initProjEnv() {
PJ_CONTEXT* C; PJ_CONTEXT* C;
@ -286,8 +286,8 @@ int main(int argc, char* argv[])
std::cout << "program start:\t" << getCurrentTimeString() << endl;; std::cout << "program start:\t" << getCurrentTimeString() << endl;;
int mode = 1; int mode = 1;
GDALAllRegister(); GDALAllRegister();
if (argc == 1) { // 测试参数 if (argc == 1) { // 测试参数
// 算法说明 // 算法说明
std::cout << "========================== description ================================================" << endl; std::cout << "========================== description ================================================" << endl;
std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl; std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl;
std::cout << "mode 1: Preprocess\n "; std::cout << "mode 1: Preprocess\n ";
@ -315,7 +315,7 @@ 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;
@ -327,7 +327,7 @@ int main(int argc, char* argv[])
else if (mode == 1) { else if (mode == 1) {
PreProcess(argc, argv); PreProcess(argc, argv);
} }
else if (mode == 2) { // RPC 计算模块 else if (mode == 2) { // RPC 计算模块
calIncident_localIncident_angle(argc, argv); calIncident_localIncident_angle(argc, argv);
} }
else if (mode == 3) { else if (mode == 3) {
@ -356,16 +356,16 @@ int main(int argc, char* argv[])
} }
} }
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
std::cout << "program over：\t" << getCurrentTimeString() << endl; std::cout << "program over：\t" << getCurrentTimeString() << endl;
} }
// 运行程序: Ctrl + F5 或调试 >“开始执行(不调试)”菜单 // 运行程序: Ctrl + F5 或调试 >“开始执行(不调试)”菜单
// 调试程序: F5 或调试 >“开始调试”菜单 // 调试程序: F5 或调试 >“开始调试”菜单
// 入门使用技巧: // 入门使用技巧:
// 1. 使用解决方案资源管理器窗口添加/管理文件 // 1. 使用解决方案资源管理器窗口添加/管理文件
// 2. 使用团队资源管理器窗口连接到源代码管理 // 2. 使用团队资源管理器窗口连接到源代码管理
// 3. 使用输出窗口查看生成输出和其他消息 // 3. 使用输出窗口查看生成输出和其他消息
// 4. 使用错误列表窗口查看错误 // 4. 使用错误列表窗口查看错误
// 5. 转到“项目”>“添加新项”以创建新的代码文件，或转到“项目”>“添加现有项”以将现有代码文件添加到项目 // 5. 转到“项目”>“添加新项”以创建新的代码文件，或转到“项目”>“添加现有项”以将现有代码文件添加到项目
// 6. 将来，若要再次打开此项目，请转到“文件”>“打开”>“项目”并选择 .sln 文件 // 6. 将来，若要再次打开此项目，请转到“文件”>“打开”>“项目”并选择 .sln 文件

BIN
SIMOrthoProgram.rar Normal file
View File

Binary file not shown.

22
SateOrbit.cpp
View File

@ -20,7 +20,7 @@ OrbitPoly::OrbitPoly()
OrbitPoly::OrbitPoly(int polynum, Eigen::MatrixX<double> polySatellitePara, double SatelliteModelStartTime) OrbitPoly::OrbitPoly(int polynum, Eigen::MatrixX<double> polySatellitePara, double SatelliteModelStartTime)
{ {
if (polySatellitePara.rows() != polynum||polySatellitePara.cols()!=6) { if (polySatellitePara.rows() != polynum||polySatellitePara.cols()!=6) {
throw exception("轨道模型参数形状不一致"); throw exception("?????????????????");
} }
this->polySatellitePara = polySatellitePara; this->polySatellitePara = polySatellitePara;
this->SatelliteModelStartTime = SatelliteModelStartTime; this->SatelliteModelStartTime = SatelliteModelStartTime;
@ -42,27 +42,27 @@ Eigen::MatrixX<double> OrbitPoly::SatelliteSpacePoint(long double satellitetime)
if (this->polySatellitePara.rows() != polynum || this->polySatellitePara.cols() != 6) { if (this->polySatellitePara.rows() != polynum || this->polySatellitePara.cols() != 6) {
throw exception("the size of satellitetime has error!! row: p1,p2,p3,p4 col: x,y,z,vx,vy,vz "); throw exception("the size of satellitetime has error!! row: p1,p2,p3,p4 col: x,y,z,vx,vy,vz ");
} }
// 更新轨道时间 // ?????????
double satellitetime2 =double( satellitetime - this->SatelliteModelStartTime); double satellitetime2 =double( satellitetime - this->SatelliteModelStartTime);
Eigen::MatrixX<double> satetime(1, polynum); Eigen::MatrixX<double> satetime(1, polynum);
for (int i = 0; i < polynum; i++) { for (int i = 0; i < polynum; i++) {
satetime(0, i) = pow(satellitetime2, i); satetime(0, i) = pow(satellitetime2, i);
} }
// 计算 // ????
Eigen::MatrixX<double> satellitePoints(1, 6); Eigen::MatrixX<double> satellitePoints(1, 6);
satellitePoints = satetime * polySatellitePara; satellitePoints = satetime * polySatellitePara;
return satellitePoints; return satellitePoints;
} }
/// <summary> /// <summary>
/// 获取指定时刻的卫星轨道坐标 /// ????????????????????
/// </summary> /// </summary>
/// <param name="satellitetime">卫星时刻</param> /// <param name="satellitetime">???????</param>
/// <param name="SatelliteModelStartTime">模型起算时间</param> /// <param name="SatelliteModelStartTime">??????????</param>
/// <param name="polySatellitePara">模型参数[x1,y1,z1,vx1,vy1,vz1; x2,y2,z2,vx2,vy2,vz2; ..... ]</param> /// <param name="polySatellitePara">??????[x1,y1,z1,vx1,vy1,vz1; x2,y2,z2,vx2,vy2,vz2; ..... ]</param>
/// <param name="polynum">模型参数数量</param> /// <param name="polynum">??????????</param>
/// <returns>卫星坐标</returns> /// <returns>????????</returns>
Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double>& satellitetime, double SatelliteModelStartTime, Eigen::MatrixX<double>& polySatellitePara, int polynum) Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double>& satellitetime, double SatelliteModelStartTime, Eigen::MatrixX<double>& polySatellitePara, int polynum)
{ {
if (satellitetime.cols() != 1) { if (satellitetime.cols() != 1) {
@ -71,7 +71,7 @@ Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double>& satellitetim
if (polySatellitePara.rows() != polynum || polySatellitePara.cols() != 6) { if (polySatellitePara.rows() != polynum || polySatellitePara.cols() != 6) {
throw exception("the size of satellitetime has error!! row: p1,p2,p3,p4 col: x,y,z,vx,vy,vz "); throw exception("the size of satellitetime has error!! row: p1,p2,p3,p4 col: x,y,z,vx,vy,vz ");
} }
// 更新轨道时间 // ?????????
int satellitetime_num = satellitetime.rows(); int satellitetime_num = satellitetime.rows();
satellitetime = satellitetime.array() - SatelliteModelStartTime; satellitetime = satellitetime.array() - SatelliteModelStartTime;
Eigen::MatrixX<double> satelliteTime(satellitetime_num, polynum); Eigen::MatrixX<double> satelliteTime(satellitetime_num, polynum);
@ -79,7 +79,7 @@ Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double>& satellitetim
satelliteTime.col(i) = satellitetime.array().pow(i); satelliteTime.col(i) = satellitetime.array().pow(i);
} }
// 计算 // ????
Eigen::MatrixX<double> satellitePoints(satellitetime_num, 6); Eigen::MatrixX<double> satellitePoints(satellitetime_num, 6);
satellitePoints = satelliteTime * polySatellitePara; satellitePoints = satelliteTime * polySatellitePara;
return satellitePoints; return satellitePoints;

18
SateOrbit.h
View File

@ -1,12 +1,12 @@
#pragma once #pragma once
/// ///
/// 计算卫星轨道 /// 计算卫星轨道
/// ///
//#define EIGEN_USE_MKL_ALL //#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2 //#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h> //#include <mkl.h>
// 本地方法 // 本地方法
#include "baseTool.h" #include "baseTool.h"
#include <iostream> #include <iostream>
@ -23,7 +23,7 @@ using namespace Eigen;
/// <summary> /// <summary>
/// 多项式轨道模型 /// 多项式轨道模型
/// </summary> /// </summary>
class OrbitPoly { class OrbitPoly {
public: public:
@ -43,11 +43,11 @@ public:
/// <summary> /// <summary>
/// 获取指定时刻的卫星轨道坐标 /// 获取指定时刻的卫星轨道坐标
/// </summary> /// </summary>
/// <param name="satellitetime">卫星时刻</param> /// <param name="satellitetime">卫星时刻</param>
/// <param name="SatelliteModelStartTime">模型起算时间</param> /// <param name="SatelliteModelStartTime">模型起算时间</param>
/// <param name="polySatellitePara">模型参数[x1,y1,z1,vx1,vy1,vz1; x2,y2,z2,vx2,vy2,vz2; ..... ]</param> /// <param name="polySatellitePara">模型参数[x1,y1,z1,vx1,vy1,vz1; x2,y2,z2,vx2,vy2,vz2; ..... ]</param>
/// <param name="polynum">模型参数数量</param> /// <param name="polynum">模型参数数量</param>
/// <returns>卫星坐标</returns> /// <returns>卫星坐标</returns>
Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double> &satellitetime, double SatelliteModelStartTime, Eigen::MatrixX<double>& polySatellitePara, int polynum = 4); Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double> &satellitetime, double SatelliteModelStartTime, Eigen::MatrixX<double>& polySatellitePara, int polynum = 4);

86
baseTool.cpp
View File

@ -125,7 +125,7 @@ Landpoint crossProduct(const Landpoint& a, const Landpoint& b) {
/// n21=n2xn1 /// n21=n2xn1
/// n=(n21+n32+n43+n41)*0.25; /// n=(n21+n32+n43+n41)*0.25;
/// </summary> /// </summary>
/// <param name="p0">Ŀ<EFBFBD><EFBFBD>ￄ1<EFBFBD>7</param> /// <param name="p0">目锟斤拷锟?1锟?7</param>
/// <param name="p1"><3E><>Χ<EFBFBD><CEA7>1</param> /// <param name="p1"><3E><>Χ<EFBFBD><CEA7>1</param>
/// <param name="p2"><3E><>Χ<EFBFBD><CEA7>2</param> /// <param name="p2"><3E><>Χ<EFBFBD><CEA7>2</param>
/// <param name="p3"><3E><>Χ<EFBFBD><CEA7>3</param> /// <param name="p3"><3E><>Χ<EFBFBD><CEA7>3</param>
@ -195,8 +195,8 @@ Landpoint getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landp
/// p31 p32 p33 p34 /// p31 p32 p33 p34
/// p41 p42 p43 p44 /// p41 p42 p43 p44
/// </summary> /// </summary>
/// <param name="u"><EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7</param> /// <param name="u">锟斤拷元锟斤拷锟斤拷锟斤拷锟叫★拷锟斤拷锟斤拷锟?1锟?7</param>
/// <param name="v"><EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7</param> /// <param name="v">锟斤拷元锟斤拷锟斤拷锟斤拷锟叫★拷锟斤拷锟斤拷锟?1锟?7</param>
/// <param name="img">4x4 <20><>ֵ<EFBFBD><D6B5><EFBFBD><EFBFBD></param> /// <param name="img">4x4 <20><>ֵ<EFBFBD><D6B5><EFBFBD><EFBFBD></param>
/// <returns></returns> /// <returns></returns>
complex<double> Cubic_Convolution_interpolation(double u, double v, Eigen::MatrixX<complex<double>> img) complex<double> Cubic_Convolution_interpolation(double u, double v, Eigen::MatrixX<complex<double>> img)
@ -341,17 +341,17 @@ Landpoint XYZ2LLA(const Landpoint& XYZ) {
} }
/// <summary> /// <summary>
/// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ<EFBFBD><EFBFBD><EFBFBD>ȡӰ<EFBFBD>ￄ1<EFBFBD>7 /// 锟斤拷锟斤拷图锟斤拷锟饺∮帮拷锟?1锟?7
/// </summary> /// </summary>
/// <param name="dem_path"></param> /// <param name="dem_path"></param>
gdalImage::gdalImage(string raster_path) gdalImage::gdalImage(string raster_path)
{ {
omp_lock_t lock; omp_lock_t lock;
omp_init_lock(&lock); // <20><>ʼ<EFBFBD><CABC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> omp_init_lock(&lock); // <20><>ʼ<EFBFBD><CABC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
omp_set_lock(&lock); //<EFBFBD><EFBFBD>û<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 omp_set_lock(&lock); //锟斤拷没锟斤拷锟斤拷锟?1锟?7
this->img_path = raster_path; this->img_path = raster_path;
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 GDALAllRegister();// 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
// <20><>DEMӰ<4D><D3B0> // <20><>DEMӰ<4D><D3B0>
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.c_str(), GA_ReadOnly));//<2F><>ֻ<EFBFBD><D6BB><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD>Ӱ<EFBFBD><D3B0> GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.c_str(), GA_ReadOnly));//<2F><>ֻ<EFBFBD><D6BB><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD>Ӱ<EFBFBD><D3B0>
this->width = rasterDataset->GetRasterXSize(); this->width = rasterDataset->GetRasterXSize();
@ -406,10 +406,10 @@ void gdalImage::setData(Eigen::MatrixXd data)
Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids = 1) Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids = 1)
{ {
omp_lock_t lock; omp_lock_t lock;
omp_init_lock(&lock); // <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> omp_init_lock(&lock);
omp_set_lock(&lock); //<2F><>û<EFBFBD><C3BB><EFBFBD><EFBFBD>ￄ1<EFBF84>7 omp_set_lock(&lock);
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 GDALAllRegister();
// <EFBFBD><EFBFBD>DEMӰ<EFBFBD><EFBFBD>
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_ReadOnly));//<2F><>ֻ<EFBFBD><D6BB><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD>Ӱ<EFBFBD><D3B0> GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_ReadOnly));//<2F><>ֻ<EFBFBD><D6BB><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD>Ӱ<EFBFBD><D3B0>
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
@ -473,15 +473,13 @@ Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count,
void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0, int band_ids = 1) void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0, int band_ids = 1)
{ {
omp_lock_t lock; omp_lock_t lock;
omp_init_lock(&lock); // <20><>ʼ<EFBFBD><CABC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> omp_init_lock(&lock);
omp_set_lock(&lock); //<2F><>û<EFBFBD><C3BB><EFBFBD><EFBFBD>ￄ1<EFBF84>7 omp_set_lock(&lock);
// <20>Ϸ<EFBFBD><CFB7>Լ<EFBFBD>ￄ1<EFBF84>7
// д<><D0B4><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) { if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
string tip = "д<EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD>Χ<EFBFBD><EFBFBD>Խ<EFBFBD><EFBFBD>Ŀ<EFBFBD>귶Χ" + this->img_path; string tip = "file path: " + this->img_path;
throw exception(tip.c_str()); throw exception(tip.c_str());
} }
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 GDALAllRegister();
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = nullptr; GDALDataset* poDstDS = nullptr;
if (boost::filesystem::exists(this->img_path)) { if (boost::filesystem::exists(this->img_path)) {
@ -491,7 +489,7 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
poDstDS = poDriver->Create(this->img_path.c_str(), this->width, this->height, this->band_num, GDT_Float32, NULL); // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> poDstDS = poDriver->Create(this->img_path.c_str(), this->width, this->height, this->band_num, GDT_Float32, NULL); // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
poDstDS->SetProjection(this->projection.c_str()); poDstDS->SetProjection(this->projection.c_str());
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
double gt_ptr[6]; double gt_ptr[6];
for (int i = 0; i < this->gt.rows(); i++) { for (int i = 0; i < this->gt.rows(); i++) {
for (int j = 0; j < this->gt.cols(); j++) { for (int j = 0; j < this->gt.cols(); j++) {
@ -533,7 +531,7 @@ void gdalImage::setNoDataValue(double nodatavalue = -9999, int band_ids = 1)
{ {
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 GDALAllRegister();// 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
//GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); //GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.c_str(), GA_Update)); GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.c_str(), GA_Update));
poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue); poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue);
@ -656,7 +654,7 @@ GDALRPCInfo gdalImage::getRPC()
//<2F><>Ԫ<EFBFBD><D4AA><EFBFBD><EFBFBD><EFBFBD>л<EFBFBD>ȡRPC<50><43>Ϣ //<2F><>Ԫ<EFBFBD><D4AA><EFBFBD><EFBFBD><EFBFBD>л<EFBFBD>ȡRPC<50><43>Ϣ
char** papszRPC = pDS->GetMetadata("RPC"); char** papszRPC = pDS->GetMetadata("RPC");
//<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD>RPC<EFBFBD><EFBFBD>Ϣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɽṹ<EFBFBD>ￄ1<EFBFBD>7 //锟斤拷锟斤拷取锟斤拷RPC锟斤拷息锟斤拷锟斤拷山峁癸拷锟?1锟?7
GDALRPCInfo oInfo; GDALRPCInfo oInfo;
GDALExtractRPCInfo(papszRPC, &oInfo); GDALExtractRPCInfo(papszRPC, &oInfo);
@ -691,7 +689,7 @@ Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points)
Eigen::MatrixXd gdalImage::getHist(int bandids) Eigen::MatrixXd gdalImage::getHist(int bandids)
{ {
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 GDALAllRegister();// 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
// <20><>DEMӰ<4D><D3B0> // <20><>DEMӰ<4D><D3B0>
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_ReadOnly));//<2F><>ֻ<EFBFBD><D6BB><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD>Ӱ<EFBFBD><D3B0> GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_ReadOnly));//<2F><>ֻ<EFBFBD><D6BB><EFBFBD><EFBFBD>ʽ<EFBFBD><CABD>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD>Ӱ<EFBFBD><D3B0>
@ -723,7 +721,7 @@ gdalImage CreategdalImage(string img_path, int height, int width, int band_num,
if (boost::filesystem::exists(img_path.c_str())) { if (boost::filesystem::exists(img_path.c_str())) {
boost::filesystem::remove(img_path.c_str()); boost::filesystem::remove(img_path.c_str());
} }
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 GDALAllRegister();// 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = poDriver->Create(img_path.c_str(), width, height, band_num, GDT_Float32, NULL); // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> GDALDataset* poDstDS = poDriver->Create(img_path.c_str(), width, height, band_num, GDT_Float32, NULL); // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
if (need_gt) { if (need_gt) {
@ -766,9 +764,9 @@ void clipGdalImage(string in_path, string out_path, DemBox box, double pixelinte
/*** /***
* ң<EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD><EFBFBD>ز<EFBFBD><EFBFBD><EFBFBD> (Ҫ<EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͶӰ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߲<EFBFBD>̈́1<EFBFBD>7) * 遥锟斤拷影锟斤拷锟截诧拷锟斤拷 (要锟斤拷影锟斤拷锟斤拷锟斤拷锟酵队帮拷锟斤拷锟斤拷锟斤拷卟锟酵?1锟?7)
* @param pszSrcFile <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD><EFBFBD><EFBFBD>·<EFBFBD><EFBFBD> * @param pszSrcFile <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD><EFBFBD><EFBFBD>·<EFBFBD><EFBFBD>
* @param pszOutFile д<EFBFBD><EFBFBD>Ľ<EFBFBD><EFBFBD>ͼ<EFBFBD><EFBFBD><EFBFBD>·<EFBFBD>ￄ1<EFBFBD>7 * @param pszOutFile 写锟斤拷慕锟斤拷图锟斤拷锟铰凤拷锟?1锟?7
* @param eResample <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>֣<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>μ<EFBFBD>GDALResampleAlg<EFBFBD><EFBFBD><EFBFBD>壬Ĭ<EFBFBD><EFBFBD>Ϊ˫<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڲ<EFBFBD> * @param eResample <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>֣<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>μ<EFBFBD>GDALResampleAlg<EFBFBD><EFBFBD><EFBFBD>壬Ĭ<EFBFBD><EFBFBD>Ϊ˫<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڲ<EFBFBD>
* @param gt Xת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȣ<EFBFBD>Ĭ<EFBFBD>ϴ<EFBFBD>СΪ1.0<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>1ͼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD>1<EFBFBD><EFBFBD>ʾͼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD>ϵ<EFBFBD><EFBFBD>ڲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ<EFBFBD><EFBFBD>Ŀ<EFBFBD><EFBFBD>ȺͲ<EFBFBD><EFBFBD><EFBFBD>ǰͼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȵı<EFBFBD> * @param gt Xת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȣ<EFBFBD>Ĭ<EFBFBD>ϴ<EFBFBD>СΪ1.0<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>1ͼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD>1<EFBFBD><EFBFBD>ʾͼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD>ϵ<EFBFBD><EFBFBD>ڲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ<EFBFBD><EFBFBD>Ŀ<EFBFBD><EFBFBD>ȺͲ<EFBFBD><EFBFBD><EFBFBD>ǰͼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȵı<EFBFBD>
* @param new_width * @param new_width
@ -802,7 +800,7 @@ int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int
char* pszSrcWKT = NULL; char* pszSrcWKT = NULL;
pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef()); pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef());
//<EFBFBD><EFBFBD><EFBFBD>û<EFBFBD><EFBFBD>ͶӰ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD>ￄ1<EFBFBD>7 //锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
if (strlen(pszSrcWKT) <= 0) if (strlen(pszSrcWKT) <= 0)
{ {
OGRSpatialReference oSRS; OGRSpatialReference oSRS;
@ -815,7 +813,7 @@ int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int
void* hTransformArg; void* hTransformArg;
hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT, FALSE, 0.0, 1); hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT, FALSE, 0.0, 1);
std::cout << "no proj " << endl; std::cout << "no proj " << endl;
//(û<EFBFBD><EFBFBD>ͶӰ<EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߲<EFBFBD>ͨ<EFBFBD>ￄ1<EFBFBD>7) //(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7)
if (hTransformArg == NULL) if (hTransformArg == NULL)
{ {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
@ -832,7 +830,7 @@ int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int
//GDALDestroyGenImgProjTransformer(hTransformArg); //GDALDestroyGenImgProjTransformer(hTransformArg);
//<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݼￄ1<EFBFBD>7 //锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷菁锟?1锟?7
GDALDataset* pDDst = pDriver->Create(pszOutFile, new_width, new_height, nBandCount, dataType, NULL); GDALDataset* pDDst = pDriver->Create(pszOutFile, new_width, new_height, nBandCount, dataType, NULL);
if (pDDst == NULL) if (pDDst == NULL)
{ {
@ -1074,19 +1072,19 @@ WGS84_J2000::~WGS84_J2000()
} }
/// <summary> /// <summary>
/// step1 WGS 84 ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Э<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD>ￄ1<EFBFBD>7 /// step1 WGS 84 转锟斤拷锟斤拷协锟斤拷锟斤拷锟斤拷锟斤拷锟较碉拷锟?1锟?7
/// <EFBFBD><EFBFBD>ؾ<EFBFBD>γ<EFBFBD>ȸ߳<EFBFBD>BLH<EFBFBD><EFBFBD>ع<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>ת<EFBFBD>ￄ1<EFBFBD>7 BLH-- > XG /// 锟斤拷鼐锟轿筹拷雀叱锟紹LH锟斤拷毓锟斤拷锟斤拷锟较碉拷锟阶?拷锟?1锟?7 BLH-- > XG
/// ECEF, Earth Centered Earth Fixed; <EFBFBD>ع<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ <20>ο<EFBFBD>ƽ<EFBFBD>棺ƽ<E6A3BA><C6BD><EFBFBD><EFBFBD>棬<EFBFBD><E6A3AC><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><C4A3><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>CIO<49><4F><EFBFBD><EFBFBD><EFBFBD>ߴ<EFBFBD>ֱ<EFBFBD><D6B1>ƽ<EFBFBD>档 /// ECEF, Earth Centered Earth Fixed; 锟截癸拷锟斤拷锟斤拷系 锟轿匡拷平锟芥：平锟斤拷锟斤拷妫?拷锟斤拷锟斤拷锟斤拷模锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷CIO锟斤拷锟斤拷锟竭达拷直锟斤拷平锟芥。
/// x<><78>Ϊ<EFBFBD>ο<EFBFBD>ƽ<EFBFBD><C6BD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD>潻<EFBFBD>ߣ<EFBFBD>z<EFBFBD><7A>Ϊ<EFBFBD><CEAA><EFBFBD><EFBFBD>ָ<EFBFBD><D6B8>CIO<49>㡣 /// x<><78>Ϊ<EFBFBD>ο<EFBFBD>ƽ<EFBFBD><C6BD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD>潻<EFBFBD>ߣ<EFBFBD>z<EFBFBD><7A>Ϊ<EFBFBD><CEAA><EFBFBD><EFBFBD>ָ<EFBFBD><D6B8>CIO<49>㡣
/// </summary> /// </summary>
/// <param name="LBH_deg_m">B L H<EFBFBD>ֱ<EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD>γ<EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD>ؾ<EFBFBD><EFBFBD>Ⱥʹ<EFBFBD>ظ߳ￄ1<EFBFBD>7 <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ΪN * 3<><33><EFBFBD><EFBFBD></param> /// <param name="LBH_deg_m">B L H锟街憋拷为锟斤拷锟轿筹拷取锟斤拷锟截撅拷锟饺和达拷馗叱锟?1锟?7 锟斤拷锟斤拷锟斤拷锟轿狽 * 3锟斤拷锟斤拷</param>
/// <returns>XYZ Ϊ<>ع<EFBFBD><D8B9><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><CFB5> x y z<><7A><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ΪN * 3<><33><EFBFBD><EFBFBD></returns> /// <returns>XYZ Ϊ<>ع<EFBFBD><D8B9><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><CFB5> x y z<><7A><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ΪN * 3<><33><EFBFBD><EFBFBD></returns>
Eigen::MatrixXd WGS84_J2000::WGS84TECEF(Eigen::MatrixXd LBH_deg_m) Eigen::MatrixXd WGS84_J2000::WGS84TECEF(Eigen::MatrixXd LBH_deg_m)
{ {
return LLA2XYZ(LBH_deg_m); return LLA2XYZ(LBH_deg_m);
} }
/// <summary> /// <summary>
/// step 2 Э<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>τ1<EFBFBD>7 ת<><D7AA>Ϊ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ /// step 2 协锟斤拷锟斤拷锟斤拷锟斤拷锟较?1锟?7 转锟斤拷为瞬时锟斤拷锟斤拷锟斤拷锟斤拷系
/// <20><><EFBFBD><EFBFBD>Ҫ<EFBFBD>漰<EFBFBD><E6BCB0>ѯ <20><><EFBFBD><EFBFBD>ʱ<EFBFBD>̵<EFBFBD> xp,yp , <20><><EFBFBD>Բ<EFBFBD>ѯEOP<4F>ļ<EFBFBD><C4BC><EFBFBD>ã<EFBFBD><C3A3><EFBFBD><EFBFBD>ص<EFBFBD>ַ<EFBFBD><D6B7><EFBFBD><EFBFBD>http://celestrak.com/spacedata/ /// <20><><EFBFBD><EFBFBD>Ҫ<EFBFBD>漰<EFBFBD><E6BCB0>ѯ <20><><EFBFBD><EFBFBD>ʱ<EFBFBD>̵<EFBFBD> xp,yp , <20><><EFBFBD>Բ<EFBFBD>ѯEOP<4F>ļ<EFBFBD><C4BC><EFBFBD>ã<EFBFBD><C3A3><EFBFBD><EFBFBD>ص<EFBFBD>ַ<EFBFBD><D6B7><EFBFBD><EFBFBD>http://celestrak.com/spacedata/
/// earthFixedXYZ=ordinateSingleRotate('x',yp)*ordinateSingleRotate('y',xp)*earthFixedXYZ; /// earthFixedXYZ=ordinateSingleRotate('x',yp)*ordinateSingleRotate('y',xp)*earthFixedXYZ;
/// </summary> /// </summary>
@ -1128,17 +1126,17 @@ Eigen::MatrixXd WGS84_J2000::ordinateSingleRotate(int axis, double angle_deg)
/// UNTITLED <20><><EFBFBD>㵱<EFBFBD>غ<EFBFBD><D8BA><EFBFBD>ʱ,<2C><><EFBFBD><EFBFBD>ֵ<EFBFBD><D6B5>ʱ<EFBFBD><CAB1>Ϊ<EFBFBD><CEAA>λ /// UNTITLED <20><><EFBFBD>㵱<EFBFBD>غ<EFBFBD><D8BA><EFBFBD>ʱ,<2C><><EFBFBD><EFBFBD>ֵ<EFBFBD><D6B5>ʱ<EFBFBD><CAB1>Ϊ<EFBFBD><CEAA>λ
/// %dAT <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> /// %dAT <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
/// % TAI<41><49><EFBFBD><EFBFBD>ԭ<EFBFBD><D4AD>ʱ<EFBFBD><CAB1> <20><>ʱ<EFBFBD><CAB1><EFBFBD><EFBFBD>׼ TAI = UTC + dAT;% <20><><EFBFBD><EFBFBD>ԭ<EFBFBD><D4AD>ʱ<EFBFBD><CAB1> /// % TAI<41><49><EFBFBD><EFBFBD>ԭ<EFBFBD><D4AD>ʱ<EFBFBD><CAB1> <20><>ʱ<EFBFBD><CAB1><EFBFBD><EFBFBD>׼ TAI = UTC + dAT;% <20><><EFBFBD><EFBFBD>ԭ<EFBFBD><D4AD>ʱ<EFBFBD><CAB1>
/// % TT <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ TT = TAI + 32.184 <20><>2014<31><34><EFBFBD>ʱ<EFBFBD>ￄ1<EFBF84>7; /// % TT 锟斤拷锟斤拷时 TT = TAI + 32.184 锟斤拷2014锟斤拷锟绞憋拷锟?1锟?7;
/// % TDT <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1> /// % TDT <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1>
/// % ET <20><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1> /// % ET <20><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1>
/// % <20><><EFBFBD><EFBFBD>ʱ = <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1> = <20><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1> /// % <20><><EFBFBD><EFBFBD>ʱ = <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1> = <20><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1>
/// %<25><><EFBFBD><EFBFBD>ʱ=<3D><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1>=<3D><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1> /// %<25><><EFBFBD><EFBFBD>ʱ=<3D><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1>=<3D><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1>
/// </summary> /// </summary>
/// <param name="UTC"><3E><>ʽΪy m d <20><><EFBFBD><EFBFBD>d<EFBFBD><64><EFBFBD><EFBFBD>ֵΪС<CEAA><D0A1><EFBFBD><EFBFBD><EFBFBD><EFBFBD>h m s <20><><EFBFBD><EFBFBD>sec/3600+min/60+h)/24ת<34><D7AA><EFBFBD><EFBFBD>С<EFBFBD><D0A1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ۼӵ<DBBC>day <20><></param> /// <param name="UTC"><3E><>ʽΪy m d <20><><EFBFBD><EFBFBD>d<EFBFBD><64><EFBFBD><EFBFBD>ֵΪС<CEAA><D0A1><EFBFBD><EFBFBD><EFBFBD><EFBFBD>h m s <20><><EFBFBD><EFBFBD>sec/3600+min/60+h)/24ת<34><D7AA><EFBFBD><EFBFBD>С<EFBFBD><D0A1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ۼӵ<DBBC>day <20><></param>
/// <param name="dUT1">dUT1 Ϊut1-utc <20>IJ<C4B2><EEA3AC>ֵ<EFBFBD><D6B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>1<EFBFBD>룬<EFBFBD><EBA3AC>iers<72>ɻ<EFBFBD><C9BB><EFBFBD><EFBFBD>ք1<D684>7 0</param> /// <param name="dUT1">dUT1 为ut1-utc 锟侥差，锟斤拷值锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷1锟诫，锟斤拷iers锟缴伙拷锟斤拷锟街?1锟?7 0</param>
/// <param name="dAT"><3E><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> 37</param> /// <param name="dAT"><3E><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> 37</param>
/// <param name="gst_deg"><3E><><EFBFBD><EFBFBD>ԭ<EFBFBD><D4AD>ʱ<EFBFBD><CAB1> <20><>ʱ<EFBFBD><CAB1><EFBFBD><EFBFBD>׼ TAI=UTC+dAT; %<25><><EFBFBD><EFBFBD>ԭ<EFBFBD><D4AD>ʱ<EFBFBD><CAB1></param> /// <param name="gst_deg"><3E><><EFBFBD><EFBFBD>ԭ<EFBFBD><D4AD>ʱ<EFBFBD><CAB1> <20><>ʱ<EFBFBD><CAB1><EFBFBD><EFBFBD>׼ TAI=UTC+dAT; %<25><><EFBFBD><EFBFBD>ԭ<EFBFBD><D4AD>ʱ<EFBFBD><CAB1></param>
/// <param name="JDTDB"><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ TT = TAI+32.184 <20><>2014<31><34><EFBFBD>ʱ<EFBFBD>ￄ1<EFBF84>7;</param> /// <param name="JDTDB">锟斤拷锟斤拷时 TT = TAI+32.184 锟斤拷2014锟斤拷锟绞憋拷锟?1锟?7;</param>
/// <returns></returns> /// <returns></returns>
int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& gst_deg, double& JDTDB) int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& gst_deg, double& JDTDB)
{ {
@ -1156,7 +1154,7 @@ int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& g
//% JDTT = YMD2JD(UTC(1), UTC(2), UTC(3)) + dUT1 / 86400; //% JDTT = YMD2JD(UTC(1), UTC(2), UTC(3)) + dUT1 / 86400;
//% <EFBFBD><EFBFBD><EFBFBD>ȼ<EFBFBD><EFBFBD><EFBFBD>TDB, TDB<44><42><EFBFBD><EFBFBD><EFBFBD>Ķ<EFBFBD><C4B6><EFBFBD>ѧʱ<D1A7><CAB1><EFBFBD><EFBFBD>̫<EFBFBD><CCAB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ǵ<EFBFBD><C7B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>е<EFBFBD>ʱ<EFBFBD><CAB1>߶ￄ1<EFBF84>7 //% 锟斤拷锟饺硷拷锟斤拷TDB, TDB锟斤拷锟斤拷锟侥讹拷锟斤拷学时锟斤拷锟斤拷太锟斤拷锟斤拷锟斤拷锟斤拷锟角碉拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟叫碉拷时锟斤拷叨锟?1锟?7
double T = (JDTT - J2000) / 36525; double T = (JDTT - J2000) / 36525;
double temp = 0.001657 * sin(628.3076 * T + 6.2401) double temp = 0.001657 * sin(628.3076 * T + 6.2401)
+ 0.000022 * sin(575.3385 * T + 4.2970) + 0.000022 * sin(575.3385 * T + 4.2970)
@ -1167,7 +1165,7 @@ int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& g
+ 0.00001 * T * sin(628.3076 * T + 4.2490); + 0.00001 * T * sin(628.3076 * T + 4.2490);
JDTDB = JDTT + temp / 86400; JDTDB = JDTT + temp / 86400;
//% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>UT2, UT2<54><32><EFBFBD><EFBFBD>UT1<54>Ļ<EFBFBD><C4BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><D4BC>ڱ仯<DAB1><E4BBAF>õ<EFBFBD><C3B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD>ￄ1<EFBF84>7 //% 锟斤拷锟斤拷锟斤拷锟経T2, UT2锟斤拷锟斤拷UT1锟侥伙拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟皆硷拷锟节变化锟斤拷玫锟斤拷锟斤拷锟斤拷锟绞憋拷锟?1锟?7
// %% <20><><EFBFBD><EFBFBD>UT1<54><31><EFBFBD><EFBFBD>Tb, TbΪ<62>Ա<EFBFBD><D4B1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><CEAA>λ<EFBFBD><CEBB><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ԪB2000.0<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ff // %% <20><><EFBFBD><EFBFBD>UT1<54><31><EFBFBD><EFBFBD>Tb, TbΪ<62>Ա<EFBFBD><D4B1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><CEAA>λ<EFBFBD><CEBB><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ԪB2000.0<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ff
// % B2000 = 2451544.033; // % B2000 = 2451544.033;
//% Tb = (YMD2JD(UT1(1), UT1(2), UT1(3)) - B2000) / 365.2422; //% Tb = (YMD2JD(UT1(1), UT1(2), UT1(3)) - B2000) / 365.2422;
@ -1182,14 +1180,14 @@ int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& g
double T4 = T3 * T; double T4 = T3 * T;
double T5 = T4 * T; double T5 = T4 * T;
double Du = JDUT1 - J2000; double Du = JDUT1 - J2000;
double thita = 0.7790572732640 + 1.00273781191135448 * Du;//% <EFBFBD><EFBFBD>J2000<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ȧ<EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 double thita = 0.7790572732640 + 1.00273781191135448 * Du;//% 锟斤拷J2000锟斤拷锟斤拷锟阶?拷锟斤拷锟饺︼拷锟斤拷锟?1锟?7
double GMST_deg = (-0.0000000368 * T5 - 0.000029956 * T4 - 0.00000044 * T3 + 1.3915817 * T2 + 4612.156534 * T + 0.014506) / 3600 + (thita - floor(thita)) * 360; double GMST_deg = (-0.0000000368 * T5 - 0.000029956 * T4 - 0.00000044 * T3 + 1.3915817 * T2 + 4612.156534 * T + 0.014506) / 3600 + (thita - floor(thita)) * 360;
double epthilongA_deg, dertaPthi_deg, dertaEpthilong_deg, epthilong_deg; double epthilongA_deg, dertaPthi_deg, dertaEpthilong_deg, epthilong_deg;
WGS84_J2000::nutationInLongitudeCaculate(JDTDB, epthilongA_deg, dertaPthi_deg, dertaEpthilong_deg, epthilong_deg); WGS84_J2000::nutationInLongitudeCaculate(JDTDB, epthilongA_deg, dertaPthi_deg, dertaEpthilong_deg, epthilong_deg);
//[epthilongA_deg, dertaPthi_deg] = nutationInLongitudeCaculate(JDTDB); //[epthilongA_deg, dertaPthi_deg] = nutationInLongitudeCaculate(JDTDB);
//% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ྭ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>eclipticObliquitygama<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD> //% 锟斤拷锟斤拷嗑?拷露锟斤拷锟斤拷锟斤拷锟斤拷锟絜clipticObliquitygama锟斤拷锟斤拷位为锟斤拷锟斤拷
//% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7 ̫<><CCAB>ƽ<EFBFBD><C6BD><EFBFBD><EFBFBD>ￄ1<EFBF84>7 <20><><EFBFBD><EFBFBD>γ<EFBFBD>ȷ<EFBFBD><C8B7><EFBFBD> <20><><EFBFBD><EFBFBD>ƽ<EFBFBD>Ǿ<EFBFBD>(<28><><EFBFBD><EFBFBD>ƽ<EFBFBD>ƾ<EFBFBD> - ̫<><CCAB>ƽ<EFBFBD>ƾ<EFBFBD><C6BE><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD>ƾ<EFBFBD> <20><>λΪ<CEBB><CEAA><EFBFBD><EFBFBD> //% 锟斤拷锟斤拷锟斤拷锟斤拷平锟斤拷锟斤拷锟?1锟?7 太锟斤拷平锟斤拷锟斤拷锟?1锟?7 锟斤拷锟斤拷纬锟饺凤拷锟斤拷 锟斤拷锟斤拷平锟角撅拷(锟斤拷锟斤拷平锟狡撅拷 - 太锟斤拷平锟狡撅拷锟斤拷 锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷平锟狡撅拷 锟斤拷位为锟斤拷锟斤拷
double F_deg = (0.00000417 * T4 - 0.001037 * T3 - 12.7512 * T2 + 1739527262.8478 * T + 335779.526232) / 3600; double F_deg = (0.00000417 * T4 - 0.001037 * T3 - 12.7512 * T2 + 1739527262.8478 * T + 335779.526232) / 3600;
F_deg = fmod(F_deg, 360); F_deg = fmod(F_deg, 360);
double D_deg = (-0.00003169 * T4 + 0.006593 * T3 - 6.3706 * T2 + 1602961601.2090 * T + 1072260.70369) / 3600; double D_deg = (-0.00003169 * T4 + 0.006593 * T3 - 6.3706 * T2 + 1602961601.2090 * T + 1072260.70369) / 3600;
@ -1217,7 +1215,7 @@ int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& g
/// <20><><EFBFBD><EFBFBD>ƽ<EFBFBD>Ƴཻ<C6B3>ǣ<EFBFBD><C7A3>ƾ<EFBFBD><C6BE><EFBFBD><C2B6><EFBFBD><EFBFBD>ͽ<EFBFBD><CDBD><EFBFBD><EFBFBD><EFBFBD><C2B6><EFBFBD>˲ʱ<CBB2>Ƴཻ<C6B3>ǡ<EFBFBD> /// <20><><EFBFBD><EFBFBD>ƽ<EFBFBD>Ƴཻ<C6B3>ǣ<EFBFBD><C7A3>ƾ<EFBFBD><C6BE><EFBFBD><C2B6><EFBFBD><EFBFBD>ͽ<EFBFBD><CDBD><EFBFBD><EFBFBD><EFBFBD><C2B6><EFBFBD>˲ʱ<CBB2>Ƴཻ<C6B3>ǡ<EFBFBD>
/// </summary> /// </summary>
/// <param name="JD"><3E><><EFBFBD><EFBFBD><EFBFBD><EFBFBD></param> /// <param name="JD"><3E><><EFBFBD><EFBFBD><EFBFBD><EFBFBD></param>
/// <param name="accuracy"><EFBFBD><EFBFBD>ʾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ľ<EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD>ￄ1<EFBFBD>7 <20><>ֵΪ<D6B5><CEAA>normal<61><6C> <20>͡<EFBFBD>high<67><68> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>'N'<27>͡<EFBFBD>H'<27><>ʾһ<CABE>㾫<EFBFBD>Ⱥ͸߾<CDB8><DFBE><EFBFBD> <20><>Ĭ<EFBFBD><C4AC>Ϊ<EFBFBD>߾<EFBFBD><DFBE>ȼ<EFBFBD><C8BC><EFBFBD></param> /// <param name="accuracy">锟斤拷示锟斤拷锟斤拷木锟斤拷锟揭?拷锟?1锟?7 锟斤拷值为锟斤拷normal锟斤拷 锟酵★拷high锟斤拷 锟斤拷锟斤拷锟斤拷'N'锟酵★拷H'锟斤拷示一锟姐精锟饺和高撅拷锟斤拷 锟斤拷默锟斤拷为锟竭撅拷锟饺硷拷锟斤拷</param>
/// <param name="epthilongA_deg">ƽ<>Ƴཻ<C6B3><E0BDBB></param> /// <param name="epthilongA_deg">ƽ<>Ƴཻ<C6B3><E0BDBB></param>
/// <param name="dertaPthi_deg"><3E>ƾ<EFBFBD><C6BE><EFBFBD></param> /// <param name="dertaPthi_deg"><3E>ƾ<EFBFBD><C6BE><EFBFBD></param>
/// <param name="dertaEpthilong_deg"><3E>ͽ<EFBFBD><CDBD><EFBFBD><EFBFBD><EFBFBD></param> /// <param name="dertaEpthilong_deg"><3E>ͽ<EFBFBD><CDBD><EFBFBD><EFBFBD><EFBFBD></param>
@ -1266,7 +1264,7 @@ int WGS84_J2000::nutationInLongitudeCaculate(double JD, double& epthilongA_deg,
} }
/// <summary> /// <summary>
/// zA<EFBFBD><EFBFBD>thitaA<EFBFBD><EFBFBD>zetaAΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7, <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> /// zA锟斤拷thitaA锟斤拷zetaA为锟斤拷锟斤拷锟斤拷锟?1锟?7, 锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷
/// </summary> /// </summary>
/// <param name="JDTDB"></param> /// <param name="JDTDB"></param>
/// <param name="zetaA"></param> /// <param name="zetaA"></param>
@ -1318,13 +1316,13 @@ double WGS84_J2000::YMD2JD(double y, double m, double d)
/// <returns></returns> /// <returns></returns>
Eigen::MatrixXd WGS84_J2000::WGS842J2000(Eigen::MatrixXd BLH_deg_m, Eigen::MatrixXd UTC, double xp, double yp, double dut1, double dat) Eigen::MatrixXd WGS84_J2000::WGS842J2000(Eigen::MatrixXd BLH_deg_m, Eigen::MatrixXd UTC, double xp, double yp, double dut1, double dat)
{ {
//step 1 step1 WGS 84 ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Э<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD>ￄ1<EFBFBD>7 //step 1 step1 WGS 84 转锟斤拷锟斤拷协锟斤拷锟斤拷锟斤拷锟斤拷锟较碉拷锟?1锟?7
Eigen::MatrixXd earthFixedXYZ = WGS84_J2000::WGS84TECEF(BLH_deg_m).transpose(); Eigen::MatrixXd earthFixedXYZ = WGS84_J2000::WGS84TECEF(BLH_deg_m).transpose();
//step 2 Э<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>τ1<EFBFBD>7 ת<><D7AA>Ϊ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> //step 2 协锟斤拷锟斤拷锟斤拷锟斤拷锟较?1锟?7 转锟斤拷为瞬时锟斤拷锟斤拷锟斤拷锟斤拷
// <20><><EFBFBD><EFBFBD>Ҫ<EFBFBD>漰<EFBFBD><E6BCB0>ѯ <20><><EFBFBD><EFBFBD>ʱ<EFBFBD>̵<EFBFBD> xp,yp , <20><><EFBFBD>Բ<EFBFBD>ѯEOP<4F>ļ<EFBFBD><C4BC><EFBFBD>ã<EFBFBD><C3A3><EFBFBD><EFBFBD>ص<EFBFBD>ַ<EFBFBD><D6B7><EFBFBD><EFBFBD>http://celestrak.com/spacedata/ // <20><><EFBFBD><EFBFBD>Ҫ<EFBFBD>漰<EFBFBD><E6BCB0>ѯ <20><><EFBFBD><EFBFBD>ʱ<EFBFBD>̵<EFBFBD> xp,yp , <20><><EFBFBD>Բ<EFBFBD>ѯEOP<4F>ļ<EFBFBD><C4BC><EFBFBD>ã<EFBFBD><C3A3><EFBFBD><EFBFBD>ص<EFBFBD>ַ<EFBFBD><D6B7><EFBFBD><EFBFBD>http://celestrak.com/spacedata/
earthFixedXYZ = ordinateSingleRotate(1, yp) * ordinateSingleRotate(2, xp) * earthFixedXYZ; earthFixedXYZ = ordinateSingleRotate(1, yp) * ordinateSingleRotate(2, xp) * earthFixedXYZ;
//step 3 ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ ת<><D7AA>Ϊ ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ //step 3 ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ ת<><D7AA>Ϊ ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ
// <EFBFBD>ò<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD>漰<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>κ<EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD>ǵļ<EFBFBD><EFBFBD>㣬<EFBFBD><EFBFBD><EFBFBD>ڸ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>κ<EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD> <20><><EFBFBD>㷽<EFBFBD><E3B7BD> <20>ܶ࣬<DCB6><E0A3AC><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBF84>7 һ<><D2BB><EFBFBD><EFBFBD>Ϊ<EFBFBD><CEAA>ȷ<EFBFBD>ļ<EFBFBD><C4BC>㷽<EFBFBD><E3B7BD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>dut1 <20><>dat<61><74><EFBFBD><EFBFBD><EFBFBD><EFBFBD>EOP<4F>ļ<EFBFBD><C4BC><EFBFBD><EFBFBD>С<EFBFBD>EOP<4F>ļ<EFBFBD><C4BC><EFBFBD><EFBFBD>ص<EFBFBD>ַ<EFBFBD><D6B7><EFBFBD><EFBFBD> http://celestrak.org/spacedata/ // 锟矫诧拷锟斤拷锟斤拷要锟芥及锟斤拷锟斤拷锟斤拷锟斤拷锟轿猴拷锟斤拷时锟角的硷拷锟姐，锟斤拷锟节革拷锟斤拷锟斤拷锟轿猴拷锟斤拷时锟斤拷 锟斤拷锟姐方锟斤拷 锟杰多，锟斤拷锟斤拷锟斤拷锟?1锟?7 一锟斤拷锟斤拷为锟斤拷确锟侥硷拷锟姐方锟斤拷锟斤拷锟斤拷锟斤拷dut1 锟斤拷dat锟斤拷锟斤拷锟斤拷EOP锟侥硷拷锟斤拷锟叫★拷EOP锟侥硷拷锟斤拷锟截碉拷址锟斤拷锟斤拷 http://celestrak.org/spacedata/
double gst_deg, JDTDB; double gst_deg, JDTDB;
WGS84_J2000::utc2gst(UTC, dut1, dat, gst_deg, JDTDB); WGS84_J2000::utc2gst(UTC, dut1, dat, gst_deg, JDTDB);
Eigen::MatrixXd xyz = ordinateSingleRotate(3, -1 * gst_deg) * earthFixedXYZ; Eigen::MatrixXd xyz = ordinateSingleRotate(3, -1 * gst_deg) * earthFixedXYZ;

153
baseTool.h
View File

@ -1,6 +1,6 @@
#pragma once #pragma once
/// ///
/// 基本类、基本函数 /// 基本类、基本函数
/// ///
//#define EIGEN_USE_MKL_ALL //#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2 //#define EIGEN_VECTORIZE_SSE4_2
@ -27,54 +27,61 @@
using namespace std; using namespace std;
using namespace Eigen; using namespace Eigen;
#define PI_180 180/3.141592653589793238462643383279; #define PI_180 180/3.141592653589793238462643383279
#define T180_PI 3.141592653589793238462643383279/180; #define T180_PI 3.141592653589793238462643383279/180
#define Radians2Degrees(Radians) Radians*PI_180 #define Radians2Degrees(Radians) Radians*PI_180
#define Degrees2Radians(Degrees) Degrees*T180_PI #define Degrees2Radians(Degrees) Degrees*T180_PI
const double PI = 3.141592653589793238462643383279; const double PI = 3.141592653589793238462643383279;
const double epsilon = 0.000000000000001; const double epsilon = 0.000000000000001;
const double pi = 3.14159265358979323846; const double pi = 3.14159265358979323846;
const double d2r = pi / 180; const double d2r = pi / 180;
const double r2d = 180 / pi; const double r2d = 180 / pi;
const double a = 6378137.0; //椭球长半轴 const double a = 6378137.0; //椭球长半轴
const double ae = 6378137.0; //椭球长半轴 const double ae = 6378137.0; //椭球长半轴
const double ee= 0.0818191910428;// 第一偏心率 const double ee= 0.0818191910428;// 第一偏心率
const double f_inverse = 298.257223563; //扁率倒数 const double f_inverse = 298.257223563; //扁率倒数
const double b = a - a / f_inverse; const double b = a - a / f_inverse;
const double eSquare = (a * a - b * b) / (a * a); const double eSquare = (a * a - b * b) / (a * a);
const double e = sqrt(eSquare); const double e = sqrt(eSquare);
const double earth_Re = 6378136.49; const double earth_Re = 6378136.49;
const double earth_Rp = (1 - 1 / f_inverse) * earth_Re; const double earth_Rp = (1 - 1 / f_inverse) * earth_Re;
const double earth_We = 0.000072292115; const double earth_We = 0.000072292115;
///////////////////////////////////// 运行时间打印 ////////////////////////////////////////////////////////// ///////////////////////////////////// 运行时间打印 //////////////////////////////////////////////////////////
std::string getCurrentTimeString(); std::string getCurrentTimeString();
std::string getCurrentShortTimeString(); std::string getCurrentShortTimeString();
/////////////////////////////// 基本图像类 ////////////////////////////////////////////////////////// /////////////////////////////// 基本图像类 //////////////////////////////////////////////////////////
/// <summary> /// <summary>
/// 三维向量，坐标表达 /// 三维向量，坐标表达
/// </summary> /// </summary>
struct Landpoint // 点 SAR影像的像素坐标； struct Landpoint // 点 SAR影像的像素坐标；
{ {
/// <summary> /// <summary>
/// 经度x /// 经度x
/// </summary> /// </summary>
double lon; // 经度x lon pixel_col double lon; // 经度x lon pixel_col
/// <summary> /// <summary>
/// 纬度y /// 纬度y
/// </summary> /// </summary>
double lat; // 纬度y lat pixel_row double lat; // 纬度y lat pixel_row
/// <summary> /// <summary>
/// 高度z /// 高度z
/// </summary> /// </summary>
double ati; // 高程z ati pixel_time double ati; // 高程z ati pixel_time
}; };
struct Point_3d { struct Point_3d {
double x; double x;
@ -91,15 +98,15 @@ bool operator ==(const Landpoint& p1, const Landpoint& p2);
Landpoint operator *(const Landpoint& p, double scale); Landpoint operator *(const Landpoint& p, double scale);
/// <summary> /// <summary>
/// 向量A,B的夹角,角度 /// 向量A,B的夹角,角度
/// </summary> /// </summary>
/// <param name="a"></param> /// <param name="a"></param>
/// <param name="b"></param> /// <param name="b"></param>
/// <returns>角度制 0-360度，逆时针</returns> /// <returns>角度制 0-360度，逆时针</returns>
double getAngle(const Landpoint& a, const Landpoint& b); double getAngle(const Landpoint& a, const Landpoint& b);
/// <summary> /// <summary>
/// 点乘 /// 点乘
/// </summary> /// </summary>
/// <param name="p1"></param> /// <param name="p1"></param>
/// <param name="p2"></param> /// <param name="p2"></param>
@ -113,20 +120,20 @@ Landpoint crossProduct(const Landpoint& a, const Landpoint& b);
struct DemBox { struct DemBox {
double min_lat; //纬度 double min_lat; //纬度
double min_lon;//经度 double min_lon;//经度
double max_lat;//纬度 double max_lat;//纬度
double max_lon;//经度 double max_lon;//经度
}; };
/// <summary> /// <summary>
/// gdalImage图像操作类 /// gdalImage图像操作类
/// </summary> /// </summary>
class gdalImage class gdalImage
{ {
public: // 方法 public: // 方法
gdalImage(string raster_path); gdalImage(string raster_path);
~gdalImage(); ~gdalImage();
void setHeight(int); void setHeight(int);
@ -148,15 +155,15 @@ public: // 方法
Eigen::MatrixXd getHist(int bandids); Eigen::MatrixXd getHist(int bandids);
public: public:
string img_path; // 图像文件 string img_path; // 图像文件
int height; // 高 int height; // 高
int width; // 宽 int width; // 宽
int band_num;// 波段数 int band_num;// 波段数
int start_row;// int start_row;//
int start_col;// int start_col;//
int data_band_ids; int data_band_ids;
Eigen::MatrixXd gt; // 变换矩阵 Eigen::MatrixXd gt; // 变换矩阵
Eigen::MatrixXd inv_gt; // 逆变换矩阵 Eigen::MatrixXd inv_gt; // 逆变换矩阵
Eigen::MatrixXd data; Eigen::MatrixXd data;
string projection; string projection;
}; };
@ -171,45 +178,45 @@ int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResa
/////////////////////////////// 基本图像类 结束 ////////////////////////////////////////////////////////// /////////////////////////////// 基本图像类 结束 //////////////////////////////////////////////////////////
string Convert(float Num); string Convert(float Num);
std::string JoinPath(const std::string& path, const std::string& filename); std::string JoinPath(const std::string& path, const std::string& filename);
////////////////////////////// 坐标部分基本方法 ////////////////////////////////////////// ////////////////////////////// 坐标部分基本方法 //////////////////////////////////////////
/// <summary> /// <summary>
/// 将经纬度转换为地固参心坐标系 /// 将经纬度转换为地固参心坐标系
/// </summary> /// </summary>
/// <param name="XYZP">经纬度点--degree</param> /// <param name="XYZP">经纬度点--degree</param>
/// <returns>投影坐标系点</returns> /// <returns>投影坐标系点</returns>
Landpoint LLA2XYZ(const Landpoint& LLA); Landpoint LLA2XYZ(const Landpoint& LLA);
Eigen::MatrixXd LLA2XYZ(Eigen::MatrixXd landpoint); Eigen::MatrixXd LLA2XYZ(Eigen::MatrixXd landpoint);
/// <summary> /// <summary>
/// 将地固参心坐标系转换为经纬度 /// 将地固参心坐标系转换为经纬度
/// </summary> /// </summary>
/// <param name="XYZ">固参心坐标系</param> /// <param name="XYZ">固参心坐标系</param>
/// <returns>经纬度--degree</returns> /// <returns>经纬度--degree</returns>
Landpoint XYZ2LLA(const Landpoint& XYZ); Landpoint XYZ2LLA(const Landpoint& XYZ);
////////////////////////////// 坐标部分基本方法 ////////////////////////////////////////// ////////////////////////////// 坐标部分基本方法 //////////////////////////////////////////
/// <summary> /// <summary>
/// 计算地表坡度向量 /// 计算地表坡度向量
/// </summary> /// </summary>
/// <param name="p0">固参心坐标系</param> /// <param name="p0">固参心坐标系</param>
/// <param name="p1">固参心坐标系</param> /// <param name="p1">固参心坐标系</param>
/// <param name="p2">固参心坐标系</param> /// <param name="p2">固参心坐标系</param>
/// <param name="p3">固参心坐标系</param> /// <param name="p3">固参心坐标系</param>
/// <param name="p4">固参心坐标系</param> /// <param name="p4">固参心坐标系</param>
/// <returns>向量角度</returns> /// <returns>向量角度</returns>
//Landpoint getSlopeVector(Landpoint& p0, Landpoint& p1, Landpoint& p2, Landpoint& p3, Landpoint& p4, Eigen::MatrixXd UTC, double xp, double yp, double dut1, double dat); //Landpoint getSlopeVector(Landpoint& p0, Landpoint& p1, Landpoint& p2, Landpoint& p3, Landpoint& p4, Eigen::MatrixXd UTC, double xp, double yp, double dut1, double dat);
Landpoint getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landpoint& p2, const Landpoint& p3, const Landpoint& p4); Landpoint getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landpoint& p2, const Landpoint& p3, const Landpoint& p4);
////////////////////////////// 插值 //////////////////////////////////////////// ////////////////////////////// 插值 ////////////////////////////////////////////
complex<double> Cubic_Convolution_interpolation(double u,double v,Eigen::MatrixX<complex<double>> img); complex<double> Cubic_Convolution_interpolation(double u,double v,Eigen::MatrixX<complex<double>> img);
@ -234,8 +241,8 @@ inline double operator/(Point_3d a, Point_3d b) {
}; };
inline bool onSegment(Point_3d Pi, Point_3d Pj, Point_3d Q) inline bool onSegment(Point_3d Pi, Point_3d Pj, Point_3d Q)
{ {
if ((Q.x - Pi.x) * (Pj.y - Pi.y) == (Pj.x - Pi.x) * (Q.y - Pi.y) //叉乘 if ((Q.x - Pi.x) * (Pj.y - Pi.y) == (Pj.x - Pi.x) * (Q.y - Pi.y) //叉乘
//保证Q点坐标在pi,pj之间 //保证Q点坐标在pi,pj之间
&& min(Pi.x, Pj.x) <= Q.x && Q.x <= max(Pi.x, Pj.x) && min(Pi.x, Pj.x) <= Q.x && Q.x <= max(Pi.x, Pj.x)
&& min(Pi.y, Pj.y) <= Q.y && Q.y <= max(Pi.y, Pj.y)) && min(Pi.y, Pj.y) <= Q.y && Q.y <= max(Pi.y, Pj.y))
return true; return true;
@ -269,7 +276,7 @@ inline Point_3d invBilinear(Point_3d p, Point_3d a, Point_3d b, Point_3d c, Poin
{ {
float w = k1 * k1 - 4.0 * k0 * k2; float w = k1 * k1 - 4.0 * k0 * k2;
if (w < 0.0){ if (w < 0.0){
// 可能在边界上 // 可能在边界上
if (onSegment(a, b, p)) { if (onSegment(a, b, p)) {
Point_3d tt = b - a; Point_3d tt = b - a;
Point_3d ttpa = p - a; Point_3d ttpa = p - a;
@ -328,21 +335,21 @@ inline Point_3d invBilinear(Point_3d p, Point_3d a, Point_3d b, Point_3d c, Poin
// //
// WGS84 到J2000 坐标系的变换 // WGS84 到J2000 坐标系的变换
// 参考网址：https://blog.csdn.net/hit5067/article/details/116894616 // 参考网址：https://blog.csdn.net/hit5067/article/details/116894616
// 资料网址：http://celestrak.org/spacedata/ // 资料网址：http://celestrak.org/spacedata/
// 参数文件： // 参数文件：
// a. Earth Orientation Parameter 文件： http://celestrak.org/spacedata/EOP-Last5Years.csv // a. Earth Orientation Parameter 文件： http://celestrak.org/spacedata/EOP-Last5Years.csv
// b. Space Weather Data 文件： http://celestrak.org/spacedata/SW-Last5Years.csv // b. Space Weather Data 文件： http://celestrak.org/spacedata/SW-Last5Years.csv
// 备注：上述文件是自2017年-五年内 // 备注：上述文件是自2017年-五年内
/** /**
在wgs84 坐标系转到J2000 坐标系 主要 涉及到坐标的相互转换。一般给定的WGS坐标为 给定时刻的 t ,BLH 在wgs84 坐标系转到J2000 坐标系 主要 涉及到坐标的相互转换。一般给定的WGS坐标为 给定时刻的 t ,BLH
转换步骤： 转换步骤：
step 1: WGS 84 转换到协议地球坐标系 step 1: WGS 84 转换到协议地球坐标系
step 2: 协议地球坐标系 转换为瞬时地球坐标系 step 2: 协议地球坐标系 转换为瞬时地球坐标系
step 3: 瞬时地球坐标系 转换为 瞬时真天球坐标系 step 3: 瞬时地球坐标系 转换为 瞬时真天球坐标系
step 4: 瞬时真天球坐标系 转到瞬时平天球 坐标系 step 4: 瞬时真天球坐标系 转到瞬时平天球 坐标系
step 5: 瞬时平天球坐标系转换为协议天球坐标系（J2000） step 5: 瞬时平天球坐标系转换为协议天球坐标系（J2000）
**/ **/
@ -362,26 +369,26 @@ public:
~WGS84_J2000(); ~WGS84_J2000();
public: public:
// step1 WGS 84 转换到协议地球坐标系。 // step1 WGS 84 转换到协议地球坐标系。
static Eigen::MatrixXd WGS84TECEF(Eigen::MatrixXd WGS84_Lon_lat_ait); static Eigen::MatrixXd WGS84TECEF(Eigen::MatrixXd WGS84_Lon_lat_ait);
//step 2 协议地球坐标系 转换为瞬时地球坐标系 //step 2 协议地球坐标系 转换为瞬时地球坐标系
static Eigen::MatrixXd ordinateSingleRotate(int axis, double angle_deg); static Eigen::MatrixXd ordinateSingleRotate(int axis, double angle_deg);
// step 3 瞬时地球坐标系 转换为 瞬时真天球坐标系 // step 3 瞬时地球坐标系 转换为 瞬时真天球坐标系
// xyz= ordinateSingleRotate('z',-gst_deg)*earthFixedXYZ; // xyz= ordinateSingleRotate('z',-gst_deg)*earthFixedXYZ;
static int utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& gst_deg, double& JDTDB); static int utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& gst_deg, double& JDTDB);
// step 4 瞬时真天球坐标系 转到瞬时平天球 坐标系 // step 4 瞬时真天球坐标系 转到瞬时平天球 坐标系
static int nutationInLongitudeCaculate(double JD, double& epthilongA_deg, double& dertaPthi_deg, double& dertaEpthilong_deg, double& epthilong_deg); static int nutationInLongitudeCaculate(double JD, double& epthilongA_deg, double& dertaPthi_deg, double& dertaEpthilong_deg, double& epthilong_deg);
// step5 瞬时平天球坐标系转换为协议天球坐标系（J2000）函数中 JDTDB 为给定时刻 的地球动力学时对应的儒略日，其计算方法由步骤三中的函数给出。 // step5 瞬时平天球坐标系转换为协议天球坐标系（J2000）函数中 JDTDB 为给定时刻 的地球动力学时对应的儒略日，其计算方法由步骤三中的函数给出。
// xyz=ordinateSingleRotate('Z',zetaA)*ordinateSingleRotate('y',-thitaA)*ordinateSingleRotate('z',zA)*xyz; // xyz=ordinateSingleRotate('Z',zetaA)*ordinateSingleRotate('y',-thitaA)*ordinateSingleRotate('z',zA)*xyz;
static int precessionAngle(double JDTDB, double& zetaA, double& thitaA, double& zA); static int precessionAngle(double JDTDB, double& zetaA, double& thitaA, double& zA);
// YMD2JD 同时 YMD2JD函数为 年月日转换为儒略日，具体说明 见公元纪年法（儒略历-格里高历）转儒略日_${王小贱}的博客-CSDN博客_年积日计算公式 // YMD2JD 同时 YMD2JD函数为 年月日转换为儒略日，具体说明 见公元纪年法（儒略历-格里高历）转儒略日_${王小贱}的博客-CSDN博客_年积日计算公式
static double YMD2JD(double y, double m, double d); static double YMD2JD(double y, double m, double d);
static Eigen::MatrixXd WGS842J2000(Eigen::MatrixXd BLH_deg_m, Eigen::MatrixXd UTC, double xp, double yp, double dut1, double dat); static Eigen::MatrixXd WGS842J2000(Eigen::MatrixXd BLH_deg_m, Eigen::MatrixXd UTC, double xp, double yp, double dut1, double dat);
static Landpoint WGS842J2000(Landpoint LBH_deg_m, Eigen::MatrixXd UTC, double xp, double yp, double dut1, double dat); static Landpoint WGS842J2000(Landpoint LBH_deg_m, Eigen::MatrixXd UTC, double xp, double yp, double dut1, double dat);
public: public:
static std::string EOP_File_Path; static std::string EOP_File_Path;
static std::string Space_Weather_Data; static std::string Space_Weather_Data;
// IAU2000模型有77项11列 // IAU2000模型有77项11列
static Eigen::Matrix<double, 77, 11> IAU2000ModelParams; static Eigen::Matrix<double, 77, 11> IAU2000ModelParams;
}; };

2
interpolation.h
View File

@ -1,6 +1,6 @@
#pragma once #pragma once
/** /**
专门用于插值计算的类库 专门用于插值计算的类库
*/ */
#include <stdio.h> #include <stdio.h>

4
resource.h
View File

@ -1,6 +1,6 @@
//{{NO_DEPENDENCIES}} //{{NO_DEPENDENCIES}}
// Microsoft Visual C++ 生成的包含文件。 // Microsoft Visual C++ 生成的包含文件。
// 供 SIMOrthoProgram.rc 使用 // 供 SIMOrthoProgram.rc 使用
// //
#define IDR_PROJ1 101 #define IDR_PROJ1 101

698
simptsn.cpp
View File

File diff suppressed because it is too large Load Diff

140
simptsn.h
View File

@ -1,11 +1,11 @@
#pragma once #pragma once
/// ///
/// 仿真成像算法 /// 仿真成像算法
/// ///
//#define EIGEN_USE_MKL_ALL //#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2 //#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h> //#include <mkl.h>
// 本地方法 // 本地方法
#include "baseTool.h" #include "baseTool.h"
#include <iostream> #include <iostream>
#include <Eigen/Core> #include <Eigen/Core>
@ -18,11 +18,11 @@
using namespace std; using namespace std;
using namespace Eigen; using namespace Eigen;
///////////// ptsn 算法 ///////////////////// ///////////// ptsn 算法 /////////////////////
class PSTNAlgorithm { class PSTNAlgorithm {
public: public:
// 初始化与析构函数 // 初始化与析构函数
PSTNAlgorithm(); PSTNAlgorithm();
PSTNAlgorithm(string parameterPath); PSTNAlgorithm(string parameterPath);
~PSTNAlgorithm(); ~PSTNAlgorithm();
@ -30,53 +30,53 @@ public:
int dem2lat_lon(std::string dem_path, std::string lon_path, std::string lat_path); 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); int dem2SAR_RC(std::string dem_path, std::string sim_rc_path);
// 模拟算法函数 // 模拟算法函数
Eigen::MatrixX<double> calNumericalDopplerValue(Eigen::MatrixX<double> R);// 数值模拟法计算多普勒频移值 Eigen::MatrixX<double> calNumericalDopplerValue(Eigen::MatrixX<double> R);// 数值模拟法计算多普勒频移值
double calNumericalDopplerValue(double R);// 数值模拟法计算多普勒频移值 double calNumericalDopplerValue(double R);// 数值模拟法计算多普勒频移值
Eigen::MatrixX<double> calTheoryDopplerValue(Eigen::MatrixX<double> R, Eigen::MatrixX<double> R_s1, Eigen::MatrixX<double> V_s1);//根据理论模型计算多普勒频移值 Eigen::MatrixX<double> calTheoryDopplerValue(Eigen::MatrixX<double> R, Eigen::MatrixX<double> R_s1, Eigen::MatrixX<double> V_s1);//根据理论模型计算多普勒频移值
double calTheoryDopplerValue(double R, Eigen::MatrixX<double> R_s1, Eigen::MatrixX<double> V_s1); double calTheoryDopplerValue(double R, Eigen::MatrixX<double> R_s1, Eigen::MatrixX<double> V_s1);
Eigen::MatrixX<double> PSTN(Eigen::MatrixX<double> TargetPostion, double ave_dem, double dt, bool isCol ); // 输入DEM Eigen::MatrixX<double> PSTN(Eigen::MatrixX<double> TargetPostion, double ave_dem, double dt, bool isCol ); // 输入DEM
//Eigen::MatrixXd WGS842J2000(Eigen::MatrixXd blh); //Eigen::MatrixXd WGS842J2000(Eigen::MatrixXd blh);
//Landpoint WGS842J2000(Landpoint p); //Landpoint WGS842J2000(Landpoint p);
public: // WGS84 到 J2000 之间的变换参数 public: // WGS84 到 J2000 之间的变换参数
Eigen::MatrixXd UTC; Eigen::MatrixXd UTC;
double Xp = 0.204071; double Xp = 0.204071;
double Yp = 0.318663; double Yp = 0.318663;
double Dut1 = 0.0366742; double Dut1 = 0.0366742;
double Dat = 37; double Dat = 37;
public: public:
// std::string dem_path; // 模拟用的DEM影像 // std::string dem_path; // 模拟用的DEM影像
// std::string lon_path;// 经度影像 // std::string lon_path;// 经度影像
// std::string lat_path;// 纬度影像 // std::string lat_path;// 纬度影像
// std::string simsar_path;// 仿真影像坐标 // std::string simsar_path;// 仿真影像坐标
int height; // 影像的高 int height; // 影像的高
int width; int width;
// 入射角 // 入射角
double near_in_angle;// 近入射角 double near_in_angle;// 近入射角
double far_in_angle;// 远入射角 double far_in_angle;// 远入射角
// SAR的成像参数 // SAR的成像参数
double widthspace;// 距离向分辨率 double widthspace;// 距离向分辨率
double R0;//近斜距 double R0;//近斜距
double LightSpeed; // 光速 double LightSpeed; // 光速
double lamda;// 波长 double lamda;// 波长
double refrange;// 参考斜距 double refrange;// 参考斜距
double imgStartTime;// 成像起始时间 double imgStartTime;// 成像起始时间
double PRF;// 脉冲重复率 double PRF;// 脉冲重复率
int doppler_paramenter_number;//多普勒系数个数 int doppler_paramenter_number;//多普勒系数个数
MatrixX<double> doppler_paras;// 多普勒系数 MatrixX<double> doppler_paras;// 多普勒系数
OrbitPoly orbit; // 轨道模型 OrbitPoly orbit; // 轨道模型
}; };
/// <summary> /// <summary>
/// 获取局地入射角，角度值 /// 获取局地入射角，角度值
/// </summary> /// </summary>
/// <param name="satepoint"></param> /// <param name="satepoint"></param>
/// <param name="landpoint"></param> /// <param name="landpoint"></param>
@ -86,7 +86,7 @@ double getlocalIncAngle(Landpoint& satepoint, Landpoint& landpoint, Landpoint& s
/// <summary> /// <summary>
/// 侧视入射角,角度值 /// 侧视入射角,角度值
/// </summary> /// </summary>
/// <param name="satepoint"></param> /// <param name="satepoint"></param>
/// <param name="landpoint"></param> /// <param name="landpoint"></param>
@ -95,12 +95,12 @@ double getIncAngle(Landpoint& satepoint, Landpoint& landpoint);
/// <summary> /// <summary>
/// ASF计算方法 /// ASF计算方法
/// </summary> /// </summary>
class ASFOrthClass { class ASFOrthClass {
public: public:
Eigen::MatrixXd Satellite2ECR(Eigen::Vector3d Rsc, Eigen::Vector3d Vsc); // 根据卫星坐标计算变换矩阵 M Eigen::MatrixXd Satellite2ECR(Eigen::Vector3d Rsc, Eigen::Vector3d Vsc); // 根据卫星坐标计算变换矩阵 M
Eigen::Vector3d UnitVectorOfSatelliteAndTarget(double alpha, double beta, Eigen::MatrixXd M); Eigen::Vector3d UnitVectorOfSatelliteAndTarget(double alpha, double beta, Eigen::MatrixXd M);
double GetLookFromRangeYaw(double R, double alpha, double beta, Eigen::Vector3d SatellitePosition, Eigen::Vector3d SatelliteVelocity, double R_threshold, double H = 0); double GetLookFromRangeYaw(double R, double alpha, double beta, Eigen::Vector3d SatellitePosition, Eigen::Vector3d SatelliteVelocity, double R_threshold, double H = 0);
@ -114,7 +114,7 @@ public:
}; };
/// <summary> /// <summary>
/// 仿真处理流程 /// 仿真处理流程
/// </summary> /// </summary>
class simProcess { class simProcess {
public: public:
@ -123,10 +123,10 @@ public:
int parameters(int argc, char* argv[]); int parameters(int argc, char* argv[]);
int InitSimulationSAR(std::string paras_path,std::string workspace_path,std::string out_dir_path,std::string in_dem_path,std::string in_sar_path); // 间接定位法 int InitSimulationSAR(std::string paras_path,std::string workspace_path,std::string out_dir_path,std::string in_dem_path,std::string in_sar_path); // 间接定位法
int CreateSARDEM(); int CreateSARDEM();
int dem2SAR(); // 切换行号 int dem2SAR(); // 切换行号
int SARIncidentAngle(); int SARIncidentAngle();
int SARSimulationWGS(); int SARSimulationWGS();
int SARSimulation(); int SARSimulation();
@ -134,30 +134,30 @@ public:
int ReflectTable_WGS2Range(); int ReflectTable_WGS2Range();
int InitRPCIncAngle(std::string paras_path, std::string workspace_path, std::string out_dir_path, std::string in_dem_path,std::string in_rpc_lon_lat_path,std::string out_inc_angle_path, std::string out_local_inc_angle_path, std::string out_inc_angle_geo_path, std::string out_local_inc_angle_geo_path); int InitRPCIncAngle(std::string paras_path, std::string workspace_path, std::string out_dir_path, std::string in_dem_path,std::string in_rpc_lon_lat_path,std::string out_inc_angle_path, std::string out_local_inc_angle_path, std::string out_inc_angle_geo_path, std::string out_local_inc_angle_geo_path);
int dem2SAR_row(); // 切换行号 int dem2SAR_row(); // 切换行号
int SARIncidentAngle_RPC(); int SARIncidentAngle_RPC();
int createRPCDEM(); int createRPCDEM();
// 格网离散插值 // 格网离散插值
int Scatter2Grid(std::string lon_lat_path, std::string data_tiff, std::string grid_path, double space); int Scatter2Grid(std::string lon_lat_path, std::string data_tiff, std::string grid_path, double space);
// ASF 模块计算 // ASF 模块计算
int InitASFSAR(std::string paras_path, std::string workspace_path, std::string out_dir_path, std::string in_dem_path); int InitASFSAR(std::string paras_path, std::string workspace_path, std::string out_dir_path, std::string in_dem_path);
int ASF(std::string in_dem_path, std::string out_latlon_path, ASFOrthClass asfclass, PSTNAlgorithm PSTN); int ASF(std::string in_dem_path, std::string out_latlon_path, ASFOrthClass asfclass, PSTNAlgorithm PSTN);
//int dem2simSAR(std::string dem_path, std::string out_sim_path, PSTNAlgorithm PSTN); // 根据DEM生成输出仿真图像 //int dem2simSAR(std::string dem_path, std::string out_sim_path, PSTNAlgorithm PSTN); // 根据DEM生成输出仿真图像
int sim_SAR(std::string dem_path, std::string sim_rc_path, std::string sim_sar_path, PSTNAlgorithm PSTN); int sim_SAR(std::string dem_path, std::string sim_rc_path, std::string sim_sar_path, PSTNAlgorithm PSTN);
int dem2aspect_slope(std::string dem_path, std::string aspect_path, std::string slope_path); int dem2aspect_slope(std::string dem_path, std::string aspect_path, std::string slope_path);
int sim_sum_SAR(std::string sim_dem_path,std::string sim_rc_path, std::string sim_orth_path, std::string sim_sum_path,PSTNAlgorithm pstn); int sim_sum_SAR(std::string sim_dem_path,std::string sim_rc_path, std::string sim_orth_path, std::string sim_sum_path,PSTNAlgorithm pstn);
void interpolation_GTC_sar_sigma(std::string in_rc_wgs84_path, std::string in_ori_sar_path, std::string out_orth_sar_path, PSTNAlgorithm pstn); void interpolation_GTC_sar_sigma(std::string in_rc_wgs84_path, std::string in_ori_sar_path, std::string out_orth_sar_path, PSTNAlgorithm pstn);
// 原始影像强度图 // 原始影像强度图
int ori_sar_power(std::string ori_sar_path, std::string out_sar_power_path); int ori_sar_power(std::string ori_sar_path, std::string out_sar_power_path);
// 匹配偏移模型构建 // 匹配偏移模型构建
int createImageMatchModel(std::string ori_sar_rsc_path, std::string ori_sar_rsc_jpg_path, std::string sim_sar_path, std::string sim_sar_jpg_path, std::string matchmodel_path); int createImageMatchModel(std::string ori_sar_rsc_path, std::string ori_sar_rsc_jpg_path, std::string sim_sar_path, std::string sim_sar_jpg_path, std::string matchmodel_path);
// 纠正影像 // 纠正影像
int correctOrth_rc(std::string sim_rc_path,ImageMatch matchmodel); int correctOrth_rc(std::string sim_rc_path,ImageMatch matchmodel);
int correct_ati(std::string orth_rc_path, std::string dem_path, std::string out_inclocal_path, std::string out_dem_path, std::string out_count_path, PSTNAlgorithm PSTN); int correct_ati(std::string orth_rc_path, std::string dem_path, std::string out_inclocal_path, std::string out_dem_path, std::string out_count_path, PSTNAlgorithm PSTN);
@ -175,58 +175,58 @@ public:
void CorrectionFromSLC2Geo(string in_lon_lat_path, string in_radar_path, string out_radar_path, int in_band_ids); void CorrectionFromSLC2Geo(string in_lon_lat_path, string in_radar_path, string out_radar_path, int in_band_ids);
// 内部参数 // 内部参数
public: public:
ImageMatch matchmodel; ImageMatch matchmodel;
bool isMatchModel; bool isMatchModel;
std::string workSpace_path;// 工作空间路径 std::string workSpace_path;// 工作空间路径
std::string outSpace_path;// 输出工作空间路径 std::string outSpace_path;// 输出工作空间路径
PSTNAlgorithm pstn; // 参数组件 PSTNAlgorithm pstn; // 参数组件
///// RPC 局地入射角 /////////////////////////////////// ///// RPC 局地入射角 ///////////////////////////////////
std::string in_rpc_lon_lat_path; std::string in_rpc_lon_lat_path;
std::string out_dir_path; std::string out_dir_path;
std::string workspace_path; std::string workspace_path;
// 临时文件 // 临时文件
std::string rpc_wgs_path; std::string rpc_wgs_path;
std::string ori_sim_count_tiff;// 映射角文件 std::string ori_sim_count_tiff;// 映射角文件
// 输出文件 // 输出文件
std::string out_inc_angle_rpc_path; std::string out_inc_angle_rpc_path;
std::string out_local_inc_angle_rpc_path; std::string out_local_inc_angle_rpc_path;
// 输出文件 // 输出文件
std::string out_inc_angle_geo_path; std::string out_inc_angle_geo_path;
std::string out_local_inc_angle_geo_path; std::string out_local_inc_angle_geo_path;
///// 正向仿真方法 ///////// ///// 正向仿真方法 /////////
// 临时产品文件 // 临时产品文件
std::string in_dem_path; // 输入DEM std::string in_dem_path; // 输入DEM
std::string dem_path; // 重采样之后DEM std::string dem_path; // 重采样之后DEM
std::string dem_r_path;// 行号影像 std::string dem_r_path;// 行号影像
std::string dem_rc_path; // 行列号 std::string dem_rc_path; // 行列号
std::string sar_sim_wgs_path; std::string sar_sim_wgs_path;
std::string sar_sim_path; std::string sar_sim_path;
std::string in_sar_path; std::string in_sar_path;
std::string sar_jpg_path; std::string sar_jpg_path;
std::string sar_power_path; std::string sar_power_path;
// 最终产品 // 最终产品
std::string out_dem_rc_path; // 经纬度与行列号变换文件 std::string out_dem_rc_path; // 经纬度与行列号变换文件
std::string out_dem_slantRange_path; // 斜距文件 std::string out_dem_slantRange_path; // 斜距文件
std::string out_plant_slantRange_path;// 平面斜距文件 std::string out_plant_slantRange_path;// 平面斜距文件
//// ASF 方法 ///////////// //// ASF 方法 /////////////
std::string out_lon_path;// 经度 std::string out_lon_path;// 经度
std::string out_lat_path;// 纬度 std::string out_lat_path;// 纬度
std::string out_slantRange_path;// 斜距文件 std::string out_slantRange_path;// 斜距文件
/// 共同文件 /// /// 共同文件 ///
std::string out_localIncidenct_path;// 计算局地入射角 std::string out_localIncidenct_path;// 计算局地入射角
std::string out_incidence_path; // 入射角文件 std::string out_incidence_path; // 入射角文件
std::string out_ori_sim_tiff;// 映射角文件' std::string out_ori_sim_tiff;// 映射角文件'
std::string out_sim_ori_tiff; std::string out_sim_ori_tiff;
}; };

8
test_moudel.cpp
View File

@ -88,12 +88,12 @@ int test_mode_0()
int test_mode_1(string rootpath) int test_mode_1(string rootpath)
{ {
// .\baseTool\x64\Release\SIMOrthoProgram.exe 1 D:\MicroWorkspace\C-SAR\Ortho\Temporary\unpack\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_AHV_L10003923848\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_HH_L10003923848.tiff // .\baseTool\x64\Release\SIMOrthoProgram.exe 1 D:\MicroWorkspace\C-SAR\Ortho\Temporary\unpack\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_AHV_L10003923848\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_HH_L10003923848.tiff
//<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> //???????
std::string root_path = rootpath;// "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary2"; std::string root_path = rootpath;// "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary2";
std::string workspace = "D:\\MicroSAR\\C-SAR\\MicroWorkspace\\Ortho\\Temporary"; std::string workspace = "D:\\MicroSAR\\C-SAR\\MicroWorkspace\\Ortho\\Temporary";
std::string parameter_path =workspace+ "\\package\\orth_para.txt"; // std::string parameter_path =workspace+ "\\package\\orth_para.txt"; //
std::string dem_path = workspace + "\\TestDEM\\mergedDEM.tif"; // std::string dem_path = workspace + "\\TestDEM\\pindem.tif"; //
std::string in_sar_path = workspace + "\\unpack\\GF3_KAS_FSI_020253_E110.8_N25.5_20200614_L1A_HHHV_L10004871459\\GF3_KAS_FSI_020253_E110.8_N25.5_20200614_L1A_HH_L10004871459.tiff"; // std::string in_sar_path = workspace + "\\unpack\\GF3_KAS_FSI_020253_E110.8_N25.5_20200614_L1A_HHHV_L10004871459\\GF3_KAS_FSI_020253_E110.8_N25.5_20200614_L1A_HH_L10004871459.tiff"; //
std::string work_path = workspace ; // std::string work_path = workspace ; //
std::string taget_path = workspace + "\\package"; // std::string taget_path = workspace + "\\package"; //
@ -124,7 +124,7 @@ int test_ASF(string rootpath) {
std::string out_GEC_lon_lat_path = root_path + "\\output\\GEC_lon_lat.tif"; std::string out_GEC_lon_lat_path = root_path + "\\output\\GEC_lon_lat.tif";
std::string out_clip_DEM_path = root_path + "\\output\\Orth_dem.tif"; std::string out_clip_DEM_path = root_path + "\\output\\Orth_dem.tif";
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
std::cout << "Ԥ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>DEM <20><>Χ" << 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::cout << "out_GEC_dem_path out_GTC_rc_path out_GEC_lon_lat_path out_clip_dem_path" << endl; std::cout << "out_GEC_dem_path out_GTC_rc_path out_GEC_lon_lat_path out_clip_dem_path" << endl;
std::cout << "==========================================================================" << endl; std::cout << "==========================================================================" << endl;
@ -277,6 +277,6 @@ int test_LLA2XYZ_XYZ2LLA()
int test_vector_operator() int test_vector_operator()
{ {
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD> // ????
return 0; return 0;
} }