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

入所定稿版本，功能包括：ASF、间接定位法、插值法

master
czh 2022-12-07 17:12:34 +08:00
commit 6872c62bce
28 changed files with 7567 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
.gitattributes vendored Normal file
View File

@ -0,0 +1,63 @@
###############################################################################
# Set default behavior to automatically normalize line endings.
###############################################################################
* text=auto
###############################################################################
# Set default behavior for command prompt diff.
#
# This is need for earlier builds of msysgit that does not have it on by
# default for csharp files.
# Note: This is only used by command line
###############################################################################
#*.cs diff=csharp
###############################################################################
# Set the merge driver for project and solution files
#
# Merging from the command prompt will add diff markers to the files if there
# are conflicts (Merging from VS is not affected by the settings below, in VS
# the diff markers are never inserted). Diff markers may cause the following
# file extensions to fail to load in VS. An alternative would be to treat
# these files as binary and thus will always conflict and require user
# intervention with every merge. To do so, just uncomment the entries below
###############################################################################
#*.sln merge=binary
#*.csproj merge=binary
#*.vbproj merge=binary
#*.vcxproj merge=binary
#*.vcproj merge=binary
#*.dbproj merge=binary
#*.fsproj merge=binary
#*.lsproj merge=binary
#*.wixproj merge=binary
#*.modelproj merge=binary
#*.sqlproj merge=binary
#*.wwaproj merge=binary
###############################################################################
# behavior for image files
#
# image files are treated as binary by default.
###############################################################################
#*.jpg binary
#*.png binary
#*.gif binary
###############################################################################
# diff behavior for common document formats
#
# Convert binary document formats to text before diffing them. This feature
# is only available from the command line. Turn it on by uncommenting the
# entries below.
###############################################################################
#*.doc diff=astextplain
#*.DOC diff=astextplain
#*.docx diff=astextplain
#*.DOCX diff=astextplain
#*.dot diff=astextplain
#*.DOT diff=astextplain
#*.pdf diff=astextplain
#*.PDF diff=astextplain
#*.rtf diff=astextplain
#*.RTF diff=astextplain

56
.gitignore vendored Normal file
View File

@ -0,0 +1,56 @@
# Prerequisites
*.d
# Compiled Object files
*.slo
*.lo
*.o
*.obj
# Precompiled Headers
*.gch
*.pch
# Compiled Dynamic libraries
*.so
*.dylib
*.dll
# Fortran module files
*.mod
*.smod
# Compiled Static libraries
*.lai
*.la
*.a
*.lib
# Executables
*.exe
*.out
*.app
*.dll
x64/
x64/*
.vs/
.vs/*
/x64/*
/.vs/*
./x64/*
./.vs/*
./x64/*
/x64/*
*.ipch
*.db
*.pdb
*.tlog
*.log
*.pdb
*.db
*.tiff
*.tif
*.jpg
Temporary*/

510
ImageMatch.cpp Normal file
View File

@ -0,0 +1,510 @@
#pragma once
//#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
// 本地方法
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#include <omp.h>
#include "baseTool.h"
#include "simptsn.h"
#include "SateOrbit.h"
#include <gdal_utils.h>
#include <opencv2/core.hpp>
#include <opencv2/opencv.hpp>
#include <opencv2/core/eigen.hpp>
#include "ImageMatch.h"
#include <opencv2/core/eigen.hpp>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace Eigen;
using namespace cv;
/// <summary>
/// cpnvert gdal to jpg
/// </summary>
/// <param name="gdal_path"></param>
/// <param name="jpg_path"></param>
/// <param name="band_ids"></param>
/// <returns></returns>
int ImageMatch::gdal2JPG(std::string gdal_path, std::string jpg_path, int band_ids)
{
std::cout << "convert gdal to jpg , beigining : \t" << getCurrentTimeString() << endl;
gdalImage gdalimg(gdal_path);
cv::Mat img(gdalimg.height, gdalimg.width, CV_8U);
{
int start_ids = 0;
int line_invert = int(80000000 / gdalimg.width);
Eigen::MatrixXd temp(line_invert, gdalimg.width);
double min_value = 0, temp_min = 0;
double max_value = 0, temp_max = 0;
// 线性拉伸2%
Eigen::MatrixXd hist = gdalimg.getHist(band_ids);
int count = 0;
int sum_count = gdalimg.height * gdalimg.width;
int rows = hist.rows();
for (int i = 0; i < rows; i++) {
if ((count+ hist(i, 1)) / sum_count > 0.01) {
min_value = hist(i, 0);
break;
}
count = count + hist(i, 1);
}
count = 0;
for (int i = rows - 1; i >= 0; i--) {
if ((count + hist(i, 1)) / sum_count > 0.01) {
max_value = hist(i, 0);
break;
}
count = count + hist(i, 1);
}
// 重新缩放最大值，最小值
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;
temp = gdalimg.getData(start_ids, 0, line_invert, gdalimg.width, band_ids);
temp = (temp.array() - min_value) * 250 / (max_value - min_value);
for (int i = 0; i < line_invert; i++) {
for (int j = 0; j < gdalimg.width; j++) {
if (temp(i, j) < 0) { temp(i, j) = 0; }
if (temp(i, j) > 250) { temp(i, j) = 250; }
uchar tempvalue = uchar(temp(i, j));;
img.at<uchar>(i + start_ids, j) = tempvalue;
}
}
start_ids = start_ids + line_invert;
} while (start_ids < gdalimg.height);
cv::Mat result1(gdalimg.height, gdalimg.width, CV_8U);
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); //在这个填入你要的图片质量
cv::Mat out = result1;
//cv::resize(img, out, cv::Size(img.cols, img.rows), cv::INTER_AREA);
cv::imwrite(jpg_path, out, compression_params);
}
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 << "=========================================\n";
return 0;
}
/// <summary>
/// 获取模型匹配点
/// </summary>
/// <param name="ori_power_path"></param>
/// <param name="sim_sum_path"></param>
/// <returns></returns>
Eigen::MatrixXd ImageMatch::ImageMatch_ori_sim(std::string ori_power_path, std::string sim_sum_path, int roughSize, int PreciseSize, int scale, int searchSize, int roughStep,int preciseStep)
{
std::cout << "match point between ori and sim , overing : \t" << getCurrentTimeString() << endl;
std::cout << "match parameters : \n";
std::cout << "rough match templet size : \t" << roughSize << endl;
std::cout << "Precise match templet size : \t" << PreciseSize << endl;
std::cout << "Precise match wait search size : \t" << searchSize << endl;
std::cout << "Precise match scale : \t" << scale << endl;
std::cout << "Precise match step : \t" << preciseStep << endl;
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);
cv::Mat sim = openJPG(sim_sum_path);
int sim_rows = sim.rows;
int sim_cols = sim.cols;
int BlockCount = 0;
for (int i = 0; i < sim_rows; i=i+ roughSize) {
for (int j = 0; j < sim_cols; j=j+ roughSize) {
if (i + 2 * roughSize < sim_rows && j + 2 * roughSize < sim_cols) {
BlockCount++;
}
}
}
Eigen::MatrixXd sim_block(BlockCount,3);
BlockCount = 0;
for (int i = 0; i < sim_rows; i = i + roughSize) {
for (int j = 0; j < sim_cols; j = j + roughSize) {
if (i + 2 * roughSize < sim_rows && j + 2 * roughSize < sim_cols) {
sim_block(BlockCount, 0) = i;
sim_block(BlockCount, 1) = j;
BlockCount++;
}
}
}
Eigen::MatrixXd Tempmatchpoints(BlockCount, 6);
Tempmatchpoints = Tempmatchpoints.array() * 0;
int count = 0;
#pragma omp parallel for num_threads(8)
for (int ii = 0; ii < BlockCount; ii++) {
int i = sim_block(ii, 0);
int j = sim_block(ii, 1);
cv::Mat templet_mat = sim(Rect(j, i, 2 * roughSize, 2 * roughSize));
cv::Scalar mean1;
cv::Mat stddevMat;
cv::meanStdDev(templet_mat, mean1, stddevMat);
double minvalue = 0;
double maxvalue = 0;
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;
int len_j = search_j + 4 * searchSize;
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);//用于检测矩阵中最大值和最小值的位置
double offset_j = maxLoc.x + search_j;
double offset_i = maxLoc.y + search_i;
//cv::Scalar mean1;
//cv::Mat stddevMat;
//cv::meanStdDev(templet_mat, mean1, stddevMat);
Tempmatchpoints(ii, 0) = double(j); //sim_x
Tempmatchpoints(ii, 1) = double(i); //sim_y
Tempmatchpoints(ii, 2) = offset_j; // sim_x
Tempmatchpoints(ii, 3) = offset_i; // sim_y
Tempmatchpoints(ii, 4) = tempmaxVal; // maxVal
Tempmatchpoints(ii, 5) = sig; // maxVal
std::cout << 100.0*count/ BlockCount<<"\t:\t"<< count << " / " << BlockCount<<"\t" << Tempmatchpoints(ii, 0) << "\t" << Tempmatchpoints(ii, 1) << "\t" << Tempmatchpoints(ii, 2)
<< "\t" << Tempmatchpoints(ii, 3) << "\t" << Tempmatchpoints(ii, 4) << "\t" << Tempmatchpoints(ii, 5) << endl;
count++;
}
//BlockCount = 0;
count = 0;
for (int ii = 0; ii < BlockCount; ii++) {
if (Tempmatchpoints(ii, 4) > 0.7) {
std::cout << Tempmatchpoints(ii, 0) << "\t" << Tempmatchpoints(ii, 1) << "\t" << Tempmatchpoints(ii, 2)
<< "\t" << Tempmatchpoints(ii, 3) << "\t" << Tempmatchpoints(ii, 4) << "\t" << Tempmatchpoints(ii, 5) << endl;
//BlockCount++;
count++;
}
}
Eigen::MatrixXd matchpoints(count, 5);
count = 0;
for (int ii = 0; ii < BlockCount; ii++) {
if (Tempmatchpoints(ii, 4) > 0.7) {
matchpoints(count, 0)=Tempmatchpoints(ii, 0); //sim_x
matchpoints(count, 1) = Tempmatchpoints(ii, 1); //sim_y
matchpoints(count, 2) = Tempmatchpoints(ii, 2);// sim_x
matchpoints(count, 3) = Tempmatchpoints(ii, 3);// sim_y
matchpoints(count, 4) = Tempmatchpoints(ii, 4); // maxVal
count++;
}
}
/*
// step 1: 粗匹配，分块均匀匹配
std::cout << "rough match , begining : \t" << getCurrentTimeString() << endl;
double offset_x = 0, offset_y = 0;
cv::Mat mask = sim;
cv::Mat temp_ori = ori;
int offsetcount = 0;
Eigen::MatrixXd matchpoints;
{
int row_count = ori.rows;
int col_count = ori.cols;
int count = 0;
for (int i = roughSize; i < row_count; i = i + roughSize + 500) { // y
for (int j = roughSize; j < col_count; j = j + roughSize + 200) { //x
count = count + 1;
}
}
matchpoints= Eigen::MatrixXd::Zero(count, 5);//ori_x,ori_y,sim_x,sim_y,maxval
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); // 初始化互斥锁
count = 0;
int search_count = 0;
#pragma omp parallel for num_threads(8)
for (int i = roughSize; i < row_count; i = i + roughSize+ 500) { // y
double tempminVal = 100, tempmaxVal = 0;
cv::Mat templet_mat;
cv::Mat result;
for (int j = roughSize; j < col_count; j = j + roughSize+ 200) { //x
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); //获得互斥器
Point minLoc;
Point maxLoc;
Point matchLoc;
cv::minMaxLoc(result, &tempminVal, &tempmaxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
if (tempmaxVal >= maxVal) {
offset_x = maxLoc.x - (j - roughSize);
offset_y = maxLoc.y - (i - roughSize);
//maxVal = tempmaxVal;
matchpoints(count, 0) = (j - roughSize); //ori_x
matchpoints(count, 1) = (i - roughSize); //ori_y
double temp= maxLoc.x;
matchpoints(count, 2) = temp; // sim_x
temp = maxLoc.y;
matchpoints(count, 3) = temp; // sim_y
matchpoints(count, 4) = tempmaxVal; // maxVal
count = count + 1;
offsetcount += 1;
}
search_count = search_count + 1;
std::cout << j - roughSize << "\t" << i - roughSize << "\t" << maxLoc.x << "\t" << maxLoc.y << "\t" << tempmaxVal << "\t" << search_count << "\t" << matchpoints.rows() << endl;
omp_unset_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;
std::cout << "offset X : \t" << offset_x << endl;
std::cout << "offset Y : \t" << offset_y << endl;
std::cout << "maxVal : \t" << maxVal << endl;
}
std::cout << "rough match out : \t" << getCurrentTimeString() << endl;
// step1.1: 粗匹配绘制结果
std::string rough_math_path = sim_sum_path;
boost::algorithm::replace_last(rough_math_path, ".", "_ori_in_sim.");
std::cout << "ori in sim :\t" << rough_math_path << endl;
cv::Mat ori_in_sim = sim;
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); //在这个填入你要的图片质量
cv::Mat out = ori_in_sim;
cv::imwrite(rough_math_path, out, compression_params);
if (offsetcount == 0) { // 表示全局粗匹配失败，无法进行精校准，考虑直接定位法
std::cout << "there are not effective point in rought match \t" << endl;
return Eigen::MatrixXd(0, 5);
}
Eigen::MatrixXd matchpointstemp(offsetcount,5);
for (int i = 0; i < matchpoints.rows(); i++) {
if (matchpoints(i, 4) > 0.7) {
for (int j = 0; j < matchpoints.cols(); j++) {
matchpointstemp(i, j) = matchpoints(i, j);
}
}
}
matchpoints = matchpointstemp;
*/
//std::cout << "rough match , overing : \t" << getCurrentTimeString() << endl;
//// step 2: 精匹配，
//std::cout << "Precise match , begining : \t" << getCurrentTimeString() << endl;
//Eigen::MatrixXd matchpoints;
//{
// int row_count = ori.rows;
// int col_count = ori.cols;
// int count = 0;
// for (int i = PreciseSize; i < row_count; i = i + preciseStep) { // y
// for (int j = PreciseSize; j < col_count; j = j + preciseStep) { //x
// count = count + 1;
// }
// }
// matchpoints(count, 5);//ori_x,ori_y,sim_x,sim_y,maxval
// double templeta_size = PreciseSize * scale;
// double search_size = searchSize * scale;
// count = 0;
// double ori_start_x, ori_start_y;
// double sim_start_x, sim_start_y;
// row_count = row_count - PreciseSize;
// col_count = col_count - PreciseSize;
// omp_lock_t 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;
// cv::Mat search_mat;
// cv::Mat resample_templet_mat;
// cv::Mat resample_search_mat;
// 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));
// resample_templet_mat = resampledMat(templet_mat, templeta_size);
// resample_search_mat = resampledMat(search_mat, search_size);
// matchTemplate(sim, templet_mat, result, cv::TM_CCORR_NORMED);
// // 通过函数 minMaxLoc 定位最匹配的位置;
// cv::Point minLoc; cv::Point maxLoc;
// cv::Point matchLoc;
// cv::minMaxLoc(result, &minVal, &maxVal, &minLoc, &maxLoc);//用于检测矩阵中最大值和最小值的位置
// if (maxVal > 0.7) {
// omp_set_lock(&lock); //获得互斥器
//
// matchpoints(count, 0) = ori_start_x; //ori_x
// matchpoints(count, 1) = ori_start_y; //ori_y
// matchpoints(count, 2) = sim_start_x + maxLoc.x * 1.0 / scale; // sim_x
// matchpoints(count, 3) = sim_start_y + maxLoc.y * 1.0 / scale; // sim_y
// matchpoints(count, 4) = maxVal; // maxVal
// count = count + 1;
//
// omp_unset_lock(&lock); //释放互斥器
// }
// }
// }
// omp_destroy_lock(&lock); //销毁互斥器
//}
//
std::cout << "Precise match , ending : \t" << getCurrentTimeString() << endl;
std::cout << "=======================================================================" << endl;
std::cout << "match point result: \t" << getCurrentTimeString() << endl;
std::cout << "=======================================================================" << endl;
{
std::cout << "ori_x\tori_y\tsim_x\tsim_y\tmaxval \t" << endl;
int count = matchpoints.rows();
for (int i = 0; i < count; i++) {
std::cout << matchpoints(i, 0) << "\t" << matchpoints(i, 1) << "\t" << matchpoints(i, 2) << "\t" << matchpoints(i, 3) << "\t" << matchpoints(i, 4) << endl;
}
}
std::cout << "=======================================================================" << endl;
return matchpoints;
}
Eigen::MatrixXd ImageMatch::CreateMatchModel(Eigen::MatrixXd offsetXY_matrix)
{
// ori_x,oir_y,to_x,to_y,maxval
// 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
temp.col(1) = offsetXY_matrix.col(3).array(); // r
temp.col(2) = offsetXY_matrix.col(2).array(); // c
temp.col(3) = offsetXY_matrix.col(3).array().pow(2);//r2
temp.col(4) = offsetXY_matrix.col(2).array().pow(2);//c2
temp.col(5) = offsetXY_matrix.col(2).array()* offsetXY_matrix.col(3).array();//r*c
Eigen::MatrixXd matchmodel(2, 6);
Eigen::MatrixXd tempx= temp.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(offset_x);//x c
Eigen::MatrixXd tempy = temp.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(offset_y);//y r
matchmodel.row(1) = tempx;
matchmodel.row(0) = tempy;
return matchmodel;
}
Eigen::MatrixXd ImageMatch::correctMatchModel(Eigen::MatrixXd r, Eigen::MatrixXd c)
{
Eigen::MatrixXd a0 = Eigen::MatrixXd::Ones(1,r.cols());
Eigen::MatrixXd a1 = r.array();
Eigen::MatrixXd a2 = c.array();
Eigen::MatrixXd a3 = r.array().pow(2);
Eigen::MatrixXd a4 = c.array().pow(2);
Eigen::MatrixXd a5 = r.array() * c.array();
Eigen::MatrixXd offset(2, r.cols());//r,c
offset.row(0) = r.array() + this->match_model(0, 0) * a0.array() + this->match_model(0, 1) * a1.array() + this->match_model(0, 2) * a2.array() + this->match_model(0, 3) * a3.array() + this->match_model(0, 4) * a4.array() + this->match_model(0, 5) * a5.array();
offset.row(1) = c.array() + this->match_model(1, 0) * a0.array() + this->match_model(1, 1) * a1.array() + this->match_model(1, 2) * a2.array() + this->match_model(1, 3) * a3.array() + this->match_model(1, 4) * a4.array() + this->match_model(1, 5) * a5.array();
return offset;
}
int ImageMatch::outMatchModel(std::string matchmodel_path)
{
ofstream fout(matchmodel_path, ios::trunc);
fout << "model:" << endl;
fout << this->match_model(0, 0) << " " << this->match_model(0, 1) << " " << this->match_model(0, 2) << " " << this->match_model(0, 3) << " " << this->match_model(0, 4) << " " << this->match_model(0, 5) << endl;
fout << this->match_model(1, 0) << " " << this->match_model(1, 1) << " " << this->match_model(1, 2) << " " << this->match_model(1, 3) << " " << this->match_model(1, 4) << " " << this->match_model(1, 5) << endl;
fout << "model_points:" << endl;
std::cout << "ori_x\tori_y\tsim_x\tsim_y\tmaxval \t" << endl;
int count = this->offsetXY_matrix.rows();
for (int i = 0; i < count; i++) {
std::cout << this->offsetXY_matrix(i, 0) << "\t" << this->offsetXY_matrix(i, 1) << "\t" << this->offsetXY_matrix(i, 2) << "\t" << this->offsetXY_matrix(i, 3) << "\t" << this->offsetXY_matrix(i, 4) << endl;
}
fout.close();
return 0;
}
/// <summary>
/// 读取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新出现的空指针常量
{
throw new exception("图片文件不存在");
}
return image;
}
cv::Mat resampledMat(cv::Mat& image, int targetSize, int interpolation)
{
cv::Mat out;
cv::resize(image, out, cv::Size(targetSize, targetSize), interpolation);
return out;
}

39
ImageMatch.h Normal file
View File

@ -0,0 +1,39 @@
#pragma once
//#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
// ±¾µØ·½·¨
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#include <boost/filesystem.hpp>
#include <omp.h>
#include "baseTool.h"
#include "simptsn.h"
#include "SateOrbit.h"
#include <gdal_utils.h>
#include <opencv2/core.hpp>
#include <opencv2/opencv.hpp>
#include <opencv2/core/eigen.hpp>
using namespace std;
using namespace Eigen;
class ImageMatch
{
public:
int gdal2JPG(std::string gdal_path,std::string jpg_path,int band_ids);
Eigen::MatrixXd ImageMatch_ori_sim(std::string ori_power_path, std::string sim_sum_path, int roughSize=500, int Precise=300,int scale=5, int searchSize=1000,int roughStep=400 ,int preciseStep=300);
Eigen::MatrixXd CreateMatchModel(Eigen::MatrixXd offsetXY_matrix);
Eigen::MatrixXd correctMatchModel(Eigen::MatrixXd r, Eigen::MatrixXd c);
int outMatchModel(std::string matchmodel_path);
//²ÎÊý
Eigen::MatrixXd offsetXY_matrix;
Eigen::MatrixXd match_model;
};
cv::Mat openJPG(std::string jpg_path);
cv::Mat resampledMat(cv::Mat& templet, int targetSize, int interpolation = cv::INTER_AREA);

21
LICENSE.txt Normal file
View File

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

1
OctreeNode.cpp Normal file
View File

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

264
OctreeNode.h Normal file
View File

@ -0,0 +1,264 @@
#pragma once
#include <iostream>
using namespace std;
//定义八叉树节点类
template<class T>
struct OctreeNode
{
T data; //节点数据
T xmin, xmax; //节点坐标，即六面体个顶点的坐标
T ymin, ymax;
T zmin, zmax;
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 //节点类
(T nodeValue = T(),
T xminValue = T(), T xmaxValue = T(),
T yminValue = T(), T ymaxValue = T(),
T zminValue = T(), T zmaxValue = T(),
OctreeNode<T>* top_left_front_Node = NULL,
OctreeNode<T>* top_left_back_Node = NULL,
OctreeNode<T>* top_right_front_Node = NULL,
OctreeNode<T>* top_right_back_Node = NULL,
OctreeNode<T>* bottom_left_front_Node = NULL,
OctreeNode<T>* bottom_left_back_Node = NULL,
OctreeNode<T>* bottom_right_front_Node = NULL,
OctreeNode<T>* bottom_right_back_Node = NULL)
:data(nodeValue),
xmin(xminValue), xmax(xmaxValue),
ymin(yminValue), ymax(ymaxValue),
zmin(zminValue), zmax(zmaxValue),
top_left_front(top_left_front_Node),
top_left_back(top_left_back_Node),
top_right_front(top_right_front_Node),
top_right_back(top_right_back_Node),
bottom_left_front(bottom_left_front_Node),
bottom_left_back(bottom_left_back_Node),
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
if (maxdepth >= 0)
{
root = new OctreeNode<T>();
//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 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);
createOctree(root->top_right_back, maxdepth, xmax - xm, xmax, ymin, ymax - ym, zmax - zm, zmax);
createOctree(root->bottom_left_front, maxdepth, xmin, xmax - xm, ymax - ym, ymax, zmin, zmax - zm);
createOctree(root->bottom_left_back, maxdepth, xmin, xmax - xm, ymin, ymax - ym, zmin, zmax - zm);
createOctree(root->bottom_right_front, maxdepth, xmax - xm, xmax, ymax - ym, ymax, zmin, zmax - zm);
createOctree(root->bottom_right_back, maxdepth, xmax - xm, xmax, ymin, ymax - ym, zmin, zmax - zm);
}
}
int i = 1;
//先序遍历八叉树
template <class T>
void preOrder(OctreeNode<T>*& p)
{
if (p)
{
//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;
i += 1;
cout << endl;
preOrder(p->top_left_front);
preOrder(p->top_left_back);
preOrder(p->top_right_front);
preOrder(p->top_right_back);
preOrder(p->bottom_left_front);
preOrder(p->bottom_left_back);
preOrder(p->bottom_right_front);
preOrder(p->bottom_right_back);
cout << endl;
}
}
//求八叉树的深度
template<class T>
int depth(OctreeNode<T>*& p)
{
if (p == NULL)
return -1;
int h = depth(p->top_left_front);
return h + 1;
}
template<class T>
int num(OctreeNode<T>*& p)
{
if (p == NULL)
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;
if (1 == num)
result = 1;
else
{
for (int i = 1; i < num; i++)
result = 2 * result;
}
return result;
}
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;
int tmaxdepth = 0;
double txm = 1, tym = 1, tzm = 1;
double xm = (p->xmax - p->xmin) / 2;
double ym = (p->ymax - p->ymin) / 2;
double zm = (p->ymax - p->ymin) / 2;
times++;
if (x > xmax || x<xmin || y>ymax || y<ymin || z>zmax || z < zmin)
{
//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 << "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;
return 1;
}
else if (x < (p->xmax - xm) && y < (p->ymax - ym) && z < (p->zmax - zm))
{
find(p->bottom_left_back, x, y, z);
}
else if (x < (p->xmax - xm) && y<(p->ymax - ym) && z>(p->zmax - zm))
{
find(p->top_left_back, x, y, z);
}
else if (x > (p->xmax - xm) && y < (p->ymax - ym) && z < (p->zmax - zm))
{
find(p->bottom_right_back, x, y, z);
}
else if (x > (p->xmax - xm) && y<(p->ymax - ym) && z>(p->zmax - zm))
{
find(p->top_right_back, x, y, z);
}
else if (x<(p->xmax - xm) && y>(p->ymax - ym) && z < (p->zmax - zm))
{
find(p->bottom_left_front, x, y, z);
}
else if (x<(p->xmax - xm) && y>(p->ymax - ym) && z > (p->zmax - zm))
{
find(p->top_left_front, x, y, z);
}
else if (x > (p->xmax - xm) && y > (p->ymax - ym) && z < (p->zmax - zm))
{
find(p->bottom_right_front, x, y, z);
}
else if (x > (p->xmax - xm) && y > (p->ymax - ym) && z > (p->zmax - zm))
{
find(p->top_right_front, x, y, z);
}
}
//main函数
/*
int main()
{
OctreeNode<double>* rootNode = NULL;
int choiced = 0;
cout << "系统开始前请先创建八叉树" << endl;
cout << "请输入最大递归深度：" << endl;
cin >> maxdepth;
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)
{
tmaxdepth = cal(maxdepth);
txm = (xmax - xmin) / tmaxdepth;
tym = (ymax - ymin) / tmaxdepth;
tzm = (zmax - zmin) / tmaxdepth;
createOctree(rootNode, maxdepth, xmin, xmax, ymin, ymax, zmin, zmax);
}
while (true)
{
system("cls");
cout << "请选择操作：\n";
cout << "\t1.计算空间中区域的个数\n";
cout << "\t2.先序遍历八叉树\n";
cout << "\t3.查看树深度\n";
cout << "\t4.查找节点 \n";
cout << "\t0.退出\n";
cin >> choiced;
if (choiced == 0)
return 0;
if (choiced == 1)
{
system("cls");
cout << "空间区域个数" << endl;
cout << num(rootNode);
}
if (choiced == 2)
{
system("cls");
cout << "先序遍历八叉树结果：/n";
i = 1;
preOrder(rootNode);
cout << endl;
system("pause");
}
if (choiced == 3)
{
system("cls");
int dep = depth(rootNode);
cout << "此八叉树的深度为" << dep + 1 << endl;
system("pause");
}
if (choiced == 4)
{
system("cls");
cout << "请输入您希望查找的点的坐标，顺序如下：x,y,z\n";
double x, y, z;
cin >> x >> y >> z;
times = 0;
cout << endl << "开始搜寻该点……" << endl;
find(rootNode, x, y, z);
system("pause");
}
else
{
system("cls");
cout << "\n\n错误选择！\n";
system("pause");
}
}
*/

1
README.md Normal file
View File

@ -0,0 +1 @@
# SIMOrthoProgram

25
RPC_Correct.cpp Normal file
View File

@ -0,0 +1,25 @@
#pragma once
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#include <boost/filesystem.hpp>
#include <omp.h>
#include "baseTool.h"
#include "simptsn.h"
#include "SateOrbit.h"
#include "ImageMatch.h"
#include <gdal_utils.h>
#include <armadillo>
#include <proj.h>
#include "gdal_priv.h"
#include "gdal_alg.h"
#include "RPC_Correct.h"
using namespace std;
using namespace Eigen;
using namespace arma;

35
RPC_Correct.h Normal file
View File

@ -0,0 +1,35 @@
#pragma once
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#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"
#include "SateOrbit.h"
#include "ImageMatch.h"
using namespace std;
using namespace Eigen;
using namespace arma;
using namespace std;
using namespace Eigen;
using namespace arma;
// רÃÅÓÃÓÚ½âÎöRPC
class RPC_Correct
{
};

BIN
SIMOrthoProgram.aps Normal file
View File

Binary file not shown.

368
SIMOrthoProgram.cpp Normal file
View File

@ -0,0 +1,368 @@
// SIMOrthoProgram.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//
//#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
//#include <mkl.h>
#include <stdlib.h>
#include <direct.h>
// gdal
#include <proj.h>
#include <string>'
#include "gdal_priv.h"
#include "ogr_geometry.h"
#include "gdalwarper.h"
#include "baseTool.h"
#include "simptsn.h"
#include "test_moudel.h"
#include <Windows.h>
using namespace std;
using namespace Eigen;
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
//输入参数
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::cout << "==========================================================================" << endl;
std::cout << "in parameters:========================================================" << endl;
std::cout << "parameters file path:\t" << parameter_path << endl;
std::cout << "input dem image(WGS84)" << dem_path << endl;
std::cout << "the sar image:\n" << in_sar_path << endl;
std::cout << "the work path for outputing temp file :\t" << work_path << endl;
std::cout << "the out file for finnal file:\t" << taget_path << endl;
simProcess process;
std::cout << "==========================================================================" << endl;
process.InitSimulationSAR(parameter_path, work_path, taget_path, dem_path, in_sar_path);
std::cout << "==========================================================================" << endl;
}
//mode 2
void calIncident_localIncident_angle(int argc, char* argv[]) {
std::cout << "mode 2: get incident angle and local incident angle by rc_wgs84 and dem and statellite model:\n";
std::cout << "SIMOrthoProgram.exe 2 in_parameter_path in_dem_path in_rc_wgs84_path out_incident_angle_path out_local_incident_angle_path";
std::string parameter_path = argv[2];
std::string dem_path = argv[3];
std::string in_rc_wgs84_path = argv[4];
std::string out_incident_angle_path = argv[5];
std::string out_local_incident_angle_path = argv[6];
simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
std::cout << "==========================================================================" << endl;
process.pstn = pstn;
process.calcalIncident_localIncident_angle(dem_path, in_rc_wgs84_path, out_incident_angle_path, out_local_incident_angle_path, pstn);
}
// mode 3
void calInterpolation_cubic_Wgs84_rc_sar(int argc, char* argv[]) {
std::cout << "mode 3: interpolation(cubic convolution) orth sar value by rc_wgs84 and ori_sar image and model:\n ";
std::cout << "SIMOrthoProgram.exe 3 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path";
std::string parameter_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\orth_para.txt"; argv[2];
std::string in_rc_wgs84_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\dem_rc.tiff"; argv[3];
std::string in_ori_sar_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\unpack\\GF3_MYN_QPSI_011437_E99.2_N28.6_20181012_L1A_AHV_L10003514912\\GF3_MYN_QPSI_011437_E99.2_N28.6_20181012_L1A_VV_L10003514912.tiff"; argv[4];
std::string out_orth_sar_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\GF3_MYN_QPSI_011437_E99.2_N28.6_20181012_L1A_VV_L10003514912_GTC_rpc_geo.tif"; argv[5];
parameter_path = argv[2];
in_rc_wgs84_path = argv[3];
in_ori_sar_path = argv[4];
out_orth_sar_path = argv[5];
simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
process.pstn = pstn;
std::cout << "==========================================================================" << endl;
process.interpolation_GTC_sar(in_rc_wgs84_path, in_ori_sar_path, out_orth_sar_path, pstn);
}
// mode 4 处理 RPC的入射角
void getRPC_Incident_localIncident_angle(int argc, char* argv[]) {
std::cout << "mode 4: get RPC incident and local incident angle sar model:";
std::cout << "SIMOrthoProgram.exe 4 in_parameter_path in_dem_path in_rpc_rc_path out_rpc_dem_path out_incident_angle_path out_local_incident_angle_path";
std::string parameter_path = argv[2];
std::string in_dem_path = argv[3];
std::string in_rpc_rc_path = argv[4];
std::string out_rpc_dem_path = argv[5];
std::string out_incident_angle_path = argv[6];
std::string out_local_incident_angle_path = argv[7];
simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
std::cout << "==========================================================================" << endl;
process.CreateRPC_DEM(in_rpc_rc_path, in_dem_path, out_rpc_dem_path);
process.calcalIncident_localIncident_angle(out_rpc_dem_path, in_rpc_rc_path, out_incident_angle_path, out_local_incident_angle_path, pstn);
}
//mode 5
void cal_ori_2_power_tiff(int argc, char* argv[]) {
std::cout << "mode 5: convert ori tiff to power tiff:";
std::cout << "SIMOrthoProgram.exe 5 in_ori_path out_power_path";
std::string in_ori_path = argv[2];
std::string out_power_path = argv[3];
simProcess process;
process.ori_sar_power(in_ori_path, out_power_path);
}
// mode 6
void cal_GEC_Incident_localIncident_angle(int argc, char* argv[]) {
std::cout << "mode 6: get gec incident and local incident angle sar model:";
std::cout << "SIMOrthoProgram.exe 6 in_parameter_path in_dem_path in_gec_lon_lat_path out_incident_angle_path out_local_incident_angle_path";
std::string parameter_path = argv[2];
std::string dem_path = argv[3];
std::string in_gec_lon_lat_path = argv[4];
std::string out_incident_angle_path = argv[5];
std::string out_local_incident_angle_path = argv[6];
simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
std::cout << "==========================================================================" << endl;
process.calGEC_Incident_localIncident_angle(dem_path, in_gec_lon_lat_path, out_incident_angle_path, out_local_incident_angle_path, pstn);
}
// mode 7
void RPC_inangle(int argc, char* argv[]) {
std::cout << "mode 7: get rpc incident and local incident angle sar model:";
std::cout << "SIMOrthoProgram.exe 7 in_parameter_path in_dem_path in_gec_lon_lat_path work_path taget_path out_incident_angle_path out_local_incident_angle_path";
std::string parameter_path = argv[2];
std::string dem_path = argv[3];
std::string in_gec_lon_lat_path = argv[4];
std::string work_path = argv[5];
std::string taget_path = argv[6];
std::string out_incident_angle_path = argv[7];
std::string out_local_incident_angle_path = argv[8];
std::string out_incident_angle_geo_path = argv[9];
std::string out_local_incident_angle_geo_path = argv[10];
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.InitRPCIncAngle(parameter_path, work_path, taget_path, dem_path, in_gec_lon_lat_path, out_incident_angle_path, out_local_incident_angle_path, out_incident_angle_geo_path, out_local_incident_angle_geo_path);
std::cout << "==========================================================================" << endl;
}
// mode 3
void calInterpolation_cubic_Wgs84_rc_sar_sigma(int argc, char* argv[]) {
std::cout << "mode 3: interpolation(cubic convolution) orth sar value by rc_wgs84 and ori_sar image and model:\n ";
std::cout << "SIMOrthoProgram.exe 3 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path";
std::string parameter_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\orth_para.txt"; argv[2];
std::string in_rc_wgs84_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\dem_rc.tiff"; argv[3];
std::string in_ori_sar_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\unpack\\GF3_MYN_QPSI_011437_E99.2_N28.6_20181012_L1A_AHV_L10003514912\\GF3_MYN_QPSI_011437_E99.2_N28.6_20181012_L1A_VV_L10003514912.tiff"; argv[4];
std::string out_orth_sar_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\GF3_MYN_QPSI_011437_E99.2_N28.6_20181012_L1A_VV_L10003514912_GTC_rpc_geo.tif"; argv[5];
parameter_path = argv[2];
in_rc_wgs84_path = argv[3];
in_ori_sar_path = argv[4];
out_orth_sar_path = argv[5];
simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
process.pstn = pstn;
std::cout << "==========================================================================" << endl;
process.interpolation_GTC_sar_sigma(in_rc_wgs84_path, in_ori_sar_path, out_orth_sar_path, pstn);
}
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];
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);
std::cout << "==========================================================================" << endl;
}
void Scatter2Grid_lon_lat(int argc, char* argv[]) {
std::cout << "mode 10";
std::cout << "SIMOrthoProgram.exe 10 lon_lat_path data_tiff grid_path space";
std::string lon_lat_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\RD_ori_sim.tif";
std::string data_tiff = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1B_HH_L10003923848_db.tif";
std::string grid_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary\\package\\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1B_HH_L10003923848__db.tif";
double space = 5;
lon_lat_path = argv[2];
data_tiff = argv[3];
grid_path = argv[4];
space = stod(argv[5]);
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.Scatter2Grid(lon_lat_path, data_tiff, grid_path, space);
std::cout << "==========================================================================" << endl;
}
string GetExePath()
{
char szFilePath[MAX_PATH + 1] = { 0 };
GetModuleFileNameA(NULL, szFilePath, MAX_PATH);
/*
strrchr:函数功能：查找一个字符c在另一个字符串str中末次出现的位置（也就是从str的右侧开始查找字符c首次出现的位置），
并返回这个位置的地址。如果未能找到指定字符，那么函数将返回NULL。
使用这个地址返回从最后一个字符c到str末尾的字符串。
*/
(strrchr(szFilePath, '\\'))[0] = 0; // 删除文件名，只获得路径字串//
string path = szFilePath;
return path;
}
/// <summary>
/// 初始化
/// </summary>
void initProjEnv() {
PJ_CONTEXT* C;
C = proj_context_create();
std::cout << "========================== init PROJ ================================================" << endl;
string exepath = GetExePath();
char buffer[10240];
int i;
for (i = 0; i < exepath.length(); i++) {
buffer[i] = exepath[i];
}
buffer[i] = '\0';
const char* proj_share_path = buffer;
proj_context_set_search_paths(C, 1, &proj_share_path);
char* buf = nullptr;
size_t sz = 0;
if (_dupenv_s(&buf, &sz, "PROJ_LIB") == 0 && buf != nullptr)
{
printf("PROJ_LIB = %s\n", buf);
std::string newEnv = "PROJ_LIB=" + std::string(buffer);
_putenv(newEnv.c_str());
}
else {
std::string newEnv = "PROJ_LIB=" + std::string(buffer);
_putenv(newEnv.c_str());
}
if (_dupenv_s(&buf, &sz, "PROJ_LIB") == 0 && buf != nullptr)
{
std::cout << "after PROJ_LIB = " << buf << endl;
}
free(buf);
std::cout << "========================================================================================" << endl;
}
int main(int argc, char* argv[])
{
initProjEnv();
//WGS84_J2000();
cout << "test\t" << acos(-1) << endl;
cout << getAngle(Landpoint{ -3421843,5089485,3630606 }, Landpoint{ -2609414,4763328,3332879 }) << endl;;
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;
GDALAllRegister();
if (argc == 1) { // 测试参数
// 算法说明
std::cout << "========================== description ================================================" << endl;
std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl;
std::cout << "mode 1: Preprocess\n ";
std::cout << "SIMOrthoProgram.exe 1 in_parameter_path in_dem_path in_ori_sar_path in_work_path in_taget_path out_GEC_dem_path out_GTC_rc_path out_GEC_lon_lat_path out_clip_dem_path" << endl;
std::cout << "mode 2: get incident angle and local incident angle by rc_wgs84 and dem and statellite model:\n";
std::cout << "SIMOrthoProgram.exe 2 in_parameter_path in_dem_path in_rc_wgs84_path out_incident_angle_path out_local_incident_angle_path";
std::cout << "mode 3: interpolation(cubic convolution) orth sar value by rc_wgs84 and ori_sar image and model:\n ";
std::cout << "SIMOrthoProgram.exe 3 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path";
std::cout << "mode 4: get RPC incident and local incident angle sar model:";
std::cout << "SIMOrthoProgram.exe 4 in_parameter_path in_dem_path in_rpc_rc_path out_rpc_dem_path out_incident_angle_path out_local_incident_angle_path";
std::cout << "mode 5: interpolation(cubic convolution) orth sar value by gec_lon_lat and dem and ori_sar image and sar model:";
std::cout << "SIMOrthoProgram.exe 5 in_parameter_path in_gec_lon_lat_path in_dem_path in_sar_path out_orth_sar_path";
std::cout << "mode 6: get gec incident and local incident angle sar model:";
std::cout << "SIMOrthoProgram.exe 6 in_parameter_path in_dem_path in_gec_lon_lat_path out_incident_angle_path out_local_incident_angle_path";
if (mode == 10) {
Scatter2Grid_lon_lat(argc, argv);
}
else {
test_main(mode, "D:\\MicroSAR\\C-SAR\\Ortho\\Ortho\\Temporary");
}
//calInterpolation_cubic_Wgs84_rc_sar(argc, argv);
}
else if (argc > 1) { // 预处理模块
std::cout << "=============================description V2.0 =============================================" << endl;
std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl;
std::cout << "algorithm moudel:.exe [modeparamert] {otherParaments}" << endl;
mode = stoi(argv[1]);
if (mode == 0) {
test_main(mode, argv[2]);
}
else if (mode == 1) {
PreProcess(argc, argv);
}
else if (mode == 2) { // RPC 计算模块
calIncident_localIncident_angle(argc, argv);
}
else if (mode == 3) {
calInterpolation_cubic_Wgs84_rc_sar(argc, argv);
}
else if (mode == 4) {
getRPC_Incident_localIncident_angle(argc, argv);
}
else if (mode == 5) {
cal_ori_2_power_tiff(argc, argv);
}
else if (mode == 6) {
cal_GEC_Incident_localIncident_angle(argc, argv);
}
else if (mode == 7) {
RPC_inangle(argc, argv);
}
else if (mode == 8) {
createRPC_lon_lat(argc, argv);
}
else if (mode == 9) {
calInterpolation_cubic_Wgs84_rc_sar_sigma(argc, argv);
}
else if (mode == 10) {
Scatter2Grid_lon_lat(argc, argv);
}
}
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
std::cout << "program over：\t" << getCurrentTimeString() << endl;
}
// 运行程序: Ctrl + F5 或调试 >“开始执行(不调试)”菜单
// 调试程序: F5 或调试 >“开始调试”菜单
// 入门使用技巧:
// 1. 使用解决方案资源管理器窗口添加/管理文件
// 2. 使用团队资源管理器窗口连接到源代码管理
// 3. 使用输出窗口查看生成输出和其他消息
// 4. 使用错误列表窗口查看错误
// 5. 转到“项目”>“添加新项”以创建新的代码文件，或转到“项目”>“添加现有项”以将现有代码文件添加到项目
// 6. 将来，若要再次打开此项目，请转到“文件”>“打开”>“项目”并选择 .sln 文件

69
SIMOrthoProgram.rc Normal file
View File

@ -0,0 +1,69 @@
// Microsoft Visual C++ generated resource script.
//
#include "resource.h"
#define APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE 2 resource.
//
#include "winres.h"
/////////////////////////////////////////////////////////////////////////////
#undef APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
// 中文(简体，中国) resources
#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_CHS)
LANGUAGE LANG_CHINESE, SUBLANG_CHINESE_SIMPLIFIED
#pragma code_page(936)
#ifdef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// TEXTINCLUDE
//
1 TEXTINCLUDE
BEGIN
"resource.h\0"
END
2 TEXTINCLUDE
BEGIN
"#include ""winres.h""\r\n"
"\0"
END
3 TEXTINCLUDE
BEGIN
"\r\n"
"\0"
END
#endif // APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// proj
//
IDR_PROJ1 proj "proj.db"
#endif // 中文(简体，中国) resources
/////////////////////////////////////////////////////////////////////////////
#ifndef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE 3 resource.
//
/////////////////////////////////////////////////////////////////////////////
#endif // not APSTUDIO_INVOKED

31
SIMOrthoProgram.sln Normal file
View File

@ -0,0 +1,31 @@

Microsoft Visual Studio Solution File, Format Version 12.00
# Visual Studio Version 17
VisualStudioVersion = 17.2.32602.215
MinimumVisualStudioVersion = 10.0.40219.1
Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "SIMOrthoProgram", "SIMOrthoProgram.vcxproj", "{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}"
EndProject
Global
GlobalSection(SolutionConfigurationPlatforms) = preSolution
Debug|x64 = Debug|x64
Debug|x86 = Debug|x86
Release|x64 = Release|x64
Release|x86 = Release|x86
EndGlobalSection
GlobalSection(ProjectConfigurationPlatforms) = postSolution
{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}.Debug|x64.ActiveCfg = Debug|x64
{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}.Debug|x64.Build.0 = Debug|x64
{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}.Debug|x86.ActiveCfg = Debug|Win32
{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}.Debug|x86.Build.0 = Debug|Win32
{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}.Release|x64.ActiveCfg = Release|x64
{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}.Release|x64.Build.0 = Release|x64
{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}.Release|x86.ActiveCfg = Release|Win32
{7722B0A9-572B-4E32-AFBE-4DC88898F8EE}.Release|x86.Build.0 = Release|Win32
EndGlobalSection
GlobalSection(SolutionProperties) = preSolution
HideSolutionNode = FALSE
EndGlobalSection
GlobalSection(ExtensibilityGlobals) = postSolution
SolutionGuid = {CE6D0980-E92A-4188-91CE-EEE5749F908C}
EndGlobalSection
EndGlobal

196
SIMOrthoProgram.vcxproj Normal file
View File

@ -0,0 +1,196 @@
<?xml version="1.0" encoding="utf-8"?>
<Project DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
<ItemGroup Label="ProjectConfigurations">
<ProjectConfiguration Include="Debug|Win32">
<Configuration>Debug</Configuration>
<Platform>Win32</Platform>
</ProjectConfiguration>
<ProjectConfiguration Include="Release|Win32">
<Configuration>Release</Configuration>
<Platform>Win32</Platform>
</ProjectConfiguration>
<ProjectConfiguration Include="Debug|x64">
<Configuration>Debug</Configuration>
<Platform>x64</Platform>
</ProjectConfiguration>
<ProjectConfiguration Include="Release|x64">
<Configuration>Release</Configuration>
<Platform>x64</Platform>
</ProjectConfiguration>
</ItemGroup>
<PropertyGroup Label="Globals">
<VCProjectVersion>16.0</VCProjectVersion>
<Keyword>Win32Proj</Keyword>
<ProjectGuid>{7722b0a9-572b-4e32-afbe-4dc88898f8ee}</ProjectGuid>
<RootNamespace>SIMOrthoProgram</RootNamespace>
</PropertyGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'" Label="Configuration">
<ConfigurationType>Application</ConfigurationType>
<UseDebugLibraries>true</UseDebugLibraries>
<PlatformToolset>v143</PlatformToolset>
<CharacterSet>Unicode</CharacterSet>
</PropertyGroup>
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'" Label="Configuration">
<ConfigurationType>Application</ConfigurationType>
<UseDebugLibraries>false</UseDebugLibraries>
<PlatformToolset>v143</PlatformToolset>
<WholeProgramOptimization>true</WholeProgramOptimization>
<CharacterSet>Unicode</CharacterSet>
</PropertyGroup>
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
<ConfigurationType>Application</ConfigurationType>
<UseDebugLibraries>false</UseDebugLibraries>
<PlatformToolset>v143</PlatformToolset>
<CharacterSet>Unicode</CharacterSet>
<UseInteloneMKL>Sequential</UseInteloneMKL>
</PropertyGroup>
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
<ConfigurationType>Application</ConfigurationType>
<UseDebugLibraries>false</UseDebugLibraries>
<PlatformToolset>v143</PlatformToolset>
<WholeProgramOptimization>true</WholeProgramOptimization>
<CharacterSet>Unicode</CharacterSet>
<UseInteloneMKL>Parallel</UseInteloneMKL>
<UseIntelMPI>false</UseIntelMPI>
</PropertyGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
<ImportGroup Label="ExtensionSettings">
</ImportGroup>
<ImportGroup Label="Shared">
</ImportGroup>
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
</ImportGroup>
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
</ImportGroup>
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
</ImportGroup>
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
</ImportGroup>
<PropertyGroup Label="UserMacros" />
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
<LibraryPath>C:\ProgramData\Miniconda3\envs\orth\Library\lib;C:\Program Files (x86)\Windows Kits\10\Lib\10.0.19041.0\ucrt\x64;$(oneMKLOmpLibDir);$(LibraryPath)</LibraryPath>
<CopyLocalProjectReference>true</CopyLocalProjectReference>
<CopyCppRuntimeToOutputDir>false</CopyCppRuntimeToOutputDir>
<ExternalIncludePath>C:\ProgramData\Miniconda3\envs\orth\Library\include;$(ExternalIncludePath)</ExternalIncludePath>
</PropertyGroup>
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
<CopyLocalProjectReference>true</CopyLocalProjectReference>
<CopyLocalDebugSymbols>true</CopyLocalDebugSymbols>
<CopyCppRuntimeToOutputDir>true</CopyCppRuntimeToOutputDir>
<IncludePath>$(IncludePath)</IncludePath>
<ExecutablePath>$(ExecutablePath)</ExecutablePath>
<ExternalIncludePath>$(VC_IncludePath);$(WindowsSDK_IncludePath);</ExternalIncludePath>
<LibraryPath>$(LibraryPath)</LibraryPath>
<ReferencePath>$(VC_ReferencesPath_x64);</ReferencePath>
<LibraryWPath>$(WindowsSDK_MetadataPath);</LibraryWPath>
</PropertyGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
<ClCompile>
<WarningLevel>Level3</WarningLevel>
<SDLCheck>true</SDLCheck>
<PreprocessorDefinitions>WIN32;_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
<ConformanceMode>true</ConformanceMode>
</ClCompile>
<Link>
<SubSystem>Console</SubSystem>
<GenerateDebugInformation>true</GenerateDebugInformation>
</Link>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">
<ClCompile>
<WarningLevel>Level3</WarningLevel>
<FunctionLevelLinking>true</FunctionLevelLinking>
<IntrinsicFunctions>true</IntrinsicFunctions>
<SDLCheck>true</SDLCheck>
<PreprocessorDefinitions>WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
<ConformanceMode>true</ConformanceMode>
</ClCompile>
<Link>
<SubSystem>Console</SubSystem>
<EnableCOMDATFolding>true</EnableCOMDATFolding>
<OptimizeReferences>true</OptimizeReferences>
<GenerateDebugInformation>true</GenerateDebugInformation>
</Link>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
<ClCompile>
<WarningLevel>Level3</WarningLevel>
<SDLCheck>true</SDLCheck>
<PreprocessorDefinitions>_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
<ConformanceMode>true</ConformanceMode>
<OpenMPSupport>true</OpenMPSupport>
<LanguageStandard>Default</LanguageStandard>
<LanguageStandard_C>stdc11</LanguageStandard_C>
<RuntimeLibrary>MultiThreaded</RuntimeLibrary>
</ClCompile>
<Link>
<SubSystem>Console</SubSystem>
<GenerateDebugInformation>true</GenerateDebugInformation>
<AdditionalDependencies>gdal_i.lib;%(AdditionalDependencies)</AdditionalDependencies>
</Link>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
<ClCompile>
<WarningLevel>Level3</WarningLevel>
<FunctionLevelLinking>
</FunctionLevelLinking>
<IntrinsicFunctions>false</IntrinsicFunctions>
<SDLCheck>true</SDLCheck>
<PreprocessorDefinitions>_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
<ConformanceMode>true</ConformanceMode>
<RuntimeLibrary>MultiThreaded</RuntimeLibrary>
<OpenMPSupport>true</OpenMPSupport>
<LanguageStandard>Default</LanguageStandard>
<LanguageStandard_C>Default</LanguageStandard_C>
<Optimization>MaxSpeed</Optimization>
<FavorSizeOrSpeed>Neither</FavorSizeOrSpeed>
<EnableFiberSafeOptimizations>false</EnableFiberSafeOptimizations>
<WholeProgramOptimization>false</WholeProgramOptimization>
<MultiProcessorCompilation>true</MultiProcessorCompilation>
<AdditionalOptions>/bigobj %(AdditionalOptions)</AdditionalOptions>
</ClCompile>
<Link>
<SubSystem>Console</SubSystem>
<EnableCOMDATFolding>true</EnableCOMDATFolding>
<OptimizeReferences>true</OptimizeReferences>
<GenerateDebugInformation>true</GenerateDebugInformation>
<StackCommitSize>
</StackCommitSize>
<AdditionalDependencies>%(AdditionalDependencies)</AdditionalDependencies>
</Link>
</ItemDefinitionGroup>
<ItemGroup>
<ClCompile Include="baseTool.cpp" />
<ClCompile Include="ImageMatch.cpp" />
<ClCompile Include="OctreeNode.cpp" />
<ClCompile Include="RPC_Correct.cpp" />
<ClCompile Include="simptsn.cpp" />
<ClCompile Include="SateOrbit.cpp" />
<ClCompile Include="SIMOrthoProgram.cpp" />
<ClCompile Include="test_moudel.cpp" />
</ItemGroup>
<ItemGroup>
<ClInclude Include="baseTool.h" />
<ClInclude Include="ImageMatch.h" />
<ClInclude Include="OctreeNode.h" />
<ClInclude Include="resource.h" />
<ClInclude Include="RPC_Correct.h" />
<ClInclude Include="simptsn.h" />
<ClInclude Include="SateOrbit.h" />
<ClInclude Include="test_moudel.h" />
</ItemGroup>
<ItemGroup>
<ResourceCompile Include="SIMOrthoProgram.rc" />
</ItemGroup>
<ItemGroup>
<None Include="proj.db" />
</ItemGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
<ImportGroup Label="ExtensionTargets">
</ImportGroup>
</Project>

77
SIMOrthoProgram.vcxproj.filters Normal file
View File

@ -0,0 +1,77 @@
<?xml version="1.0" encoding="utf-8"?>
<Project ToolsVersion="4.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
<ItemGroup>
<Filter Include="源文件">
<UniqueIdentifier>{4FC737F1-C7A5-4376-A066-2A32D752A2FF}</UniqueIdentifier>
<Extensions>cpp;c;cc;cxx;c++;cppm;ixx;def;odl;idl;hpj;bat;asm;asmx</Extensions>
</Filter>
<Filter Include="头文件">
<UniqueIdentifier>{93995380-89BD-4b04-88EB-625FBE52EBFB}</UniqueIdentifier>
<Extensions>h;hh;hpp;hxx;h++;hm;inl;inc;ipp;xsd</Extensions>
</Filter>
<Filter Include="资源文件">
<UniqueIdentifier>{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}</UniqueIdentifier>
<Extensions>rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms</Extensions>
</Filter>
</ItemGroup>
<ItemGroup>
<ClCompile Include="SIMOrthoProgram.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="baseTool.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="simptsn.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="SateOrbit.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="ImageMatch.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="test_moudel.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="RPC_Correct.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="OctreeNode.cpp">
<Filter>源文件</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<ClInclude Include="baseTool.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="simptsn.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="SateOrbit.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="ImageMatch.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="resource.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="test_moudel.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="RPC_Correct.h">
<Filter>头文件</Filter>
</ClInclude>
<ClInclude Include="OctreeNode.h">
<Filter>头文件</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ResourceCompile Include="SIMOrthoProgram.rc">
<Filter>资源文件</Filter>
</ResourceCompile>
</ItemGroup>
<ItemGroup>
<None Include="proj.db" />
</ItemGroup>
</Project>

4
SIMOrthoProgram.vcxproj.user Normal file
View File

@ -0,0 +1,4 @@
<?xml version="1.0" encoding="utf-8"?>
<Project ToolsVersion="Current" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
<PropertyGroup />
</Project>

89
SateOrbit.cpp Normal file
View File

@ -0,0 +1,89 @@
#include "SateOrbit.h"
//#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
//#include <mkl.h>
#include <armadillo>
using namespace std;
using namespace Eigen;
using namespace arma;
OrbitPoly::OrbitPoly()
{
}
OrbitPoly::OrbitPoly(int polynum, Eigen::MatrixX<double> polySatellitePara, double SatelliteModelStartTime)
{
if (polySatellitePara.rows() != polynum||polySatellitePara.cols()!=6) {
throw exception("轨道模型参数形状不一致");
}
this->polySatellitePara = polySatellitePara;
this->SatelliteModelStartTime = SatelliteModelStartTime;
this->polynum = polynum;
}
OrbitPoly::~OrbitPoly()
{
}
Eigen::MatrixX<double> OrbitPoly::SatelliteSpacePoint(Eigen::MatrixX<double> satellitetime)
{
Eigen::MatrixX<double> result= SatelliteSpacePoints(satellitetime, this->SatelliteModelStartTime, this->polySatellitePara, this->polynum);
return result;
}
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>
Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double>& satellitetime, double SatelliteModelStartTime, Eigen::MatrixX<double>& polySatellitePara, int polynum)
{
if (satellitetime.cols() != 1) {
throw exception("the size of satellitetime has error!!");
}
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);
for (int i = 0; i < polynum; i++) {
satelliteTime.col(i) = satellitetime.array().pow(i);
}
// 计算
Eigen::MatrixX<double> satellitePoints(satellitetime_num, 6);
satellitePoints = satelliteTime * polySatellitePara;
return satellitePoints;
}

53
SateOrbit.h Normal file
View File

@ -0,0 +1,53 @@
#pragma once
///
/// 计算卫星轨道
///
//#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
// 本地方法
#include "baseTool.h"
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
//#include <mkl.h>
//#include <armadillo>
using namespace std;
using namespace Eigen;
//using namespace arma;
/// <summary>
/// 多项式轨道模型
/// </summary>
class OrbitPoly {
public:
//OrbitPoly(std::string orbitModelPath);
OrbitPoly();
OrbitPoly(int polynum, Eigen::MatrixX<double> polySatellitePara, double SatelliteModelStartTime);
~OrbitPoly();
Eigen::MatrixX<double> SatelliteSpacePoint(Eigen::MatrixX<double> satellitetime);
Eigen::MatrixX<double> SatelliteSpacePoint(long double satellitetime);
public:
int polynum;
Eigen::MatrixX<double> polySatellitePara;
double SatelliteModelStartTime;
};
/// <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>
Eigen::MatrixX<double> SatelliteSpacePoints(Eigen::MatrixX<double> &satellitetime, double SatelliteModelStartTime, Eigen::MatrixX<double>& polySatellitePara, int polynum = 4);

19
WGS84_J2000.cpp Normal file
View File

@ -0,0 +1,19 @@
#pragma once
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#include <math.h>
#include <string>
#include <omp.h>
#include <complex>
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <ogrsf_frmts.h>
#include <opencv2/opencv.hpp>
#include <opencv2/core/eigen.hpp>
#include <opencv2/features2d.hpp>
#include <fstream>
#include "baseTool.h"
#include"WGS84_J2000.h"

10
WGS84_J2000.h Normal file
View File

@ -0,0 +1,10 @@
#pragma once
#include <stdlib.h>
#include <vector>
#include <iostream>
#include <omp.h>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#include "baseTool.h"
#include <string>

1348
baseTool.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

388
baseTool.h Normal file
View File

@ -0,0 +1,388 @@
#pragma once
///
/// 基本类、基本函数
///
//#define EIGEN_USE_MKL_ALL
//#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
//#include <mkl.h>
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#include <string>
#include <omp.h>
#include <complex>
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <ogrsf_frmts.h>
#include <opencv2/opencv.hpp>
#include <opencv2/core/eigen.hpp>
#include <opencv2/features2d.hpp>
#include <fstream>
using namespace std;
using namespace Eigen;
#define PI_180 180/3.141592653589793238462643383279;
#define T180_PI 3.141592653589793238462643383279/180;
#define Radians2Degrees(Radians) Radians*PI_180
#define Degrees2Radians(Degrees) Degrees*T180_PI
const double PI = 3.141592653589793238462643383279;
const double epsilon = 0.000000000000001;
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 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影像的像素坐标；
{
/// <summary>
/// 经度x
/// </summary>
double lon; // 经度x lon pixel_col
/// <summary>
/// 纬度y
/// </summary>
double lat; // 纬度y lat pixel_row
/// <summary>
/// 高度z
/// </summary>
double ati; // 高程z ati pixel_time
};
struct Point_3d {
double x;
double y;
double z;
};
Landpoint operator +(const Landpoint& p1, const Landpoint& p2);
Landpoint operator -(const Landpoint& p1, const Landpoint& p2);
bool operator ==(const Landpoint& p1, const Landpoint& p2);
Landpoint operator *(const Landpoint& p, double scale);
/// <summary>
/// 向量A,B的夹角,角度
/// </summary>
/// <param name="a"></param>
/// <param name="b"></param>
/// <returns>角度制 0-360度，逆时针</returns>
double getAngle(const Landpoint& a, const Landpoint& b);
/// <summary>
/// 点乘
/// </summary>
/// <param name="p1"></param>
/// <param name="p2"></param>
/// <returns></returns>
double dot(const Landpoint& p1, const Landpoint& p2);
double getlength(const Landpoint& p1);
Landpoint crossProduct(const Landpoint& a, const Landpoint& b);
struct DemBox {
double min_lat; //纬度
double min_lon;//经度
double max_lat;//纬度
double max_lon;//经度
};
/// <summary>
/// gdalImage图像操作类
/// </summary>
class gdalImage
{
public: // 方法
gdalImage(string raster_path);
~gdalImage();
void setHeight(int);
void setWidth(int);
void setTranslationMatrix(Eigen::MatrixXd gt);
void setData(Eigen::MatrixXd);
Eigen::MatrixXd getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
void saveImage(Eigen::MatrixXd,int start_row,int start_col, int band_ids);
void saveImage();
void setNoDataValue(double nodatavalue,int band_ids);
int InitInv_gt();
Landpoint getRow_Col(double lon, double lat);
Landpoint getLandPoint(double i, double j, double ati);
double mean(int bandids=1);
double max(int bandids=1);
double min(int bandids=1);
GDALRPCInfo getRPC();
Eigen::MatrixXd getLandPoint(Eigen::MatrixXd points);
Eigen::MatrixXd getHist(int bandids);
public:
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 data;
string projection;
};
gdalImage CreategdalImage(string img_path, int height, int width, int band_num, Eigen::MatrixXd gt, std::string projection, bool need_gt=true);
void clipGdalImage(string in_path, string out_path, DemBox box, double pixelinterval);
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample);
int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample=GRIORA_Bilinear);
/////////////////////////////// 基本图像类 结束 //////////////////////////////////////////////////////////
string Convert(float Num);
std::string JoinPath(const std::string& path, const std::string& filename);
////////////////////////////// 坐标部分基本方法 //////////////////////////////////////////
/// <summary>
/// 将经纬度转换为地固参心坐标系
/// </summary>
/// <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>
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>
//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);
complex<double> Cubic_kernel_weight(double s);
double Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12, Landpoint p22);
inline float cross2d(Point_3d a, Point_3d b) { return a.x * b.y - a.y * b.x; }
inline Point_3d operator-(Point_3d a, Point_3d b) {
return Point_3d{ a.x - b.x, a.y - b.y, a.z - b.z };
};
inline Point_3d operator+(Point_3d a, Point_3d b) {
return Point_3d{ a.x + b.x, a.y +b.y, a.z + b.z };
};
inline double operator/(Point_3d a, Point_3d b) {
return sqrt(pow(a.x,2)+ pow(a.y, 2))/sqrt(pow(b.x, 2)+ pow(b.y, 2));
};
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之间
&& 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;
else
return false;
}
inline Point_3d invBilinear(Point_3d p, Point_3d a, Point_3d b, Point_3d c, Point_3d d)
{
Point_3d res ;
Point_3d e = b - a;
Point_3d f = d - a;
Point_3d g = a - b + c - d;
Point_3d h = p - a;
double k2 = cross2d(g, f);
double k1 = cross2d(e, f) + cross2d(h, g);
double k0 = cross2d(h, e);
double u, v;
// if edges are parallel, this is a linear equation
if (abs(k2) < 0.001)
{
v = -k0 / k1;
u = (h.x - f.x *v) / (e.x + g.x *v);
p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v;
return p;
}
// otherwise, it's a quadratic
else
{
float w = k1 * k1 - 4.0 * k0 * k2;
if (w < 0.0){
// 可能在边界上
if (onSegment(a, b, p)) {
Point_3d tt = b - a;
Point_3d ttpa = p - a;
double scater=ttpa / tt;
if(scater<0||scater>1){ return { -9999,-9999,-9999 }; }
p.z = a.z + scater * tt.z;
return p;
}
else if (onSegment(b, c, p)) {
Point_3d tt = c-b;
Point_3d ttpa = p - b;
double scater = ttpa / tt;
if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; }
p.z = b.z + scater * tt.z;
return p;
}
else if (onSegment(c, d, p)) {
Point_3d tt = d-c;
Point_3d ttpa = p - c;
double scater = ttpa / tt;
if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; }
p.z = c.z + scater * tt.z;
return p;
}
else if (onSegment(d, a, p)) {
Point_3d tt = a-d;
Point_3d ttpa = p - d;
double scater = ttpa / tt;
if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; }
p.z = d.z + scater * tt.z;
return p;
}
return { -9999,-9999,-9999 };
}
else {
w = sqrt(w);
float ik2 = 0.5 / k2;
float v = (-k1 - w) * ik2;
float u = (h.x - f.x * v) / (e.x + g.x * v);
if (u < 0.0 || u>1.0 || v < 0.0 || v>1.0)
{
v = (-k1 + w) * ik2;
u = (h.x - f.x * v) / (e.x + g.x * v);
}
p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v;
return p;
}
}
p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v;
return p;
}
//
// 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）
**/
inline double sind(double degree) {
return sin(degree * d2r);
}
inline double cosd(double d) {
return cos(d * d2r);
}
/*
class WGS84_J2000
{
public:
WGS84_J2000();
~WGS84_J2000();
public:
// step1 WGS 84 转换到协议地球坐标系。
static Eigen::MatrixXd WGS84TECEF(Eigen::MatrixXd WGS84_Lon_lat_ait);
//step 2 协议地球坐标系 转换为瞬时地球坐标系
static Eigen::MatrixXd ordinateSingleRotate(int axis, double angle_deg);
// 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 瞬时真天球坐标系 转到瞬时平天球 坐标系
static int nutationInLongitudeCaculate(double JD, double& epthilongA_deg, double& dertaPthi_deg, double& dertaEpthilong_deg, double& epthilong_deg);
// 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博客_年积日计算公式
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列
static Eigen::Matrix<double, 77, 11> IAU2000ModelParams;
};
*/

16
resource.h Normal file
View File

@ -0,0 +1,16 @@
//{{NO_DEPENDENCIES}}
// Microsoft Visual C++ 生成的包含文件。
// 供 SIMOrthoProgram.rc 使用
//
#define IDR_PROJ1 101
// Next default values for new objects
//
#ifdef APSTUDIO_INVOKED
#ifndef APSTUDIO_READONLY_SYMBOLS
#define _APS_NEXT_RESOURCE_VALUE 102
#define _APS_NEXT_COMMAND_VALUE 40001
#define _APS_NEXT_CONTROL_VALUE 1001
#define _APS_NEXT_SYMED_VALUE 101
#endif
#endif

3313
simptsn.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

232
simptsn.h Normal file
View File

@ -0,0 +1,232 @@
#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>
//#include <mkl.h>
#include <omp.h>
#include "SateOrbit.h"
#include "ImageMatch.h"
using namespace std;
using namespace Eigen;
using namespace arma;
///////////// 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);//根据理论模型计算多普勒频移值
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::MatrixXd WGS842J2000(Eigen::MatrixXd blh);
//Landpoint WGS842J2000(Landpoint p);
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 simsar_path;// 仿真影像坐标
int height; // 影像的高
int width;
// 入射角
double near_in_angle;// 近入射角
double far_in_angle;// 远入射角
// SAR的成像参数
double widthspace;// 距离向分辨率
double R0;//近斜距
double LightSpeed; // 光速
double lamda;// 波长
double refrange;// 参考斜距
double imgStartTime;// 成像起始时间
double PRF;// 脉冲重复率
int doppler_paramenter_number;//多普勒系数个数
MatrixX<double> doppler_paras;// 多普勒系数
OrbitPoly orbit; // 轨道模型
};
/// <summary>
/// 获取局地入射角，角度值
/// </summary>
/// <param name="satepoint"></param>
/// <param name="landpoint"></param>
/// <param name="slopeVector"></param>
/// <returns></returns>
double getlocalIncAngle(Landpoint& satepoint, Landpoint& landpoint, Landpoint& slopeVector);
/// <summary>
/// 侧视入射角,角度值
/// </summary>
/// <param name="satepoint"></param>
/// <param name="landpoint"></param>
/// <returns></returns>
double getIncAngle(Landpoint& satepoint, Landpoint& landpoint);
/// <summary>
/// ASF计算方法
/// </summary>
class ASFOrthClass {
public:
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);
Eigen::Vector3d GetXYZByBetaAlpha(double alpha,double beta, Eigen::Vector3d SatellitePosition,double R, Eigen::MatrixXd M);
double FD(double alpha, double beta, Eigen::Vector3d SatelliteVelocity, Eigen::Vector3d TargetVelocity,double R,double lamda, Eigen::MatrixXd M);
double FR(double alpha, double beta, Eigen::Vector3d SatellitePosition, Eigen::MatrixXd M, double R_threshold, double H = 0);
Eigen::Vector3d ASF(double R, Eigen::Vector3d SatellitePosition, Eigen::Vector3d SatelliteVelocity, Eigen::Vector3d TargetPosition, PSTNAlgorithm PSTN, double R_threshold, double H = 0, double alpha0 = 0, double beta0 = 0);
};
/// <summary>
/// 仿真处理流程
/// </summary>
class simProcess {
public:
simProcess();
~simProcess();
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 CreateSARDEM();
int dem2SAR(); // 切换行号
int SARIncidentAngle();
int SARSimulationWGS();
int SARSimulation();
int in_sar_power();
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 SARIncidentAngle_RPC();
int createRPCDEM();
// 格网离散插值
int Scatter2Grid(std::string lon_lat_path, std::string data_tiff, std::string grid_path, double space);
//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);
void interpolation_GTC_sar(std::string in_rc_wgs84_path, std::string in_ori_sar_path, std::string out_orth_sar_path, PSTNAlgorithm pstn);
// 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);
// 内部参数
public:
ImageMatch matchmodel;
bool isMatchModel;
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 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 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_slantRange_path; // 斜距文件
std::string out_plant_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;// 映射角文件
};
bool PtInRect(Point_3d pCur, Point_3d pLeftTop, Point_3d pRightTop, Point_3d pRightBelow, Point_3d pLeftBelow);

298
test_moudel.cpp Normal file
View File

@ -0,0 +1,298 @@
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
//#include <mkl.h>
#include <stdlib.h>
#include <direct.h>
// gdal
#include <proj.h>
#include <string>
#include "gdal_priv.h"
#include "ogr_geometry.h"
#include "gdalwarper.h"
#include "baseTool.h"
#include "simptsn.h"
using namespace std;
using namespace Eigen;
#include "test_moudel.h"
int test_main(int mode, string rootpath)
{
switch (mode)
{
case 0:
{
test_mode_0();
break;
};
case 1:
{
test_mode_1(rootpath);
break;
}
case 7: {
test_mode_7();
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";
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);
std::cout << "==========================================================================" << endl;
break;
}
default:
test_ASF(rootpath);
break;
}
return 0;
}
int test_mode_7() {
//
std::cout << "mode 7: get rpc incident and local incident angle sar model:";
std::cout << "SIMOrthoProgram.exe 7 in_parameter_path in_dem_path in_gec_lon_lat_path work_path taget_path out_incident_angle_path out_local_incident_angle_path";
std::string parameter_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary\\package\\orth_para.txt";
std::string dem_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary\\TestDEM\\mergedDEM.tif";
std::string in_gec_lon_lat_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary\\package\\RPC_ori_sim.tif";
std::string work_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary";
std::string taget_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary\\package";
std::string out_incident_angle_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary\\package\\inc_angle.tiff";
std::string out_local_incident_angle_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary\\package\\local_inc_angle.tiff";
std::string out_incident_angle_geo_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary\\package\\inc_angle_geo.tiff";
std::string out_local_incident_angle_geo_path = "D:\\MicroSAR\\C-SAR\\Ortho\\Temporary\\package\\local_inc_angle_geo.tiff";
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.InitRPCIncAngle(parameter_path, work_path, taget_path, dem_path, in_gec_lon_lat_path, out_incident_angle_path, out_local_incident_angle_path, out_incident_angle_geo_path, out_local_incident_angle_geo_path);
std::cout << "==========================================================================" << endl;
return 0;
}
int test_mode_0()
{
std::cout << "State:\tTEST MODE" << endl;
test_LLA2XYZ_XYZ2LLA();
return 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
//输入参数
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::cout << "==========================================================================" << endl;
std::cout << "in parameters:========================================================" << endl;
std::cout << "parameters file path:\t" << parameter_path << endl;
std::cout << "input dem image(WGS84)" << dem_path << endl;
std::cout << "the sar image:\n" << in_sar_path << endl;
std::cout << "the work path for outputing temp file :\t" << work_path << endl;
std::cout << "the out file for finnal file:\t" << taget_path << endl;
simProcess process;
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 work_path = root_path + "\\TestSim";
std::string taget_path = root_path + "\\output";
std::string out_GEC_dem_path = root_path + "\\output\\GEC_dem.tif";
std::string out_GTC_rc_path = root_path + "\\output\\GTC_rc_wgs84.tif";
std::string out_GEC_lon_lat_path = root_path + "\\output\\GEC_lon_lat.tif";
std::string out_clip_DEM_path = root_path + "\\output\\Orth_dem.tif";
std::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::cout << "out_GEC_dem_path out_GTC_rc_path out_GEC_lon_lat_path out_clip_dem_path" << endl;
std::cout << "==========================================================================" << endl;
std::cout << "in parameters:========================================================" << endl;
std::cout << "parameters file path:\t" << parameter_path << endl;
std::cout << "input dem image(WGS84)" << dem_path << endl;
std::cout << "the reference sar image for match simulation sar image:\n" << ori_sar_path << endl;
std::cout << "the work path for outputing temp file :\t" << work_path << endl;
std::cout << "the out file for finnal file:\t" << taget_path << endl;
simProcess process;
ASFOrthClass ASFClass;
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
std::string dem_slope_path = JoinPath(work_path, "dem_slope.tif");
std::string dem_aspect_path = JoinPath(work_path, "dem_aspect.tif");
std::string sim_sar_path = JoinPath(work_path, "sim_sar_sum.tif"); // sim count dem
std::string ori_sar_rsc_path = JoinPath(work_path, "ori_sar_power.tif"); // ori_power
std::string ori_sar_rsc_jpg_path = JoinPath(work_path, "ori_sar_power.jpg");
std::string sim_sar_jpg_path = JoinPath(work_path, "sim_sar_sum.jpg");
std::string matchmodel_path = JoinPath(work_path, "matchmodel.txt"); // matchmodel
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);
}
if (boost::filesystem::exists(boost::filesystem::path(out_sar_rc_path))) {
boost::filesystem::remove(boost::filesystem::path(out_sar_rc_path));
}
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);
return 0;
}
int test_mode_2()
{
std::string parameter_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\orth_para.txt";
std::string dem_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\RPC_Ortho\\RPC_DEM.tiff";
std::string in_rc_wgs84_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\RPC_Ortho\\RPC_sar_rc.tiff";
std::string out_incident_angle_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\RPC_Ortho\\RPC_incident_angle.tiff";
std::string out_local_incident_angle_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\RPC_Ortho\\RPC_local_incident_angle.tiff";
std::cout << "mode 2: get incident angle and local incident angle by rc_wgs84 and dem and statellite model:\n";
std::cout << "SIMOrthoProgram.exe 2 in_parameter_path in_dem_path in_rc_wgs84_path out_incident_angle_path out_local_incident_angle_path";
simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
std::cout << "==========================================================================" << endl;
process.calcalIncident_localIncident_angle(dem_path, in_rc_wgs84_path, out_incident_angle_path, out_local_incident_angle_path, pstn);
return 0;
}
int test_mode_3()
{
std::string parameter_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\orth_para.txt";
std::string in_rc_wgs84_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\output\\GTC_rc_wgs84.tiff";
std::string in_ori_sar_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\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 out_orth_sar_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\output\\GTC\\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_HH_L10003923848.tiff";
std::cout << "mode 3: interpolation(cubic convolution) orth sar value by rc_wgs84 and ori_sar image and model:\n ";
std::cout << "SIMOrthoProgram.exe 3 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path";
simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
std::cout << "==========================================================================" << endl;
process.interpolation_GTC_sar(in_rc_wgs84_path, in_ori_sar_path, out_orth_sar_path, pstn);
return 0;
}
int test_mode_4()
{
std::string parameter_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\orth_para.txt";
std::string in_dem_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\output\\Orth_dem.tiff";
std::string in_rpc_rc_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\output\\\RPC_Ortho\\RPC_sar_rc.tiff";
std::string out_rpc_dem_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\output\\RPC_Ortho\\RPC_DEM.tiff";
std::string out_incident_angle_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\output\\RPC_Ortho\\RPC_incident_angle.tiff";
std::string out_local_incident_angle_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\output\\RPC_Ortho\\RPC_local_incident_angle.tiff";
std::cout << "mode 4: get RPC incident and local incident angle sar model:";
std::cout << "SIMOrthoProgram.exe 4 in_parameter_path in_dem_path in_rpc_rc_path out_rpc_dem_path out_incident_angle_path out_local_incident_angle_path";
simProcess process;
std::cout << "==========================================================================" << endl;
PSTNAlgorithm pstn(parameter_path);
std::cout << "==========================================================================" << endl;
process.CreateRPC_DEM(in_rpc_rc_path, in_dem_path, out_rpc_dem_path);
process.calcalIncident_localIncident_angle(out_rpc_dem_path, in_rpc_rc_path, out_incident_angle_path, out_local_incident_angle_path, pstn);
return 0;
}
int test_mode_5()
{
std::string in_ori_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\output\\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_HH_L10003923848_GTC_geo_wp.tif";
std::string out_power_path = "D:\\MicroSAR\\C-SAR\\SIMOrthoProgram\\Temporary\\package\\GF3_SAY_QPSI_013952_E118.9_N31.5_20190404_L1A_HH_L10003923848_GTC_geo_OrthoResult.tif";
std::cout << "mode 5: convert ori tiff to power tiff:";
std::cout << "SIMOrthoProgram.exe 5 in_ori_path out_power_path";
simProcess process;
process.ori_sar_power(in_ori_path, out_power_path);
return 0;
}
int test_mode_6()
{
return 0;
}
int test_LLA2XYZ_XYZ2LLA()
{
std::cout << "===========================================" << endl;
std::cout << "TEST LLA2XYZ and XYZ2LLA " << endl;
bool pass = false;
Landpoint lla_p{ 30,40,500 };
Landpoint xyz_p = LLA2XYZ(lla_p);
Landpoint lla_p2 = XYZ2LLA(xyz_p);
std::cout << "LLA:\t" << lla_p.lon << "," << lla_p.lat << "," << lla_p.ati << endl;
std::cout << "XYZ:\t" << xyz_p.lon << "," << xyz_p.lat << "," << xyz_p.ati << endl;
std::cout<<"LLA2:\t"<< lla_p2.lon << "," << lla_p2.lat << "," << lla_p2.ati << endl;
pass = (lla_p2.lon-lla_p.lon==0) && (lla_p2.lat- lla_p.lat==0) && (lla_p2.ati-lla_p.ati==0);
std::cout << "TEST LLA2XYZ and XYZ2LLA passed:" << pass << endl;
std::cout << "===========================================" << endl;
return pass?1:0;
}
int test_vector_operator()
{
// 创建
return 0;
}

41
test_moudel.h Normal file
View File

@ -0,0 +1,41 @@
#pragma once
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
//#include <mkl.h>
#include <stdlib.h>
#include <direct.h>
// gdal
#include <proj.h>
#include <string>
#include "gdal_priv.h"
#include "ogr_geometry.h"
#include "gdalwarper.h"
#include "baseTool.h"
#include "simptsn.h"
using namespace std;
using namespace Eigen;
class test_moudel
{
};
int test_main(int mode, string rootpath);
int test_mode_0();
int test_mode_1(string rootpath);
int test_ASF(string rootpath);
int test_mode_2();
int test_mode_3();
int test_mode_4();
int test_mode_5();
int test_mode_6();
int test_mode_7();
int test_LLA2XYZ_XYZ2LLA();
int test_vector_operator();