KJ135
/
RasterProcessTool
mirror of http://172.16.0.12:5000/LAMPSARToolSoftware/RasterProcessTool.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

补充了完整的流程

Release
陈增辉 2024-11-25 14:26:46 +08:00
parent 81a8326364
commit 3d251929e4
20 changed files with 654 additions and 319 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

303
GF3ProcessToolbox/GF3CalibrationAndGeocodingClass.cpp
View File

@ -2,6 +2,11 @@
#include "SatelliteGF3xmlParser.h"
#include <QRegularExpression>
#include <QMessageBox>
#include "SatelliteOribtModel.h"
#include "boost/asio.hpp"
#include <boost/thread.hpp>
#include <thread>
#include <qprogressdialog.h>
ErrorCode GF3CalibrationRaster(QString inRasterPath, QString outRasterPath, double Qualifyvalue, double calibrationConst)
{
@ -13,7 +18,7 @@ ErrorCode GF3CalibrationRaster(QString inRasterPath, QString outRasterPath, doub
gdalImageComplex outraster(outRasterPath);
long blocklines = Memory1GB * 2 /8/ imgraster.width;
long blocklines = Memory1GB * 2 / 8 / imgraster.width;
blocklines = blocklines < 100 ? 100 : blocklines;
Eigen::MatrixXd imgArrb1 = Eigen::MatrixXd::Zero(blocklines, imgraster.width);
Eigen::MatrixXd imgArrb2 = Eigen::MatrixXd::Zero(blocklines, imgraster.width);
@ -42,7 +47,7 @@ ErrorCode ImportGF3L1ARasterToWorkspace(QString inMetaxmlPath, QString inRasterP
gf3xml.loadFile(inMetaxmlPath);
QString filename = getFileNameWidthoutExtend(inRasterPath);
SARSimulationImageL1Dataset l1dataset;
l1dataset.OpenOrNew(outWorkDirPath, filename,gf3xml.height,gf3xml.width);
l1dataset.OpenOrNew(outWorkDirPath, filename, gf3xml.height, gf3xml.width);
l1dataset.setSateAntPos(SatellitePos2SatelliteAntPos(gf3xml.antposs));
l1dataset.setCenterFreq(gf3xml.RadarCenterFrequency);
l1dataset.setcolCount(gf3xml.width);
@ -57,7 +62,7 @@ ErrorCode ImportGF3L1ARasterToWorkspace(QString inMetaxmlPath, QString inRasterP
l1dataset.setIncidenceAngleFarRange(gf3xml.incidenceAngleFarRange);
l1dataset.setIncidenceAngleNearRangeet(gf3xml.incidenceAngleNearRange);
l1dataset.setDopplerParams(gf3xml.d0, gf3xml.d1, gf3xml.d2, gf3xml.d3, gf3xml.d4);
l1dataset.setDopplerCenterCoff(gf3xml.r0, gf3xml.r1,gf3xml.r2, gf3xml.r3, gf3xml.r4);
l1dataset.setDopplerCenterCoff(gf3xml.r0, gf3xml.r1, gf3xml.r2, gf3xml.r3, gf3xml.r4);
l1dataset.setLatitudeCenter(gf3xml.latitude_center);
l1dataset.setLongitudeCenter(gf3xml.longitude_center);
l1dataset.setLatitudeTopLeft(gf3xml.latitude_topLeft);
@ -75,12 +80,12 @@ ErrorCode ImportGF3L1ARasterToWorkspace(QString inMetaxmlPath, QString inRasterP
QString outRasterpath = l1dataset.getImageRasterPath();
ErrorCode errorcode = ErrorCode::SUCCESS;
qDebug() << "导入数据:\t"<<inRasterPath;
qDebug() << "导入数据:\t" << inRasterPath;
switch (polsartype)
{
case POLARHH:
qDebug() << "极化：HH";
errorcode=GF3CalibrationRaster(inRasterPath, outRasterpath, gf3xml.HH_QualifyValue, gf3xml.HH_CalibrationConst);
errorcode = GF3CalibrationRaster(inRasterPath, outRasterpath, gf3xml.HH_QualifyValue, gf3xml.HH_CalibrationConst);
break;
case POLARHV:
qDebug() << "极化：HH";
@ -177,16 +182,16 @@ ErrorCode ImportGF3L1AProcess(QString inMetaxmlPath, QString outWorkDirPath)
switch (poltype)
{
case POLARHH:
errorcode=ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARHH);
errorcode = ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARHH);
break;
case POLARHV:
errorcode=ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARHV);
errorcode = ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARHV);
break;
case POLARVH:
errorcode=ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARVH);
errorcode = ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARVH);
break;
case POLARVV:
errorcode=ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARVV);
errorcode = ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARVV);
break;
case POLARUNKOWN:
break;
@ -207,7 +212,7 @@ ErrorCode ImportGF3L1AProcess(QString inMetaxmlPath, QString outWorkDirPath)
ErrorCode Complex2AmpRaster(QString inComplexPath, QString outRasterPath)
{
gdalImageComplex inimg(inComplexPath);
gdalImage ampimg= CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);
gdalImage ampimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);
long blocklines = Memory1GB * 2 / 8 / inimg.width;
blocklines = blocklines < 100 ? 100 : blocklines;
@ -265,14 +270,14 @@ ErrorCode Complex2dBRaster(QString inComplexPath, QString outRasterPath)
imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 2);
imgArrb1 = imgArr.array().abs().log10()*20.0;
imgArrb1 = imgArr.array().abs().log10() * 20.0;
ampimg.saveImage(imgArrb1, startrow, 0, 1);
}
qDebug() << "影像写入到：" << outRasterPath;
return ErrorCode::SUCCESS;
}
ErrorCode ResampleDEM(QString indemPath, QString outdemPath,double gridx,double gridy)
ErrorCode ResampleDEM(QString indemPath, QString outdemPath, double gridx, double gridy)
{
double gridlat = gridy / 110000.0;
double gridlon = gridx / 100000.0;
@ -283,12 +288,55 @@ ErrorCode ResampleDEM(QString indemPath, QString outdemPath,double gridx,double
return ErrorCode::SUCCESS;
}
else {
QMessageBox::warning(nullptr,u8"警告",u8"请输入WGS84坐标的DEM影像");
QMessageBox::warning(nullptr, u8"警告", u8"请输入WGS84坐标的DEM影像");
return ErrorCode::FAIL;
}
}
ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath,QString outworkdir, QString outlooktablePath)
ErrorCode RD_PSTN(double& refrange, double& lamda, double& timeR, double& R, double& tx, double& ty, double& tz, double& slopex, double& slopey, double& slopez, PolyfitSatelliteOribtModel& polyfitmodel, SatelliteOribtNode& node, double& d0, double& d1, double& d2, double& d3, double& d4)
{
double dt = 1e-6;
double inct = 0;
bool antfalg = false;
double timeR2 = 0;
//return ((R_s1.array() * V_s1.array()).rowwise().sum().array() * (-2) / (R * this->lamda))[0];
//double t = (R - this->refrange) * (1e6 / LIGHTSPEED);
//return this->doppler_paras(0) + this->doppler_paras(1) * t + this->doppler_paras(2) * t * t + this->doppler_paras(3) * t * t * t + this->doppler_paras(4) * t * t * t * t;
double R1 = 0;
double R2 = 0;
double dplerTheory1 = 0;
double dplerTheory2 = 0;
double dplerNumber1 = 0;
double dplerNumber2 = 0;
double dplerR = 0;
// Rst=Rs-Rt;
for (int i = 0; i < 50; i++) {
polyfitmodel.getSatelliteOribtNode(timeR, node, antfalg);
R1 = std::sqrt(std::pow(node.Px - tx, 2) + std::pow(node.Py - ty, 2) + std::pow(node.Pz - tz, 2));
dplerTheory1 = (-2 / lamda) * (((node.Px - tx) * (node.Vx + EARTHWE * ty) + (node.Py - ty) * (node.Vy - EARTHWE * tz) + (node.Pz - tz) * (node.Vz - 0)) / R1);
dplerR = (R1 - refrange) * 1e6 / LIGHTSPEED;
dplerNumber1 = d0 + dplerR * d1 + std::pow(dplerR, 2) * d2 + std::pow(dplerR, 3) * d3 + std::pow(dplerR, 4) * d4;
timeR2 = timeR + dt;
polyfitmodel.getSatelliteOribtNode(timeR2, node, antfalg);
R2 = std::sqrt(std::pow(node.Px - tx, 2) + std::pow(node.Py - ty, 2) + std::pow(node.Pz - tz, 2));
dplerTheory2 = (-2 / lamda) * (((node.Px - tx) * (node.Vx + EARTHWE * ty) + (node.Py - ty) * (node.Vy - EARTHWE * tz) + (node.Pz - tz) * (node.Vz - 0)) / R2);
dplerR = (R2 - refrange) * 1e6 / LIGHTSPEED;
dplerNumber2 = d0 + dplerR * d1 + std::pow(dplerR, 2) * d2 + std::pow(dplerR, 3) * d3 + std::pow(dplerR, 4) * d4;
inct = dt * (dplerTheory2 - dplerNumber1) / (dplerTheory1 - dplerTheory2);
if (std::abs(dplerTheory1 - dplerTheory2) < 1e-9 || std::abs(inct) < dt) {
break;
}
timeR = timeR - inct;
}
R = R1; // 斜距
return ErrorCode::SUCCESS;
}
ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath, QString outworkdir, QString outlooktablePath,QString outLocalIncidenceAnglePath,bool localincAngleFlag)
{
QString indemfilename = getFileNameWidthoutExtend(indemPath);
@ -298,36 +346,239 @@ ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath,QString outw
gdalImage demxyz(resampleDEMXYZPath);// X,Y,Z
gdalImage demslope(resampleDEMSLOPEPath);// X,Y,Z
gdalImage rasterRC = CreategdalImage(outlooktablePath, demxyz.height, demxyz.width,2, demxyz.gt, demxyz.projection, true, true);// X,Y,Z
gdalImage rasterRC = CreategdalImage(outlooktablePath, demxyz.height, demxyz.width, 2, demxyz.gt, demxyz.projection, true, true);// X,Y,Z
gdalImage localincangleRaster;
if (localincAngleFlag) {
localincangleRaster = CreategdalImage(outLocalIncidenceAnglePath, demxyz.height, demxyz.width, 2, demxyz.gt, demxyz.projection, true, true);// X,Y,Z
}
SARSimulationImageL1Dataset l1dataset;
l1dataset.Open(inxmlPath);
// RD 相关参数
QVector<double> dopplers = l1dataset.getDopplerParams();
QVector<double> dopplers = l1dataset.getDopplerParams(); // 多普勒参数
double d0 = dopplers[0];
double d1 = dopplers[1];
double d2 = dopplers[2];
double d3 = dopplers[3];
double d4 = dopplers[4];
double fs = l1dataset.getFs();
double prf = (l1dataset.getEndImageTime()-l1dataset.getStartImageTime())/(l1dataset.getrowCount()-1);
double fs = l1dataset.getFs();// Fs采样
double prf = (l1dataset.getEndImageTime() - l1dataset.getStartImageTime()) / (l1dataset.getrowCount() - 1);// PRF 采样
double nearRange = l1dataset.getNearRange();
double imagestarttime = l1dataset.getStartImageTime();
double refrange = l1dataset.getRefRange();
double lamda = (LIGHTSPEED*1e-6)/ l1dataset.getCenterFreq();
// 构建轨道模型
PolyfitSatelliteOribtModel polyfitmodel;
QVector < SatelliteAntPos > antposes = l1dataset.getXmlSateAntPos();
polyfitmodel.setSatelliteOribtStartTime(imagestarttime);
for (long i = 0; i < antposes.size(); i++) {
SatelliteOribtNode node;
node.time = antposes[i].time;
node.Px = antposes[i].Px;
node.Py = antposes[i].Py;
node.Pz = antposes[i].Pz;
node.Vx = antposes[i].Vx;
node.Vy = antposes[i].Vy;
node.Vz = antposes[i].Vz;
polyfitmodel.addOribtNode(node);
}
polyfitmodel.polyFit(3, false);
qDebug() << "-----------------------------------";
qDebug() << "satellite polyfit model params:\n";
qDebug() << polyfitmodel.getSatelliteOribtModelParamsString();
qDebug() << "-----------------------------------";
// 开始计算查找表
//1.计算分块
long cpucore_num = std::thread::hardware_concurrency();
long blockline = Memory1MB * 500 / 8 / cpucore_num / 4 / l1dataset.getcolCount();
blockline = blockline < 50 ? 50 : blockline;
//2.迭代计算
long colcount = l1dataset.getcolCount();
long rowcount = l1dataset.getrowCount();
long startRId = 0;
long processNumber = 0;
QProgressDialog progressDialog(u8"查找表生成中", u8"终止", 0, rowcount);
progressDialog.setWindowTitle(u8"查找表生成中");
progressDialog.setWindowModality(Qt::WindowModal);
progressDialog.setAutoClose(true);
progressDialog.setValue(0);
progressDialog.setMaximum(rowcount);
progressDialog.setMinimum(0);
progressDialog.show();
omp_lock_t lock;
omp_init_lock(&lock);
#pragma omp parallel for num_threads(cpucore_num-1)
for (startRId = 0; startRId < rowcount; startRId = startRId + blockline) {
Eigen::MatrixXd sar_r = rasterRC.getData(startRId, 0, blockline, colcount, 1);
Eigen::MatrixXd sar_c = rasterRC.getData(startRId, 0, blockline, colcount, 2);
Eigen::MatrixXd dem_x = demxyz.getData(startRId, 0, blockline, colcount, 1);
Eigen::MatrixXd dem_y = demxyz.getData(startRId, 0, blockline, colcount, 2);
Eigen::MatrixXd dem_z = demxyz.getData(startRId, 0, blockline, colcount, 3);
// 逐像素迭代计算
double timeR = 0;
long blockrows = sar_r.rows();
long blockcols = sar_r.cols();
double tx = 0;
double ty = 0;
double tz = 0;
double R = 0;
double slopex=0, slopey=0, slopez=0;
SatelliteOribtNode node{0,0,0,0,0,0,0,0};
bool antflag = false;
for (long i = 0; i < blockrows; i++) {
for (long j = 0; j < blockcols; j++) {
tx = dem_x(i, j);
ty = dem_y(i, j);
tz = dem_z(i, j);
if (RD_PSTN(refrange,lamda, timeR,R,tx,ty,tz,slopex,slopey,slopez,polyfitmodel,node,d0,d1,d2,d3,d4) == ErrorCode::SUCCESS) {
sar_r(i, j) = timeR * prf;
sar_c(i, j) = ((R - nearRange) / 2 / LIGHTSPEED) * fs;
}
else {
sar_r(i, j) = -1;
sar_c(i, j) = -1;
}
}
}
Eigen::MatrixXd Angle_Arr = Eigen::MatrixXd::Zero(dem_x.rows(),dem_x.cols()).array()+181;
if (localincAngleFlag) {
Eigen::MatrixXd demslope_x = demslope.getData(startRId, 0, blockline, colcount, 1);
Eigen::MatrixXd demslope_y = demslope.getData(startRId, 0, blockline, colcount, 2);
Eigen::MatrixXd demslope_z = demslope.getData(startRId, 0, blockline, colcount, 3);
Eigen::MatrixXd Angle_Arr = localincangleRaster.getData(startRId, 0, blockline, colcount, 1);
double Rst_x = 0, Rst_y = 0, Rst_z = 0, localangle = 0;
double slopeR = 0;
for (long i = 0; i < blockrows; i++) {
for (long j = 0; j < blockcols; j++) {
timeR = sar_r(i, j) / prf;
slopex = demslope_x(i, j);
slopey = demslope_y(i, j);
slopez = demslope_z(i, j);
tx = dem_x(i, j);
ty = dem_y(i, j);
tz = dem_z(i, j);
polyfitmodel.getSatelliteOribtNode(timeR, node, antflag);
Rst_x = node.Px - tx;
Rst_y = node.Py - ty;
Rst_z = node.Pz - tz;
R = std::sqrt(Rst_x*Rst_x+Rst_y*Rst_y+Rst_z*Rst_z);
slopeR = std::sqrt(slopex*slopex + slopey * slopey + slopez * slopez);
localangle = ((slopex * Rst_x) + (slopey * Rst_y) + (slopez * Rst_z));
localangle = std::acos(localangle / R/slopeR)*r2d;
Angle_Arr(i, j) = localangle;
}
}
}
// 保存结果
omp_set_lock(&lock);
rasterRC.saveImage(sar_r, startRId, 0, 1);
rasterRC.saveImage(sar_c, startRId, 0, 2);
if (localincAngleFlag) {
localincangleRaster.saveImage(Angle_Arr, startRId, 0, 1);
}
else {}
processNumber = processNumber + blockrows;
std::cout << "\rprocess bar:\t" << processNumber * 100.0 / rowcount << " % " << "\t\t\t";
if (progressDialog.maximum() <= processNumber) {
processNumber = progressDialog.maximum() - 1;
}
progressDialog.setValue(processNumber);
omp_unset_lock(&lock);
}
progressDialog.close();
return ErrorCode::SUCCESS;
}
ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath, QString outworkdir, double gridx, double gridy)
ErrorCode GF3OrthSLC( QString inRasterPath, QString inlooktablePath, QString outRasterPath)
{
gdalImage slcRaster(inRasterPath);//
gdalImage looktableRaster(inlooktablePath);//
gdalImage outRaster(outRasterPath);//
if (outRaster.height != looktableRaster.height || outRaster.width != looktableRaster.width) {
qDebug() << "look table size is not same as outRaster size"<< looktableRaster.height <<"!="<<outRaster.height<<" "<<looktableRaster.width<<"!="<<outRaster.width;
return ErrorCode::FAIL;
}
Eigen::MatrixXd slcImg = slcRaster.getData(0, 0, slcRaster.height, slcRaster.width, 1);
Eigen::MatrixXd sar_r = looktableRaster.getData(0, 0, looktableRaster.height, looktableRaster.width, 1);
Eigen::MatrixXd sar_c = looktableRaster.getData(0, 0, looktableRaster.height, looktableRaster.width, 1);
Eigen::MatrixXd outImg = outRaster.getData(0, 0, outRaster.height, outRaster.width, 1);
long lookRows = sar_r.rows();
long lookCols = sar_r.cols();
long slcRows = slcImg.rows();
long slcCols = slcImg.cols();
long lastr = 0;
long nextr = 0;
long lastc = 0;
long nextc = 0;
double Bileanervalue = 0;
// 插值
Landpoint p0{ 0,0,0 }, p11{ 0,0,0 }, p21{ 0,0,0 }, p12{ 0,0,0 }, p22{ 0,0,0 }, p{ 0,0,0 };
QProgressDialog progressDialog(u8"正射图像生成中", u8"终止", 0, lookRows);
progressDialog.setWindowTitle(u8"正射图像生成中");
progressDialog.setWindowModality(Qt::WindowModal);
progressDialog.setAutoClose(true);
progressDialog.setValue(0);
progressDialog.setMaximum(lookRows);
progressDialog.setMinimum(0);
progressDialog.show();
for (long i = 0; i < lookRows; i++) {
for (long j = 0; j < lookCols; j++) {
p0 = {sar_r(i,j),sar_c(i,j),0};
lastr = std::floor(sar_r(i, j));
nextr = std::ceil(sar_r(i, j));
lastc = std::floor(sar_c(i, j));
nextc = std::ceil(sar_c(i, j));
if (lastr < 0 || lastc < 0 || nextr >= slcRows || nextc >= slcCols) {
continue;
}
p11 = { sar_r(i,j),sar_c(i,j),0 };
p11 = Landpoint{ 0,0,slcImg(lastr,lastc) };
p21 = Landpoint{ 0,1,slcImg(nextr,lastc) };
p12 = Landpoint{ 1,0,slcImg(lastr,nextc) };
p22 = Landpoint{ 1,1,slcImg(nextr,nextc) };
Bileanervalue = Bilinear_interpolation(p0, p11, p21, p12, p22);
outImg(i, j) = Bileanervalue;
}
progressDialog.setValue(i);
}
outRaster.saveImage(outImg, 0, 0, 1);
progressDialog.close();
return ErrorCode::SUCCESS;
}
ErrorCode GF3RDProcess(QString inxmlPath, QString indemPath, QString outworkdir, double gridx, double gridy)
{
SARSimulationImageL1Dataset l1dataset;
l1dataset.Open(inxmlPath);
DemBox box = l1dataset.getExtend();
double dlon = 0.1*(box.max_lon - box.min_lon);
double dlon = 0.1 * (box.max_lon - box.min_lon);
double dlat = 0.1 * (box.max_lat - box.min_lat);
double minlat = box.min_lat - dlat;
@ -338,18 +589,20 @@ ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath, QString out
// 裁剪DEM
QString filename = getFileNameWidthoutExtend(l1dataset.getxmlFilePath());
filename = filename.right(5);
QString clipDEMPath = JoinPath(outworkdir, getFileNameWidthoutExtend(indemPath) + filename+"_clip.tif");
QString clipDEMPath = JoinPath(outworkdir, getFileNameWidthoutExtend(indemPath) + filename + "_clip.tif");
cropRasterByLatLon(indemPath.toUtf8().constData(), clipDEMPath.toUtf8().constData(), minlon, maxlon, minlat, maxlat);
QString resampleDEMPath = JoinPath(outworkdir, getFileNameWidthoutExtend(indemPath) + filename + "_clip_resample.tif");
resampleRaster(clipDEMPath.toUtf8().constData(), resampleDEMPath.toUtf8().constData(), gridx, gridy);
QString outlooktablePath= JoinPath(outworkdir, getFileNameWidthoutExtend(l1dataset.getxmlFilePath()) + "_looktable.tif");
if (GF3RDCreateLookTable(inxmlPath, resampleDEMPath, outworkdir, outlooktablePath) != ErrorCode::SUCCESS) {
qDebug() << "查找表生成错误：\t"+ getFileNameWidthoutExtend(inxmlPath);
QString outlooktablePath = JoinPath(outworkdir, getFileNameWidthoutExtend(l1dataset.getxmlFilePath()) + "_looktable.tif");
QString outlocalAnglePath = JoinPath(outworkdir, getFileNameWidthoutExtend(l1dataset.getxmlFilePath()) + "_localAngle.tif");
QString outOrthPath=JoinPath(outworkdir, getFileNameWidthoutExtend(l1dataset.getxmlFilePath()) + "_orth.tif");
if (GF3RDCreateLookTable(inxmlPath, resampleDEMPath, outworkdir, outlooktablePath, outlocalAnglePath,true) != ErrorCode::SUCCESS) {
qDebug() << "查找表生成错误：\t" + getFileNameWidthoutExtend(inxmlPath);
return ErrorCode::FAIL;
}
return ErrorCode::SUCCESS;
}

21
GF3ProcessToolbox/GF3CalibrationAndGeocodingClass.h
View File

@ -1,10 +1,12 @@
#pragma once
#include "../BaseLibraryCPP/BaseConstVariable.h"
#include "../BaseLibraryCPP/ImageOperatorBase.h"
#include "../BaseLibraryCPP/LogInfoCls.h"
#include "../BaseLibraryCPP/SARSimulationImageL1.h"
#include "BaseConstVariable.h"
#include "ImageOperatorBase.h"
#include "LogInfoCls.h"
#include "SARSimulationImageL1.h"
#include "SatelliteOribtModel.h"
#include <QString>
// 数据定标
ErrorCode GF3CalibrationRaster(QString inRasterPath, QString outRasterPath, double Qualifyvalue, double calibrationConst);
@ -23,12 +25,17 @@ ErrorCode Complex2dBRaster(QString inComplexPath, QString outRasterPath);
ErrorCode ResampleDEM(QString indemPath, QString outdemPath, double gridx, double gridy);
// RD Ëã·¨Àà
ErrorCode RD_PSTN(double& refrange,double& lamda, double& timeR, double& R, double& tx, double& ty, double& tz, double& slopex, double& slopey, double& slopez, PolyfitSatelliteOribtModel& polyfitmodel, SatelliteOribtNode& node,double& d0,double& d1, double& d2, double& d3, double& d4);
//创建查找表
ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath, QString outworkdir, QString outlooktablePath);
ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath, QString outworkdir, QString outlooktablePath, QString outLocalIncidenceAnglePath, bool localincAngleFlag=false);
ErrorCode GF3OrthSLC( QString inRasterPath, QString inlooktablePath, QString outRasterPath);
// 正射处理流程
ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath, QString outworkdir, double gridx, double gridy);
ErrorCode GF3RDProcess(QString inxmlPath, QString indemPath, QString outworkdir, double gridx, double gridy);

4
GF3ProcessToolbox/QComplex2AmpPhase.cpp
View File

@ -1,8 +1,8 @@
#include "QComplex2AmpPhase.h"
#include <QtWidgets>
#include <QFileDialog>
#include "../BaseLibraryCPP/FileOperator.h"
#include "../BaseLibraryCPP/BaseTool.h"
#include "FileOperator.h"
#include "BaseTool.h"
#include "GF3CalibrationAndGeocodingClass.h"
QComplex2AmpPhase::QComplex2AmpPhase(QWidget *parent)

2
GF3ProcessToolbox/SateOrbit.h
View File

@ -7,7 +7,7 @@
//#define EIGEN_VECTORIZE_SSE4_2
//#include <mkl.h>
// 本地方法
#include "../BaseLibraryCPP/BaseTool.h"
#include "BaseTool.h"
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>

2
GF3ProcessToolbox/SatelliteGF3xmlParser.h
View File

@ -3,7 +3,7 @@
#include <QFile>
#include <QDomDocument>
#include <QDebug>
#include "../BaseLibraryCPP/SARSimulationImageL1.h"
#include "SARSimulationImageL1.h"

2
GF3ProcessToolbox/WGS84_J2000.cpp
View File

@ -15,5 +15,5 @@
#include <opencv2/core/eigen.hpp>
#include <opencv2/features2d.hpp>
#include <fstream>
#include "../BaseLibraryCPP/BaseTool.h"
#include "BaseTool.h"
#include"WGS84_J2000.h"

2
GF3ProcessToolbox/WGS84_J2000.h
View File

@ -6,5 +6,5 @@
#include <Eigen/Core>
#include <Eigen/Dense>
#include <time.h>
#include "../BaseLibraryCPP/BaseTool.h"
#include "BaseTool.h"
#include <string>

38
GF3ProcessToolbox/simptsn.cpp
View File

@ -14,7 +14,7 @@
#include <time.h>
#include <boost/filesystem.hpp>
#include <omp.h>
#include "../BaseLibraryCPP/BaseTool.h"
#include "BaseTool.h"
#include "simptsn.h"
#include "SateOrbit.h"
#include <gdal_utils.h>
@ -28,7 +28,7 @@
#include <vector>
#include <gdalgrid.h>
#include <boost/filesystem.hpp>
#include "../BaseLibraryCPP/ImageOperatorBase.h"
#include "ImageOperatorBase.h"
#include <opencv2/opencv.hpp>
@ -88,7 +88,7 @@ PSTNAlgorithm::PSTNAlgorithm(QString infile_path)
getline(infile, buf); this->PRF = stod(buf); qDebug() << "PRF\t" << this->PRF ;
getline(infile, buf); this->refrange = stod(buf); qDebug() << "refrange\t" << this->refrange ;
getline(infile, buf); this->fs = stod(buf); qDebug() << "fs\t" << this->fs ;
this->widthspace = (1 / this->fs) * LIGHTSPEED;
getline(infile, buf); this->doppler_paramenter_number = stoi(buf); qDebug() << "doppler_paramenter_number\t" << this->doppler_paramenter_number ;
// 读取多普勒系数
@ -275,8 +275,8 @@ Eigen::MatrixX<double> PSTNAlgorithm::PSTN(Eigen::MatrixX<double> landpoints, do
satepoint = satestate(Eigen::all, { 0,1,2 });
sate_land = satepoint - landpoint;
R1 = (sate_land.array().pow(2).array().rowwise().sum().array().sqrt())[0];;
//SARPoint(j, 2) = getIncAngle(satePoints, landPoint);
//SARPoint(j, 2) = (R1 - this->R0 - (this->refrange - this->R0) / 3) / this->widthspace;
SARPoint(j, 2) = getIncAngle(satePoints, landPoint);
SARPoint(j, 2) = (R1 - this->R0 - (this->refrange - this->R0) / 3) / this->widthspace;
SARPoint(j, 2) = getRangeColumn(R1, this->R0, this->fs, this->lamda);// (R1 - this->R0) / this->widthspace; //
}
//qDebug() <<"SARPoints(200, 0):\t" << SARPoint(200, 0) << "\t" << SARPoint(200, 1) << "\t" << SARPoint(200, 2) << "\t" ;
@ -1861,7 +1861,7 @@ int simProcess::correct_ati(QString orth_rc_path, QString dem_path, QString out_
gdalImage dem_img(dem_path);
gdalImage orth_rc(orth_rc_path);
//gdalImage ori_inclocal = CreategdalImage(out_inclocal_path, this->pstn.height, this->pstn.width, 1, orth_rc.gt, orth_rc.projection, false);// 注意这里保留仿真结果
gdalImage ori_inclocal = CreategdalImage(out_inclocal_path, this->pstn.height, this->pstn.width, 1, orth_rc.gt, orth_rc.projection, false);// 注意这里保留仿真结果
gdalImage ori_dem = CreategdalImage(out_dem_path, this->pstn.height, this->pstn.width, 1, orth_rc.gt, orth_rc.projection, false);// 注意这里保留仿真结果
gdalImage ori_count = CreategdalImage(out_count_path, this->pstn.height, this->pstn.width, 1, orth_rc.gt, orth_rc.projection, false);// 注意这里保留仿真
@ -1873,11 +1873,11 @@ int simProcess::correct_ati(QString orth_rc_path, QString dem_path, QString out_
int start_ids = 0;
do {
Eigen::MatrixXd sar_dem = ori_dem.getData(start_ids, 0, line_invert, this->pstn.width, 1);
//Eigen::MatrixXd sar_inclocal = ori_inclocal.getData(start_ids, 0, line_invert, this->pstn.width, 1);
Eigen::MatrixXd sar_inclocal = ori_inclocal.getData(start_ids, 0, line_invert, this->pstn.width, 1);
Eigen::MatrixXd sar_count = ori_count.getData(start_ids, 0, line_invert, this->pstn.width, 1);
sar_dem = sar_dem.array() * 0;
//sar_inclocal = sar_inclocal.array() * 0;
sar_inclocal = sar_inclocal.array() * 0;
sar_count = sar_count.array() * 0;
ori_dem.saveImage(sar_dem, start_ids, 0, 1);
@ -1931,7 +1931,7 @@ int simProcess::correct_ati(QString orth_rc_path, QString dem_path, QString out_
sar_dem = ori_dem.getData(r_min, 0, r_max - r_min + 1, this->pstn.width, 1).cast<double>();
sar_count = ori_count.getData(r_min, 0, r_max - r_min + 1, this->pstn.width, 1).cast<double>();
//sar_inclocal = ori_inclocal.getData(r_min, 0, r_max - r_min + 1, this->pstn.width, 1).cast<double>(); // 列
sar_inclocal = ori_inclocal.getData(r_min, 0, r_max - r_min + 1, this->pstn.width, 1).cast<double>(); // 列
rowcount = dem.rows();
colcount = dem.cols();
@ -1958,7 +1958,7 @@ int simProcess::correct_ati(QString orth_rc_path, QString dem_path, QString out_
r = sar_r(i, j);
c = sar_c(i, j);
if (r >= 0 && r < this->pstn.height && c >= 0 && c < this->pstn.width) {
//sar_inclocal(r - r_min, c) += localincangle;
sar_inclocal(r - r_min, c) += localincangle;
sar_count(r - r_min, c) += 1;
sar_dem(r - r_min, c) += dem(i, j);
}
@ -1981,11 +1981,11 @@ int simProcess::correct_ati(QString orth_rc_path, QString dem_path, QString out_
int start_ids = 0;
do {
Eigen::MatrixXd sar_dem = ori_dem.getData(start_ids, 0, line_invert, this->pstn.width, 1);
//Eigen::MatrixXd sar_inclocal = ori_inclocal.getData(start_ids, 0, line_invert, this->pstn.width, 1);
Eigen::MatrixXd sar_inclocal = ori_inclocal.getData(start_ids, 0, line_invert, this->pstn.width, 1);
Eigen::MatrixXd sar_count = ori_count.getData(start_ids, 0, line_invert, this->pstn.width, 1);
sar_dem = sar_dem.array() / sar_count.array();
//sar_inclocal = sar_inclocal.array() / sar_count.array();
sar_inclocal = sar_inclocal.array() / sar_count.array();
ori_dem.saveImage(sar_dem, start_ids, 0, 1);
//ori_inclocal.saveImage(sar_inclocal, start_ids, 0, 1);
@ -2611,17 +2611,17 @@ void simProcess::interpolation_GTC_sar_sigma(QString in_rc_wgs84_path, QString i
int start_ids = 0;
do {
Eigen::MatrixXd sar_a = sar_img.getData(start_ids, 0, line_invert, sar_rc.width, 1);
//Eigen::MatrixXd sar_b = sar_img.getData(start_ids, 0, line_invert, sar_rc.width, 2);
Eigen::MatrixXd sar_b = sar_img.getData(start_ids, 0, line_invert, sar_rc.width, 2);
sar_a = sar_a.array() * 0 - 9999;
//sar_b = sar_b.array() * 0 - 9999;
sar_b = sar_b.array() * 0 - 9999;
sar_img.saveImage(sar_a, start_ids, 0, 1);
//sar_img.saveImage(sar_a, start_ids, 0, 2);
sar_img.saveImage(sar_a, start_ids, 0, 2);
start_ids = start_ids + line_invert;
} while (start_ids < sar_rc.height);
sar_img.setNoDataValue(-9999, 1);
//sar_img.setNoDataValue(-9999, 2);
sar_img.setNoDataValue(-9999, 2);
}
// 正式计算插值
@ -2672,7 +2672,7 @@ void simProcess::interpolation_GTC_sar_sigma(QString in_rc_wgs84_path, QString i
ori_a = Eigen::MatrixXd(1, 1);
//ori_b = Eigen::MatrixXd(1, 1); // 释放内存
Eigen::MatrixXd sar_a = sar_img.getData(start_ids, 0, line_invert, sar_rc.width, 1);
//Eigen::MatrixXd sar_b = sar_img.getData(start_ids, 0, line_invert, sar_rc.width, 2);
Eigen::MatrixXd sar_b = sar_img.getData(start_ids, 0, line_invert, sar_rc.width, 2);
sar_r = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 1);
sar_c = sar_rc.getData(start_ids, 0, line_invert, sar_rc.width, 2);
row_count = sar_r.rows();
@ -2909,11 +2909,11 @@ void simProcess::interpolation_GTC_sar_sigma(QString in_rc_wgs84_path, QString i
complex<double> interpolation_value = Cubic_Convolution_interpolation(r - r1, c - c1, img_block);
//
sar_a(i, j) = interpolation_value.real();
//sar_b(i, j) = interpolation_value.imag();
sar_b(i, j) = interpolation_value.imag();
}
}
sar_img.saveImage(sar_a, start_ids, 0, 1);
//sar_img.saveImage(sar_b, start_ids, 0, 2);
sar_img.saveImage(sar_b, start_ids, 0, 2);
start_ids = start_ids + line_invert;
} while (start_ids < sar_rc.height);
}

4
GF3ProcessToolbox/simptsn.h
View File

@ -7,7 +7,7 @@
//#include <mkl.h>
// 本地方法
#include "../BaseLibraryCPP/BaseTool.h"
#include "BaseTool.h"
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
@ -55,6 +55,8 @@ public:
int height; // 影像的高
int width;
int widthspace;
// 入射角
double near_in_angle;// 近入射角
double far_in_angle;// 远入射角

2
Imageshow/ImageShowDialogClass.h
View File

@ -3,7 +3,7 @@
#ifndef IMAGESHOWDIALOGCLASS_H
#define IMAGESHOWDIALOGCLASS_H
#include "../BaseLibraryCPP/BaseConstVariable.h"
#include "BaseConstVariable.h"
#include <QDialog>
#include <QLabel>
#include "qcustomplot.h"

2
QMergeRasterProcessDialog.cpp
View File

@ -1,6 +1,6 @@
#include "QMergeRasterProcessDialog.h"
#include "QToolProcessBarDialog.h"
#include "BaseLibraryCPP/ImageOperatorBase.h"
#include "ImageOperatorBase.h"
#include <QString>
#include <QVector>
#include <QFileDialog>

10
QSimulationRTPCGUI.cpp Normal file
View File

@ -0,0 +1,10 @@
#include "QSimulationRTPCGUI.h"
QSimulationRTPCGUI::QSimulationRTPCGUI(QWidget *parent)
: QDialog(parent)
{
ui.setupUi(this);
}
QSimulationRTPCGUI::~QSimulationRTPCGUI()
{}

16
QSimulationRTPCGUI.h Normal file
View File

@ -0,0 +1,16 @@
#pragma once
#include <QDialog>
#include "ui_QSimulationRTPCGUI.h"
class QSimulationRTPCGUI : public QDialog
{
Q_OBJECT
public:
QSimulationRTPCGUI(QWidget *parent = nullptr);
~QSimulationRTPCGUI();
private:
Ui::QSimulationRTPCGUIClass ui;
};

20
QSimulationRTPCGUI.ui Normal file
View File

@ -0,0 +1,20 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>QSimulationRTPCGUIClass</class>
<widget class="QDialog" name="QSimulationRTPCGUIClass">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>600</width>
<height>400</height>
</rect>
</property>
<property name="windowTitle">
<string>RTPC回波仿真</string>
</property>
</widget>
<layoutdefault spacing="6" margin="11"/>
<resources/>
<connections/>
</ui>

2
QToolProcessBarDialog.h
View File

@ -2,7 +2,7 @@
#include <QDialog>
#include "ui_QToolProcessBarDialog.h"
#include "BaseLibraryCPP/ImageOperatorBase.h"
#include "ImageOperatorBase.h"
class QToolProcessBarDialog : public QDialog, public ShowProessAbstract
{
Q_OBJECT

5
RasterProcessTool.vcxproj
View File

@ -63,8 +63,10 @@
</ImportGroup>
<PropertyGroup Label="UserMacros" />
<PropertyGroup Condition="'$(Configuration)|$(Platform)' == 'Debug|x64'">
<IncludePath>.\SimulationSAR;.\GF3ProcessToolbox;.\BaseLibraryCPP;$(IncludePath)</IncludePath>
</PropertyGroup>
<PropertyGroup Condition="'$(Configuration)|$(Platform)' == 'Release|x64'">
<IncludePath>.\SimulationSAR;.\GF3ProcessToolbox;.\BaseLibraryCPP;$(oneMKLIncludeDir);$(IncludePath)</IncludePath>
</PropertyGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
<ClCompile>
@ -106,6 +108,7 @@
</Link>
</ItemDefinitionGroup>
<ItemGroup>
<ClCompile Include="QSimulationRTPCGUI.cpp" />
<ClCompile Include="SimulationSAR\RTPCProcessCls.cpp" />
<ClCompile Include="SimulationSAR\SARSatelliteSimulationAbstractCls.cpp" />
<ClCompile Include="SimulationSAR\SARSimulationTaskSetting.cpp" />
@ -119,6 +122,7 @@
<QtUic Include="Imageshow\ImageShowDialogClass.ui" />
<QtUic Include="Imageshow\qcustomplot.ui" />
<QtUic Include="QMergeRasterProcessDialog.ui" />
<QtUic Include="QSimulationRTPCGUI.ui" />
<QtUic Include="QToolProcessBarDialog.ui" />
<QtUic Include="RasterProcessTool.ui" />
<QtMoc Include="RasterProcessTool.h" />
@ -155,6 +159,7 @@
<ClInclude Include="BaseLibraryCPP\EchoDataFormat.h" />
<ClInclude Include="BaseLibraryCPP\FileOperator.h" />
<ClInclude Include="BaseLibraryCPP\GeoOperator.h" />
<QtMoc Include="QSimulationRTPCGUI.h" />
<ClInclude Include="SimulationSAR\RTPCProcessCls.h" />
<ClInclude Include="SimulationSAR\SARSatelliteSimulationAbstractCls.h" />
<ClInclude Include="SimulationSAR\SARSimulationTaskSetting.h" />

9
RasterProcessTool.vcxproj.filters
View File

@ -133,6 +133,9 @@
<ClCompile Include="SimulationSAR\SigmaDatabase.cpp">
<Filter>SimulationSAR</Filter>
</ClCompile>
<ClCompile Include="QSimulationRTPCGUI.cpp">
<Filter>Source Files</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<ClInclude Include="BaseLibraryCPP\BaseConstVariable.h">
@ -218,6 +221,9 @@
<QtMoc Include="GF3ProcessToolbox\QRDOrthProcessClass.h">
<Filter>GF3ProcessToolbox</Filter>
</QtMoc>
<QtMoc Include="QSimulationRTPCGUI.h">
<Filter>Header Files</Filter>
</QtMoc>
</ItemGroup>
<ItemGroup>
<QtUic Include="QMergeRasterProcessDialog.ui">
@ -241,5 +247,8 @@
<QtUic Include="GF3ProcessToolbox\QRDOrthProcessClass.ui">
<Filter>GF3ProcessToolbox</Filter>
</QtUic>
<QtUic Include="QSimulationRTPCGUI.ui">
<Filter>Form Files</Filter>
</QtUic>
</ItemGroup>
</Project>

3
SimulationSAR/SARSimulationTaskSetting.cpp
View File

@ -263,8 +263,6 @@ std::shared_ptr<AbstractSARSatelliteModel> ReadSimulationSettingsXML(QString xml
std::shared_ptr<SARSimulationTaskSetting> taskSetting = std::make_shared<SARSimulationTaskSetting>();
//SARSimulationTaskSetting taskSetting;
double starttimestamp = convertToMilliseconds(startTime.text().trimmed().toStdString());
double endtimestamp = convertToMilliseconds(endTime.text().trimmed().toStdString());
@ -291,7 +289,6 @@ std::shared_ptr<AbstractSARSatelliteModel> ReadSimulationSettingsXML(QString xml
}
taskSetting->setCenterLookAngle(centerLookAngle.text().toDouble());
taskSetting->setSARImageStartTime(starttimestamp); // 成像开始时间
taskSetting->setSARImageEndTime(endtimestamp); // 成像终止时间

49
SimulationSAR/SatelliteOribtModel.cpp
View File

@ -37,6 +37,39 @@ PolyfitSatelliteOribtModel::~PolyfitSatelliteOribtModel()
}
QString PolyfitSatelliteOribtModel::getSatelliteOribtModelParamsString()
{
QString result = "";
result += this->polyfitPx.size()-1 + "\n";
result += "----------- poly Position X -------------------------\n";
for (long i = 0; i < this->polyfitPx.size(); i++) { // 求解平方
result += QString::number(this->polyfitPx[i], 'e', 6) + "\n";
}
result += "----------- poly Position Y -------------------------\n";
for (long i = 0; i < this->polyfitPx.size(); i++) { // 求解平方
result += QString::number(this->polyfitPy[i], 'e', 6) + "\n";
}
result += "----------- poly Position Z -------------------------\n";
for (long i = 0; i < this->polyfitPx.size(); i++) { // 求解平方
result += QString::number(this->polyfitPz[i], 'e', 6) + "\n";
}
result += "----------- poly Position Vector X ------------------\n";
for (long i = 0; i < this->polyfitPx.size(); i++) { // 求解平方
result += QString::number(this->polyfitVx[i], 'e', 6) + "\n";
}
result += "----------- poly Position Vector Y ------------------\n";
for (long i = 0; i < this->polyfitPx.size(); i++) { // 求解平方
result += QString::number(this->polyfitVy[i], 'e', 6) + "\n";
}
result += "----------- poly Position Vector Z ------------------\n";
for (long i = 0; i < this->polyfitPx.size(); i++) { // 求解平方
result += QString::number(this->polyfitVz[i], 'e', 6) + "\n";
}
result+= "------------------------------------------------------\n";
return result;
}
SatelliteOribtNode PolyfitSatelliteOribtModel::getSatelliteOribtNode(double& timeFromStartTime, bool& antAzAngleFlag)
{
// 位置、速度
@ -236,22 +269,6 @@ void PolyfitSatelliteOribtModel::addOribtNode(SatelliteOribtNode node)
}
ErrorCode polynomial_fit(const std::vector<double>& x, const std::vector<double>& y, int degree, std::vector<double>& out_factor, double& out_chisq) {
int xyLength = x.size();
double* xdata = new double[xyLength];
double* ydata = new double[xyLength];
for (long i = 0; i < xyLength; i++) {
xdata[i] = x[i];
ydata[i] = y[i];
}
ErrorCode state= polyfit(xdata, ydata,xyLength, degree, out_factor, out_chisq);
delete[] xdata;
delete[] ydata; // ÊÍ·ÅÄÚ´æ
return state;
}
ErrorCode PolyfitSatelliteOribtModel::polyFit(int polynum/*=3*/, bool input_timeFromReferenceTime)
{

5
SimulationSAR/SatelliteOribtModel.h
View File

@ -20,7 +20,8 @@ class PolyfitSatelliteOribtModel:public AbstractSatelliteOribtModel
public:
PolyfitSatelliteOribtModel();
~PolyfitSatelliteOribtModel() override;
public:
QString getSatelliteOribtModelParamsString();
public: // 卫星轨道节点
virtual SatelliteOribtNode getSatelliteOribtNode(double& timeFromStartTime, bool& antAzAngleFlag) override; // 获取轨道节点
@ -100,8 +101,6 @@ private:
};
ErrorCode polynomial_fit(const std::vector<double>& x, const std::vector<double>& y, int degree, std::vector<double>& out_factor, double& out_chisq);
std::shared_ptr<AbstractSatelliteOribtModel> CreataPolyfitSatelliteOribtModel(std::vector<SatelliteOribtNode> &nodes, long double startTime,int polynum=3);