RasterProcessTool/LAMPSARProcessProgram/ToolBox/SimulationSAR/RFPCProcessCls.cpp

932 lines
36 KiB
C++

#include "stdafx.h"
#include "RFPCProcessCls.h"
#include "BaseConstVariable.h"
#include "SARSatelliteSimulationAbstractCls.h"
#include "SARSimulationTaskSetting.h"
#include "SatelliteOribtModel.h"
#include <QDebug>
#include "ImageOperatorBase.h"
#include "GeoOperator.h"
#include "EchoDataFormat.h"
#include <QDir>
#include <QDatetime>
#include <omp.h>
#include <QProgressDialog>
#include <QMessageBox>
#ifdef DEBUGSHOWDIALOG
#include "ImageShowDialogClass.h"
#endif
#ifdef __CUDANVCC___
#include "GPUTool.cuh"
#include "GPURFPC.cuh"
#include <cuda_runtime.h>
#include <cublas_v2.h>
#endif // __CUDANVCC___
//#include <Imageshow/ImageShowDialogClass.h>
CUDA_AntSate_PtrList* malloc_AntSate_PtrList(long PRFCount)
{
CUDA_AntSate_PtrList* antlist = (CUDA_AntSate_PtrList*)malloc(sizeof(CUDA_AntSate_PtrList));
antlist->h_antpx = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antpy = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antpz = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antvx = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antvy = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antvz = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antdirectx = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antdirecty = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antdirectz = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antXaxisX = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antXaxisY = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antXaxisZ = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antYaxisX = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antYaxisY = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antYaxisZ = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antZaxisX = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antZaxisY = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->h_antZaxisZ = (double*)mallocCUDAHost(sizeof(double) * PRFCount);
antlist->d_antpx = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antpy = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antpz = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antvx = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antvy = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antvz = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antdirectx = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antdirecty = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antdirectz = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antXaxisX = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antXaxisY = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antXaxisZ = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antYaxisX = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antYaxisY = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antYaxisZ = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antZaxisX = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antZaxisY = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->d_antZaxisZ = (double*)mallocCUDADevice(sizeof(double) * PRFCount);
antlist->PRF_len = PRFCount;
return antlist;
}
void Free_AntSate_PtrList(CUDA_AntSate_PtrList* antlist)
{
FreeCUDAHost(antlist->h_antpx);
FreeCUDAHost(antlist->h_antpy);
FreeCUDAHost(antlist->h_antpz);
FreeCUDAHost(antlist->h_antvx);
FreeCUDAHost(antlist->h_antvy);
FreeCUDAHost(antlist->h_antvz);
FreeCUDAHost(antlist->h_antdirectx);
FreeCUDAHost(antlist->h_antdirecty);
FreeCUDAHost(antlist->h_antdirectz);
FreeCUDAHost(antlist->h_antXaxisX);
FreeCUDAHost(antlist->h_antXaxisY);
FreeCUDAHost(antlist->h_antXaxisZ);
FreeCUDAHost(antlist->h_antYaxisX);
FreeCUDAHost(antlist->h_antYaxisY);
FreeCUDAHost(antlist->h_antYaxisZ);
FreeCUDAHost(antlist->h_antZaxisX);
FreeCUDAHost(antlist->h_antZaxisY);
FreeCUDAHost(antlist->h_antZaxisZ);
FreeCUDADevice(antlist->d_antpx);
FreeCUDADevice(antlist->d_antpy);
FreeCUDADevice(antlist->d_antpz);
FreeCUDADevice(antlist->d_antvx);
FreeCUDADevice(antlist->d_antvy);
FreeCUDADevice(antlist->d_antvz);
FreeCUDADevice(antlist->d_antdirectx);
FreeCUDADevice(antlist->d_antdirecty);
FreeCUDADevice(antlist->d_antdirectz);
FreeCUDADevice(antlist->d_antXaxisX);
FreeCUDADevice(antlist->d_antXaxisY);
FreeCUDADevice(antlist->d_antXaxisZ);
FreeCUDADevice(antlist->d_antYaxisX);
FreeCUDADevice(antlist->d_antYaxisY);
FreeCUDADevice(antlist->d_antYaxisZ);
FreeCUDADevice(antlist->d_antZaxisX);
FreeCUDADevice(antlist->d_antZaxisY);
FreeCUDADevice(antlist->d_antZaxisZ);
antlist->h_antpx = nullptr;
antlist->h_antpy = nullptr;
antlist->h_antpz = nullptr;
antlist->h_antvx = nullptr;
antlist->h_antvy = nullptr;
antlist->h_antvz = nullptr;
antlist->h_antdirectx = nullptr;
antlist->h_antdirecty = nullptr;
antlist->h_antdirectz = nullptr;
antlist->h_antXaxisX = nullptr;
antlist->h_antXaxisY = nullptr;
antlist->h_antXaxisZ = nullptr;
antlist->h_antYaxisX = nullptr;
antlist->h_antYaxisY = nullptr;
antlist->h_antYaxisZ = nullptr;
antlist->h_antZaxisX = nullptr;
antlist->h_antZaxisY = nullptr;
antlist->h_antZaxisZ = nullptr;
antlist->d_antpx = nullptr;
antlist->d_antpy = nullptr;
antlist->d_antpz = nullptr;
antlist->d_antvx = nullptr;
antlist->d_antvy = nullptr;
antlist->d_antvz = nullptr;
antlist->d_antdirectx = nullptr;
antlist->d_antdirecty = nullptr;
antlist->d_antdirectz = nullptr;
antlist->d_antXaxisX = nullptr;
antlist->d_antXaxisY = nullptr;
antlist->d_antXaxisZ = nullptr;
antlist->d_antYaxisX = nullptr;
antlist->d_antYaxisY = nullptr;
antlist->d_antYaxisZ = nullptr;
antlist->d_antZaxisX = nullptr;
antlist->d_antZaxisY = nullptr;
antlist->d_antZaxisZ = nullptr;
free(antlist);
antlist = nullptr;
}
void COPY_AntStation_FROM_HOST_GPU(std::shared_ptr<SatelliteOribtNode[]> sateOirbtNodes,
std::shared_ptr<CUDA_AntSate_PtrList> gpupptr,
long startPID,
long PRF_len)
{
assert(gpupptr->PRF_len <= PRF_len);
long prfid = 0;
for (long tempprfid = 0; tempprfid < PRF_len; tempprfid++) {
prfid = tempprfid + startPID;
gpupptr->h_antpx[tempprfid] = sateOirbtNodes[prfid].Px;
gpupptr->h_antpy[tempprfid] = sateOirbtNodes[prfid].Py;
gpupptr->h_antpz[tempprfid] = sateOirbtNodes[prfid].Pz;
gpupptr->h_antvx[tempprfid] = sateOirbtNodes[prfid].Vx;
gpupptr->h_antvy[tempprfid] = sateOirbtNodes[prfid].Vy;
gpupptr->h_antvz[tempprfid] = sateOirbtNodes[prfid].Vz; //6
gpupptr->h_antdirectx[tempprfid] = sateOirbtNodes[prfid].AntDirecX;
gpupptr->h_antdirecty[tempprfid] = sateOirbtNodes[prfid].AntDirecY;
gpupptr->h_antdirectz[tempprfid] = sateOirbtNodes[prfid].AntDirecZ;
gpupptr->h_antXaxisX[tempprfid] = sateOirbtNodes[prfid].AntXaxisX;
gpupptr->h_antXaxisY[tempprfid] = sateOirbtNodes[prfid].AntXaxisY;
gpupptr->h_antXaxisZ[tempprfid] = sateOirbtNodes[prfid].AntXaxisZ;//12
gpupptr->h_antYaxisX[tempprfid] = sateOirbtNodes[prfid].AntYaxisX;
gpupptr->h_antYaxisY[tempprfid] = sateOirbtNodes[prfid].AntYaxisY;
gpupptr->h_antYaxisZ[tempprfid] = sateOirbtNodes[prfid].AntYaxisZ;//15
gpupptr->h_antZaxisX[tempprfid] = sateOirbtNodes[prfid].AntZaxisX;
gpupptr->h_antZaxisY[tempprfid] = sateOirbtNodes[prfid].AntZaxisY;
gpupptr->h_antZaxisZ[tempprfid] = sateOirbtNodes[prfid].AntZaxisZ;//18
}
HostToDevice(gpupptr->h_antpx, gpupptr->d_antpx, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antpy, gpupptr->d_antpy, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antpz, gpupptr->d_antpz, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antvx, gpupptr->d_antvx, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antvy, gpupptr->d_antvy, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antvz, gpupptr->d_antvz, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antdirectx, gpupptr->d_antdirectx, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antdirecty, gpupptr->d_antdirecty, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antdirectz, gpupptr->d_antdirectz, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antXaxisX, gpupptr->d_antXaxisX, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antXaxisY, gpupptr->d_antXaxisY, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antXaxisZ, gpupptr->d_antXaxisZ, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antYaxisX, gpupptr->d_antYaxisX, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antYaxisY, gpupptr->d_antYaxisY, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antYaxisZ, gpupptr->d_antYaxisZ, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antZaxisX, gpupptr->d_antZaxisX, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antZaxisY, gpupptr->d_antZaxisY, sizeof(double) * PRF_len);
HostToDevice(gpupptr->h_antZaxisZ, gpupptr->d_antZaxisZ, sizeof(double) * PRF_len);
}
RFPCProcessCls::RFPCProcessCls()
{
this->PluseCount = 0;
this->PlusePoint = 0;
this->TaskSetting = nullptr;
this->EchoSimulationData = nullptr;
this->LandCoverPath = "";
this->OutEchoPath = "";
this->LandCoverPath.clear();
this->OutEchoPath.clear();
this->SigmaDatabasePtr = std::shared_ptr<SigmaDatabase>(new SigmaDatabase);
}
RFPCProcessCls::~RFPCProcessCls()
{
}
void RFPCProcessCls::setTaskSetting(std::shared_ptr < AbstractSARSatelliteModel> TaskSetting)
{
this->TaskSetting = std::shared_ptr < AbstractSARSatelliteModel>(TaskSetting);
qDebug() << "RFPCProcessCls::setTaskSetting";
}
void RFPCProcessCls::setEchoSimulationDataSetting(std::shared_ptr<EchoL0Dataset> EchoSimulationData)
{
this->EchoSimulationData = std::shared_ptr<EchoL0Dataset>(EchoSimulationData);
qDebug() << "RFPCProcessCls::setEchoSimulationDataSetting";
}
void RFPCProcessCls::setTaskFileName(QString EchoFileName)
{
this->TaskFileName = EchoFileName;
}
void RFPCProcessCls::setDEMTiffPath(QString DEMTiffPath)
{
this->demxyzPath = DEMTiffPath;
}
void RFPCProcessCls::setSloperPath(QString InSloperPath)
{
this->demsloperPath = InSloperPath;
}
void RFPCProcessCls::setLandCoverPath(QString LandCoverPath)
{
this->LandCoverPath = LandCoverPath;
}
void RFPCProcessCls::setOutEchoPath(QString OutEchoPath)
{
this->OutEchoPath = OutEchoPath;
}
ErrorCode RFPCProcessCls::Process(long num_thread)
{
// RFPC 算法
qDebug() << u8"params init ....";
ErrorCode stateCode = this->InitParams();
if (stateCode != ErrorCode::SUCCESS) {
return stateCode;
}
else {}
qDebug() << "RFPCMainProcess";
stateCode = this->InitEchoMaskArray();
if (stateCode != ErrorCode::SUCCESS) {
return stateCode;
}
else {}
qDebug() << "InitEchoMaskArray";
//stateCode = this->RFPCMainProcess(num_thread);
// 初始化回波
this->EchoSimulationData->initEchoArr(std::complex<double>(0, 0));
stateCode = this->RFPCMainProcess_GPU();
if (stateCode != ErrorCode::SUCCESS) {
return stateCode;
}
else {}
return ErrorCode::SUCCESS;
}
ErrorCode RFPCProcessCls::InitParams()
{
if (nullptr == this->TaskSetting || this->demxyzPath.isEmpty() ||
this->LandCoverPath.isEmpty() || this->demsloperPath.isEmpty()) {
return ErrorCode::RFPC_PARAMSISEMPTY;
}
else {
}
// 归一化绝对路径
this->OutEchoPath = QDir(this->OutEchoPath).absolutePath();
// 回波大小
double imgStart_end = this->TaskSetting->getSARImageEndTime() - this->TaskSetting->getSARImageStartTime();
this->PluseCount = ceil(imgStart_end * this->TaskSetting->getPRF());
double rangeTimeSample = (this->TaskSetting->getFarRange() - this->TaskSetting->getNearRange()) * 2.0 / LIGHTSPEED;
this->PlusePoint = ceil(rangeTimeSample * this->TaskSetting->getFs());
// 初始化回波存放位置
qDebug() << "--------------Echo Data Setting ---------------------------------------";
this->EchoSimulationData = std::shared_ptr<EchoL0Dataset>(new EchoL0Dataset);
this->EchoSimulationData->setCenterFreq(this->TaskSetting->getCenterFreq());
this->EchoSimulationData->setNearRange(this->TaskSetting->getNearRange());
this->EchoSimulationData->setFarRange(this->TaskSetting->getFarRange());
this->EchoSimulationData->setFs(this->TaskSetting->getFs());
this->EchoSimulationData->setBandwidth(this->TaskSetting->getBandWidth());
this->EchoSimulationData->setCenterAngle(this->TaskSetting->getCenterLookAngle());
this->EchoSimulationData->setLookSide(this->TaskSetting->getIsRightLook() ? "R" : "L");
this->EchoSimulationData->OpenOrNew(OutEchoPath, TaskFileName, PluseCount, PlusePoint);
QString tmpfolderPath = QDir(OutEchoPath).filePath("tmp");
if (QDir(tmpfolderPath).exists() == false) {
QDir(OutEchoPath).mkpath(tmpfolderPath);
}
this->tmpfolderPath = tmpfolderPath;
return ErrorCode::SUCCESS;
}
ErrorCode RFPCProcessCls::InitEchoMaskArray()
{
QString name = this->EchoSimulationData->getSimulationTaskName();
this->OutEchoMaskPath = JoinPath(this->OutEchoPath, name + "_echomask.bin");
Eigen::MatrixXd gt(2, 3);
gt(0, 0) = 0;
gt(0, 1) = 1;
gt(0, 2) = 0;
gt(1, 0) = 0;
gt(1, 1) = 0;
gt(1, 2) = 1;
gdalImage echomaskImg = CreategdalImage(this->OutEchoMaskPath,
this->EchoSimulationData->getPluseCount(),
this->EchoSimulationData->getPlusePoints(),
1,
gt, "",
false, true, true);
long cols = this->EchoSimulationData->getPlusePoints();
long rows = this->EchoSimulationData->getPluseCount();
long blocksize = Memory1GB / 8 / this->EchoSimulationData->getPlusePoints() * 4;
for (long startid = 0; startid < rows; startid = startid + blocksize) {
Eigen::MatrixXd data = echomaskImg.getData(startid, 0, blocksize, cols, 1);
data = data.array() * 0;
echomaskImg.saveImage(data, startid, 0, 1);
}
return ErrorCode::SUCCESS;
}
std::shared_ptr<SatelliteOribtNode[]> RFPCProcessCls::getSatelliteOribtNodes(double prf_time, double dt, bool antflag, long double imageStarttime)
{
std::shared_ptr<SatelliteOribtNode[]> sateOirbtNodes(new SatelliteOribtNode[this->PluseCount], delArrPtr);
{ // 姿态计算不同
qDebug() << "Ant position finished started !!!";
// 计算姿态
std::shared_ptr<double> antpos = this->EchoSimulationData->getAntPos();
double dAt = 1e-6;
double prf_time_dt = 0;
Landpoint InP{ 0,0,0 }, outP{ 0,0,0 };
for (long prf_id = 0; prf_id < this->PluseCount; prf_id++) {
prf_time = dt * prf_id;
prf_time_dt = prf_time + dAt;
SatelliteOribtNode sateOirbtNode;
SatelliteOribtNode sateOirbtNode_dAt;
this->TaskSetting->getSatelliteOribtNode(prf_time, sateOirbtNode, antflag);
this->TaskSetting->getSatelliteOribtNode(prf_time_dt, sateOirbtNode_dAt, antflag);
sateOirbtNode.AVx = (sateOirbtNode_dAt.Vx - sateOirbtNode.Vx) / dAt; // 加速度
sateOirbtNode.AVy = (sateOirbtNode_dAt.Vy - sateOirbtNode.Vy) / dAt;
sateOirbtNode.AVz = (sateOirbtNode_dAt.Vz - sateOirbtNode.Vz) / dAt;
InP.lon = sateOirbtNode.Px;
InP.lat = sateOirbtNode.Py;
InP.ati = sateOirbtNode.Pz;
outP = XYZ2LLA(InP);
antpos.get()[prf_id * 19 + 0] = prf_time + imageStarttime;
antpos.get()[prf_id * 19 + 1] = sateOirbtNode.Px;
antpos.get()[prf_id * 19 + 2] = sateOirbtNode.Py;
antpos.get()[prf_id * 19 + 3] = sateOirbtNode.Pz;
antpos.get()[prf_id * 19 + 4] = sateOirbtNode.Vx;
antpos.get()[prf_id * 19 + 5] = sateOirbtNode.Vy;
antpos.get()[prf_id * 19 + 6] = sateOirbtNode.Vz;
antpos.get()[prf_id * 19 + 7] = sateOirbtNode.AntDirecX;
antpos.get()[prf_id * 19 + 8] = sateOirbtNode.AntDirecY;
antpos.get()[prf_id * 19 + 9] = sateOirbtNode.AntDirecZ;
antpos.get()[prf_id * 19 + 10] = sateOirbtNode.AVx;
antpos.get()[prf_id * 19 + 11] = sateOirbtNode.AVy;
antpos.get()[prf_id * 19 + 12] = sateOirbtNode.AVz;
antpos.get()[prf_id * 19 + 13] = sateOirbtNode.zeroDopplerDirectX;
antpos.get()[prf_id * 19 + 14] = sateOirbtNode.zeroDopplerDirectY;
antpos.get()[prf_id * 19 + 15] = sateOirbtNode.zeroDopplerDirectZ;
antpos.get()[prf_id * 19 + 16] = outP.lon;
antpos.get()[prf_id * 19 + 17] = outP.lat;
antpos.get()[prf_id * 19 + 18] = outP.ati;
sateOirbtNodes[prf_id] = sateOirbtNode;
}
this->EchoSimulationData->saveAntPos(antpos);
antpos.reset();
qDebug() << "Ant position finished sucessfully !!!";
}
return sateOirbtNodes;
}
void RFPCProcessMain(long num_thread,
QString TansformPatternFilePath, QString ReceivePatternFilePath,
QString simulationtaskName, QString OutEchoPath,
QString GPSXmlPath, QString TaskXmlPath,QString demTiffPath, QString sloperPath, QString LandCoverPath)
{
std::shared_ptr < AbstractSARSatelliteModel> task = ReadSimulationSettingsXML(TaskXmlPath);
if (nullptr == task)
{
return;
}
else {
// 打印参数
qDebug() << "--------------Task Seting ---------------------------------------";
qDebug() << "SARImageStartTime: " << task->getSARImageStartTime();
qDebug() << "SARImageEndTime: " << task->getSARImageEndTime();
qDebug() << "BandWidth: " << task->getBandWidth();
qDebug() << "CenterFreq: " << task->getCenterFreq();
qDebug() << "PRF: " << task->getPRF();
qDebug() << "Lamda: " << task->getCenterLamda();
qDebug() << "Fs: " << task->getFs();
qDebug() << "POLAR: " << task->getPolarType();
qDebug() << "NearRange: " << task->getNearRange();
qDebug() << "FarRange: " << task->getFarRange();
qDebug() << (task->getFarRange() - task->getNearRange()) * 2 / LIGHTSPEED * task->getFs();
qDebug() << "\n\n";
}
// 1.2 设置天线方向图
std::vector<RadiationPatternGainPoint> TansformPatternGainpoints = ReadGainFile(TansformPatternFilePath);
std::shared_ptr<AbstractRadiationPattern> TansformPatternGainPtr = CreateAbstractRadiationPattern(TansformPatternGainpoints);
std::vector<RadiationPatternGainPoint> ReceivePatternGainpoints = ReadGainFile(ReceivePatternFilePath);
std::shared_ptr<AbstractRadiationPattern> ReceivePatternGainPtr = CreateAbstractRadiationPattern(ReceivePatternGainpoints);
task->setTransformRadiationPattern(TansformPatternGainPtr);
task->setReceiveRadiationPattern(ReceivePatternGainPtr);
//2. 读取GPS节点
std::vector<SatelliteOribtNode> nodes;
ErrorCode stateCode = ReadSateGPSPointsXML(GPSXmlPath, nodes);
if (stateCode != ErrorCode::SUCCESS)
{
qWarning() << QString::fromStdString(errorCode2errInfo(stateCode));
return;
}
else {}
std::shared_ptr<AbstractSatelliteOribtModel> SatelliteOribtModel = CreataPolyfitSatelliteOribtModel(nodes, task->getSARImageStartTime(), 3); // 以成像开始时间作为 时间参考起点
SatelliteOribtModel->setbeamAngle(task->getCenterLookAngle(), task->getIsRightLook()); // 设置天线方向图
if (nullptr == SatelliteOribtModel)
{
return;
}
else {
task->setSatelliteOribtModel(SatelliteOribtModel);
}
qDebug() << "-------------- RFPC init ---------------------------------------";
RFPCProcessCls RFPC;
RFPC.setTaskSetting(task); //qDebug() << "setTaskSetting";
RFPC.setTaskFileName(simulationtaskName); //qDebug() << "setTaskFileName";
RFPC.setDEMTiffPath(demTiffPath); //qDebug() << "setDEMTiffPath";
RFPC.setSloperPath(sloperPath); //qDebug() << "setDEMTiffPath";
RFPC.setLandCoverPath(LandCoverPath); //qDebug() << "setLandCoverPath";
RFPC.setOutEchoPath(OutEchoPath); //qDebug() << "setOutEchoPath";
qDebug() << "-------------- RFPC start---------------------------------------";
RFPC.Process(num_thread); // 处理程序
qDebug() << "-------------- RFPC end---------------------------------------";
}
ErrorCode RFPCProcessCls::RFPCMainProcess_GPU() {
/** 内存分配***************************************************/
long TargetMemoryMB = 500;
/** 参数区域***************************************************/
QVector<double> freqlist = this->TaskSetting->getFreqList();
long freqnum = freqlist.count();
float f0 = float(freqlist[0] / 1e9);
float dfreq = float((freqlist[1] - freqlist[0]) / 1e9);
#if (defined __PRFDEBUG__) && (defined __PRFDEBUG_PRFINF__)
double* h_freqPtr = (double*)mallocCUDAHost(sizeof(double) * freqnum);
for (long fid = 0; fid < freqnum; fid++) {
h_freqPtr[fid] = (f0 + dfreq * fid) * 1e9;
}
testOutAmpArr("freqlist.bin", h_freqPtr, freqnum, 1);
#endif
long PRFCount = this->EchoSimulationData->getPluseCount();
double NearRange = this->EchoSimulationData->getNearRange(); // 近斜距
double FarRange = this->EchoSimulationData->getFarRange();
double Pt = this->TaskSetting->getPt() * this->TaskSetting->getGri();// 发射电压 1v
double lamda = this->TaskSetting->getCenterLamda(); // 波长
double refphaseRange = this->TaskSetting->getRefphaseRange(); // 参考相位斜距
double prf_time = 0;
double dt = 1 / this->TaskSetting->getPRF();// 获取每次脉冲的时间间隔
bool antflag = true; // 计算天线方向图
long double imageStarttime = this->TaskSetting->getSARImageStartTime();
this->EchoSimulationData->getAntPos();
std::shared_ptr<SatelliteOribtNode[]> sateOirbtNodes = this->getSatelliteOribtNodes(prf_time, dt, antflag, imageStarttime);
/** 天线方向图***************************************************/
std::shared_ptr<AbstractRadiationPattern> TransformPattern = this->TaskSetting->getTransformRadiationPattern(); // 发射天线方向图
std::shared_ptr<AbstractRadiationPattern> ReceivePattern = this->TaskSetting->getReceiveRadiationPattern(); // 接收天线方向图
POLARTYPEENUM polartype = this->TaskSetting->getPolarType();
PatternImageDesc TantPatternDesc = {};
double* h_TantPattern = nullptr;
double* d_TantPattern = nullptr;
double maxTransAntPatternValue = 0;
{
// 处理发射天线方向图
double Tminphi = TransformPattern->getMinPhi();
double Tmaxphi = TransformPattern->getMaxPhi();
double Tmintheta = TransformPattern->getMinTheta();
double Tmaxtheta = TransformPattern->getMaxTheta();
long Tphinum = TransformPattern->getPhis().size();
long Tthetanum = TransformPattern->getThetas().size();
double TstartTheta = Tmintheta;
double TstartPhi = Tminphi;
double Tdtheta = (Tmaxtheta - Tmintheta) / (Tthetanum - 1);
double Tdphi = (Tmaxphi - Tminphi) / (Tphinum - 1);
h_TantPattern = (double*)mallocCUDAHost(sizeof(double) * Tthetanum * Tphinum);
d_TantPattern = (double*)mallocCUDADevice(sizeof(double) * Tthetanum * Tphinum);
for (long i = 0; i < Tthetanum; i++) {
for (long j = Tphinum - 1; j >= 0; j--) {
//h_TantPattern[i * Tphinum + j] = TransformPattern->getGainLearThetaPhi(TstartTheta + i * Tdtheta, TstartPhi + j * Tdphi);
h_TantPattern[i * Tphinum + j] = TransformPattern->getGain(TstartTheta + i * Tdtheta, TstartPhi + j * Tdphi);
}
}
testOutAntPatternTrans("TransPattern.bin", h_TantPattern, TstartTheta, Tdtheta, TstartPhi, Tdphi, Tthetanum, Tphinum);
maxTransAntPatternValue = powf(10.0, h_TantPattern[0] / 10);
for (long i = 0; i < Tthetanum; i++) {
for (long j = 0; j < Tphinum; j++) {
h_TantPattern[i * Tphinum + j] = powf(10.0, h_TantPattern[i * Tphinum + j] / 10); // 转换为线性值
if (maxTransAntPatternValue < h_TantPattern[i * Tphinum + j]) {
maxTransAntPatternValue = h_TantPattern[i * Tphinum + j];
}
}
}
HostToDevice(h_TantPattern, d_TantPattern, sizeof(double) * Tthetanum * Tphinum);
TantPatternDesc.startTheta = TstartTheta;
TantPatternDesc.startPhi = TstartPhi;
TantPatternDesc.dtheta = Tdtheta;
TantPatternDesc.dphi = Tdphi;
TantPatternDesc.phinum = Tphinum;
TantPatternDesc.thetanum = Tthetanum;
}
PatternImageDesc RantPatternDesc = {};
double* h_RantPattern = nullptr;
double* d_RantPattern = nullptr;
double maxReceiveAntPatternValue = 0;
{
// 处理接收天线方向图
double Rminphi = ReceivePattern->getMinPhi();
double Rmaxphi = ReceivePattern->getMaxPhi();
double Rmintheta = ReceivePattern->getMinTheta();
double Rmaxtheta = ReceivePattern->getMaxTheta();
long Rphinum = ReceivePattern->getPhis().size();
long Rthetanum = ReceivePattern->getThetas().size();
double RstartTheta = Rmintheta;
double RstartPhi = Rminphi;
double Rdtheta = (Rmaxtheta - Rmintheta) / (Rthetanum - 1);
double Rdphi = (Rmaxphi - Rminphi) / (Rphinum - 1);
h_RantPattern = (double*)mallocCUDAHost(sizeof(double) * Rthetanum * Rphinum);
d_RantPattern = (double*)mallocCUDADevice(sizeof(double) * Rthetanum * Rphinum);
for (long i = 0; i < Rthetanum; i++) {
for (long j = 0; j < Rphinum; j++) {
//h_RantPattern[i * Rphinum + j] = ReceivePattern->getGainLearThetaPhi(RstartTheta + i * Rdtheta, RstartPhi + j * Rdphi);
h_RantPattern[i * Rphinum + j] = ReceivePattern->getGain(RstartTheta + i * Rdtheta, RstartPhi + j * Rdphi);
}
}
testOutAntPatternTrans("ReceivePattern.bin", h_RantPattern, Rmintheta, Rdtheta, RstartPhi, Rdphi, Rthetanum, Rphinum);
maxReceiveAntPatternValue = powf(10.0, h_RantPattern[0] / 10);
for (long i = 0; i < Rthetanum; i++) {
for (long j = 0; j < Rphinum; j++) {
h_RantPattern[i * Rphinum + j] = powf(10.0, h_RantPattern[i * Rphinum + j] / 10);
if (maxReceiveAntPatternValue < h_RantPattern[i * Rphinum + j]) {
maxReceiveAntPatternValue = h_RantPattern[i * Rphinum + j];
}
}
}
HostToDevice(h_RantPattern, d_RantPattern, sizeof(double) * Rthetanum * Rphinum);
RantPatternDesc.startTheta = RstartTheta;
RantPatternDesc.startPhi = RstartPhi;
RantPatternDesc.dtheta = Rdtheta;
RantPatternDesc.dphi = Rdphi;
RantPatternDesc.phinum = Rphinum;
RantPatternDesc.thetanum = Rthetanum;
}
/** 坐标区域点***************************************************/
gdalImage demxyz(this->demxyzPath);// 地面点坐标
gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型
gdalImage demsloperxyz(this->demsloperPath);// 地面坡向
long demRow = demxyz.height;
long demCol = demxyz.width;
//处理地表覆盖
QMap<long, long> clamap;
long clamapid = 0;
long startline = 0;
{
long blokline = getBlockRows(2e4, demCol, sizeof(double),demRow);
for (startline = 0; startline < demRow; startline = startline + blokline) {
Eigen::MatrixXd clsland = demlandcls.getData(startline, 0, blokline, demlandcls.width, 1);
long clsrows = clsland.rows();
long clscols = clsland.cols();
long clsid = 0;
for (long ii = 0; ii < clsrows; ii++) {
for (long jj = 0; jj < clscols; jj++) {
clsid = clsland(ii, jj);
if (clamap.contains(clsid)) {}
else {
clamap.insert(clsid, clamapid);
clamapid = clamapid + 1;
}
}
}
}
std::cout << "class id recoding" << std::endl;
for (long id : clamap.keys()) {
std::cout << id << " -> " << clamap[id] << std::endl;
}
}
CUDASigmaParam* h_clsSigmaParam = (CUDASigmaParam*)mallocCUDAHost(sizeof(CUDASigmaParam) * clamapid);
CUDASigmaParam* d_clsSigmaParam = (CUDASigmaParam*)mallocCUDADevice(sizeof(CUDASigmaParam) * clamapid);
{
std::map<long, SigmaParam> tempSigmaParam = this->SigmaDatabasePtr->getsigmaParams(polartype);
for (long id : clamap.keys()) {
SigmaParam tempp = tempSigmaParam[id];
h_clsSigmaParam[clamap[id]].p1 = tempp.p1;
h_clsSigmaParam[clamap[id]].p2 = tempp.p2;
h_clsSigmaParam[clamap[id]].p3 = tempp.p3;
h_clsSigmaParam[clamap[id]].p4 = tempp.p4;
h_clsSigmaParam[clamap[id]].p5 = tempp.p5;
h_clsSigmaParam[clamap[id]].p6 = tempp.p6;
}
// 打印日志
std::cout << "sigma params:" << std::endl;
std::cout << "classid:\tp1\tp2\tp3\tp4\tp5\tp6" << std::endl;
for (long ii = 0; ii < clamapid; ii++) {
std::cout << ii << ":\t" << h_clsSigmaParam[ii].p1;
std::cout << "\t" << h_clsSigmaParam[ii].p2;
std::cout << "\t" << h_clsSigmaParam[ii].p3;
std::cout << "\t" << h_clsSigmaParam[ii].p4;
std::cout << "\t" << h_clsSigmaParam[ii].p5;
std::cout << "\t" << h_clsSigmaParam[ii].p6 << std::endl;
}
std::cout << "";
}
HostToDevice(h_clsSigmaParam, d_clsSigmaParam, sizeof(CUDASigmaParam) * clamapid);
qDebug() << "CUDA class Proces finished!!!";
// 处理地面坐标
long blockline = getBlockRows(TargetMemoryMB, demCol, sizeof(double), demRow);
double* h_dem_x = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
double* h_dem_y = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
double* h_dem_z = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
double* h_demsloper_x = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
double* h_demsloper_y = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
double* h_demsloper_z = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
long* h_demcls = (long*)mallocCUDAHost(sizeof(long) * blockline * demCol);
double* d_dem_x = (double*)mallocCUDADevice(sizeof(double) * blockline * demCol);
double* d_dem_y = (double*)mallocCUDADevice(sizeof(double) * blockline * demCol);
double* d_dem_z = (double*)mallocCUDADevice(sizeof(double) * blockline * demCol);
double* d_demsloper_x = (double*)mallocCUDADevice(sizeof(double) * blockline * demCol);
double* d_demsloper_y = (double*)mallocCUDADevice(sizeof(double) * blockline * demCol);
double* d_demsloper_z = (double*)mallocCUDADevice(sizeof(double) * blockline * demCol);
long* d_demcls = (long*) mallocCUDADevice(sizeof(long) * blockline * demCol);
/** 处理回波***************************************************/
long echo_block_rows = getBlockRows(1000, freqnum, sizeof(float)*2, PRFCount);
float* h_echo_block_real = (float*)mallocCUDAHost(sizeof(float) * echo_block_rows * freqnum);
float* h_echo_block_imag = (float*)mallocCUDAHost(sizeof(float) * echo_block_rows * freqnum);
float* d_echo_block_real = (float*)mallocCUDADevice(sizeof(float) * echo_block_rows * freqnum);
float* d_echo_block_imag = (float*)mallocCUDADevice(sizeof(float) * echo_block_rows * freqnum);
float* d_temp_R = (float*)mallocCUDADevice(sizeof(float) * echo_block_rows * SHAREMEMORY_FLOAT_HALF); //2GB 距离
float* d_temp_amp = (float*)mallocCUDADevice(sizeof(float) * echo_block_rows * SHAREMEMORY_FLOAT_HALF);//2GB 强度
/** 主流程处理 ***************************************************/
qDebug() << "CUDA Main Proces";
for (long sprfid = 0; sprfid < PRFCount; sprfid = sprfid + echo_block_rows) {
long PRF_len = (sprfid + echo_block_rows) < PRFCount ? echo_block_rows : (PRFCount - sprfid);
qDebug() << "Start PRF: " << sprfid << "\t-\t" << sprfid + PRF_len << "\t:copy ant list host -> GPU";
std::shared_ptr< CUDA_AntSate_PtrList> antptrlist(malloc_AntSate_PtrList(PRF_len), Free_AntSate_PtrList);
COPY_AntStation_FROM_HOST_GPU(sateOirbtNodes, antptrlist, sprfid, PRF_len);
qDebug() << "Start PRF: " << sprfid << "\t-\t" << sprfid + PRF_len << "\t:copy echo data list host -> GPU";
std::shared_ptr<std::complex<double>> echo_temp = this->EchoSimulationData->getEchoArr(sprfid, PRF_len);
for (long ii = 0; ii < PRF_len; ii++) {
for (long jj = 0; jj < freqnum; jj++) {
h_echo_block_real[ii * freqnum + jj]=echo_temp.get()[ii * freqnum + jj].real();
h_echo_block_imag[ii * freqnum + jj]=echo_temp.get()[ii * freqnum + jj].imag();
}
}
HostToDevice(h_echo_block_real, d_echo_block_real, sizeof(float) * PRF_len* freqnum);
HostToDevice(h_echo_block_imag, d_echo_block_imag, sizeof(float) * PRF_len* freqnum);
for (startline = 0; startline < demRow; startline = startline + blockline) {
Eigen::MatrixXd dem_x = demxyz.getData(startline, 0, blockline, demCol, 1); // 地面坐标
Eigen::MatrixXd dem_y = demxyz.getData(startline, 0, blockline, demCol, 2);
Eigen::MatrixXd dem_z = demxyz.getData(startline, 0, blockline, demCol, 3);
Eigen::MatrixXd demsloper_x = demsloperxyz.getData(startline, 0, blockline, demCol, 1);
Eigen::MatrixXd demsloper_y = demsloperxyz.getData(startline, 0, blockline, demCol, 2);
Eigen::MatrixXd demsloper_z = demsloperxyz.getData(startline, 0, blockline, demCol, 3);
Eigen::MatrixXd landcover = demlandcls.getData(startline, 0, blockline, demCol, 1);
long temp_dem_row = dem_x.rows();
long temp_dem_col = dem_x.cols();
long temp_dem_count = dem_x.count();
// 更新数据格式
for (long i = 0; i < temp_dem_row; i++) {
for (long j = 0; j < temp_dem_col; j++) {
h_dem_x[i * temp_dem_col + j] = double(dem_x(i, j));
h_dem_y[i * temp_dem_col + j] = double(dem_y(i, j));
h_dem_z[i * temp_dem_col + j] = double(dem_z(i, j));
h_demsloper_x[i * temp_dem_col + j] = double(demsloper_x(i, j));
h_demsloper_y[i * temp_dem_col + j] = double(demsloper_y(i, j));
h_demsloper_z[i * temp_dem_col + j] = double(demsloper_z(i, j));
h_demcls[i * temp_dem_col + j] = clamap[long(landcover(i, j))];
}
}
qDebug() << "Start PRF: " << sprfid << "\t-\t" << sprfid + PRF_len << "\t:copy target data ("<< startline<<" - "<< startline + blockline << ") host -> GPU";
HostToDevice(h_dem_x, d_dem_x , sizeof(double) * blockline * demCol);
HostToDevice(h_dem_y, d_dem_y , sizeof(double) * blockline * demCol);
HostToDevice(h_dem_z, d_dem_z , sizeof(double) * blockline * demCol);
HostToDevice(h_demsloper_x, d_demsloper_x , sizeof(double) * blockline * demCol);
HostToDevice(h_demsloper_y, d_demsloper_y , sizeof(double) * blockline * demCol);
HostToDevice(h_demsloper_z, d_demsloper_z , sizeof(double) * blockline * demCol);
HostToDevice(h_demcls, d_demcls ,sizeof(long)* blockline* demCol);
// 分块处理
qDebug() << "Start PRF: " << sprfid << "\t-\t" << sprfid + PRF_len << "\t:GPU Computer target data (" << startline << "-" << startline + blockline << ")";
CUDA_RFPC_MainProcess(
antptrlist->d_antpx, antptrlist->d_antpy, antptrlist->d_antpz,
antptrlist->d_antXaxisX, antptrlist->d_antXaxisY, antptrlist->d_antXaxisZ, // 天线坐标系的X轴
antptrlist->d_antYaxisX, antptrlist->d_antYaxisY, antptrlist->d_antYaxisZ,// 天线坐标系的Y轴
antptrlist->d_antZaxisX, antptrlist->d_antZaxisY, antptrlist->d_antZaxisZ,// 天线坐标系的Z轴
antptrlist->d_antdirectx, antptrlist->d_antdirecty, antptrlist->d_antdirectz,// 天线的指向
PRF_len, freqnum,
f0,dfreq,
Pt,
refphaseRange,
// 天线方向图
d_TantPattern,
TantPatternDesc.startTheta, TantPatternDesc.startPhi, TantPatternDesc.dtheta, TantPatternDesc.dphi, TantPatternDesc.thetanum, TantPatternDesc.phinum,
d_RantPattern,
RantPatternDesc.startTheta, RantPatternDesc.startPhi, RantPatternDesc.dtheta, RantPatternDesc.dphi, RantPatternDesc.thetanum, RantPatternDesc.phinum,
maxTransAntPatternValue, maxReceiveAntPatternValue,
NearRange, FarRange, // 近斜据
d_dem_x, d_dem_y, d_dem_z, d_demcls, temp_dem_count, // 地面坐标
d_demsloper_x, d_demsloper_y, d_demsloper_z, // 地表坡度矢量
d_clsSigmaParam, clamapid,
d_echo_block_real, d_echo_block_imag,// 输出回波
d_temp_R, d_temp_amp
);
PRINT("dem : %d ~ %d / %d , echo: %d ~ %d / %d \n", startline, startline+ temp_dem_row, demRow, sprfid, sprfid+ PRF_len, PRFCount);
}
#if (defined __PRFDEBUG__) && (defined __PRFDEBUG_PRFINF__)
float* h_temp_R = (float*)mallocCUDAHost(sizeof(float) * echo_block_rows * SHAREMEMORY_FLOAT_HALF); //2GB 距离
float* h_temp_amp = (float*)mallocCUDAHost(sizeof(float) * echo_block_rows * SHAREMEMORY_FLOAT_HALF);//2GB 强度
DeviceToHost(h_temp_R, d_temp_R, sizeof(float) * echo_block_rows * SHAREMEMORY_FLOAT_HALF);
DeviceToHost(h_temp_amp, d_temp_amp, sizeof(float) * echo_block_rows * SHAREMEMORY_FLOAT_HALF);
testOutAmpArr("temp_R.bin", h_temp_R, echo_block_rows, SHAREMEMORY_FLOAT_HALF);
testOutAmpArr("temp_Amp.bin", h_temp_amp, echo_block_rows, SHAREMEMORY_FLOAT_HALF);
FreeCUDAHost(h_temp_R);
FreeCUDAHost(h_temp_amp);
#endif
DeviceToHost(h_echo_block_real, d_echo_block_real, sizeof(float) * PRF_len * freqnum);
DeviceToHost(h_echo_block_imag, d_echo_block_imag, sizeof(float) * PRF_len * freqnum);
for (long ii = 0; ii < PRF_len; ii++) {
for (long jj = 0; jj < freqnum; jj++) {
echo_temp.get()[ii * freqnum + jj].real(h_echo_block_real[ii * freqnum + jj]);
echo_temp.get()[ii * freqnum + jj].imag(h_echo_block_imag[ii * freqnum + jj]);
}
}
this->EchoSimulationData->saveEchoArr(echo_temp, sprfid, PRF_len);
}
/** 内存释放***************************************************/
FreeCUDAHost(h_TantPattern);
FreeCUDAHost(h_RantPattern);
FreeCUDADevice(d_TantPattern);
FreeCUDADevice(d_RantPattern);
FreeCUDAHost(h_dem_x);
FreeCUDAHost(h_dem_y);
FreeCUDAHost(h_dem_z);
FreeCUDAHost(h_demsloper_x);
FreeCUDAHost(h_demsloper_y);
FreeCUDAHost(h_demsloper_z);
FreeCUDAHost(h_demcls);
FreeCUDAHost(h_echo_block_real);
FreeCUDAHost(h_echo_block_imag);
FreeCUDADevice(d_dem_x);
FreeCUDADevice(d_dem_y);
FreeCUDADevice(d_dem_z);
FreeCUDADevice(d_demsloper_x);
FreeCUDADevice(d_demsloper_y);
FreeCUDADevice(d_demsloper_z);
FreeCUDADevice(d_demcls);
FreeCUDADevice(d_echo_block_real);
FreeCUDADevice(d_echo_block_imag);
FreeCUDADevice(d_temp_R);
FreeCUDADevice(d_temp_amp);
return ErrorCode::SUCCESS;
}