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

调整了 正射计算模块修改了RPC，增加了DEM传入接口

master
剑古敛锋 2023-02-14 10:55:14 +08:00
parent 929e67e347
commit cb3db54958
20 changed files with 1311 additions and 792 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_VECTORIZE_SSE4_2
//#include <mkl.h>
// 本地方法
// 本地方法
#include <iostream>
#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);
double min_value = 0, temp_min = 0;
double max_value = 0, temp_max = 0;
// 线性拉伸2%
// 线性拉伸2%
Eigen::MatrixXd hist = gdalimg.getHist(band_ids);
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);
}
// 重新缩放最大值，最小值
std::cout << "min value ：\t" << min_value << "\n";
std::cout << "max value ：\t" << max_value << "\n";
// 重新缩放最大值，最小值
std::cout << "min value ：\t" << min_value << "\n";
std::cout << "max value ：\t" << max_value << "\n";
start_ids = 0;
do {
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;
//bilateralFilter(img, result1, 10, 80, 50);
vector<int> compression_params;
compression_params.push_back(cv::IMWRITE_JPEG_QUALITY); //选择jpeg
compression_params.push_back(100); //在这个填入你要的图片质量
compression_params.push_back(cv::IMWRITE_JPEG_QUALITY); //选择jpeg
compression_params.push_back(100); //在这个填入你要的图片质量
cv::Mat out = result1;
//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 << "=========================================\n";
std::cout << " convert gdal to jpg :\n";
std::cout << "input gdal img file path ：\t" << gdal_path << "\n";
std::cout << "input gdal img band ：\t" << band_ids << "\n";
std::cout << "out jpg file path ：\t" << jpg_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 << "out jpg file path ：\t" << jpg_path << "\n";
std::cout << "=========================================\n";
return 0;
}
@ -110,7 +110,7 @@ int ImageMatch::gdal2JPG(std::string gdal_path, std::string jpg_path, int band_i
/// <summary>
/// 获取模型匹配点
/// 获取模型匹配点
/// </summary>
/// <param name="ori_power_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 sim image path : \t" << sim_sum_path << endl;
//读取影像
//读取影像
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);
double minvalue = 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);
if (sig >1) {
//continue;
}
// 构建搜索域
// 构建搜索域
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 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_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 result;
matchTemplate(serch, templet_mat, result, cv::TM_CCOEFF_NORMED);
Point minLoc;
Point maxLoc;
Point matchLoc;
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_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;
double offset_x = 0, offset_y = 0;
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;
omp_lock_t lock;
std::cout << "ori_x\tori_y\tsim_x\tsim_y\tmaxval \t" << endl;
omp_init_lock(&lock); // 初始化互斥锁
omp_init_lock(&lock); // 初始化互斥锁
count = 0;
int search_count = 0;
#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));
matchTemplate(sim, templet_mat, result, cv::TM_CCOEFF_NORMED);
//normalize(result, result, 1, 0, cv::NORM_MINMAX);
// 通过函数 minMaxLoc 定位最匹配的位置;
omp_set_lock(&lock); //获得互斥器
// 通过函数 minMaxLoc 定位最匹配的位置;
omp_set_lock(&lock); //获得互斥器
Point minLoc;
Point maxLoc;
Point matchLoc;
cv::minMaxLoc(result, &tempminVal, &tempmaxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
cv::minMaxLoc(result, &tempminVal, &tempmaxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
if (tempmaxVal >= maxVal) {
offset_x = maxLoc.x - (j - 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;
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_y = offset_y * 1.0 / offsetcount;
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 << "rough match out : \t" << getCurrentTimeString() << endl;
// step1.1: 粗匹配绘制结果
// step1.1: 粗匹配绘制结果
std::string rough_math_path = sim_sum_path;
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);
vector<int> compression_params;
compression_params.push_back(cv::IMWRITE_JPEG_QUALITY); //选择jpeg
compression_params.push_back(100); //在这个填入你要的图片质量
compression_params.push_back(cv::IMWRITE_JPEG_QUALITY); //选择jpeg
compression_params.push_back(100); //在这个填入你要的图片质量
cv::Mat out = ori_in_sim;
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;
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;
//// step 2: 精匹配，
//// step 2: 精匹配，
//std::cout << "Precise match , begining : \t" << getCurrentTimeString() << endl;
//Eigen::MatrixXd matchpoints;
//{
@ -366,8 +366,8 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
// row_count = row_count - PreciseSize;
// col_count = col_count - PreciseSize;
// omp_lock_t lock;
// omp_init_lock(&lock); // 初始化互斥锁
// // 以搜索范围为核心计算
// omp_init_lock(&lock); // 初始化互斥锁
// // 以搜索范围为核心计算
// #pragma omp parallel for num_threads(8)
// for (int i = searchSize + offset_y; i < row_count; i = i + preciseStep) { // y
// cv::Mat templet_mat;
@ -377,13 +377,13 @@ Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::
// cv::Mat result;
// double minVal = 100, maxVal = 0;
// for (int j = searchSize + offset_x; j < col_count; j = j + preciseStep) { //x
// // 计算 起始点
// // 计算 起始点
// sim_start_x = j - searchSize;
// sim_start_y = i - searchSize;
// ori_start_x = (j - searchSize / 2) - PreciseSize / 2;
// ori_start_y = (i - searchSize / 2) - PreciseSize / 2;
// // 匹配模板，待匹配模板
// // 匹配模板，待匹配模板
// templet_mat = ori(Rect(ori_start_x, ori_start_y, PreciseSize, PreciseSize));
// 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);
// matchTemplate(sim, templet_mat, result, cv::TM_CCORR_NORMED);
// // 通过函数 minMaxLoc 定位最匹配的位置;
// // 通过函数 minMaxLoc 定位最匹配的位置;
// cv::Point minLoc; cv::Point maxLoc;
// cv::Point matchLoc;
// cv::minMaxLoc(result, &minVal, &maxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
// cv::minMaxLoc(result, &minVal, &maxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
// if (maxVal > 0.7) {
// omp_set_lock(&lock); //获得互斥器
// omp_set_lock(&lock); //获得互斥器
//
// matchpoints(count, 0) = ori_start_x; //ori_x
// 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
// 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
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 temp(offset_x.rows(), 6);
temp.col(0) = temp.col(0).array()*0+1; //1
@ -488,16 +488,16 @@ int ImageMatch::outMatchModel(std::string matchmodel_path)
}
/// <summary>
/// 读取jpg 文件
/// 读取jpg 文件
/// </summary>
/// <param name="jpg_path"></param>
/// <returns></returns>
cv::Mat openJPG(std::string 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;
}

4
ImageMatch.h
View File

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

78
OctreeNode.h
View File

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

2
RPC_Correct.cpp
View File

@ -10,7 +10,6 @@
#include "SateOrbit.h"
#include "ImageMatch.h"
#include <gdal_utils.h>
#include <armadillo>
#include <proj.h>
#include "gdal_priv.h"
#include "gdal_alg.h"
@ -19,7 +18,6 @@
using namespace std;
using namespace Eigen;
using namespace arma;

8
RPC_Correct.h
View File

@ -6,11 +6,9 @@
#include <boost/filesystem.hpp>
#include <omp.h>
#include <gdal_utils.h>
#include <armadillo>
#include <proj.h>
#include "gdal_priv.h"
#include "gdal_alg.h"
#include <armadillo>
//#include <mkl.h>
#include "baseTool.h"
#include "simptsn.h"
@ -18,13 +16,9 @@
#include "ImageMatch.h"
using namespace std;
using namespace Eigen;
using namespace arma;
using namespace std;
using namespace Eigen;
using namespace arma;
// 专门用于解析RPC
// 专门用于解析RPC

21
SIMOrthoProgram.cpp
View File

@ -24,7 +24,7 @@
using namespace std;
using namespace Eigen;
//mode 1
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
@ -32,11 +32,12 @@ void PreProcess(int argc, char* argv[])
std::cout << "==========================================================================" << 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::string parameter_path = argv[2];
std::string dem_path = argv[3];
std::string in_sar_path = argv[4];
std::string work_path = argv[5];
std::string taget_path = argv[6];
std::string parameter_path = argv[2]; // 参数文件
std::string dem_path = argv[3]; // dem 文件
std::string in_sar_path = argv[4]; // 输入SAR文件
std::string work_path = argv[5]; // 目标空间文件
std::string taget_path = argv[6]; // 输出坐标映射文件
//std::string Incident_path = argv[7];// 输出入射角文件
std::cout << "==========================================================================" << endl;
std::cout << "in parameters:========================================================" << endl;
@ -190,11 +191,13 @@ void createRPC_lon_lat(int argc, char* argv[]) {
std::cout << "mode 8";
std::cout << "SIMOrthoProgram.exe 8 in_rpc_tiff out_lon_lat_path";
std::string in_rpc_tiff = argv[2];
std::string out_lon_lat_path = argv[3];
std::string in_dem_tiff = argv[3];
std::string out_lon_lat_path = argv[4];
simProcess process;
//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::cout << "==========================================================================" << endl;
process.CreateRPC_refrenceTable(in_rpc_tiff, out_lon_lat_path);
process.CreateRPC_refrenceTable(in_rpc_tiff, in_dem_tiff,out_lon_lat_path);
std::cout << "==========================================================================" << endl;
}
@ -281,7 +284,7 @@ int main(int argc, char* argv[])
Landpoint p2 = { -3421843,5089485,3630606 }; Landpoint p1 = { -2609414,4763328,3332879 };
cout << getIncAngle(p2, p1) << endl;;
std::cout << "program start:\t" << getCurrentTimeString() << endl;;
int mode = 10;
int mode = 1;
GDALAllRegister();
if (argc == 1) { // 测试参数
// 算法说明

4
SIMOrthoProgram.rc
View File

@ -13,7 +13,7 @@
#undef APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
// 中文(简体，中国) resources
// 中文(简体，中国) resources
#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_CHS)
LANGUAGE LANG_CHINESE, SUBLANG_CHINESE_SIMPLIFIED
@ -52,7 +52,7 @@ END
IDR_PROJ1 proj "proj.db"
#endif // 中文(简体，中国) resources
#endif // 中文(简体，中国) resources
/////////////////////////////////////////////////////////////////////////////

3
SIMOrthoProgram.vcxproj
View File

@ -167,6 +167,7 @@
<ItemGroup>
<ClCompile Include="baseTool.cpp" />
<ClCompile Include="ImageMatch.cpp" />
<ClCompile Include="interpolation.cpp" />
<ClCompile Include="OctreeNode.cpp" />
<ClCompile Include="RPC_Correct.cpp" />
<ClCompile Include="simptsn.cpp" />
@ -177,6 +178,8 @@
<ItemGroup>
<ClInclude Include="baseTool.h" />
<ClInclude Include="ImageMatch.h" />
<ClInclude Include="interpolation.h" />
<ClInclude Include="linterp.h" />
<ClInclude Include="OctreeNode.h" />
<ClInclude Include="resource.h" />
<ClInclude Include="RPC_Correct.h" />

6
SIMOrthoProgram.vcxproj.filters
View File

@ -39,6 +39,9 @@
<ClCompile Include="OctreeNode.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="interpolation.cpp">
<Filter>源文件</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<ClInclude Include="baseTool.h">
@ -65,6 +68,9 @@
<ClInclude Include="OctreeNode.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="interpolation.h">
<Filter>头文件</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ResourceCompile Include="SIMOrthoProgram.rc">

25
SateOrbit.cpp
View File

@ -7,12 +7,9 @@
#include <Eigen/Dense>
#include <time.h>
//#include <mkl.h>
#include <armadillo>
using namespace std;
using namespace Eigen;
using namespace arma;
@ -23,7 +20,7 @@ OrbitPoly::OrbitPoly()
OrbitPoly::OrbitPoly(int polynum, Eigen::MatrixX<double> polySatellitePara, double SatelliteModelStartTime)
{
if (polySatellitePara.rows() != polynum||polySatellitePara.cols()!=6) {
throw exception("轨道模型参数形状不一致");
throw exception("轨道模型参数形状不一致");
}
this->polySatellitePara = polySatellitePara;
this->SatelliteModelStartTime = SatelliteModelStartTime;
@ -45,27 +42,27 @@ Eigen::MatrixX<double> OrbitPoly::SatelliteSpacePoint(long double satellitetime)
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 ");
}
// 更新轨道时间
// 更新轨道时间
double satellitetime2 =double( satellitetime - this->SatelliteModelStartTime);
Eigen::MatrixX<double> satetime(1, polynum);
for (int i = 0; i < polynum; i++) {
satetime(0, i) = pow(satellitetime2, i);
}
// 计算
// 计算
Eigen::MatrixX<double> satellitePoints(1, 6);
satellitePoints = satetime * polySatellitePara;
return satellitePoints;
}
/// <summary>
/// 获取指定时刻的卫星轨道坐标
/// 获取指定时刻的卫星轨道坐标
/// </summary>
/// <param name="satellitetime">卫星时刻</param>
/// <param name="SatelliteModelStartTime">模型起算时间</param>
/// <param name="polySatellitePara">模型参数[x1,y1,z1,vx1,vy1,vz1; x2,y2,z2,vx2,vy2,vz2; ..... ]</param>
/// <param name="polynum">模型参数数量</param>
/// <returns>卫星坐标</returns>
/// <param name="satellitetime">卫星时刻</param>
/// <param name="SatelliteModelStartTime">模型起算时间</param>
/// <param name="polySatellitePara">模型参数[x1,y1,z1,vx1,vy1,vz1; x2,y2,z2,vx2,vy2,vz2; ..... ]</param>
/// <param name="polynum">模型参数数量</param>
/// <returns>卫星坐标</returns>
Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double>& satellitetime, double SatelliteModelStartTime, Eigen::MatrixX<double>& polySatellitePara, int polynum)
{
if (satellitetime.cols() != 1) {
@ -74,7 +71,7 @@ Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double>& satellitetim
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 ");
}
// 更新轨道时间
// 更新轨道时间
int satellitetime_num = satellitetime.rows();
satellitetime = satellitetime.array() - SatelliteModelStartTime;
Eigen::MatrixX<double> satelliteTime(satellitetime_num, polynum);
@ -82,7 +79,7 @@ Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double>& satellitetim
satelliteTime.col(i) = satellitetime.array().pow(i);
}
// 计算
// 计算
Eigen::MatrixX<double> satellitePoints(satellitetime_num, 6);
satellitePoints = satelliteTime * polySatellitePara;
return satellitePoints;

18
SateOrbit.h
View File

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

304
baseTool.cpp
View File

@ -1,6 +1,6 @@
#pragma once
///
/// 基本类、基本函数
/// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ࡢ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
///
//#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2
@ -38,7 +38,7 @@ using namespace cv;
/**
* @brief GetCurrentTime 获取当前时间
* @brief GetCurrentTime <EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD>ǰʱ<EFBFBD><EFBFBD>
* @return
*/
std::string getCurrentTimeString() {
@ -117,20 +117,20 @@ Landpoint crossProduct(const Landpoint& a, const Landpoint& b) {
}
/// <summary>
/// 计算地表坡度向量，输入点都是WGS84坐标系，且p1为北点，
/// 北 N
/// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ر<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>㶼<EFBFBD><EFBFBD>WGS84<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>p1Ϊ<EFBFBD><EFBFBD><EFBFBD>㣬
/// <EFBFBD><EFBFBD> N
/// W p1
/// p4 p0 p2 E
/// p3 S
/// n21=n2xn1
/// n=(n21+n32+n43+n41)*0.25;
/// </summary>
/// <param name="p0">目标点</param>
/// <param name="p1">周围点1</param>
/// <param name="p2">周围点2</param>
/// <param name="p3">周围点3</param>
/// <param name="p4">周围点4</param>
/// <returns>向量角度</returns>
/// <param name="p0">Ŀ<EFBFBD><EFBFBD>ￄ1<EFBFBD>7</param>
/// <param name="p1"><EFBFBD><EFBFBD>Χ<EFBFBD><EFBFBD>1</param>
/// <param name="p2"><EFBFBD><EFBFBD>Χ<EFBFBD><EFBFBD>2</param>
/// <param name="p3"><EFBFBD><EFBFBD>Χ<EFBFBD><EFBFBD>3</param>
/// <param name="p4"><EFBFBD><EFBFBD>Χ<EFBFBD><EFBFBD>4</param>
/// <returns><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƕ<EFBFBD></returns>
Landpoint getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landpoint& p2, const Landpoint& p3, const Landpoint& p4) {
Landpoint n0 = LLA2XYZ(p0),
@ -139,11 +139,11 @@ Landpoint getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landp
n3 = LLA2XYZ(p3),
n4 = LLA2XYZ(p4);
Landpoint n01 = n1 - n0, n02 = n2 - n0, n03 = n3 - n0, n04 = n4 - n0;
// 以n01为正北方向
// <EFBFBD><EFBFBD>n01Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Landpoint np01 = p1-p0, np02 = p2-p0, np03 = p3-p0, np04 = p4-p0;
double a2 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np02.lon,np02.lat,0 });// 01->02 逆时针
double a3 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np03.lon,np03.lat,0 });// 01->03 逆时针
double a4 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np04.lon,np04.lat,0 });// 01->04 逆时针
double a2 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np02.lon,np02.lat,0 });// 01->02 <EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>
double a3 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np03.lon,np03.lat,0 });// 01->03 <EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>
double a4 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np04.lon,np04.lat,0 });// 01->04 <EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>
//std::cout << a2 << "\t" << a3 << "\t" << a4 << endl;
a2 = 360 - a2;
a3 = 360 - a3;
@ -188,24 +188,24 @@ Landpoint getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landp
}
/// <summary>
/// 根据定义计算大小
/// <EFBFBD><EFBFBD><EFBFBD>ݶ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С
/// p11 p12 p13 p14
/// p21 p22 p23 p24
/// p(2+u,2+v)
/// p31 p32 p33 p34
/// p41 p42 p43 p44
/// </summary>
/// <param name="u">像元行坐标的小数部分</param>
/// <param name="v">像元行坐标的小数部分</param>
/// <param name="img">4x4 插值区域</param>
/// <param name="u"><EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>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="img">4x4 <EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD></param>
/// <returns></returns>
complex<double> Cubic_Convolution_interpolation(double u, double v, Eigen::MatrixX<complex<double>> img)
{
if (img.rows() != 4 || img.cols() != 4) {
throw exception("the size of img's block is not right");
}
// 计算掩模版
Eigen::MatrixX<complex<double>> wrc(1, 4);// 使用 complex<double> 类型主要原因是为了获取值
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD>
Eigen::MatrixX<complex<double>> wrc(1, 4);// ʹ<EFBFBD><EFBFBD> complex<double> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫԭ<D2AA><D4AD><EFBFBD><EFBFBD>Ϊ<EFBFBD>˻<EFBFBD>ȡֵ
Eigen::MatrixX<complex<double>> wcr(4, 1);//
for (int i = 0; i < 4; i++) {
wrc(0, i) = Cubic_kernel_weight(u+1-i); // u+1,u,u-1,u-2
@ -261,10 +261,10 @@ double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpo
/// <summary>
/// 将经纬度转换为地固参心坐标系
/// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>γ<EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD>Ϊ<EFBFBD>ع̲<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ
/// </summary>
/// <param name="XYZP">经纬度点--degree</param>
/// <returns>投影坐标系点</returns>
/// <param name="XYZP"><EFBFBD><EFBFBD>γ<EFBFBD>ȵ<EFBFBD>--degree</param>
/// <returns>ͶӰ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD></returns>
Landpoint LLA2XYZ(const Landpoint& LLA) {
double L = LLA.lon * d2r;
double B = LLA.lat * d2r;
@ -284,7 +284,7 @@ Landpoint LLA2XYZ(const Landpoint& LLA) {
}
/// <summary>
/// 矩阵（n,3) lon,lat,ati
/// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>n,3) lon,lat,ati
/// </summary>
/// <param name="landpoint"></param>
/// <returns></returns>
@ -309,13 +309,13 @@ Eigen::MatrixXd LLA2XYZ(Eigen::MatrixXd landpoint)
/// <summary>
/// 将地固参心坐标系转换为经纬度
/// <EFBFBD><EFBFBD><EFBFBD>ع̲<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD>γ<EFBFBD><EFBFBD>
/// </summary>
/// <param name="XYZ">固参心坐标系</param>
/// <returns>经纬度--degree</returns>
/// <param name="XYZ"><EFBFBD>̲<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ</param>
/// <returns><EFBFBD><EFBFBD>γ<EFBFBD><EFBFBD>--degree</returns>
Landpoint XYZ2LLA(const Landpoint& XYZ) {
double tmpX = XYZ.lon;// 经度 x
double temY = XYZ.lat;// 纬度 y
double tmpX = XYZ.lon;// <EFBFBD><EFBFBD><EFBFBD><EFBFBD> x
double temY = XYZ.lat;// γ<EFBFBD><EFBFBD> y
double temZ = XYZ.ati;
double curB = 0;
@ -341,42 +341,42 @@ Landpoint XYZ2LLA(const Landpoint& XYZ) {
}
/// <summary>
/// 根据图像获取影像
/// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ<EFBFBD><EFBFBD><EFBFBD>ȡӰ<EFBFBD>ￄ1<EFBFBD>7
/// </summary>
/// <param name="dem_path"></param>
gdalImage::gdalImage(string raster_path)
{
omp_lock_t lock;
omp_init_lock(&lock); // 初始化互斥锁
omp_set_lock(&lock); //获得互斥器
omp_init_lock(&lock); // <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
omp_set_lock(&lock); //<EFBFBD><EFBFBD>û<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
this->img_path = raster_path;
GDALAllRegister();// 注册格式的驱动
// 打开DEM影像
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.c_str(), GA_ReadOnly));//以只读方式读取地形影像
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
// <EFBFBD><EFBFBD>DEMӰ<EFBFBD><EFBFBD>
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.c_str(), GA_ReadOnly));//<EFBFBD><EFBFBD>ֻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD>
this->width = rasterDataset->GetRasterXSize();
this->height = rasterDataset->GetRasterYSize();
this->band_num = rasterDataset->GetRasterCount();
double* gt = new double[6];
// 获得仿射矩阵
// <EFBFBD><EFBFBD>÷<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
rasterDataset->GetGeoTransform(gt);
this->gt = Eigen::MatrixXd(2, 3);
this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5];
this->projection = rasterDataset->GetProjectionRef();
// 逆仿射矩阵
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//double* inv_gt = new double[6];;
//GDALInvGeoTransform(gt, inv_gt); // 逆仿射矩阵
// 构建投影
//GDALInvGeoTransform(gt, inv_gt); // <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͶӰ
GDALFlushCache((GDALDatasetH)rasterDataset);
GDALClose((GDALDatasetH)rasterDataset);
rasterDataset = NULL;// 指针置空
rasterDataset = NULL;// ָ<EFBFBD><EFBFBD><EFBFBD>ÿ<EFBFBD>
this->InitInv_gt();
delete[] gt;
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); //释放互斥器
omp_destroy_lock(&lock); //销毁互斥器
omp_unset_lock(&lock); //<EFBFBD>ͷŻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
omp_destroy_lock(&lock); //<EFBFBD><EFBFBD><EFBFBD>ٻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
}
gdalImage::~gdalImage()
@ -406,11 +406,11 @@ 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)
{
omp_lock_t lock;
omp_init_lock(&lock); // 初始化互斥锁
omp_set_lock(&lock); //获得互斥器
GDALAllRegister();// 注册格式的驱动
// 打开DEM影像
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_ReadOnly));//以只读方式读取地形影像
omp_init_lock(&lock); // <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
omp_set_lock(&lock); //<EFBFBD><EFBFBD>û<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
// <EFBFBD><EFBFBD>DEMӰ<EFBFBD><EFBFBD>
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_ReadOnly));//<EFBFBD><EFBFBD>ֻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD>
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids);
@ -456,15 +456,15 @@ Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count,
delete[] temp;
}
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); //释放互斥器
omp_destroy_lock(&lock); //销毁互斥器
omp_unset_lock(&lock); //<EFBFBD>ͷŻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
omp_destroy_lock(&lock); //<EFBFBD><EFBFBD><EFBFBD>ٻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return datamatrix;
}
/// <summary>
/// 写入遥感影像
/// д<EFBFBD><EFBFBD>ң<EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD>
/// </summary>
/// <param name="data"></param>
/// <param name="start_row"></param>
@ -473,25 +473,25 @@ 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)
{
omp_lock_t lock;
omp_init_lock(&lock); // 初始化互斥锁
omp_set_lock(&lock); //获得互斥器
// 合法性检测
// 写入数据
omp_init_lock(&lock); // <EFBFBD><EFBFBD>ʼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
omp_set_lock(&lock); //<EFBFBD><EFBFBD>û<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
// <EFBFBD>Ϸ<EFBFBD><EFBFBD>Լ<EFBFBD>ￄ1<EFBFBD>7
// д<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
string tip = "写入影像范围超越了目标范围" + this->img_path;
string tip = "д<EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD>Χ<EFBFBD><EFBFBD>Խ<EFBFBD><EFBFBD>Ŀ<EFBFBD>귶Χ" + this->img_path;
throw exception(tip.c_str());
}
GDALAllRegister();// 注册格式的驱动
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = nullptr;
if (boost::filesystem::exists(this->img_path)) {
poDstDS = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_Update));
}
else {
poDstDS = poDriver->Create(this->img_path.c_str(), this->width, this->height, this->band_num, GDT_Float32, NULL); // 输出结果
poDstDS = poDriver->Create(this->img_path.c_str(), this->width, this->height, this->band_num, GDT_Float32, NULL); // <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
poDstDS->SetProjection(this->projection.c_str());
// 设置转换矩阵
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
double gt_ptr[6];
for (int i = 0; i < this->gt.rows(); i++) {
for (int j = 0; j < this->gt.cols(); j++) {
@ -520,8 +520,8 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
GDALClose((GDALDatasetH)poDstDS);
//GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
delete[] databuffer;
omp_unset_lock(&lock); //释放互斥器
omp_destroy_lock(&lock); //销毁互斥器
omp_unset_lock(&lock); //<EFBFBD>ͷŻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
omp_destroy_lock(&lock); //<EFBFBD><EFBFBD><EFBFBD>ٻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
}
void gdalImage::saveImage()
@ -533,7 +533,7 @@ void gdalImage::setNoDataValue(double nodatavalue = -9999, int band_ids = 1)
{
GDALAllRegister();// 注册格式的驱动
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
//GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.c_str(), GA_Update));
poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue);
@ -650,13 +650,13 @@ GDALRPCInfo gdalImage::getRPC()
{
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
CPLSetConfigOption("GDAL_DATA", "./data");
GDALAllRegister();//注册驱动
//打开数据
GDALAllRegister();//ע<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
GDALDataset* pDS = (GDALDataset*)GDALOpen(this->img_path.c_str(), GA_ReadOnly);
//从元数据中获取RPC信息
//<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>л<EFBFBD>ȡRPC<EFBFBD><EFBFBD>Ϣ
char** papszRPC = pDS->GetMetadata("RPC");
//将获取的RPC信息构造成结构体
//<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD>RPC<EFBFBD><EFBFBD>Ϣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɽṹ<EFBFBD>ￄ1<EFBFBD>7
GDALRPCInfo oInfo;
GDALExtractRPCInfo(papszRPC, &oInfo);
@ -672,9 +672,9 @@ Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points)
}
Eigen::MatrixXd result(points.rows(), 3);
result.col(2) = points.col(2);// 高程
result.col(2) = points.col(2);// <EFBFBD>߳<EFBFBD>
points.col(2) = points.col(2).array() * 0 + 1;
points.col(0).swap(points.col(2));// 交换
points.col(0).swap(points.col(2));// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Eigen::MatrixXd gts(3, 2);
gts.col(0) = this->gt.row(0);
gts.col(1) = this->gt.row(1);
@ -691,14 +691,14 @@ Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points)
Eigen::MatrixXd gdalImage::getHist(int bandids)
{
GDALAllRegister();// 注册格式的驱动
// 打开DEM影像
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_ReadOnly));//以只读方式读取地形影像
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
// <EFBFBD><EFBFBD>DEMӰ<EFBFBD><EFBFBD>
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.c_str(), GA_ReadOnly));//<EFBFBD><EFBFBD>ֻ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD>
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* xBand = rasterDataset->GetRasterBand(bandids);
//获取图像直方图
//<EFBFBD><EFBFBD>ȡͼ<EFBFBD><EFBFBD>ֱ<EFBFBD><EFBFBD>ͼ
double dfMin = this->min(bandids);
double dfMax = this->max(bandids);
@ -723,13 +723,13 @@ gdalImage CreategdalImage(string img_path, int height, int width, int band_num,
if (boost::filesystem::exists(img_path.c_str())) {
boost::filesystem::remove(img_path.c_str());
}
GDALAllRegister();// 注册格式的驱动
GDALAllRegister();// ע<EFBFBD><EFBFBD><EFBFBD>ʽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ￄ1<EFBFBD>7
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = poDriver->Create(img_path.c_str(), width, height, band_num, GDT_Float32, NULL); // 输出结果
GDALDataset* poDstDS = poDriver->Create(img_path.c_str(), width, height, band_num, GDT_Float32, NULL); // <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
if (need_gt) {
poDstDS->SetProjection(projection.c_str());
// 设置转换矩阵
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
double gt_ptr[6] = { 0 };
for (int i = 0; i < gt.rows(); i++) {
for (int j = 0; j < gt.cols(); j++) {
@ -750,7 +750,7 @@ gdalImage CreategdalImage(string img_path, int height, int width, int band_num,
/// <summary>
/// 裁剪DEM
/// <EFBFBD>ü<EFBFBD>DEM
/// </summary>
/// <param name="in_path"></param>
/// <param name="out_path"></param>
@ -766,17 +766,17 @@ void clipGdalImage(string in_path, string out_path, DemBox box, double pixelinte
/***
* 遥感影像重采样 (要求影像必须有投影，否则走不通)
* @param pszSrcFile 输入文件的路径
* @param pszOutFile 写入的结果图像的路径
* @param eResample 采样模式，有五种，具体参见GDALResampleAlg定义，默认为双线性内插
* @param gt X转换采样比，默认大小为1.0，大于1图像变大，小于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>̈́1<EFBFBD>7)
* @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 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 new_width
* @param new_hegiht
* @retrieve 0 成功
* @retrieve -1 打开源文件失败
* @retrieve -2 创建新文件失败
* @retrieve -3 处理过程中出错
* @retrieve 0 <EFBFBD>ɹ<EFBFBD>
* @retrieve -1 <EFBFBD><EFBFBD>Դ<EFBFBD>ļ<EFBFBD>ʧ<EFBFBD><EFBFBD>
* @retrieve -2 <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD>ʧ<EFBFBD><EFBFBD>
* @retrieve -3 <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>г<EFBFBD><EFBFBD><EFBFBD>
*/
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample)
{
@ -802,12 +802,12 @@ int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int
char* pszSrcWKT = NULL;
pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef());
//如果没有投影，人为设置一个
//<EFBFBD><EFBFBD><EFBFBD>û<EFBFBD><EFBFBD>ͶӰ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD>ￄ1<EFBFBD>7
if (strlen(pszSrcWKT) <= 0)
{
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
//oSRS.SetUTM(50, true); //北半球 东经120度
//oSRS.SetUTM(50, true); //<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD>120<32><30>
//oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&pszSrcWKT);
}
@ -815,7 +815,7 @@ int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int
void* hTransformArg;
hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT, FALSE, 0.0, 1);
std::cout << "no proj " << endl;
//(没有投影的影像到这里就走不通了)
//(û<EFBFBD><EFBFBD>ͶӰ<EFBFBD><EFBFBD>Ӱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߲<EFBFBD>ͨ<EFBFBD>ￄ1<EFBFBD>7)
if (hTransformArg == NULL)
{
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
@ -832,7 +832,7 @@ int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int
//GDALDestroyGenImgProjTransformer(hTransformArg);
//创建结果数据集
//<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݼￄ1<EFBFBD>7
GDALDataset* pDDst = pDriver->Create(pszOutFile, new_width, new_height, nBandCount, dataType, NULL);
if (pDDst == NULL)
{
@ -944,7 +944,7 @@ int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResa
/// <summary>
/// 数字转字符串
/// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD>ַ<EFBFBD><EFBFBD><EFBFBD>
/// </summary>
/// <param name="Num"></param>
/// <returns></returns>
@ -958,7 +958,7 @@ string Convert(float Num)
/// <summary>
/// 路径拼接
/// ·<EFBFBD><EFBFBD>ƴ<EFBFBD><EFBFBD>
/// </summary>
/// <param name="path"></param>
/// <param name="filename"></param>
@ -980,6 +980,9 @@ std::string JoinPath(const std::string& path, const std::string& filename)
}
}
/// <summary>
////////////////////////////////////////////////////////////////////// WGS84 to J2000 /////////////////////////////////////////////////////////////////////////////////////////
/// </summary>
@ -1071,24 +1074,24 @@ WGS84_J2000::~WGS84_J2000()
}
/// <summary>
/// step1 WGS 84 转换到协议地球坐标系。
/// 大地经纬度高程BLH与地固坐标系的转换 BLH-- > XG
/// ECEF, Earth Centered Earth Fixed; 地固坐标系 参考平面：平赤道面，即过地心，并且与地心与CIO点连线垂直的平面。
/// x轴为参考平面与格林尼治平面交线，z轴为地心指向CIO点。
/// step1 WGS 84 ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Э<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD>ￄ1<EFBFBD>7
/// <EFBFBD><EFBFBD>ؾ<EFBFBD>γ<EFBFBD>ȸ߳<EFBFBD>BLH<EFBFBD><EFBFBD>ع<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>ת<EFBFBD>ￄ1<EFBFBD>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>档
/// x<EFBFBD><EFBFBD>Ϊ<EFBFBD>ο<EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD>潻<EFBFBD>ߣ<EFBFBD>z<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ָ<EFBFBD><EFBFBD>CIO<EFBFBD>㡣
/// </summary>
/// <param name="LBH_deg_m">B L H分别为大地纬度、大地经度和大地高程 输入参数为N * 3矩阵</param>
/// <returns>XYZ 为地固坐标系的 x y z，输出参数为N * 3矩阵</returns>
/// <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>
/// <returns>XYZ Ϊ<EFBFBD>ع<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD> 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)
{
return LLA2XYZ(LBH_deg_m);
}
/// <summary>
/// step 2 协议地球坐标系 转换为瞬时地球坐标系
/// 这主要涉及查询 给定时刻的 xp,yp , 可以查询EOP文件获得，下载地址如下http://celestrak.com/spacedata/
/// step 2 Э<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>τ1<EFBFBD>7 ת<><D7AA>Ϊ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ
/// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD>漰<EFBFBD><EFBFBD>ѯ <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;
/// </summary>
/// <param name="axis">axis 表示围绕旋转的坐标轴 '1' 表示围绕 x轴逆时针旋转;'2' 表示围绕 y轴逆时针旋转;'3'表示围绕 z轴逆时针旋转</param>
/// <param name="angle_deg">angle_deg 表示旋转的角度 默认 degree</param>
/// <param name="axis">axis <EFBFBD><EFBFBD>ʾΧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> '1' <20><>ʾΧ<CABE><CEA7> x<><78><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1><EFBFBD><EFBFBD>ת;'2' <20><>ʾΧ<CABE><CEA7> y<><79><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1><EFBFBD><EFBFBD>ת;'3'<27><>ʾΧ<CABE><CEA7> z<><7A><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1><EFBFBD><EFBFBD>ת</param>
/// <param name="angle_deg">angle_deg <EFBFBD><EFBFBD>ʾ<EFBFBD><EFBFBD>ת<EFBFBD>ĽǶ<EFBFBD> Ĭ<><C4AC> degree</param>
/// <returns></returns>
Eigen::MatrixXd WGS84_J2000::ordinateSingleRotate(int axis, double angle_deg)
{
@ -1119,41 +1122,41 @@ Eigen::MatrixXd WGS84_J2000::ordinateSingleRotate(int axis, double angle_deg)
}
/// <summary>
/// step 3 瞬时地球坐标系 转换为 瞬时真天球坐标系
/// 将utc时间转换为格林尼治恒星时
/// step 3 ˲ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ ת<><D7AA>Ϊ ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ
/// <EFBFBD><EFBFBD>utcʱ<EFBFBD><EFBFBD>ת<EFBFBD><EFBFBD>Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>κ<EFBFBD><EFBFBD><EFBFBD>ʱ
/// xyz= ordinateSingleRotate('z',-gst_deg)*earthFixedXYZ;
/// UNTITLED 计算当地恒星时,返回值以时秒为单位
/// %dAT 润秒数
/// % TAI国际原子时间 该时间最准 TAI = UTC + dAT;% 国际原子时间
/// % TT 地球时 TT = TAI + 32.184 至2014年的时候;
/// % TDT 地球动力学时间
/// % ET 历书时间
/// % 地球时 = 地球动力学时间 = 历书时间
/// %地球时=地球动力学时间=历书时间
/// UNTITLED <EFBFBD><EFBFBD><EFBFBD>㵱<EFBFBD>غ<EFBFBD><EFBFBD><EFBFBD>ʱ,<2C><><EFBFBD><EFBFBD>ֵ<EFBFBD><D6B5>ʱ<EFBFBD><CAB1>Ϊ<EFBFBD><CEAA>λ
/// %dAT <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
/// % TAI<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ԭ<EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD> <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;
/// % TDT <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<EFBFBD><EFBFBD>
/// % ET <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>
/// % <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ = <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1> = <20><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1>
/// %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ=<3D><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1>=<3D><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1>
/// </summary>
/// <param name="UTC">格式为y m d 其中d的数值为小数，将h m s 按（sec/3600+min/60+h)/24转换成小数，并累加到day 上</param>
/// <param name="dUT1">dUT1 为ut1-utc 的差，数值不超过正负1秒，查iers可获得数值 0</param>
/// <param name="dAT">润秒数 37</param>
/// <param name="gst_deg">国际原子时间 该时间最准 TAI=UTC+dAT; %国际原子时间</param>
/// <param name="JDTDB">地球时 TT = TAI+32.184 至2014年的时候;</param>
/// <param name="UTC"><EFBFBD><EFBFBD>ʽΪ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="dAT"><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> 37</param>
/// <param name="gst_deg"><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ԭ<EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD> <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>
/// <returns></returns>
int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& gst_deg, double& JDTDB)
{
//function[gst_deg, JDTDB] = utc2gst(UTC, dUT1, dAT)
//% TDT 地球动力学时间
// % ET 历书时间
// % 地球时 = 地球动力学时间 = 历书时间
//% TDT <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<EFBFBD><EFBFBD>
// % ET <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD>
// % <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ = <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ѧʱ<D1A7><CAB1> = <20><><EFBFBD><EFBFBD>ʱ<EFBFBD><CAB1>
double J2000 = 2451545;
double JDutc = YMD2JD(UTC(0,0), UTC(0,1), UTC(0,2));
double JDUT1 = JDutc + dUT1 / 86400;// % UT1 为世界时，世界时由于自转的不均匀，因此与UTC时间有dUT1的差别，dUT1在各国卫星授时信号中会以0.1秒的精度给出IERS经过处理后，会以1e - 5的精度给出。
double JDUT1 = JDutc + dUT1 / 86400;// % UT1 Ϊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD>IJ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>UTCʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>dUT1<EFBFBD>IJ<EFBFBD><EFBFBD>dUT1<EFBFBD>ڸ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD>ź<EFBFBD><EFBFBD>л<EFBFBD><EFBFBD><EFBFBD>0.1<EFBFBD><EFBFBD>ľ<EFBFBD><EFBFBD>ȸ<EFBFBD><EFBFBD><EFBFBD>IERS<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>󣬻<EFBFBD><EFBFBD><EFBFBD>1e - 5<>ľ<EFBFBD><C4BE>ȸ<EFBFBD><C8B8><EFBFBD><EFBFBD><EFBFBD>
double dT = 32.184 + dAT - dUT1;
double JDTT = JDUT1 + dT / 86400;
//% JDTT = YMD2JD(UTC(1), UTC(2), UTC(3)) + dUT1 / 86400;
//% 首先计算TDB, TDB是质心动力学时，是太阳月球行星等天体星历表中的时间尺度
//% <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
double T = (JDTT - J2000) / 36525;
double temp = 0.001657 * sin(628.3076 * T + 6.2401)
+ 0.000022 * sin(575.3385 * T + 4.2970)
@ -1164,29 +1167,29 @@ int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& g
+ 0.00001 * T * sin(628.3076 * T + 4.2490);
JDTDB = JDTT + temp / 86400;
//% 下面计算UT2, UT2是在UT1的基础上修正地球周期性季节变化后得到的世界时间
// %% 根据UT1计算Tb, Tb为以贝塞尔年为单位，从历元B2000.0起算ff
//% <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
// %% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>UT1<EFBFBD><EFBFBD><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;
//% Tb = (YMD2JD(UT1(1), UT1(2), UT1(3)) - B2000) / 365.2422;
//% dTs = 0.022 * sin(2 * pi * Tb) - 0.012 * cos(2 * pi * Tb) - 0.006 * sin(4 * pi * Tb) + 0.007 * cos(4 * pi * Tb);
//% dTs = dTs / 86400;
//% UT2 = UT1 + [0 0 dTs];
//% 下面计算格林尼治平恒星时GMST; 单位为度
//% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱGMST; <20><>λΪ<CEBB><CEAA>
T = (JDTDB - J2000) / 36525;
double T2 = T * T;
double T3 = T2 * T;
double T4 = T3 * T;
double T5 = T4 * T;
double Du = JDUT1 - J2000;
double thita = 0.7790572732640 + 1.00273781191135448 * Du;//% 自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 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;
WGS84_J2000::nutationInLongitudeCaculate(JDTDB, epthilongA_deg, dertaPthi_deg, dertaEpthilong_deg, epthilong_deg);
//[epthilongA_deg, dertaPthi_deg] = nutationInLongitudeCaculate(JDTDB);
//% 计算赤经章动引起的误差eclipticObliquitygama，单位为角秒
//% 计算月球平近点角 太阳平近点角 月球纬度辐角 日月平角距(月球平黄经 - 太阳平黄经） 月球升交点平黄经 单位为角秒
//% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ྭ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>eclipticObliquitygama<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λΪ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//% <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>
double F_deg = (0.00000417 * T4 - 0.001037 * T3 - 12.7512 * T2 + 1739527262.8478 * T + 335779.526232) / 3600;
F_deg = fmod(F_deg, 360);
double D_deg = (-0.00003169 * T4 + 0.006593 * T3 - 6.3706 * T2 + 1602961601.2090 * T + 1072260.70369) / 3600;
@ -1210,15 +1213,15 @@ int WGS84_J2000::utc2gst(Eigen::MatrixXd UTC, double dUT1, double dAT, double& g
/// <summary>
/// step 4 瞬时真天球坐标系 转到瞬时平天球 坐标系
/// 计算平黄赤交角，黄经章动、和交角章动、瞬时黄赤交角。
/// step 4 ˲ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ ת<><D7AA>˲ʱƽ<CAB1><C6BD><EFBFBD><EFBFBD> <20><><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>
/// </summary>
/// <param name="JD">儒略日</param>
/// <param name="accuracy">表示计算的精度要求 数值为’normal‘ 和’high‘ 或者用'N'和’H'表示一般精度和高精度 。默认为高精度计算</param>
/// <param name="epthilongA_deg">平黄赤交角</param>
/// <param name="dertaPthi_deg">黄经章动</param>
/// <param name="dertaEpthilong_deg">和交角章动</param>
/// <param name="epthilong_deg">瞬时黄赤交角</param>
/// <param name="JD"><EFBFBD><EFBFBD><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="epthilongA_deg">ƽ<EFBFBD>Ƴཻ<EFBFBD><EFBFBD></param>
/// <param name="dertaPthi_deg"><EFBFBD>ƾ<EFBFBD><EFBFBD><EFBFBD></param>
/// <param name="dertaEpthilong_deg"><EFBFBD>ͽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD></param>
/// <param name="epthilong_deg">˲ʱ<EFBFBD>Ƴཻ<EFBFBD><EFBFBD></param>
/// <returns></returns>
int WGS84_J2000::nutationInLongitudeCaculate(double JD, double& epthilongA_deg, double& dertaPthi_deg, double& dertaEpthilong_deg, double& epthilong_deg)
{
@ -1227,9 +1230,9 @@ int WGS84_J2000::nutationInLongitudeCaculate(double JD, double& epthilongA_deg,
double T3 = T2 * T;
double T4 = T3 * T;
double T5 = T4 * T;
//% 计算月球平近点角lunarMeanAnomaly_deg(l_deg) 太阳平近点角SolarMeanAnomaly_deg(solarl_deg)
// % 月球纬度辐角lunarLatitudeAngle_deg(F_deg) 日月平角距diffLunarSolarElestialLongitude_deg(D_deg月球平黄经 - 太阳平黄经）
// % 月球升交点平黄经SolarAscendingNodeElestialLongitude_deg(Omiga_deg)
//% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>lunarMeanAnomaly_deg(l_deg) ̫<><CCAB>ƽ<EFBFBD><C6BD><EFBFBD><EFBFBD><EFBFBD>SolarMeanAnomaly_deg(solarl_deg)
// % <EFBFBD><EFBFBD><EFBFBD><EFBFBD>γ<EFBFBD>ȷ<EFBFBD><EFBFBD><EFBFBD>lunarLatitudeAngle_deg(F_deg) <20><><EFBFBD><EFBFBD>ƽ<EFBFBD>Ǿ<EFBFBD>diffLunarSolarElestialLongitude_deg(D_deg<65><67><EFBFBD><EFBFBD>ƽ<EFBFBD>ƾ<EFBFBD> - ̫<><CCAB>ƽ<EFBFBD>ƾ<EFBFBD><C6BE><EFBFBD>
// % <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD>ƾ<EFBFBD>SolarAscendingNodeElestialLongitude_deg(Omiga_deg)
double l_deg = (-0.00024470 * T4 + 0.051635 * T3 + 31.8792 * T2 + 1717915923.2178 * T + 485868.249036) / 3600;
l_deg = fmod(l_deg, 360);
double solarl_deg = (-0.00001149 * T4 + 0.000136 * T3 - 0.5532 * T2 + 129596581.0481 * T + 1287104.79305) / 3600;
@ -1244,7 +1247,7 @@ int WGS84_J2000::nutationInLongitudeCaculate(double JD, double& epthilongA_deg,
epthilongA_deg = (-0.0000000434 * T5 - 0.000000576 * T4 + 0.00200340 * T3 - 0.0001831 * T2 - 46.836769 * T + 84381.406) / 3600;
epthilongA_deg = epthilongA_deg - floor(epthilongA_deg / 360) * 360;
//% IAU2000模型有77项
//% IAU2000ģ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>77<EFBFBD><EFBFBD>
Eigen::MatrixXd elestialLonNutation = WGS84_J2000::IAU2000ModelParams;
dertaPthi_deg = -3.75e-08;
@ -1263,7 +1266,7 @@ int WGS84_J2000::nutationInLongitudeCaculate(double JD, double& epthilongA_deg,
}
/// <summary>
/// zA、thitaA、zetaA为赤道岁差角, 计算赤道岁差角
/// 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>
/// </summary>
/// <param name="JDTDB"></param>
/// <param name="zetaA"></param>
@ -1288,7 +1291,7 @@ int WGS84_J2000::precessionAngle(double JDTDB, double& zetaA, double& thitaA, do
}
/// <summary>
/// 同时 YMD2JD函数为 年月日转换为儒略日，
/// ͬʱ YMD2JD<4A><44><EFBFBD><EFBFBD>Ϊ <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ת<EFBFBD><D7AA>Ϊ<EFBFBD><CEAA><EFBFBD><EFBFBD><EFBFBD>գ<EFBFBD>
/// </summary>
/// <param name="yy"></param>
/// <param name="mm"></param>
@ -1306,7 +1309,7 @@ double WGS84_J2000::YMD2JD(double y, double m, double d)
/// <summary>
/// WGS84 ת J2000
/// </summary>
/// <param name="BLH_deg_m">WGS84 经纬度</param>
/// <param name="BLH_deg_m">WGS84 <EFBFBD><EFBFBD>γ<EFBFBD><EFBFBD></param>
/// <param name="UTC">1x3 Y,m,d</param>
/// <param name="xp">EARTH ORIENTATION PARAMETER Xp</param>
/// <param name="yp">EARTH ORIENTATION PARAMETER Yp</param>
@ -1315,21 +1318,21 @@ double WGS84_J2000::YMD2JD(double y, double m, double d)
/// <returns></returns>
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 转换到协议地球坐标系。
//step 1 step1 WGS 84 ת<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Э<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD>ￄ1<EFBFBD>7
Eigen::MatrixXd earthFixedXYZ = WGS84_J2000::WGS84TECEF(BLH_deg_m).transpose();
//step 2 协议地球坐标系 转换为瞬时地球坐标
// 这主要涉及查询 给定时刻的 xp,yp , 可以查询EOP文件获得，下载地址如下http://celestrak.com/spacedata/
//step 2 Э<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>τ1<EFBFBD>7 ת<><D7AA>Ϊ˲ʱ<CBB2><CAB1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD>漰<EFBFBD><EFBFBD>ѯ <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;
//step 3 瞬时地球坐标系 转换为 瞬时真天球坐标系
// 该步骤主要涉及到格林尼治恒星时角的计算，关于格林尼治恒星时角 计算方法 很多，下面给出 一个较为精确的计算方法。其中dut1 和dat参数在EOP文件中有。EOP文件下载地址如下 http://celestrak.org/spacedata/
//step 3 ˲ʱ<EFBFBD><EFBFBD><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/
double gst_deg, JDTDB;
WGS84_J2000::utc2gst(UTC, dut1, dat, gst_deg, JDTDB);
Eigen::MatrixXd xyz = ordinateSingleRotate(3, -1 * gst_deg) * earthFixedXYZ;
//step 4 瞬时真天球坐标系 转到瞬时平天球 坐标系
//step 4 ˲ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ ת<><D7AA>˲ʱƽ<CAB1><C6BD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD>ϵ
double epthilongA_deg, dertaPthi_deg, dertaEpthilong_deg, epthilong_deg;
WGS84_J2000::nutationInLongitudeCaculate(JDTDB, epthilongA_deg, dertaPthi_deg, dertaEpthilong_deg, epthilong_deg);
xyz = ordinateSingleRotate(1, -epthilongA_deg) * ordinateSingleRotate(3, dertaPthi_deg) * ordinateSingleRotate(1, dertaEpthilong_deg + epthilongA_deg) * xyz;
//step5 瞬时平天球坐标系转换为协议天球坐标系（J2000）
//step5 ˲ʱƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵת<EFBFBD><EFBFBD>ΪЭ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ϵ<EFBFBD><EFBFBD>J2000<EFBFBD><EFBFBD>
double zetaA, thitaA, zA;
WGS84_J2000::precessionAngle(JDTDB, zetaA, thitaA, zA);
xyz = ordinateSingleRotate(3, zetaA) * ordinateSingleRotate(2, -thitaA) * ordinateSingleRotate(3, zA) * xyz;
@ -1345,4 +1348,7 @@ Landpoint WGS84_J2000::WGS842J2000(Landpoint p, Eigen::MatrixXd UTC, double xp,
return Landpoint{ blh(0,0),blh(0,1) ,blh(0,0) };
//return Landpoint();
}
*/
*/

142
baseTool.h
View File

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

1
interpolation.cpp Normal file
View File

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

10
interpolation.h Normal file
View File

@ -0,0 +1,10 @@
#pragma once
/**
专门用于插值计算的类库
*/
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

423
linterp.h Normal file
View File

@ -0,0 +1,423 @@
//
// Copyright (c) 2012 Ronaldo Carpio
//
// Permission to use, copy, modify, distribute and sell this software
// and its documentation for any purpose is hereby granted without fee,
// provided that the above copyright notice appear in all copies and
// that both that copyright notice and this permission notice appear
// in supporting documentation. The authors make no representations
// about the suitability of this software for any purpose.
// It is provided "as is" without express or implied warranty.
//
/*
This is a C++ header-only library for N-dimensional linear interpolation on a rectangular grid. Implements two methods:
* Multilinear: Interpolate using the N-dimensional hypercube containing the point. Interpolation step is O(2^N)
* Simplicial: Interpolate using the N-dimensional simplex containing the point. Interpolation step is O(N log N), but less accurate.
Requires boost/multi_array library.
For a description of the algorithms, see:
* Weiser & Zarantonello (1988), "A Note on Piecewise Linear and Multilinear Table Interpolation in Many Dimensions", _Mathematics of Computation_ 50 (181), p. 189-196
* Davies (1996), "Multidimensional Triangulation and Interpolation for Reinforcement Learning", _Proceedings of Neural Information Processing Systems 1996_
*/
#ifndef _linterp_h
#define _linterp_h
#include <assert.h>
#include <math.h>
#include <stdarg.h>
#include <float.h>
#include <cstdarg>
#include <string>
#include <vector>
#include <array>
#include <functional>
#include <boost/multi_array.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/storage.hpp>
using std::vector;
using std::array;
typedef unsigned int uint;
typedef vector<int> iVec;
typedef vector<double> dVec;
// TODO:
// - specify behavior past grid boundaries.
// 1) clamp
// 2) return a pre-determined value (e.g. NaN)
// compile-time params:
// 1) number of dimensions
// 2) scalar type T
// 3) copy data or not (default: false). The grids will always be copied
// 4) ref count class (default: none)
// 5) continuous or not
// run-time constructor params:
// 1) f
// 2) grids
// 3) behavior outside grid: default=clamp
// 4) value to return outside grid, defaut=nan
struct EmptyClass {};
template <int N, class T, bool CopyData = true, bool Continuous = true, class ArrayRefCountT = EmptyClass, class GridRefCountT = EmptyClass>
class NDInterpolator {
public:
typedef T value_type;
typedef ArrayRefCountT array_ref_count_type;
typedef GridRefCountT grid_ref_count_type;
static const int m_N = N;
static const bool m_bCopyData = CopyData;
static const bool m_bContinuous = Continuous;
typedef boost::numeric::ublas::array_adaptor<T> grid_type;
typedef boost::const_multi_array_ref<T, N> array_type;
typedef std::unique_ptr<array_type> array_type_ptr;
array_type_ptr m_pF;
ArrayRefCountT m_ref_F; // reference count for m_pF
vector<T> m_F_copy; // if CopyData == true, this holds the copy of F
vector<grid_type> m_grid_list;
vector<GridRefCountT> m_grid_ref_list; // reference counts for grids
vector<vector<T> > m_grid_copy_list; // if CopyData == true, this holds the copies of the grids
// constructors assume that [f_begin, f_end) is a contiguous array in C-order
// non ref-counted constructor.
template <class IterT1, class IterT2, class IterT3>
NDInterpolator(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end) {
init(grids_begin, grids_len_begin, f_begin, f_end);
}
// ref-counted constructor
template <class IterT1, class IterT2, class IterT3, class RefCountIterT>
NDInterpolator(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end, ArrayRefCountT &refF, RefCountIterT grid_refs_begin) {
init_refcount(grids_begin, grids_len_begin, f_begin, f_end, refF, grid_refs_begin);
}
template <class IterT1, class IterT2, class IterT3>
void init(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end) {
set_grids(grids_begin, grids_len_begin, m_bCopyData);
set_f_array(f_begin, f_end, m_bCopyData);
}
template <class IterT1, class IterT2, class IterT3, class RefCountIterT>
void init_refcount(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end, ArrayRefCountT &refF, RefCountIterT grid_refs_begin) {
set_grids(grids_begin, grids_len_begin, m_bCopyData);
set_grids_refcount(grid_refs_begin, grid_refs_begin + N);
set_f_array(f_begin, f_end, m_bCopyData);
set_f_refcount(refF);
}
template <class IterT1, class IterT2>
void set_grids(IterT1 grids_begin, IterT2 grids_len_begin, bool bCopy) {
m_grid_list.clear();
m_grid_ref_list.clear();
m_grid_copy_list.clear();
for (int i=0; i<N; i++) {
int gridLength = grids_len_begin[i];
if (bCopy == false) {
T const *grid_ptr = &(*grids_begin[i]);
m_grid_list.push_back(grid_type(gridLength, (T*) grid_ptr )); // use the given pointer
} else {
m_grid_copy_list.push_back( vector<T>(grids_begin[i], grids_begin[i] + grids_len_begin[i]) ); // make our own copy of the grid
T *begin = &(m_grid_copy_list[i][0]);
m_grid_list.push_back(grid_type(gridLength, begin)); // use our copy
}
}
}
template <class IterT1, class RefCountIterT>
void set_grids_refcount(RefCountIterT refs_begin, RefCountIterT refs_end) {
assert(refs_end - refs_begin == N);
m_grid_ref_list.assign(refs_begin, refs_begin + N);
}
// assumes that [f_begin, f_end) is a contiguous array in C-order
template <class IterT>
void set_f_array(IterT f_begin, IterT f_end, bool bCopy) {
unsigned int nGridPoints = 1;
array<int,N> sizes;
for (unsigned int i=0; i<m_grid_list.size(); i++) {
sizes[i] = m_grid_list[i].size();
nGridPoints *= sizes[i];
}
int f_len = f_end - f_begin;
if ( (m_bContinuous && f_len != nGridPoints) || (!m_bContinuous && f_len != 2 * nGridPoints) ) {
throw std::invalid_argument("f has wrong size");
}
for (unsigned int i=0; i<m_grid_list.size(); i++) {
if (!m_bContinuous) { sizes[i] *= 2; }
}
m_F_copy.clear();
if (bCopy == false) {
m_pF.reset(new array_type(f_begin, sizes));
} else {
m_F_copy = vector<T>(f_begin, f_end);
m_pF.reset(new array_type(&m_F_copy[0], sizes));
}
}
void set_f_refcount(ArrayRefCountT &refF) {
m_ref_F = refF;
}
// -1 is before the first grid point
// N-1 (where grid.size() == N) is after the last grid point
int find_cell(int dim, T x) const {
grid_type const &grid(m_grid_list[dim]);
if (x < *(grid.begin())) return -1;
else if (x >= *(grid.end()-1)) return grid.size()-1;
else {
auto i_upper = std::upper_bound(grid.begin(), grid.end(), x);
return i_upper - grid.begin() - 1;
}
}
// return the value of f at the given cell and vertex
T get_f_val(array<int,N> const &cell_index, array<int,N> const &v_index) const {
array<int,N> f_index;
if (m_bContinuous) {
for (int i=0; i<N; i++) {
if (cell_index[i] < 0) {
f_index[i] = 0;
} else if (cell_index[i] >= m_grid_list[i].size()-1) {
f_index[i] = m_grid_list[i].size()-1;
} else {
f_index[i] = cell_index[i] + v_index[i];
}
}
} else {
for (int i=0; i<N; i++) {
if (cell_index[i] < 0) {
f_index[i] = 0;
} else if (cell_index[i] >= m_grid_list[i].size()-1) {
f_index[i] = (2*m_grid_list[i].size())-1;
} else {
f_index[i] = 1 + (2*cell_index[i]) + v_index[i];
}
}
}
return (*m_pF)(f_index);
}
T get_f_val(array<int,N> const &cell_index, int v) const {
array<int,N> v_index;
for (int dim=0; dim<N; dim++) {
v_index[dim] = (v >> (N-dim-1)) & 1; // test if the i-th bit is set
}
return get_f_val(cell_index, v_index);
}
};
template <int N, class T, bool CopyData = true, bool Continuous = true, class ArrayRefCountT = EmptyClass, class GridRefCountT = EmptyClass>
class InterpSimplex : public NDInterpolator<N,T,CopyData,Continuous,ArrayRefCountT,GridRefCountT> {
public:
typedef NDInterpolator<N,T,CopyData,Continuous,ArrayRefCountT,GridRefCountT> super;
template <class IterT1, class IterT2, class IterT3>
InterpSimplex(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end)
: super(grids_begin, grids_len_begin, f_begin, f_end)
{}
template <class IterT1, class IterT2, class IterT3, class RefCountIterT>
InterpSimplex(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end, ArrayRefCountT &refF, RefCountIterT ref_begins)
: super(grids_begin, grids_len_begin, f_begin, f_end, refF, ref_begins)
{}
template <class IterT>
T interp(IterT x_begin) const {
array<T,1> result;
array< array<T,1>, N > coord_iter;
for (int i=0; i<N; i++) {
coord_iter[i][0] = x_begin[i];
}
interp_vec(1, coord_iter.begin(), coord_iter.end(), result.begin());
return result[0];
}
template <class IterT1, class IterT2>
void interp_vec(int n, IterT1 coord_iter_begin, IterT1 coord_iter_end, IterT2 i_result) const {
assert(N == coord_iter_end - coord_iter_begin);
array<int,N> cell_index, v_index;
array<std::pair<T, int>,N> xipair;
int c;
T y, v0, v1;
//mexPrintf("%d\n", n);
for (int i=0; i<n; i++) { // for each point
for (int dim=0; dim<N; dim++) {
typename super::grid_type const &grid(super::m_grid_list[dim]);
c = this->find_cell(dim, coord_iter_begin[dim][i]);
//mexPrintf("%d\n", c);
if (c == -1) { // before first grid point
y = 1.0;
} else if (c == grid.size()-1) { // after last grid point
y = 0.0;
} else {
//mexPrintf("%f %f\n", grid[c], grid[c+1]);
y = (coord_iter_begin[dim][i] - grid[c]) / (grid[c + 1] - grid[c]);
if (y < 0.0) y=0.0;
else if (y > 1.0) y=1.0;
}
xipair[dim].first = y;
xipair[dim].second = dim;
cell_index[dim] = c;
}
// sort xi's and get the permutation
std::sort(xipair.begin(), xipair.end(), [](std::pair<T, int> const &a, std::pair<T, int> const &b) {
return (a.first < b.first);
});
// walk the vertices of the simplex determined by the permutation
for (int j=0; j<N; j++) {
v_index[j] = 1;
}
v0 = this->get_f_val(cell_index, v_index);
y = v0;
for (int j=0; j<N; j++) {
v_index[xipair[j].second]--;
v1 = this->get_f_val(cell_index, v_index);
y += (1.0 - xipair[j].first) * (v1-v0); // interpolate
v0 = v1;
}
*i_result++ = y;
}
}
};
template <int N, class T, bool CopyData = true, bool Continuous = true, class ArrayRefCountT = EmptyClass, class GridRefCountT = EmptyClass>
class InterpMultilinear : public NDInterpolator<N,T,CopyData,Continuous,ArrayRefCountT,GridRefCountT> {
public:
typedef NDInterpolator<N,T,CopyData,Continuous,ArrayRefCountT,GridRefCountT> super;
template <class IterT1, class IterT2, class IterT3>
InterpMultilinear(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end)
: super(grids_begin, grids_len_begin, f_begin, f_end)
{}
template <class IterT1, class IterT2, class IterT3, class RefCountIterT>
InterpMultilinear(IterT1 grids_begin, IterT2 grids_len_begin, IterT3 f_begin, IterT3 f_end, ArrayRefCountT &refF, RefCountIterT ref_begins)
: super(grids_begin, grids_len_begin, f_begin, f_end, refF, ref_begins)
{}
template <class IterT1, class IterT2>
static T linterp_nd_unitcube(IterT1 f_begin, IterT1 f_end, IterT2 xi_begin, IterT2 xi_end) {
int n = xi_end - xi_begin;
int f_len = f_end - f_begin;
assert(1 << n == f_len);
T sub_lower, sub_upper;
if (n == 1) {
sub_lower = f_begin[0];
sub_upper = f_begin[1];
} else {
sub_lower = linterp_nd_unitcube(f_begin, f_begin + (f_len/2), xi_begin + 1, xi_end);
sub_upper = linterp_nd_unitcube(f_begin + (f_len/2), f_end, xi_begin + 1, xi_end);
}
T result = sub_lower + (*xi_begin)*(sub_upper - sub_lower);
return result;
}
template <class IterT>
T interp(IterT x_begin) const {
array<T,1> result;
array< array<T,1>, N > coord_iter;
for (int i=0; i<N; i++) {
coord_iter[i][0] = x_begin[i];
}
interp_vec(1, coord_iter.begin(), coord_iter.end(), result.begin());
return result[0];
}
template <class IterT1, class IterT2>
void interp_vec(int n, IterT1 coord_iter_begin, IterT1 coord_iter_end, IterT2 i_result) const {
assert(N == coord_iter_end - coord_iter_begin);
array<int,N> index;
int c;
T y, xi;
vector<T> f(1 << N);
array<T,N> x;
for (int i=0; i<n; i++) { // loop over each point
for (int dim=0; dim<N; dim++) { // loop over each dimension
auto const &grid(super::m_grid_list[dim]);
xi = coord_iter_begin[dim][i];
c = this->find_cell(dim, coord_iter_begin[dim][i]);
if (c == -1) { // before first grid point
y = 1.0;
} else if (c == grid.size()-1) { // after last grid point
y = 0.0;
} else {
y = (coord_iter_begin[dim][i] - grid[c]) / (grid[c + 1] - grid[c]);
if (y < 0.0) y=0.0;
else if (y > 1.0) y=1.0;
}
index[dim] = c;
x[dim] = y;
}
// copy f values at vertices
for (int v=0; v < (1 << N); v++) { // loop over each vertex of hypercube
f[v] = this->get_f_val(index, v);
}
*i_result++ = linterp_nd_unitcube(f.begin(), f.end(), x.begin(), x.end());
}
}
};
typedef InterpSimplex<1,double> NDInterpolator_1_S;
typedef InterpSimplex<2,double> NDInterpolator_2_S;
typedef InterpSimplex<3,double> NDInterpolator_3_S;
typedef InterpSimplex<4,double> NDInterpolator_4_S;
typedef InterpSimplex<5,double> NDInterpolator_5_S;
typedef InterpMultilinear<1,double> NDInterpolator_1_ML;
typedef InterpMultilinear<2,double> NDInterpolator_2_ML;
typedef InterpMultilinear<3,double> NDInterpolator_3_ML;
typedef InterpMultilinear<4,double> NDInterpolator_4_ML;
typedef InterpMultilinear<5,double> NDInterpolator_5_ML;
// C interface
extern "C" {
void linterp_simplex_1(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult);
void linterp_simplex_2(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult);
void linterp_simplex_3(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult);
}
void inline linterp_simplex_1(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult) {
const int N=1;
size_t total_size = 1;
for (int i=0; i<N; i++) {
total_size *= grid_len_begin[i];
}
InterpSimplex<N, double, false> interp_obj(grids_begin, grid_len_begin, pF, pF + total_size);
interp_obj.interp_vec(xi_len, xi_begin, xi_begin + N, pResult);
}
void inline linterp_simplex_2(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult) {
const int N=2;
size_t total_size = 1;
for (int i=0; i<N; i++) {
total_size *= grid_len_begin[i];
}
InterpSimplex<N, double, false> interp_obj(grids_begin, grid_len_begin, pF, pF + total_size);
interp_obj.interp_vec(xi_len, xi_begin, xi_begin + N, pResult);
}
void inline linterp_simplex_3(double **grids_begin, int *grid_len_begin, double *pF, int xi_len, double **xi_begin, double *pResult) {
const int N=3;
size_t total_size = 1;
for (int i=0; i<N; i++) {
total_size *= grid_len_begin[i];
}
InterpSimplex<N, double, false> interp_obj(grids_begin, grid_len_begin, pF, pF + total_size);
interp_obj.interp_vec(xi_len, xi_begin, xi_begin + N, pResult);
}
#endif //_linterp_h

4
resource.h
View File

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

763
simptsn.cpp
View File

File diff suppressed because it is too large Load Diff

155
simptsn.h
View File

@ -1,14 +1,13 @@
#pragma once
///
/// 仿真成像算法
/// 仿真成像算法
///
//#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
// 本地方法
// 本地方法
#include "baseTool.h"
#include <iostream>
#include <armadillo>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
@ -18,68 +17,66 @@
#include "ImageMatch.h"
using namespace std;
using namespace Eigen;
using namespace arma;
///////////// ptsn 算法 /////////////////////
///////////// ptsn 算法 /////////////////////
class PSTNAlgorithm {
public:
// 初始化与析构函数
// 初始化与析构函数
PSTNAlgorithm();
PSTNAlgorithm(string parameterPath);
~PSTNAlgorithm();
int dem2lat_lon(std::string dem_path, std::string lon_path, std::string lat_path);
int dem2SAR_RC(std::string dem_path, std::string sim_rc_path);
// 模拟算法函数
Eigen::MatrixX<double> calNumericalDopplerValue(Eigen::MatrixX<double> R);// 数值模拟法计算多普勒频移值
double calNumericalDopplerValue(double R);// 数值模拟法计算多普勒频移值
Eigen::MatrixX<double> calTheoryDopplerValue(Eigen::MatrixX<double> R, Eigen::MatrixX<double> R_s1, Eigen::MatrixX<double> V_s1);//根据理论模型计算多普勒频移值
// 模拟算法函数
Eigen::MatrixX<double> calNumericalDopplerValue(Eigen::MatrixX<double> R);// 数值模拟法计算多普勒频移值
double calNumericalDopplerValue(double R);// 数值模拟法计算多普勒频移值
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);
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);
//Landpoint WGS842J2000(Landpoint p);
public: // WGS84 到 J2000 之间的变换参数
public: // WGS84 到 J2000 之间的变换参数
Eigen::MatrixXd UTC;
double Xp = 0.204071;
double Yp = 0.318663;
double Dut1 = 0.0366742;
double Dat = 37;
public:
// std::string dem_path; // 模拟用的DEM影像
// std::string lon_path;// 经度影像
// std::string lat_path;// 纬度影像
// std::string dem_path; // 模拟用的DEM影像
// std::string lon_path;// 经度影像
// std::string lat_path;// 纬度影像
// std::string simsar_path;// 仿真影像坐标
int height; // 影像的高
// std::string simsar_path;// 仿真影像坐标
int height; // 影像的高
int width;
// 入射角
double near_in_angle;// 近入射角
double far_in_angle;// 远入射角
// 入射角
double near_in_angle;// 近入射角
double far_in_angle;// 远入射角
// SAR的成像参数
double widthspace;// 距离向分辨率
double R0;//近斜距
double LightSpeed; // 光速
double lamda;// 波长
double refrange;// 参考斜距
// SAR的成像参数
double widthspace;// 距离向分辨率
double R0;//近斜距
double LightSpeed; // 光速
double lamda;// 波长
double refrange;// 参考斜距
double imgStartTime;// 成像起始时间
double PRF;// 脉冲重复率
double imgStartTime;// 成像起始时间
double PRF;// 脉冲重复率
int doppler_paramenter_number;//多普勒系数个数
MatrixX<double> doppler_paras;// 多普勒系数
int doppler_paramenter_number;//多普勒系数个数
MatrixX<double> doppler_paras;// 多普勒系数
OrbitPoly orbit; // 轨道模型
OrbitPoly orbit; // 轨道模型
};
/// <summary>
/// 获取局地入射角，角度值
/// 获取局地入射角，角度值
/// </summary>
/// <param name="satepoint"></param>
/// <param name="landpoint"></param>
@ -89,7 +86,7 @@ double getlocalIncAngle(Landpoint& satepoint, Landpoint& landpoint, Landpoint& s
/// <summary>
/// 侧视入射角,角度值
/// 侧视入射角,角度值
/// </summary>
/// <param name="satepoint"></param>
/// <param name="landpoint"></param>
@ -98,12 +95,12 @@ double getIncAngle(Landpoint& satepoint, Landpoint& landpoint);
/// <summary>
/// ASF计算方法
/// ASF计算方法
/// </summary>
class ASFOrthClass {
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);
double GetLookFromRangeYaw(double R, double alpha, double beta, Eigen::Vector3d SatellitePosition, Eigen::Vector3d SatelliteVelocity, double R_threshold, double H = 0);
@ -117,7 +114,7 @@ public:
};
/// <summary>
/// 仿真处理流程
/// 仿真处理流程
/// </summary>
class simProcess {
public:
@ -126,11 +123,10 @@ public:
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 InitASFSAR(std::string paras_path, std::string workspace_path, std::string out_dir_path, std::string in_dem_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 dem2SAR(); // 切换行号
int dem2SAR(); // 切换行号
int SARIncidentAngle();
int SARSimulationWGS();
int SARSimulation();
@ -138,33 +134,34 @@ public:
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 dem2SAR_row(); // 切换行号
int dem2SAR_row(); // 切换行号
int SARIncidentAngle_RPC();
int createRPCDEM();
// 格网离散插值
// 格网离散插值
int Scatter2Grid(std::string lon_lat_path, std::string data_tiff, std::string grid_path, double space);
// ASF 模块计算
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 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 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);
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 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 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 ASF(std::string in_dem_path, std::string out_latlon_path, std::string in_sar_rc_path, ASFOrthClass asfclass, PSTNAlgorithm PSTN);
void calcalIncident_localIncident_angle(std::string in_dem_path, std::string in_rc_wgs84_path, std::string out_incident_angle_path, std::string out_local_incident_angle_path, PSTNAlgorithm PSTN);
void calGEC_Incident_localIncident_angle(std::string in_dem_path, std::string in_gec_lon_lat_path,std::string out_incident_angle_path, std::string out_local_incident_angle_path,PSTNAlgorithm PSTN);
@ -173,60 +170,64 @@ public:
// RPC
void CreateRPC_DEM(std::string in_rpc_rc_path, std::string in_dem_path, std::string out_rpc_dem_path);
void CreateRPC_refrenceTable(string in_rpc_tiff_path, string out_rpc_lon_lat_tiff_path);
void CreateRPC_refrenceTable(string in_rpc_tiff_path, string in_merge_dem, string out_rpc_lon_lat_tiff_path);
// 内部参数
void CorrectionFromSLC2Geo(string in_lon_lat_path, string in_radar_path, string out_radar_path, int in_band_ids);
// 内部参数
public:
ImageMatch matchmodel;
bool isMatchModel;
std::string workSpace_path;// 工作空间路径
std::string outSpace_path;// 输出工作空间路径
PSTNAlgorithm pstn; // 参数组件
///// RPC 局地入射角 ///////////////////////////////////
std::string workSpace_path;// 工作空间路径
std::string outSpace_path;// 输出工作空间路径
PSTNAlgorithm pstn; // 参数组件
///// RPC 局地入射角 ///////////////////////////////////
std::string in_rpc_lon_lat_path;
std::string out_dir_path;
std::string workspace_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_local_inc_angle_rpc_path;
// 输出文件
// 输出文件
std::string out_inc_angle_geo_path;
std::string out_local_inc_angle_geo_path;
///// 正向仿真方法 /////////
///// 正向仿真方法 /////////
// 临时产品文件
std::string in_dem_path; // 输入DEM
std::string dem_path; // 重采样之后DEM
std::string dem_r_path;// 行号影像
std::string dem_rc_path; // 行列号
// 临时产品文件
std::string in_dem_path; // 输入DEM
std::string dem_path; // 重采样之后DEM
std::string dem_r_path;// 行号影像
std::string dem_rc_path; // 行列号
std::string sar_sim_wgs_path;
std::string sar_sim_path;
std::string in_sar_path;
std::string sar_jpg_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_plant_slantRange_path;// 平面斜距文件
std::string out_dem_slantRange_path; // 斜距文件
std::string out_plant_slantRange_path;// 平面斜距文件
//// ASF 方法 /////////////
std::string out_lon_path;// 经度
std::string out_lat_path;// 纬度
std::string out_slantRange_path;// 斜距文件
//// ASF 方法 /////////////
std::string out_lon_path;// 经度
std::string out_lat_path;// 纬度
std::string out_slantRange_path;// 斜距文件
/// 共同文件 ///
std::string out_localIncidenct_path;// 计算局地入射角
std::string out_incidence_path; // 入射角文件
std::string out_ori_sim_tiff;// 映射角文件
/// 共同文件 ///
std::string out_localIncidenct_path;// 计算局地入射角
std::string out_incidence_path; // 入射角文件
std::string out_ori_sim_tiff;// 映射角文件'
std::string out_sim_ori_tiff;
};
bool PtInRect(Point_3d pCur, Point_3d pLeftTop, Point_3d pRightTop, Point_3d pRightBelow, Point_3d pLeftBelow);

50
test_moudel.cpp
View File

@ -38,12 +38,13 @@ int test_main(int mode, string rootpath)
break;
}
case 8: {
std::string in_rpc_tiff = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\GF3_MYN_QPSI_011437_E99.2_N28.6_20181012_L1B_HH_L10003514912_db.tif";
std::string out_lon_lat_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\RPC_ori_sim.tif";
std::string in_rpc_tiff = "D:\\MicroSAR\\C-SAR\\MicroWorkspace\\Ortho\\Temporary\\GF3_MYN_QPSI_011437_E99.2_N28.6_20181012_L1B_HH_L10003514912_db.tif";
std::string out_lon_lat_path = "D:\\MicroSAR\\C-SAR\\MicroWorkspace\\Ortho\\Temporary\\RPC_ori_sim.tif";
std::string dem_path = "D:\\MicroSAR\\C-SAR\\MicroWorkspace\\Ortho\\Temporary\\TestDEM\\mergedDEM.tif";
simProcess process;
//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::cout << "==========================================================================" << endl;
process.CreateRPC_refrenceTable(in_rpc_tiff, out_lon_lat_path);
process.CreateRPC_refrenceTable(in_rpc_tiff, dem_path,out_lon_lat_path);
std::cout << "==========================================================================" << endl;
break;
}
@ -87,15 +88,15 @@ int test_mode_0()
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
//输入参数
//<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
std::string root_path = rootpath;// "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary2";
std::string parameter_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\orth_para.txt"; //
std::string dem_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\TestDEM\\mergedDEM.tif"; //
std::string in_sar_path = "D:\\MicroSAR\\C-SAR\\Ortho\\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::string work_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary"; //
std::string taget_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package"; //
std::string workspace = "D:\\MicroSAR\\C-SAR\\MicroWorkspace\\Ortho\\Temporary";
std::string parameter_path =workspace+ "\\package\\orth_para.txt"; //
std::string dem_path = workspace + "\\TestDEM\\mergedDEM.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 work_path = workspace ; //
std::string taget_path = workspace + "\\package"; //
std::cout << "==========================================================================" << endl;
std::cout << "in parameters:========================================================" << endl;
std::cout << "parameters file path:\t" << parameter_path << endl;
@ -107,29 +108,15 @@ int test_mode_1(string rootpath)
std::cout << "==========================================================================" << endl;
process.InitSimulationSAR(parameter_path, work_path, taget_path, dem_path, in_sar_path);
std::cout << "==========================================================================" << endl;
////// 模拟计算
//process.dem2SAR_rc(dem_path, out_simrc_path, this->dem_path, sim_sar_orth_path, pstn);
//exit(1);
//process.sim_SAR(this->dem_path, out_simrc_path, sim_sar_orth_path, pstn);
//process.sim_sum_SAR(this->dem_path, out_simrc_path, sim_sar_orth_path, sim_sar_path, pstn);
//
////// 地形结果
//process.dem2aspect_slope(this->dem_path, dem_aspect_path, dem_slope_path);
////// 原始图计算
//process.ori_sar_power(ori_sar_path, ori_sar_rsc_path);
//// 创建匹配模型
//process.createImageMatchModel(ori_sar_rsc_path, ori_sar_rsc_jpg_path, sim_sar_path, sim_sar_jpg_path, matchmodel_path);
return 0;
}
int test_ASF(string rootpath) {
std::string root_path = rootpath;// "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary2";
std::string parameter_path = root_path + "\\package\\orth_para.txt";
std::string dem_path = root_path + "\\TestDEM\\mergedDEM.tif";
std::string ori_sar_path = root_path + "\\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::string ori_sar_path = root_path + "\\unpack\\" + "GF3_KAS_FSI_020253_E110.8_N25.5_20200614_L1A_HHHV_L10004871459\\GF3_KAS_FSI_020253_E110.8_N25.5_20200614_L1A_HH_L10003923848.tiff";
std::string work_path = root_path + "\\TestSim";
std::string taget_path = root_path + "\\output";
std::string out_GEC_dem_path = root_path + "\\output\\GEC_dem.tif";
@ -137,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_clip_DEM_path = root_path + "\\output\\Orth_dem.tif";
std::cout << "==========================================================================" << endl;
std::cout << "预处理计算结果，可以计算出DEM 范围" << 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 << "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 << "==========================================================================" << endl;
@ -152,8 +139,7 @@ int test_ASF(string rootpath) {
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
std::cout << "==========================================================================" << endl;
// 拼接参数
//
dem_path = out_clip_DEM_path;// JoinPath(work_path, "clipDEM.tif");
std::string out_simrc_path = JoinPath(work_path, "dem_sar_rc.tif"); //r,c
std::string sim_sar_orth_path = JoinPath(work_path, "sim_sar_orth.tif"); // orth_sim,orth_inclocal
@ -169,13 +155,11 @@ int test_ASF(string rootpath) {
std::string out_dem_count_path = JoinPath(work_path, "dem_count.tif");
// 输出结果
std::string out_dem_path = out_GEC_dem_path;
std::string out_inclocal_path = JoinPath(taget_path, "inclocal.tif");
std::string out_sar_rc_path = out_GTC_rc_path;
std::string out_lon_lat_path = out_GEC_lon_lat_path;
// 校正正射坐标计算结果
if (process.isMatchModel) {
process.correctOrth_rc(out_simrc_path, process.matchmodel);
}
@ -184,10 +168,10 @@ int test_ASF(string rootpath) {
}
boost::filesystem::copy_file(boost::filesystem::path(out_simrc_path), boost::filesystem::path(out_sar_rc_path));
//// 计算影像对应的高程
process.correct_ati(out_sar_rc_path, dem_path, out_inclocal_path, out_dem_path, out_dem_count_path, pstn);
// ASF
process.ASF(out_dem_path, out_lon_lat_path, out_sar_rc_path, ASFClass, pstn);
process.ASF(out_dem_path, out_lon_lat_path, ASFClass, pstn);
return 0;
}
@ -293,6 +277,6 @@ int test_LLA2XYZ_XYZ2LLA()
int test_vector_operator()
{
// 创建
// <EFBFBD><EFBFBD><EFBFBD><EFBFBD>
return 0;
}