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

613 lines
23 KiB
C++
Raw Blame History

#include "GF3CalibrationAndGeocodingClass.h"
#include "SatelliteGF3xmlParser.h"
#include <QRegularExpression>
#include <QMessageBox>
#include "GF3PSTNClass.h"
#include "boost/asio.hpp"
#include <boost/thread.hpp>
#include <thread>
#include <qprogressdialog.h>
ErrorCode GF3CalibrationRaster(QString inRasterPath, QString outRasterPath, double Qualifyvalue, double calibrationConst)
{
gdalImage imgraster(inRasterPath.toUtf8().constData());
if (!isExists(outRasterPath)) {
gdalImageComplex outraster = CreategdalImageComplex(outRasterPath, imgraster.height, imgraster.width, 1, imgraster.gt, imgraster.projection, false, true);
}
gdalImageComplex outraster(outRasterPath);
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);
Eigen::MatrixXcd imgArr = Eigen::MatrixXcd::Zero(blocklines, imgraster.width);
double quayCoff = Qualifyvalue * 1.0 / 32767;
double caliCoff = std::pow(10.0, (calibrationConst * 1.0 / 20));
qDebug() << u8"¶¨±êϵÊý:\t" << quayCoff / caliCoff;
long startrow = 0;
for (startrow = 0; startrow < imgraster.height; startrow = startrow + blocklines) {
imgArrb1 = imgraster.getData(startrow, 0, blocklines, imgraster.width, 1);
imgArrb2 = imgraster.getData(startrow, 0, blocklines, imgraster.width, 2);
imgArr = outraster.getDataComplex(startrow, 0, blocklines, outraster.width, 1);
imgArr.real() = imgArrb1.array() * quayCoff / caliCoff;
imgArr.imag() = imgArrb2.array() * quayCoff / caliCoff;
outraster.saveImage(imgArr, startrow, 0, 1);
}
qDebug() << u8"Ó°ÏñдÈëµ½£º" << outRasterPath;
return ErrorCode::SUCCESS;
}
ErrorCode ImportGF3L1ARasterToWorkspace(QString inMetaxmlPath, QString inRasterPath, QString outWorkDirPath, POLARTYPEENUM polsartype)
{
SatelliteGF3xmlParser gf3xml;
gf3xml.loadFile(inMetaxmlPath);
QString filename = getFileNameWidthoutExtend(inRasterPath);
SARSimulationImageL1Dataset l1dataset;
l1dataset.OpenOrNew(outWorkDirPath, filename, gf3xml.height, gf3xml.width);
l1dataset.setSateAntPos(SatellitePos2SatelliteAntPos(gf3xml.antposs));
l1dataset.setCenterFreq(gf3xml.RadarCenterFrequency);
l1dataset.setcolCount(gf3xml.width);
l1dataset.setrowCount(gf3xml.height);
l1dataset.setNearRange(gf3xml.nearRange);
l1dataset.setRefRange(gf3xml.refRange);
l1dataset.setFs(gf3xml.eqvFs);
l1dataset.setPRF(gf3xml.eqvPRF);
l1dataset.setCenterAngle((gf3xml.incidenceAngleNearRange + gf3xml.incidenceAngleFarRange) / 2.0); // ÈëÉä½Ç
l1dataset.setDopplerParametersReferenceTime(gf3xml.DopplerParametersReferenceTime);
l1dataset.setTotalProcessedAzimuthBandWidth(gf3xml.TotalProcessedAzimuthBandWidth);
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.setLatitudeCenter(gf3xml.latitude_center);
l1dataset.setLongitudeCenter(gf3xml.longitude_center);
l1dataset.setLatitudeTopLeft(gf3xml.latitude_topLeft);
l1dataset.setLongitudeTopLeft(gf3xml.longitude_topLeft);
l1dataset.setLatitudeTopRight(gf3xml.latitude_topRight);
l1dataset.setLongitudeTopRight(gf3xml.longitude_topRight);
l1dataset.setLatitudeBottomLeft(gf3xml.latitude_bottomLeft);
l1dataset.setLongitudeBottomLeft(gf3xml.longitude_bottomLeft);
l1dataset.setLatitudeBottomRight(gf3xml.latitude_bottomRight);
l1dataset.setLongitudeBottomRight(gf3xml.longitude_bottomRight);
l1dataset.setStartImageTime(gf3xml.start);
l1dataset.setEndImageTime(gf3xml.end);
l1dataset.saveToXml();
QString outRasterpath = l1dataset.getImageRasterPath();
ErrorCode errorcode = ErrorCode::SUCCESS;
qDebug() << u8"µ¼ÈëÊý¾Ý:\t" << inRasterPath;
switch (polsartype)
{
case POLARHH:
qDebug() << u8"¼«»¯£ºHH";
errorcode = GF3CalibrationRaster(inRasterPath, outRasterpath, gf3xml.HH_QualifyValue, gf3xml.HH_CalibrationConst);
break;
case POLARHV:
qDebug() << u8"¼«»¯£ºHH";
errorcode = GF3CalibrationRaster(inRasterPath, outRasterpath, gf3xml.HV_QualifyValue, gf3xml.HV_CalibrationConst);
break;
case POLARVH:
qDebug() << u8"¼«»¯£ºVH";
errorcode = GF3CalibrationRaster(inRasterPath, outRasterpath, gf3xml.VH_QualifyValue, gf3xml.VH_CalibrationConst);
break;
case POLARVV:
qDebug() << u8"¼«»¯£ºVV";
errorcode = GF3CalibrationRaster(inRasterPath, outRasterpath, gf3xml.VV_QualifyValue, gf3xml.VV_CalibrationConst);
break;
default:
break;
}
qDebug() << u8"µ¼ÈëÊý¾Ý״̬£º\t" << QString::fromStdString(errorCode2errInfo(errorcode));
return errorcode;
}
QVector<QString> SearchGF3DataTiff(QString inMetaxmlPath)
{
// Ö¸¶¨Â·¾¶
QString xmlpath = inMetaxmlPath; // Ì滻ΪÄãµÄʵ¼Ê·¾¶
QString absPath = QFileInfo(xmlpath).absolutePath();
// »ñȡ·¾¶ËùÔÚµÄĿ¼
QDir directory(absPath);
if (!directory.exists()) {
qDebug() << u8"Directory does not exist:" << directory.absolutePath();
return QVector<QString>(0);
}
// ÉèÖùýÂËÆ÷ÒÔ½öÁгö .tif Îļþ
QStringList filters;
filters << "*.tif" << "*.TIF" << "*.tiff" << "*.TIFF";
QStringList files = directory.entryList(filters, QDir::Files);
qDebug() << u8"TIFF Files in the directory" << directory.absolutePath() << ":";
QVector<QString> filepath(0);
for (long i = 0; i < files.count(); i++) {
filepath.append(JoinPath(absPath, files[i]));
}
return filepath;
}
POLARTYPEENUM getDatasetGF3FilePolsarType(QString fileName)
{
// ÌáÈ¡¼«»¯Çé¿ö
QRegularExpression re("_([A-Z]{2})_");
QRegularExpressionMatch match = re.match(fileName);
QString polarization = "";
if (match.hasMatch()) {
polarization = match.captured(1);
polarization = polarization.toLower().replace("_", "");
qDebug() << u8"Polarization extracted:" << polarization;
if (polarization == "hh") {
return POLARTYPEENUM::POLARHH;
}
else if (polarization == "hv") {
return POLARTYPEENUM::POLARHV;
}
else if (polarization == "vh") {
return POLARTYPEENUM::POLARVH;
}
else if (polarization == "vv") {
return POLARTYPEENUM::POLARVV;
}
else {
return POLARTYPEENUM::POLARUNKOWN;
}
}
else {
qDebug() << u8"No polarization found in the path.";
return POLARTYPEENUM::POLARUNKOWN;
}
return POLARTYPEENUM::POLARUNKOWN;
}
ErrorCode ImportGF3L1AProcess(QString inMetaxmlPath, QString outWorkDirPath)
{
QVector<QString> tiffpaths = SearchGF3DataTiff(inMetaxmlPath);
ErrorCode errorcode = ErrorCode::SUCCESS;
for (long i = 0; i < tiffpaths.count(); i++) {
QString tiffpath = tiffpaths[i];
QString filenamewithoutextern = getFileNameWidthoutExtend(tiffpath);
QString filename = getFileNameFromPath(tiffpath);
POLARTYPEENUM poltype = getDatasetGF3FilePolsarType(filename);
switch (poltype)
{
case POLARHH:
errorcode = ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARHH);
break;
case POLARHV:
errorcode = ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARHV);
break;
case POLARVH:
errorcode = ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARVH);
break;
case POLARVV:
errorcode = ImportGF3L1ARasterToWorkspace(inMetaxmlPath, tiffpath, outWorkDirPath, POLARTYPEENUM::POLARVV);
break;
case POLARUNKOWN:
break;
default:
break;
}
if (errorcode == ErrorCode::SUCCESS) {
continue;
}
else {
QMessageBox::warning(nullptr, u8"´íÎó", u8"Êý¾Ýµ¼Èë´íÎó");
break;
}
}
return errorcode;
}
ErrorCode RD_PSTN(double& refrange, double& lamda, double& timeR, double& R, double& tx, double& ty, double& tz, double& slopex, double& slopey, double& slopez, GF3PolyfitSatelliteOribtModel& polyfitmodel, SatelliteOribtNode& node, double& d0, double& d1, double& d2, double& d3, double& d4)
{
double dt = 1e-4;
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 / R1) * ((node.Px - tx) * (node.Vx + EARTHWE * ty) + (node.Py - ty) * (node.Vy - EARTHWE * tx) + (node.Pz - tz) * (node.Vz - 0));
dplerTheory1 = (-2 / lamda / R1) * ((node.Px - tx) * (node.Vx + 0) + (node.Py - ty) * (node.Vy -0) + (node.Pz - tz) * (node.Vz - 0));
dplerR = (R1 - refrange) * 1e6 / LIGHTSPEED; // GF3 ¼ÆË㹫ʽ
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/R2) * ((node.Px - tx) * (node.Vx + EARTHWE * ty) + (node.Py - ty) * (node.Vy - EARTHWE * tx) + (node.Pz - tz) * (node.Vz - 0));
dplerTheory2 = (-2 / lamda/R2) * ((node.Px - tx) * (node.Vx +0) + (node.Py - ty) * (node.Vy -0) + (node.Pz - tz) * (node.Vz - 0));
//dplerR = (R2 - refrange) * 1e6 / LIGHTSPEED; // GF3 ¼ÆË㹫ʽ
//dplerNumber2 = d0 + dplerR * d1 + std::pow(dplerR, 2) * d2 + std::pow(dplerR, 3) * d3 + std::pow(dplerR, 4) * d4;
inct = dt*(dplerTheory2-dplerNumber1) / (dplerTheory2 - dplerTheory1);
if (std::abs(dplerNumber1 - dplerTheory2) < 1e-6 || std::abs(inct) < 1.0e-4) {
break;
}
inct = std::abs(inct) < 10 ?inct:inct*1e-2;
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);
QString resampleDEMXYZPath = JoinPath(outworkdir, indemfilename + "_XYZ.tif");
QString resampleDEMSLOPEPath = JoinPath(outworkdir, indemfilename + "_slopevector.tif");
DEM2XYZRasterAndSlopRaster(indemPath, resampleDEMXYZPath, resampleDEMSLOPEPath);
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 localincangleRaster;
if (localincAngleFlag) {
localincangleRaster = CreategdalImage(outLocalIncidenceAnglePath, demxyz.height, demxyz.width, 1, demxyz.gt, demxyz.projection, true, true);// X,Y,Z
}
SARSimulationImageL1Dataset l1dataset;
l1dataset.Open(inxmlPath);
// RD Ïà¹Ø²ÎÊý
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()*1e6;// Fs²ÉÑù
double prf = l1dataset.getPRF();// PRF ²ÉÑù
double nearRange = l1dataset.getNearRange();
double imagestarttime = l1dataset.getStartImageTime();
double imageendtime = l1dataset.getEndImageTime();
double refrange = l1dataset.getRefRange();
double lamda = (LIGHTSPEED*1.0e-9)/ l1dataset.getCenterFreq();
// ´òÓ¡²ÎÊý
qDebug() << "Fs:\t" << fs;
qDebug() << "prf:\t" << prf;
qDebug() << "nearRange:\t" << nearRange;
qDebug() << "imagestarttime:\t" << imagestarttime;
qDebug() << "imageendtime:\t" << imageendtime;
qDebug() << "refrange:\t" << refrange;
qDebug() << "lamda:\t" << lamda;
//´òÓ¡¶àÆÕÀÛ²ÎÊý
qDebug() << u8"-----------------------------------";
qDebug() << u8"¶àÆÕÀÕ²ÎÊý£º\n";
qDebug() << u8"d0:\t" << d0;
qDebug() << u8"d1:\t" << d1;
qDebug() << u8"d2:\t" << d2;
qDebug() << u8"d3:\t" << d3;
qDebug() << u8"d4:\t" << d4;
qDebug() << u8"-----------------------------------";
// ¹¹½¨¹ìµÀÄ£ÐÍ
GF3PolyfitSatelliteOribtModel polyfitmodel;
QVector < SatelliteAntPos > antposes = l1dataset.getXmlSateAntPos();
polyfitmodel.setSatelliteOribtStartTime((imagestarttime+imageendtime)/2);
for (long i = 0; i < antposes.size(); i++) {
if (antposes[i].time > imagestarttime - 5 && antposes[i].time < imageendtime + 5) {
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() << u8"satellite polyfit model params:\n";
qDebug() << polyfitmodel.getSatelliteOribtModelParamsString();
qDebug() << "-----------------------------------";
// ¿ªÊ¼¼ÆËã²éÕÒ±í
//1.¼ÆËã·Ö¿é
long cpucore_num = std::thread::hardware_concurrency();
long blocklineinit = Memory1MB / 8 / cpucore_num / 4 / l1dataset.getcolCount() * 8000;
blocklineinit = blocklineinit < 50 ? 50 : blocklineinit;
//2.µü´ú¼ÆËã
long colcount = rasterRC.width;//l1dataset.getcolCount();
long rowcount = rasterRC.height;//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);
for (startRId = 0; startRId < rowcount; startRId = startRId + blocklineinit) {
long blockline = blocklineinit;
if (startRId + blockline > rowcount) {
blockline = rowcount - startRId;
}
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);
long blockrows = sar_r.rows();
long blockcols = sar_r.cols();
#pragma omp parallel for
for (long i = 0; i < blockrows; i++) {
// ÖðÏñËصü´ú¼ÆËã
double timeR = 0;
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 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+ (imagestarttime + imageendtime) / 2 -imagestarttime) * prf;
sar_c(i, j) = ((R - nearRange) * 2 / LIGHTSPEED) * fs;
}
else {
sar_r(i, j) = -1;
sar_c(i, j) = -1;
}
}
}
// ±£´æ½á¹û
omp_set_lock(&lock);
rasterRC.saveImage(sar_r, startRId, 0, 1);
rasterRC.saveImage(sar_c, startRId, 0, 2);
omp_unset_lock(&lock);
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;
#pragma omp parallel for
for (long i = 0; i < blockrows; i++) {
// ÖðÏñËصü´ú¼ÆËã
double timeR = 0;
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 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);
if (localincAngleFlag) {
localincangleRaster.saveImage(Angle_Arr, startRId, 0, 1);
}
else {}
omp_unset_lock(&lock);
}
// ±£´æ½á¹û
omp_set_lock(&lock);
processNumber = processNumber + blockrows;
qDebug() << u8"\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 GF3OrthSLC( QString inRasterPath, QString inlooktablePath, QString outRasterPath)
{
gdalImage slcRaster(inRasterPath);//
gdalImage looktableRaster(inlooktablePath);//
gdalImage outRaster = CreategdalImage(outRasterPath, looktableRaster.height, looktableRaster.width, 1, looktableRaster.gt, looktableRaster.projection, true, true);// X,Y,Z
//gdalImage outRaster(outRasterPath);//
if (outRaster.height != looktableRaster.height || outRaster.width != looktableRaster.width) {
qDebug() << u8"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, 2);
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++) {
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;
}
p0 = { sar_c(i,j)-lastc, sar_r(i,j)-lastr,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 dlat = 0.1 * (box.max_lat - box.min_lat);
double minlat = box.min_lat - dlat;
double minlon = box.min_lon - dlon;
double maxlat = box.max_lat + dlat;
double maxlon = box.max_lon + dlon;
qDebug() << u8"Ó°Ïñ·¶Î§:";
qDebug() << u8"Bounding Box:";
qDebug() << u8"Latitude:" << minlat <<" - "<< maxlat;
qDebug() << u8"Longitude:" << minlon <<" - "<< maxlon;
double centerlat = (minlat + maxlat) / 2;
double centerlon = (minlon + maxlon) / 2;
long sourceespgcode = getProjectEPSGCodeByLon_Lat(centerlon, centerlat);
long demespgcode = GetEPSGFromRasterFile(indemPath);
double grid_resolution = gridx < gridy ? gridx : gridy;
double degreePerPixelX = grid_resolution / 110000.0;
double degreePerPixelY = grid_resolution / 110000.0;
bool meter2degreeflag = ConvertResolutionToDegrees(
sourceespgcode,
grid_resolution,
centerlon,
centerlat,
degreePerPixelX,
degreePerPixelY
);
if (!meter2degreeflag) {
qDebug() << u8"ת»»·Ö±æÂÊΪ¾­Î³¶Èʧ°Ü";
degreePerPixelX = grid_resolution / 110000.0;
degreePerPixelY = grid_resolution / 110000.0;
}
qDebug() << u8"DEMÓ°Ïñ·¶Î§:";
qDebug() << u8"ÊäÈë·Ö±æÂÊ£º"<<gridx<<" \t" << gridy;
qDebug() << u8"Ó°Ïñ·Ö±æÂÊ:\t" << grid_resolution;
qDebug() << u8"·Ö±æÂÊת»»Îª¾­Î³¶È:\t" << degreePerPixelX << "\t" << degreePerPixelY;
// ¼ÆËãDEMµÄ·Ö±æÂÊ
// ²Ã¼ôDEM
QString filename = getFileNameWidthoutExtend(l1dataset.getxmlFilePath());
filename = filename.right(5);
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(), degreePerPixelX, degreePerPixelY);
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() << u8"²éÕÒ±íÉú³É´íÎó£º\t" + getFileNameWidthoutExtend(inxmlPath);
return ErrorCode::FAIL;
}
return ErrorCode::SUCCESS;
}
Reference in New Issue View Git Blame Copy Permalink