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

更新仿真代码文件 

1. 修复天线坐标问题
2. 修复天线方向图 插值问题
ps: 下个版本调试整体的能量计算问题
pull/1/head
陈增辉 2024-12-20 23:38:34 +08:00
parent 2816f887a7
commit fcd4e5075b
15 changed files with 1831 additions and 819 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -355,6 +355,7 @@ healthchecksdb
# Backup folder for Package Reference Convert tool in Visual Studio 2017
MigrationBackup/
debugdata/
# Ionide (cross platform F# VS Code tools) working folder
.ionide/

BIN
BackScattering_C_SAR_V3_NOMAP.rar
View File

Binary file not shown.

3
GF3ProcessToolbox/QRDOrthProcessClass.cpp
View File

@ -50,9 +50,6 @@ void QRDOrthProcessClass::accept()
}
processdialog->close();
}
void QRDOrthProcessClass::reject()

526
README.md
View File

@ -236,3 +236,529 @@ void RTPC(float* antx, float* anty, float* antz,
ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
{
omp_set_num_threads(num_thread);// 设置openmp 线程数量
double widthSpace = LIGHTSPEED / 2 / this->TaskSetting->getFs();
double prf_time = 0;
double dt = 1 / this->TaskSetting->getPRF();// 获取每次脉冲的时间间隔
bool antflag = true; // 计算天线方向图
Landpoint LandP{ 0,0,0 };
Point3 GERpoint{ 0,0,0 };
double R = 0;
double dem_row = 0, dem_col = 0, dem_alt = 0;
long double imageStarttime = 0;
imageStarttime = this->TaskSetting->getSARImageStartTime();
//std::vector<SatelliteOribtNode> sateOirbtNodes(this->PluseCount);
std::shared_ptr<SatelliteOribtNode[]> sateOirbtNodes(new SatelliteOribtNode[this->PluseCount], delArrPtr);
{ // 姿态计算不同
// 计算姿态
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 !!!";
}
// 回波
long echoIdx = 0;
double NearRange = this->EchoSimulationData->getNearRange(); // 近斜据
double FarRange = this->EchoSimulationData->getFarRange();
double TimgNearRange = 2 * NearRange / LIGHTSPEED;
double TimgFarRange = 2 * FarRange / LIGHTSPEED;
double Fs = this->TaskSetting->getFs(); // 距离向采样率
double Pt = this->TaskSetting->getPt() * this->TaskSetting->getGri();// 发射电压 1v
//double GainAntLen = -3;// -3dB 为天线半径
long pluseCount = this->PluseCount;
double lamda = this->TaskSetting->getCenterLamda(); // 波长
// 天线方向图
std::shared_ptr<AbstractRadiationPattern> TransformPattern = this->TaskSetting->getTransformRadiationPattern(); // 发射天线方向图
std::shared_ptr<AbstractRadiationPattern> ReceivePattern = this->TaskSetting->getReceiveRadiationPattern(); // 接收天线方向图
long PlusePoint = this->EchoSimulationData->getPlusePoints();
long echoline = Memory1GB * 4 / 16 / PlusePoint;
echoline = echoline < 1000 ? 1000 : echoline;
long startecholine = 0;
for (startecholine = 0; startecholine < pluseCount; startecholine = startecholine + echoline) {
long tempecholine = echoline;
if (startecholine + tempecholine >= pluseCount) {
tempecholine = pluseCount - startecholine;
}
std::shared_ptr<std::complex<double>> echo = this->EchoSimulationData->getEchoArr(startecholine, tempecholine);
for (long i = 0; i < tempecholine * PlusePoint; i++) {
echo.get()[i] = std::complex<double>(0, 0);
}
this->EchoSimulationData->saveEchoArr(echo, startecholine, tempecholine);
}
POLARTYPEENUM polartype = this->TaskSetting->getPolarType();
#ifndef __CUDANVCC___
QMessageBox::information(this, u8"程序提示", u8"请确定安装了CUDA库");
#else
// RTPC CUDA版本
if (pluseCount * 4 * 18 > Memory1MB * 100) {
long max = Memory1MB * 100 / 4 / 20 / PluseCount;
QMessageBox::warning(nullptr, u8"仿真场景太大了", u8"当前频点数下，脉冲数量最多为：" + QString::number(max));
}
gdalImage demxyz(this->demxyzPath);// 地面点坐标
gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型
gdalImage demsloperxyz(this->demsloperPath);// 地面坡向
// 参数与分块计算
long demRow = demxyz.height;
long demCol = demxyz.width;
long blokline = 100;
// 每块 250MB*16 = 4GB
blokline = Memory1MB * 500 / 8 / demCol;
blokline = blokline < 1 ? 1 : blokline;
bool bloklineflag = false;
// 处理发射天线方向图
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);
float* h_TantPattern = (float*)mallocCUDAHost(sizeof(float) * Tthetanum * Tphinum);
float* d_TantPattern = (float*)mallocCUDADevice(sizeof(float) * Tthetanum * Tphinum);
for (long i = 0; i < Tthetanum; i++) {
for (long j = 0; j < Tphinum; j++) {
h_TantPattern[i * Tphinum + j] = TransformPattern->getGainLearThetaPhi(TstartTheta + i * Tdtheta, TstartPhi + j * Tdphi);
}
}
HostToDevice(h_TantPattern, d_TantPattern, sizeof(float) * Tthetanum * Tphinum);
// 处理接收天线方向图
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);
float* h_RantPattern = (float*)mallocCUDAHost(sizeof(float) * Rthetanum * Rphinum);
float* d_RantPattern = (float*)mallocCUDADevice(sizeof(float) * 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);
}
}
HostToDevice(h_RantPattern, d_RantPattern, sizeof(float) * Rthetanum * Rphinum);
//处理地表覆盖
QMap<long, long> clamap;
long clamapid = 0;
long startline = 0;
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;
}
}
}
}
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;
}
}
HostToDevice(h_clsSigmaParam, d_clsSigmaParam, sizeof(CUDASigmaParam) * clamapid);
// 临时变量声明
Eigen::MatrixXd dem_x = demxyz.getData(0, 0, blokline, demxyz.width, 1); // 地面坐标
long tempDemRows = dem_x.rows();
long tempDemCols = dem_x.cols();
Eigen::MatrixXd dem_y = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd dem_z = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd demsloper_x = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd demsloper_y = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd demsloper_z = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd sloperAngle = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
float* h_dem_x = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_dem_y = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_dem_z = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_x = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_y = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_z = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_angle = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
long* h_demcls = (long*)mallocCUDAHost(sizeof(long) * blokline * tempDemCols);
float* d_dem_x = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols); // 7
float* d_dem_y = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_dem_z = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_x = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_y = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_z = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_angle = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
long* d_demcls = (long*)mallocCUDADevice(sizeof(long) * blokline * tempDemCols);
// 回波
cuComplex* h_echoAmp = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * blokline * tempDemCols);
cuComplex* d_echoAmp = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * blokline * tempDemCols);
int* h_echoAmpFID = (int*)mallocCUDAHost(sizeof(int) * blokline * tempDemCols);
int* d_echoAmpFID = (int*)mallocCUDADevice(sizeof(int) * blokline * tempDemCols);
Eigen::MatrixXd landcover = Eigen::MatrixXd::Zero(blokline, tempDemCols);// 地面覆盖类型
for (startline = 0; startline < demRow; startline = startline + blokline) {
long newblokline = blokline;
if ((startline + blokline) >= demRow) {
newblokline = demRow - startline;
bloklineflag = true;
}
dem_x = demxyz.getData(startline, 0, newblokline, demxyz.width, 1); // 地面坐标
dem_y = demxyz.getData(startline, 0, newblokline, demxyz.width, 2);
dem_z = demxyz.getData(startline, 0, newblokline, demxyz.width, 3);
demsloper_x = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 1);
demsloper_y = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 2);
demsloper_z = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 3);
sloperAngle = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 4);
landcover = demlandcls.getData(startline, 0, newblokline, demlandcls.width, 1);
if (bloklineflag) {
FreeCUDAHost(h_dem_x); FreeCUDADevice(d_dem_x);
FreeCUDAHost(h_dem_y); FreeCUDADevice(d_dem_y);
FreeCUDAHost(h_dem_z); FreeCUDADevice(d_dem_z);
FreeCUDAHost(h_demsloper_x); FreeCUDADevice(d_demsloper_x);
FreeCUDAHost(h_demsloper_y); FreeCUDADevice(d_demsloper_y);
FreeCUDAHost(h_demsloper_z); FreeCUDADevice(d_demsloper_z); //6
FreeCUDAHost(h_demsloper_angle); FreeCUDADevice(d_demsloper_angle);//7
FreeCUDAHost(h_demcls); FreeCUDADevice(d_demcls);//7
FreeCUDAHost(h_echoAmp); FreeCUDADevice(d_echoAmp);//19
FreeCUDAHost(h_echoAmpFID); FreeCUDADevice(d_echoAmpFID);//19
h_dem_x = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_dem_y = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_dem_z = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_x = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_y = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_z = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_angle = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demcls = (long*)mallocCUDAHost(sizeof(long) * newblokline * tempDemCols);
d_dem_x = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols); // 7
d_dem_y = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_z = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_x = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_y = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_z = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_angle = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demcls = (long*)mallocCUDADevice(sizeof(long) * newblokline * tempDemCols);
h_echoAmp = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * newblokline * tempDemCols);;
d_echoAmp = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * newblokline * tempDemCols);;
h_echoAmpFID = (int*)mallocCUDAHost(sizeof(int) * newblokline * tempDemCols);
d_echoAmpFID = (int*)mallocCUDADevice(sizeof(int) * newblokline * tempDemCols);
}
{ // 处理 dem -> 数量
float temp_dem_x;
float temp_dem_y;
float temp_dem_z;
float temp_demsloper_x;
float temp_demsloper_y;
float temp_demsloper_z;
float temp_sloperAngle;
long temp_demclsid;
for (long i = 0; i < newblokline; i++) {
for (long j = 0; j < demxyz.width; j++) {
temp_dem_x = float(dem_x(i, j));
temp_dem_y = float(dem_y(i, j));
temp_dem_z = float(dem_z(i, j));
temp_demsloper_x = float(demsloper_x(i, j));
temp_demsloper_y = float(demsloper_y(i, j));
temp_demsloper_z = float(demsloper_z(i, j));
temp_sloperAngle = float(sloperAngle(i, j));
temp_demclsid = long(landcover(i, j));
h_dem_x[i * demxyz.width + j] = temp_dem_x;
h_dem_y[i * demxyz.width + j] = temp_dem_y;
h_dem_z[i * demxyz.width + j] = temp_dem_z;
h_demsloper_x[i * demxyz.width + j] = temp_demsloper_x;
h_demsloper_y[i * demxyz.width + j] = temp_demsloper_y;
h_demsloper_z[i * demxyz.width + j] = temp_demsloper_z;
h_demsloper_angle[i * demxyz.width + j] = temp_sloperAngle;
h_demcls[i * demxyz.width + j] = clamap[temp_demclsid];
}
}
}
HostToDevice((void*)h_dem_x, (void*)d_dem_x, sizeof(float) * newblokline * tempDemCols); // 复制 机器 -> GPU
HostToDevice((void*)h_dem_y, (void*)d_dem_y, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_z, (void*)d_dem_z, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_x, (void*)d_demsloper_x, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_y, (void*)d_demsloper_y, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_z, (void*)d_demsloper_z, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_angle, (void*)d_demsloper_angle, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demcls, (void*)d_demcls, sizeof(float) * newblokline * tempDemCols);//地表覆盖
// 临时文件声明
float antpx = 0;
float antpy = 0;
float antpz = 0;
float antvx = 0;
float antvy = 0;
float antvz = 0;
float antdirectx = 0;
float antdirecty = 0;
float antdirectz = 0;
float antXaxisX = 0;
float antXaxisY = 0;
float antXaxisZ = 0;
float antYaxisX = 0;
float antYaxisY = 0;
float antYaxisZ = 0;
float antZaxisX = 0;
float antZaxisY = 0;
float antZaxisZ = 0;
int pixelcount = newblokline * tempDemCols;
std::cout << " GPU Memory init finished!!!!" << std::endl;
long echoline = Memory1GB * 4 / 16 / PlusePoint;
echoline = echoline < 1000 ? 1000 : echoline;
long startecholine = 0;
for (startecholine = 0; startecholine < pluseCount; startecholine = startecholine + echoline) {
long tempecholine = echoline;
if (startecholine + tempecholine >= pluseCount) {
tempecholine = pluseCount - startecholine;
}
std::shared_ptr<std::complex<double>> echo = this->EchoSimulationData->getEchoArr(startecholine, tempecholine);
long prfid = 0;
for (long tempprfid = 0; tempprfid < tempecholine; tempprfid++) {
{// 计算
prfid = tempprfid + startecholine;
// 天线位置
antpx = sateOirbtNodes[prfid].Px;
antpy = sateOirbtNodes[prfid].Py;
antpz = sateOirbtNodes[prfid].Pz;
antvx = sateOirbtNodes[prfid].Vx;
antvy = sateOirbtNodes[prfid].Vy;
antvz = sateOirbtNodes[prfid].Vz; //6
antdirectx = sateOirbtNodes[prfid].AntDirecX;
antdirecty = sateOirbtNodes[prfid].AntDirecY;
antdirectz = sateOirbtNodes[prfid].AntDirecZ; // 9 天线指向
antXaxisX = sateOirbtNodes[prfid].AntXaxisX;
antXaxisY = sateOirbtNodes[prfid].AntXaxisY;
antXaxisZ = sateOirbtNodes[prfid].AntXaxisZ;//12 天线坐标系
antYaxisX = sateOirbtNodes[prfid].AntYaxisX;
antYaxisY = sateOirbtNodes[prfid].AntYaxisY;
antYaxisZ = sateOirbtNodes[prfid].AntYaxisZ;//15
antZaxisX = sateOirbtNodes[prfid].AntZaxisX;
antZaxisY = sateOirbtNodes[prfid].AntZaxisY;
antZaxisZ = sateOirbtNodes[prfid].AntZaxisZ;//18
//CUDATestHelloWorld(1, 20);
CUDA_RTPC_SiglePRF(
antpx, antpy, antpz,// 天线坐标
antXaxisX, antXaxisY, antXaxisZ, // 天线坐标系
antYaxisX, antYaxisY, antYaxisZ, //
antZaxisX, antZaxisY, antZaxisZ,
antdirectx, antdirecty, antdirectz,// 天线指向
d_dem_x, d_dem_y, d_dem_z,
d_demcls, // 地面坐标
d_demsloper_x, d_demsloper_y, d_demsloper_z, d_demsloper_angle,// 地面坡度
d_TantPattern, TstartTheta, TstartPhi, Tdtheta, Tdphi, Tthetanum, Tphinum,// 天线方向图相关
d_RantPattern, RstartTheta, RstartPhi, Rdtheta, Rdphi, Rthetanum, Rphinum,// 天线方向图相关
lamda, Fs, NearRange, Pt, PlusePoint, // 参数
d_clsSigmaParam, clamapid,// 地表覆盖类型-sigma插值对应函数-ulaby
d_echoAmp, d_echoAmpFID,
newblokline, tempDemCols);
DeviceToHost(h_echoAmpFID, d_echoAmpFID, sizeof(long) * newblokline * tempDemCols);
DeviceToHost(h_echoAmp, d_echoAmp, sizeof(long) * newblokline * tempDemCols);
for (long i = 0; i < pixelcount; i++) {
echo.get()[tempprfid * PlusePoint + h_echoAmpFID[i]] =
echo.get()[tempprfid * PlusePoint + h_echoAmpFID[i]]
+ std::complex<double>(h_echoAmp[i].x, h_echoAmp[i].y);
}
if (tempprfid % 100 == 0) {
std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "] dem:\t" << startline << "\t-\t" << startline + newblokline <<" count:\t"<< demRow << "\t:\t pluse :\t" << prfid << " / " << pluseCount << std::endl;
}
}
}
this->EchoSimulationData->saveEchoArr(echo, startecholine, tempecholine);
}
}
std::cout << std::endl;
// 地面数据释放
FreeCUDAHost(h_dem_x); FreeCUDADevice(d_dem_x);
FreeCUDAHost(h_dem_y); FreeCUDADevice(d_dem_y);
FreeCUDAHost(h_dem_z); FreeCUDADevice(d_dem_z);
FreeCUDAHost(h_demsloper_x); FreeCUDADevice(d_demsloper_x);
FreeCUDAHost(h_demsloper_y); FreeCUDADevice(d_demsloper_y);
FreeCUDAHost(h_demsloper_z); FreeCUDADevice(d_demsloper_z); //6
FreeCUDAHost(h_demsloper_angle); FreeCUDADevice(d_demsloper_angle); //7
FreeCUDAHost(h_demcls); FreeCUDADevice(d_demcls);//7
FreeCUDAHost(h_echoAmp); FreeCUDADevice(d_echoAmp);//19
FreeCUDAHost(h_echoAmpFID); FreeCUDADevice(d_echoAmpFID);//19
FreeCUDAHost(h_TantPattern); FreeCUDADevice(d_TantPattern);
FreeCUDAHost(h_RantPattern); FreeCUDADevice(d_RantPattern);
FreeCUDAHost(h_clsSigmaParam); FreeCUDADevice(d_clsSigmaParam);
#endif
this->EchoSimulationData->saveToXml();
return ErrorCode::SUCCESS;
}
Eigen::MatrixXd plusetemp = Eigen::MatrixXd::Zero(newblokline, tempDemCols);
for (long ii = 0; ii < newblokline; ii++) {
for (long jj = 0; jj < tempDemCols; jj++) {
//plusetemp(ii, jj) = h_amp[ii * tempDemCols + jj];
plusetemp(ii, jj) = std::abs(std::complex<double>(h_echoAmp[ii * tempDemCols + jj].x, h_echoAmp[ii * tempDemCols + jj].y));
}
}
std::cout << "max:" << plusetemp.maxCoeff() << std::endl;
std::cout << "min:" << plusetemp.minCoeff() << std::endl;
Eigen::MatrixXd plusetempID = Eigen::MatrixXd::Zero(newblokline, tempDemCols);
for (long ii = 0; ii < newblokline; ii++) {
for (long jj = 0; jj < tempDemCols; jj++) {
//plusetemp(ii, jj) = h_amp[ii * tempDemCols + jj];
plusetempID(ii, jj) = h_FreqID[ii * tempDemCols + jj];
}
}
std::cout << "max ID:" << plusetempID.maxCoeff() << std::endl;
std::cout << "min ID:" << plusetempID.minCoeff() << std::endl;

6
RasterProcessTool.vcxproj
View File

@ -80,6 +80,10 @@
<DebugInformationFormat>EditAndContinue</DebugInformationFormat>
<FavorSizeOrSpeed>Speed</FavorSizeOrSpeed>
</ClCompile>
<Link>
<LinkTimeCodeGeneration>Default</LinkTimeCodeGeneration>
<LargeAddressAware>true</LargeAddressAware>
</Link>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)' == 'Debug|x64'" Label="Configuration">
<ClCompile>
@ -105,7 +109,7 @@
<Link>
<SubSystem>Console</SubSystem>
<GenerateDebugInformation>DebugFull</GenerateDebugInformation>
<EnableCOMDATFolding>true</EnableCOMDATFolding>
<EnableCOMDATFolding>false</EnableCOMDATFolding>
<OptimizeReferences>true</OptimizeReferences>
</Link>
</ItemDefinitionGroup>

12
SimulationSAR/QImageSARRTPC.ui
View File

@ -7,7 +7,7 @@
<x>0</x>
<y>0</y>
<width>873</width>
<height>500</height>
<height>499</height>
</rect>
</property>
<property name="windowTitle">
@ -25,7 +25,7 @@
<x>0</x>
<y>0</y>
<width>853</width>
<height>451</height>
<height>450</height>
</rect>
</property>
<layout class="QGridLayout" name="gridLayout">
@ -77,7 +77,7 @@
</size>
</property>
<property name="text">
<string>D:/Programme/vs2022/RasterMergeTest/TestData/ant/ant_model_setting_Horn_conical1_FarField-phi.csv</string>
<string>D:/Programme/vs2022/RasterMergeTest/TestData/ant/ant_model_setting_Horn_conical1_FarField-receive.csv</string>
</property>
</widget>
</item>
@ -90,7 +90,7 @@
</size>
</property>
<property name="text">
<string>D:/Programme/vs2022/RasterMergeTest/TestData/landcover_aligned2</string>
<string>D:/Programme/vs2022/RasterMergeTest/TestData/landcover_aligned2.dat</string>
</property>
</widget>
</item>
@ -103,7 +103,7 @@
</size>
</property>
<property name="text">
<string>D:/Programme/vs2022/RasterMergeTest/TestData/outData/</string>
<string>D:/Programme/vs2022/RasterMergeTest/LAMPCAE_SCANE/</string>
</property>
</widget>
</item>
@ -142,7 +142,7 @@
</size>
</property>
<property name="text">
<string>D:/Programme/vs2022/RasterMergeTest/TestData/ant/ant_model_setting_Horn_conical1_FarField-theta.csv</string>
<string>D:/Programme/vs2022/RasterMergeTest/TestData/ant/ant_model_setting_Horn_conical1_FarField-trans.csv</string>
</property>
</widget>
</item>

658
SimulationSAR/RTPCProcessCls.cpp
View File

@ -22,6 +22,7 @@
#ifdef __CUDANVCC___
#include "GPUTool.cuh"
#endif // __CUDANVCC___
#include <Imageshow/ImageShowDialogClass.h>
@ -130,7 +131,8 @@ ErrorCode RTPCProcessCls::Process(long num_thread)
else {}
qDebug() << "RTPCMainProcess";
stateCode = this->RTPCMainProcess(num_thread);
//stateCode = this->RTPCMainProcess(num_thread);
stateCode = this->RTPCMainProcess_GPU( );
if (stateCode != ErrorCode::SUCCESS) {
return stateCode;
@ -240,7 +242,7 @@ ErrorCode RTPCProcessCls::DEMPreprocess()
line_invert = 1000;
long start_ids = 0;
long dem_rows = 0, dem_cols = 0;
for (start_ids = 1; start_ids < demds.height; start_ids = start_ids + line_invert) {
Eigen::MatrixXd demdata = demds.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
long startlineid = start_ids;
@ -297,9 +299,10 @@ ErrorCode RTPCProcessCls::DEMPreprocess()
return ErrorCode::SUCCESS;
}
ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
ErrorCode RTPCProcessCls::RTPCMainProcess_GPU( )
{
omp_set_num_threads(num_thread);// 设置openmp 线程数量
double widthSpace = LIGHTSPEED / 2 / this->TaskSetting->getFs();
double prf_time = 0;
@ -313,8 +316,10 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
long double imageStarttime = 0;
imageStarttime = this->TaskSetting->getSARImageStartTime();
//std::vector<SatelliteOribtNode> sateOirbtNodes(this->PluseCount);
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;
@ -366,9 +371,6 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
// 回波
long echoIdx = 0;
double NearRange = this->EchoSimulationData->getNearRange(); // 近斜据
double FarRange = this->EchoSimulationData->getFarRange();
@ -384,18 +386,12 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
// 天线方向图
std::shared_ptr<AbstractRadiationPattern> TransformPattern = this->TaskSetting->getTransformRadiationPattern(); // 发射天线方向图
std::shared_ptr<AbstractRadiationPattern> ReceivePattern = this->TaskSetting->getReceiveRadiationPattern(); // 接收天线方向图
std::shared_ptr<std::complex<double>> echo = this->EchoSimulationData->getEchoArr();
long PlusePoint = this->EchoSimulationData->getPlusePoints();
// 初始化 为 0
for (long i = 0; i < pluseCount * PlusePoint; i++) {
echo.get()[i] = std::complex<double>(0, 0);
}
this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
// 初始化回波
this->EchoSimulationData->initEchoArr(std::complex<double>(0, 0));
POLARTYPEENUM polartype = this->TaskSetting->getPolarType();
#ifndef __CUDANVCC___
QMessageBox::information(this, u8"程序提示", u8"请确定安装了CUDA库");
#else
@ -403,9 +399,9 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
// RTPC CUDA版本
if (pluseCount * 4 * 18 > Memory1MB * 100) {
long max = Memory1MB * 100 / 4 / 20 / PluseCount;
QMessageBox::warning(nullptr, u8"仿真场景太大了", u8"当前频点数下，脉冲数量最多为："+QString::number(max));
QMessageBox::warning(nullptr, u8"仿真场景太大了", u8"当前频点数下，脉冲数量最多为：" + QString::number(max));
}
gdalImage demxyz(this->demxyzPath);// 地面点坐标
gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型
gdalImage demsloperxyz(this->demsloperPath);// 地面坡向
@ -435,17 +431,26 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
double Tdtheta = (Tmaxtheta - Tmintheta) / (Tthetanum - 1);
double Tdphi = (Tmaxphi - Tminphi) / (Tphinum - 1);
float* h_TantPattern = (float*)mallocCUDAHost(sizeof(float) * Tthetanum * Tphinum);
float* d_TantPattern = (float*)mallocCUDADevice(sizeof(float) * Tthetanum * Tphinum);
for (long i = 0; i < Tthetanum; i++) {
for (long j = 0; j < Tphinum; j++) {
h_TantPattern[i*Tphinum+j] = TransformPattern->getGainLearThetaPhi(TstartTheta + i * Tdtheta, TstartPhi + j * Tdphi);
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);
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);
}
}
HostToDevice(h_TantPattern, d_TantPattern, sizeof(float)* Tthetanum* Tphinum);
// 处理接收天线方向图
double Rminphi = ReceivePattern->getMinPhi();
double Rmaxphi = ReceivePattern->getMaxPhi();
@ -463,13 +468,25 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
float* h_RantPattern = (float*)mallocCUDAHost(sizeof(float) * Rthetanum * Rphinum);
float* d_RantPattern = (float*)mallocCUDADevice(sizeof(float) * 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->getGainLearThetaPhi(RstartTheta + i * Rdtheta, RstartPhi + j * Rdphi);
h_RantPattern[i * Rphinum + j] = ReceivePattern->getGain(RstartTheta + i * Rdtheta, RstartPhi + j * Rdphi);
}
}
HostToDevice(h_RantPattern, d_RantPattern, sizeof(float)* Rthetanum* Rphinum);
testOutAntPatternTrans("ReceivePattern.bin", h_RantPattern, Rmintheta, Rdtheta, RstartPhi, Rdphi, Rthetanum, Rphinum);
for (long i = 0; i < Tthetanum; i++) {
for (long j = 0; j < Tphinum; j++) {
h_RantPattern[i * Tphinum + j] = powf(10.0, h_RantPattern[i * Tphinum + j] / 10);
}
}
HostToDevice(h_RantPattern, d_RantPattern, sizeof(float) * Rthetanum * Rphinum);
//处理地表覆盖
QMap<long, long> clamap;
long clamapid = 0;
@ -491,11 +508,17 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
}
}
}
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);
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;
@ -505,10 +528,28 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
h_clsSigmaParam[clamap[id]].p5 = tempp.p5;
h_clsSigmaParam[clamap[id]].p6 = tempp.p6;
}
}
HostToDevice(h_clsSigmaParam, d_clsSigmaParam, sizeof(CUDASigmaParam) * clamapid);
// 临时变量声明
// 打印日志
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);
Eigen::MatrixXd dem_x = demxyz.getData(0, 0, blokline, demxyz.width, 1); // 地面坐标
long tempDemRows = dem_x.rows();
long tempDemCols = dem_x.cols();
@ -520,33 +561,115 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
Eigen::MatrixXd demsloper_z = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd sloperAngle = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
float* h_dem_x = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_dem_y = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_dem_z = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_x = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_y = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_z = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_angle = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
long* h_demcls = (long*)mallocCUDAHost(sizeof(long) * blokline * tempDemCols);
float* h_dem_x;
float* h_dem_y;
float* h_dem_z;
float* h_demsloper_x;
float* h_demsloper_y;
float* h_demsloper_z;
float* h_demsloper_angle;
float* d_dem_x = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols); // 7
float* d_dem_y = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_dem_z = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_x = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_y = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_z = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_angle = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
long* d_demcls = (long*)mallocCUDADevice(sizeof(long) * blokline * tempDemCols);
float* d_dem_x;
float* d_dem_y;
float* d_dem_z;
float* d_demsloper_x;
float* d_demsloper_y;
float* d_demsloper_z;
float* d_demsloper_angle;
h_dem_x=(float* )mallocCUDAHost( sizeof(float) * blokline * tempDemCols);
h_dem_y=(float* )mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
h_dem_z=(float* )mallocCUDAHost( sizeof(float) * blokline * tempDemCols);
h_demsloper_x=(float* )mallocCUDAHost( sizeof(float) * blokline * tempDemCols);
h_demsloper_y=(float* )mallocCUDAHost( sizeof(float) * blokline * tempDemCols);
h_demsloper_z=(float* )mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
h_demsloper_angle= (float*)mallocCUDAHost( sizeof(float) * blokline * tempDemCols);
d_dem_x=(float* )mallocCUDADevice( sizeof(float) * blokline * tempDemCols); // 7
d_dem_y=(float* )mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_dem_z=(float* )mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_demsloper_x=(float* )mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_demsloper_y=(float* )mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_demsloper_z=(float* )mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_demsloper_angle= (float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
float* h_dem_theta; // 天线方向图
float* h_dem_phi;
float* d_dem_theta;
float* d_dem_phi;
h_dem_theta=(float* )mallocCUDAHost( sizeof(float) * blokline * tempDemCols);
h_dem_phi= (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
d_dem_theta= (float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);// 9
d_dem_phi= (float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
HostToDevice((void*)h_dem_theta, (void*)d_dem_theta, sizeof(float) * blokline * tempDemCols);
HostToDevice((void*)h_dem_phi, (void*)d_dem_phi, sizeof(float) * blokline * tempDemCols);
// 提前声明参数变量
float* h_R;// 辐射方向
float* h_localangle;//入射角
float* d_R;// 辐射方向
float* d_localangle;//入射角
h_R=(float* )mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
h_localangle= (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols); // 11
d_R= (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
d_localangle= (float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
float* h_RstX;
float* h_RstY;
float* h_RstZ;
float* d_RstX;
float* d_RstY;
float* d_RstZ;
h_RstX=(float*)mallocCUDAHost( sizeof(float) * blokline * tempDemCols);
h_RstY=(float*)mallocCUDAHost( sizeof(float) * blokline * tempDemCols);
h_RstZ=(float*)mallocCUDAHost( sizeof(float) * blokline * tempDemCols); // 14
d_RstX=(float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_RstY=(float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_RstZ=(float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
float* h_sigma0;
float* h_TransAnt;
float* h_ReciveAnt;
float* d_sigma0;
float* d_TransAnt;
float* d_ReciveAnt;
h_sigma0= (float*)mallocCUDAHost( sizeof(float)* blokline* tempDemCols);
h_TransAnt = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
h_ReciveAnt = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols); // 17
d_sigma0= (float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_TransAnt= (float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
d_ReciveAnt= (float*)mallocCUDADevice( sizeof(float) * blokline * tempDemCols);
// 回波
cuComplex* h_echoAmp = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * blokline * tempDemCols);
cuComplex* d_echoAmp = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * blokline * tempDemCols);
int* h_echoAmpFID = (int*)mallocCUDAHost(sizeof(int) * blokline * tempDemCols);
int* d_echoAmpFID = (int*)mallocCUDADevice(sizeof(int) * blokline * tempDemCols);
cuComplex* h_echoAmp;
cuComplex* d_echoAmp;
h_echoAmp=(cuComplex*)mallocCUDAHost(sizeof(cuComplex) * blokline * tempDemCols);
d_echoAmp=(cuComplex*)mallocCUDADevice( sizeof(cuComplex) * blokline * tempDemCols); //19
Eigen::MatrixXd landcover= Eigen::MatrixXd::Zero(blokline, tempDemCols);// 地面覆盖类型
long* h_FreqID;
long* d_FreqID;
h_FreqID=(long*)mallocCUDAHost( sizeof(long) * blokline * tempDemCols);
d_FreqID=(long*)mallocCUDADevice( sizeof(long) * blokline * tempDemCols); //21
// 地表覆盖类型
Eigen::MatrixXd landcover = Eigen::MatrixXd::Zero(blokline, tempDemCols);// 地面覆盖类型
long* h_demcls = (long*)mallocCUDAHost(sizeof(long) * blokline * tempDemCols);
long* d_demcls = (long*)mallocCUDADevice(sizeof(long) * blokline * tempDemCols);
float* h_amp=(float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* d_amp=(float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
for (startline = 0; startline < demRow; startline = startline + blokline) {
long newblokline = blokline;
@ -554,7 +677,6 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
newblokline = demRow - startline;
bloklineflag = true;
}
dem_x = demxyz.getData(startline, 0, newblokline, demxyz.width, 1); // 地面坐标
dem_y = demxyz.getData(startline, 0, newblokline, demxyz.width, 2);
dem_z = demxyz.getData(startline, 0, newblokline, demxyz.width, 3);
@ -562,8 +684,8 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
demsloper_y = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 2);
demsloper_z = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 3);
sloperAngle = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 4);
landcover = demlandcls.getData(startline, 0, newblokline, demlandcls.width, 1);
landcover = demlandcls.getData(startline, 0, newblokline, demlandcls.width, 1);
if (bloklineflag) {
FreeCUDAHost(h_dem_x); FreeCUDADevice(d_dem_x);
@ -573,9 +695,20 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
FreeCUDAHost(h_demsloper_y); FreeCUDADevice(d_demsloper_y);
FreeCUDAHost(h_demsloper_z); FreeCUDADevice(d_demsloper_z); //6
FreeCUDAHost(h_demsloper_angle); FreeCUDADevice(d_demsloper_angle);//7
FreeCUDAHost(h_demcls); FreeCUDADevice(d_demcls);//7
FreeCUDAHost(h_dem_theta); FreeCUDADevice(d_dem_theta);
FreeCUDAHost(h_dem_phi); FreeCUDADevice(d_dem_phi); //9
FreeCUDAHost(h_R); FreeCUDADevice(d_R);
FreeCUDAHost(h_localangle); FreeCUDADevice(d_localangle); //11
FreeCUDAHost(h_RstX); FreeCUDADevice(d_RstX);
FreeCUDAHost(h_RstY); FreeCUDADevice(d_RstY);
FreeCUDAHost(h_RstZ); FreeCUDADevice(d_RstZ); //14
FreeCUDAHost(h_sigma0); FreeCUDADevice(d_sigma0);
FreeCUDAHost(h_TransAnt); FreeCUDADevice(d_TransAnt);
FreeCUDAHost(h_ReciveAnt); FreeCUDADevice(d_ReciveAnt); //17
FreeCUDAHost(h_echoAmp); FreeCUDADevice(d_echoAmp);//19
FreeCUDAHost(h_echoAmpFID); FreeCUDADevice(d_echoAmpFID);//19
FreeCUDAHost(h_FreqID); FreeCUDADevice(d_FreqID);//20
FreeCUDAHost(h_demcls); FreeCUDADevice(d_demcls);
FreeCUDAHost(h_amp); FreeCUDADevice(d_amp);
h_dem_x = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_dem_y = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
@ -583,185 +716,234 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
h_demsloper_x = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_y = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_z = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_angle = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_angle = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
h_dem_theta = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_dem_phi = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_R = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_localangle = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_RstX = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_RstY = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_RstZ = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_sigma0 = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_TransAnt = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_ReciveAnt = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_echoAmp = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * newblokline * tempDemCols);
h_FreqID = (long*)mallocCUDAHost(sizeof(long) * newblokline * tempDemCols);
h_demcls = (long*)mallocCUDAHost(sizeof(long) * newblokline * tempDemCols);
h_amp = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
d_dem_x = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols); // 7
d_dem_y = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_z = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_x = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_y = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_z = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_angle = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_x=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_y=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_z=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_x=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_y=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_z=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);//6
d_demsloper_angle=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);//7
d_dem_theta=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_phi=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);// 9
d_R=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_localangle=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_RstX=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_RstY=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_RstZ=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_sigma0=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_TransAnt=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_ReciveAnt=(float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_echoAmp=(cuComplex*)mallocCUDADevice(sizeof(cuComplex) * newblokline * tempDemCols);
d_FreqID=(long*)mallocCUDADevice(sizeof(long) * newblokline * tempDemCols);
d_demcls = (long*)mallocCUDADevice(sizeof(long) * newblokline * tempDemCols);
h_echoAmp = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * newblokline * tempDemCols);;
d_echoAmp = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * newblokline * tempDemCols);;
h_echoAmpFID = (int*)mallocCUDAHost(sizeof(int) * newblokline * tempDemCols);
d_echoAmpFID = (int*)mallocCUDADevice(sizeof(int) * newblokline * tempDemCols);
d_amp = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
}
{ // 处理 dem -> 数量
float temp_dem_x;
float temp_dem_y;
float temp_dem_z;
float temp_demsloper_x;
float temp_demsloper_y;
float temp_demsloper_z;
float temp_sloperAngle;
long temp_demclsid;
for (long i = 0; i < newblokline; i++) {
for (long j = 0; j < demxyz.width; j++) {
temp_dem_x= float(dem_x(i, j)) ;
temp_dem_y= float(dem_y(i, j)) ;
temp_dem_z= float(dem_z(i, j)) ;
temp_demsloper_x= float(demsloper_x(i, j));
temp_demsloper_y= float(demsloper_y(i, j));
temp_demsloper_z= float(demsloper_z(i, j));
temp_sloperAngle= float(sloperAngle(i, j));
temp_demclsid = long(landcover(i,j));
h_dem_x[i * demxyz.width + j] = temp_dem_x ;
h_dem_y[i * demxyz.width + j] = temp_dem_y ;
h_dem_z[i * demxyz.width + j] = temp_dem_z ;
h_demsloper_x[i * demxyz.width + j] = temp_demsloper_x ;
h_demsloper_y[i * demxyz.width + j] = temp_demsloper_y ;
h_demsloper_z[i * demxyz.width + j] = temp_demsloper_z ;
h_demsloper_angle[i * demxyz.width + j] = temp_sloperAngle;
h_demcls[i * demxyz.width + j] = clamap[temp_demclsid];
}
//# pragma omp parallel for
for (long i = 0; i < newblokline; i++) {
for (long j = 0; j < demxyz.width; j++) {
h_dem_x[i * demxyz.width + j] = float(dem_x(i, j));
h_dem_y[i * demxyz.width + j] = float(dem_y(i, j));
h_dem_z[i * demxyz.width + j] = float(dem_z(i, j));
h_demsloper_x[i * demxyz.width + j] = float(demsloper_x(i, j));
h_demsloper_y[i * demxyz.width + j] = float(demsloper_y(i, j));
h_demsloper_z[i * demxyz.width + j] = float(demsloper_z(i, j));
h_demsloper_angle[i * demxyz.width + j] = float(sloperAngle(i, j));
h_demcls[i * demxyz.width + j] = clamap[long(landcover(i, j))];
h_amp[i * demxyz.width + j] = 0.0f;
}
}
HostToDevice((void*)h_dem_x, (void*)d_dem_x, sizeof(float) * newblokline * tempDemCols); // 复制 机器 -> GPU
HostToDevice((void*)h_dem_y, (void*)d_dem_y, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_z, (void*)d_dem_z, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_x, (void*)d_demsloper_x, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_y, (void*)d_demsloper_y, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_z, (void*)d_demsloper_z, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_angle, (void*)d_demsloper_angle, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demcls, (void*)d_demcls, sizeof(float) * newblokline * tempDemCols);//地表覆盖
HostToDevice((void*)h_demsloper_angle, (void*)d_demsloper_angle, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_theta, (void*)d_dem_theta, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_phi, (void*)d_dem_phi, sizeof(float) * newblokline* tempDemCols);
HostToDevice((void*)h_demcls, (void*)d_demcls, sizeof(long) * newblokline* tempDemCols);
HostToDevice((void*)h_amp, (void*)d_amp, sizeof(float) * newblokline* tempDemCols);
// 临时文件声明
float antpx =0;
float antpy =0;
float antpz =0;
float antvx =0;
float antvy =0;
float antvz =0;
float antdirectx =0;
float antdirecty =0;
float antdirectz =0;
float antXaxisX =0;
float antXaxisY =0;
float antXaxisZ =0;
float antYaxisX =0;
float antYaxisY =0;
float antYaxisZ =0;
float antZaxisX =0;
float antZaxisY =0;
float antZaxisZ = 0;
long pixelcount = newblokline * tempDemCols;
long echoblockline = Memory1MB * 2000 / 8 / 2 / PlusePoint;
long startprfid = 0;
for (startprfid = 0; startprfid < pluseCount; startprfid = startprfid + echoblockline) {
long templine = startprfid + echoblockline < PluseCount ? echoblockline : PluseCount - startprfid;
std::shared_ptr<std::complex<double>> echotemp = this->EchoSimulationData->getEchoArr(startprfid, templine);
for (long tempprfid = 0; tempprfid < templine; tempprfid++) {
{// 计算
long prfid = tempprfid + startprfid;
// 天线位置
float antpx = sateOirbtNodes[prfid].Px;
float antpy = sateOirbtNodes[prfid].Py;
float antpz = sateOirbtNodes[prfid].Pz;
float antvx = sateOirbtNodes[prfid].Vx;
float antvy = sateOirbtNodes[prfid].Vy;
float antvz = sateOirbtNodes[prfid].Vz; //6
float antdirectx = sateOirbtNodes[prfid].AntDirecX;
float antdirecty = sateOirbtNodes[prfid].AntDirecY;
float antdirectz = sateOirbtNodes[prfid].AntDirecZ; // 9 天线指向
float antXaxisX = sateOirbtNodes[prfid].AntXaxisX;
float antXaxisY = sateOirbtNodes[prfid].AntXaxisY;
float antXaxisZ = sateOirbtNodes[prfid].AntXaxisZ;//12 天线坐标系
float antYaxisX = sateOirbtNodes[prfid].AntYaxisX;
float antYaxisY = sateOirbtNodes[prfid].AntYaxisY;
float antYaxisZ = sateOirbtNodes[prfid].AntYaxisZ;//15
float antZaxisX = sateOirbtNodes[prfid].AntZaxisX;
float antZaxisY = sateOirbtNodes[prfid].AntZaxisY;
float antZaxisZ = sateOirbtNodes[prfid].AntZaxisZ;//18
#ifdef __PRFDEBUG__
testOutAmpArr("h_dem_x.bin", h_dem_x, newblokline, tempDemCols);
testOutAmpArr("h_dem_y.bin", h_dem_y, newblokline, tempDemCols);
testOutAmpArr("h_dem_z.bin", h_dem_z, newblokline, tempDemCols);
std::cout << "ant Position=[" << antpx << "," << antpy << "," << antpz << "]" << std::endl;
#endif // __PRFDEBUG__
int pixelcount = newblokline * tempDemCols;
std::cout << " GPU Memory init finished!!!!" << std::endl;
for (long prfid = 0; prfid < pluseCount; prfid++) {
{// 计算
// 天线位置
antpx = sateOirbtNodes[prfid].Px;
antpy = sateOirbtNodes[prfid].Py;
antpz = sateOirbtNodes[prfid].Pz;
antvx = sateOirbtNodes[prfid].Vx;
antvy = sateOirbtNodes[prfid].Vy;
antvz = sateOirbtNodes[prfid].Vz; //6
antdirectx = sateOirbtNodes[prfid].AntDirecX;
antdirecty = sateOirbtNodes[prfid].AntDirecY;
antdirectz = sateOirbtNodes[prfid].AntDirecZ; // 9 天线指向
antXaxisX = sateOirbtNodes[prfid].AntXaxisX;
antXaxisY = sateOirbtNodes[prfid].AntXaxisY;
antXaxisZ = sateOirbtNodes[prfid].AntXaxisZ;//12 天线坐标系
antYaxisX = sateOirbtNodes[prfid].AntYaxisX;
antYaxisY = sateOirbtNodes[prfid].AntYaxisY;
antYaxisZ = sateOirbtNodes[prfid].AntYaxisZ;//15
antZaxisX = sateOirbtNodes[prfid].AntZaxisX;
antZaxisY = sateOirbtNodes[prfid].AntZaxisY;
antZaxisZ = sateOirbtNodes[prfid].AntZaxisZ;//18
//CUDATestHelloWorld(1, 20);
CUDA_RTPC_SiglePRF(
antpx, antpy, antpz,// 天线坐标
antXaxisX, antXaxisY, antXaxisZ, // 天线坐标系
antYaxisX, antYaxisY, antYaxisZ, //
antZaxisX, antZaxisY, antZaxisZ,
antdirectx, antdirecty, antdirectz,// 天线指向
d_dem_x, d_dem_y, d_dem_z,
d_demcls, // 地面坐标
d_demsloper_x, d_demsloper_y, d_demsloper_z, d_demsloper_angle,// 地面坡度
d_TantPattern, TstartTheta,TstartPhi, Tdtheta, Tdphi, Tthetanum, Tphinum,// 天线方向图相关
d_RantPattern, RstartTheta, RstartPhi, Rdtheta, Rdphi, Rthetanum,Rphinum,// 天线方向图相关
lamda, Fs, NearRange, Pt, PlusePoint, // 参数
d_clsSigmaParam, clamapid,// 地表覆盖类型-sigma插值对应函数-ulaby
d_echoAmp, d_echoAmpFID,
newblokline,tempDemCols);
DeviceToHost(h_echoAmpFID, d_echoAmpFID, sizeof(long)* newblokline* tempDemCols);
DeviceToHost(h_echoAmp, d_echoAmp, sizeof(long)* newblokline* tempDemCols);
for (long i = 0; i < pixelcount; i++) {
echo.get()[prfid*PlusePoint+ h_echoAmpFID[i]] = echo.get()[prfid * PlusePoint + h_echoAmpFID[i]]
+std::complex<double>(h_echoAmp[i].x, h_echoAmp[i].y);
make_VectorA_B(antpx, antpy, antpz, d_dem_x, d_dem_y, d_dem_z, d_RstX, d_RstY, d_RstZ, pixelcount); // Rst = Rs - Rt; 地面-> 指向
Norm_Vector(d_RstX, d_RstY, d_RstZ, d_R, pixelcount); // R
cosAngle_VA_AB(d_RstX, d_RstY, d_RstZ, d_demsloper_x, d_demsloper_y, d_demsloper_z, d_localangle, pixelcount); // 局部入射角
SatelliteAntDirectNormal(d_RstX, d_RstY, d_RstZ,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antdirectx, antdirecty, antdirectz,
d_dem_theta, d_dem_phi, pixelcount);// 计算角度
#ifdef __PRFDEBUG__
DeviceToHost(h_RstX, d_RstX, sizeof(float)* pixelcount);
DeviceToHost(h_RstY, d_RstY, sizeof(float)* pixelcount);
DeviceToHost(h_RstZ, d_RstZ, sizeof(float)* pixelcount);
testOutAmpArr("h_RstX.bin", h_RstX, newblokline, tempDemCols);
testOutAmpArr("h_RstY.bin", h_RstY, newblokline, tempDemCols);
testOutAmpArr("h_RstZ.bin", h_RstZ, newblokline, tempDemCols);
#endif // __PRFDEBUG__
AntPatternInterpGain(d_dem_theta, d_dem_phi, d_TransAnt, d_TantPattern, TstartTheta, TstartPhi, Tdtheta, Tdphi, Tthetanum, Tphinum, pixelcount);
AntPatternInterpGain(d_dem_theta, d_dem_phi, d_ReciveAnt, d_RantPattern, RstartTheta, RstartPhi, Rdtheta, Rdphi, Rthetanum, Rphinum, pixelcount);
#ifdef __PRFDEBUG__
DeviceToHost(h_dem_theta, d_dem_theta, sizeof(float)* pixelcount); // 从GPU -> 主机
DeviceToHost(h_dem_phi, d_dem_phi, sizeof(float)* pixelcount);
DeviceToHost(h_localangle, d_localangle, sizeof(float)* pixelcount);
testOutAmpArr("h_localangle.bin", h_localangle, newblokline, tempDemCols);
DeviceToHost(h_TransAnt, d_TransAnt, sizeof(float)* pixelcount);
DeviceToHost(h_ReciveAnt, d_ReciveAnt, sizeof(float)* pixelcount);
testOutAmpArr("ant_theta.bin", h_dem_theta, newblokline, tempDemCols);
testOutAmpArr("ant_phi.bin", h_dem_phi, newblokline, tempDemCols);
testOutAmpArr("antPattern_Trans.bin", h_TransAnt, newblokline, tempDemCols);
testOutAmpArr("antPattern_Receive.bin", h_ReciveAnt, newblokline, tempDemCols);
#endif // __PRFDEBUG__
CUDAInterpSigma(d_demcls, d_amp, d_localangle, pixelcount, d_clsSigmaParam, clamapid);
#ifdef __PRFDEBUG__
DeviceToHost(h_amp, d_amp, sizeof(float)* pixelcount);
testOutAmpArr("amp.bin", h_amp,newblokline, tempDemCols);
testOutClsArr("h_demcls.bin", h_demcls,newblokline, tempDemCols);
#endif // __PRFDEBUG__
// 计算回波
calculationEcho(d_amp, d_TransAnt, d_ReciveAnt, d_localangle, d_R, d_demsloper_angle, NearRange, Fs, Pt, lamda, PlusePoint, d_echoAmp, d_FreqID, pixelcount);
DeviceToHost(h_echoAmp, d_echoAmp, sizeof(cuComplex) * pixelcount);
DeviceToHost(h_FreqID, d_FreqID, sizeof(long) * pixelcount);
//DeviceToHost(h_amp, d_amp, sizeof(float) * pixelcount);
#ifdef __PRFDEBUG__
float* h_echoAmp_real = (float*)mallocCUDAHost(sizeof(float) * pixelcount);
float* h_echoAmp_imag = (float*)mallocCUDAHost(sizeof(float) * pixelcount);
for (long freqi = 0; freqi < pixelcount; freqi++) {
h_echoAmp_real[freqi] = h_echoAmp[freqi].x;
h_echoAmp_imag[freqi] = h_echoAmp[freqi].y;
}
testOutAmpArr("h_echoAmp_real.bin", h_echoAmp_real, newblokline, tempDemCols);
testOutAmpArr("h_echoAmp_imag.bin", h_echoAmp_imag, newblokline, tempDemCols);
testOutClsArr("h_FreqID.bin", h_FreqID, newblokline, tempDemCols);
FreeCUDAHost(h_echoAmp_real);
FreeCUDAHost(h_echoAmp_imag);
exit(0);
#endif // __PRFDEBUG__
for (long freqi = 0; freqi < pixelcount; freqi++) {
long pluseid = h_FreqID[freqi];
echotemp.get()[tempprfid * PlusePoint + pluseid] = std::complex<double>(h_echoAmp[freqi].x, h_echoAmp[freqi].y);
}
if (prfid % 1000 == 0) {
std::cout << "[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "] dem:\t" << startline << "\t-\t" << startline + newblokline << "\t:\t pluse :\t" << prfid << " / " << pluseCount << std::endl;
}
}
//for (long i = 0; i < PlusePoint; i++) {
// std::cout << echo.get()[prfid * PlusePoint + i] << std::endl;
//}
if (prfid % 100 == 0) {
std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "] dem:\t" << startline << "\t-\t" << startline + newblokline << "\t:\t pluse :\t" << prfid << " / " << pluseCount << std::endl;
}
}
this->EchoSimulationData->saveEchoArr(echotemp, startprfid, templine);
}
}
}
std::cout << std::endl;
// 地面数据释放
FreeCUDAHost(h_dem_x); FreeCUDADevice(d_dem_x);
FreeCUDAHost(h_dem_y); FreeCUDADevice(d_dem_y);
FreeCUDAHost(h_dem_z); FreeCUDADevice(d_dem_z);
FreeCUDAHost(h_demsloper_x); FreeCUDADevice(d_demsloper_x);
FreeCUDAHost(h_demsloper_y); FreeCUDADevice(d_demsloper_y);
FreeCUDAHost(h_demsloper_z); FreeCUDADevice(d_demsloper_z); //6
FreeCUDAHost(h_demsloper_angle); FreeCUDADevice(d_demsloper_angle); //7
FreeCUDAHost(h_demcls); FreeCUDADevice(d_demcls);//7
FreeCUDAHost(h_echoAmp); FreeCUDADevice(d_echoAmp);//19
FreeCUDAHost(h_echoAmpFID); FreeCUDADevice(d_echoAmpFID);//19
FreeCUDAHost(h_dem_x); FreeCUDADevice(d_dem_x);
FreeCUDAHost(h_dem_y); FreeCUDADevice(d_dem_y);
FreeCUDAHost(h_dem_z); FreeCUDADevice(d_dem_z);
FreeCUDAHost(h_demsloper_x); FreeCUDADevice(d_demsloper_x);
FreeCUDAHost(h_demsloper_y); FreeCUDADevice(d_demsloper_y);
FreeCUDAHost(h_demsloper_z); FreeCUDADevice(d_demsloper_z); //6
FreeCUDAHost(h_demsloper_angle); FreeCUDADevice(d_demsloper_angle); //7
// 临时变量释放
FreeCUDAHost(h_dem_theta); FreeCUDADevice(d_dem_theta);
FreeCUDAHost(h_dem_phi); FreeCUDADevice(d_dem_phi);// 9
FreeCUDAHost(h_R); FreeCUDADevice(d_R);
FreeCUDAHost(h_localangle); FreeCUDADevice(h_localangle); //11
FreeCUDAHost(h_RstX); FreeCUDADevice(d_RstX);
FreeCUDAHost(h_RstY); FreeCUDADevice(d_RstY);
FreeCUDAHost(h_RstZ); FreeCUDADevice(d_RstZ); //14
FreeCUDAHost(h_sigma0); FreeCUDADevice(d_sigma0);
FreeCUDAHost(h_TransAnt); FreeCUDADevice(d_TransAnt);
FreeCUDAHost(h_ReciveAnt); FreeCUDADevice(d_ReciveAnt); //17
FreeCUDAHost(h_echoAmp); FreeCUDADevice(d_echoAmp);//19
FreeCUDAHost(h_FreqID); FreeCUDADevice(d_FreqID);//20
FreeCUDAHost(h_demcls); FreeCUDADevice(d_demcls);
FreeCUDAHost(h_amp); FreeCUDADevice(d_amp);
FreeCUDAHost(h_TantPattern); FreeCUDADevice(d_TantPattern);
FreeCUDAHost(h_RantPattern); FreeCUDADevice(d_RantPattern);
FreeCUDAHost(h_clsSigmaParam); FreeCUDADevice(d_clsSigmaParam);
#endif
this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
this->EchoSimulationData->saveToXml();
return ErrorCode::SUCCESS;
}
void RTPCProcessMain(long num_thread, QString TansformPatternFilePath, QString ReceivePatternFilePath, QString simulationtaskName, QString OutEchoPath, QString GPSXmlPath, QString TaskXmlPath, QString demTiffPath, QString LandCoverPath, QString HHSigmaPath, QString HVSigmaPath, QString VHSigmaPath, QString VVSigmaPath)
{
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);
std::shared_ptr < AbstractSARSatelliteModel> task = ReadSimulationSettingsXML(TaskXmlPath);
@ -785,6 +967,12 @@ void RTPCProcessMain(long num_thread, QString TansformPatternFilePath, QString R
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);
@ -825,3 +1013,73 @@ void RTPCProcessMain(long num_thread, QString TansformPatternFilePath, QString R
rtpc.Process(num_thread); // 处理程序
qDebug() << "-------------- RTPC end---------------------------------------";
}
void testOutAntPatternTrans(QString antpatternfilename,float* antPatternArr,
double starttheta, double deltetheta,
double startphi, double deltaphi,
long thetanum, long phinum)
{
Eigen::MatrixXd antPatternMatrix(thetanum,phinum);
for (long t = 0; t < thetanum; ++t) {
for (long p = 0; p < phinum; ++p) {
long index = t * phinum + p;
if (index < thetanum * phinum) {
antPatternMatrix(t, p) = static_cast<double>(antPatternArr[index]); // Copy to Eigen matrix
}
}
}
Eigen::MatrixXd gt(2, 3);
gt(0, 0)=startphi;//x
gt(0, 1)=deltaphi;
gt(0, 2)=0;
gt(1, 0)=starttheta;
gt(1, 1)=0;
gt(1, 2)=deltetheta;
QString antpatternfilepath = getDebugDataPath(antpatternfilename);
gdalImage ds= CreategdalImage(antpatternfilepath, thetanum, phinum, 1, gt, "", true, true, true);
ds.saveImage(antPatternMatrix, 0, 0, 1);
}
void testOutClsArr(QString filename, long* amp, long rowcount, long colcount) {
Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
for (long hii = 0; hii < rowcount; hii++) {
for (long hjj = 0; hjj < colcount; hjj++) {
h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
}
}
QString ampPath = getDebugDataPath(filename);
saveEigenMatrixXd2Bin(h_amp_img, ampPath);
std::cout << filename.toLocal8Bit().constData() << std::endl;
std::cout << "max:\t" << h_amp_img.maxCoeff() << std::endl;
std::cout << "min:\t" << h_amp_img.minCoeff() << std::endl;
}
void testOutAmpArr(QString filename, float* amp, long rowcount, long colcount)
{
Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
for (long hii = 0; hii < rowcount; hii++) {
for (long hjj = 0; hjj < colcount; hjj++) {
h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
}
}
QString ampPath = getDebugDataPath(filename);
saveEigenMatrixXd2Bin(h_amp_img, ampPath);
std::cout << filename.toLocal8Bit().constData() << std::endl;
std::cout << "max:\t" << h_amp_img.maxCoeff() << std::endl;
std::cout << "min:\t" << h_amp_img.minCoeff() << std::endl;
}

10
SimulationSAR/RTPCProcessCls.h
View File

@ -69,8 +69,8 @@ public:
private: // 处理流程
ErrorCode InitParams();// 1. 初始化参数
ErrorCode DEMPreprocess(); // 2. 裁剪DEM范围
ErrorCode RTPCMainProcess(long num_thread);
//ErrorCode RTPCMainProcess(long num_thread);
ErrorCode RTPCMainProcess_GPU();
private:
QString demxyzPath;
QString demmaskPath;
@ -79,5 +79,7 @@ private:
void RTPCProcessMain(long num_thread,QString TansformPatternFilePath,QString ReceivePatternFilePath,QString simulationtaskName, QString OutEchoPath, QString GPSXmlPath,QString TaskXmlPath,QString demTiffPath, QString LandCoverPath, QString HHSigmaPath, QString HVSigmaPath, QString VHSigmaPath, QString VVSigmaPath);
// ²âÊÔ
void testOutAntPatternTrans(QString antpatternfilename,float* antPatternArr,double starttheta,double deltetheta,double startphi,double deltaphi,long thetanum,long phinum );
void testOutAmpArr(QString filename, float* amp, long rowcount, long colcount);
void testOutClsArr(QString filename, long* amp, long rowcount, long colcount);

781
SimulationSAR/SARSatelliteSimulationAbstractCls.cpp
View File

File diff suppressed because it is too large Load Diff

6
SimulationSAR/SARSatelliteSimulationAbstractCls.h
View File

@ -101,19 +101,19 @@ struct RadiationPatternGainPoint {
/// <summary>
/// ÌìÏß·½ÏòͼµÄ»ñÈ¡
/// 注意这里存在一定的插值方法
/// 注意这里使用 双线性插值
/// </summary>
class AbstractRadiationPattern {
public:
AbstractRadiationPattern();
virtual ~AbstractRadiationPattern();
public:
virtual double getGain(double& theta, double& phi);
virtual double getGain(double theta, double phi);
virtual std::vector<RadiationPatternGainPoint> getGainList();
virtual ErrorCode getGain(double& theta, double& phi, double& GainValue);
virtual ErrorCode getGainLinear(double& theta, double& phi, double& GainValue);
double getGainLearThetaPhi(double theta, double phi);
virtual ErrorCode setGain(double& theta, double& phi, double& GainValue);
virtual ErrorCode setGain(double theta, double phi, double GainValue);
virtual ErrorCode RecontructGainMatrix(double threshold=-3);
virtual std::vector<double> getThetas();
virtual std::vector<double> getPhis();

27
SimulationSAR/SatelliteOribtModel.cpp
View File

@ -230,15 +230,32 @@ ErrorCode PolyfitSatelliteOribtModel::getAntnnaDirection(SatelliteOribtNode& nod
double nexttime = node.time + 1e-6;
SatelliteOribtNode node1 = this->getSatelliteOribtNode(nexttime, flag);
//std::cout << "getAntnnaDirection corrdination " << std::endl;
double Vx = (node1.Px - node.Px);
double Vy = (node1.Py - node.Py);
double Vz = (node1.Pz - node.Pz);
// 代码测试部分
//node.Px = 0;
//node.Py = 0;
//node.Pz = 1;
//Vx = 1, Vy = 0, Vz = 0;
// 1. 计算天线指向
Eigen::Vector3d axisZ0 = { -1 * node.Px ,-1 * node.Py,-1 * node.Pz }; // z 轴 --波位角为0时，天线指向的反方向
Eigen::Vector3d axisX0 = { (node1.Px - node.Px) , (node1.Py - node.Py), (node1.Pz - node.Pz) }; // x 轴 --飞行方向
Eigen::Vector3d axisX0 = { Vx,Vy,Vz }; // x 轴 --飞行方向
Eigen::Vector3d axisY0 = axisZ0.cross(axisX0); // y 轴 --右手定则 -- 初始坐标系
//std::cout << "axis_X0=[ " << axisX0.x() << "," << axisX0.y() << "," << axisX0.z() << "]" << std::endl;
//std::cout << "axis_Y0=[ " << axisY0.x() << "," << axisY0.y() << "," << axisY0.z() << "]" << std::endl;
//std::cout << "axis_Z0=[ " << axisZ0.x() << "," << axisZ0.y() << "," << axisZ0.z() << "]" << std::endl;
double rotateAngle = this->RightLook ? -this->beamAngle : this->beamAngle; // 旋转角度 左（逆时针）：theta , 右（顺时针）： -theta
//std::cout << "rotateAngle=" << rotateAngle << std::endl;
//std::cout << "Look side:\t" << (this->RightLook ? "right" : "left") << std::endl;
// 1.2. 根据波位角，确定卫星绕X轴-飞行轴
Eigen::Matrix3d rotateMatrixBeam = rotationMatrix(axisX0, rotateAngle*d2r); // 旋转矩阵
axisZ0=rotateMatrixBeam*axisZ0; // 旋转矩阵
@ -271,6 +288,10 @@ ErrorCode PolyfitSatelliteOribtModel::getAntnnaDirection(SatelliteOribtNode& nod
node.AntZaxisZ =axisZ0[2];
//std::cout << "axis_X=[" << axisX0.x() << "," << axisX0.y() << "," << axisX0.z() << "]" << std::endl;
//std::cout << "axis_Y=[" << axisY0.x() << "," << axisY0.y() << "," << axisY0.z() << "]" << std::endl;
//std::cout << "axis_Z=[" << axisZ0.x() << "," << axisZ0.y() << "," << axisZ0.z() << "]" << std::endl;
//std::cout << "------------------------------------" << std::endl;
return ErrorCode::SUCCESS;
}

599
SimulationSAR/TBPImageAlgCls.cpp
View File

@ -5,8 +5,80 @@
#include <QString>
#include <cmath>
#include <QProgressDialog>
#include <QMessageBox>
#include "GPUTool.cuh"
void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZPath)
{
// 创建坐标系统
long prfcount = echoL0ds->getPluseCount();
long freqcount = echoL0ds->getPlusePoints();
Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(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 xyzRaster = CreategdalImage(outPixelXYZPath, prfcount, freqcount, 3, gt, QString(""), false, true);
std::shared_ptr<double> antpos = echoL0ds->getAntPos();
double dx = LIGHTSPEED / 2 / echoL0ds->getFs();
double Rnear = echoL0ds->getNearRange();
long echocol = 1073741824 / 8 / 4 / prfcount*4;
std::cout << "echocol:\t " << echocol << std::endl;
echocol = echocol < 3000 ? 3000 : echocol;
long startcolidx = 0;
for (startcolidx = 0; startcolidx < freqcount; startcolidx = startcolidx + echocol) {
long tempechocol = echocol;
if (startcolidx + tempechocol >= freqcount) {
tempechocol = freqcount - startcolidx;
}
std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "] imgxyz :\t" << startcolidx << "\t-\t" << startcolidx + tempechocol << " / " << freqcount << std::endl;
Eigen::MatrixXd demx = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 1);
Eigen::MatrixXd demy = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 2);
Eigen::MatrixXd demz = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 3);
for (long i = 0; i < prfcount; i++) {
double Px = 0;
double Py = 0;
double Pz = 0;
double AntDirectX = 0;
double AntDirectY = 0;
double AntDirectZ = 0;
double R = 0;
double NormAnt = 0;
Px = antpos.get()[i * 19 + 1];
Py = antpos.get()[i * 19 + 2];
Pz = antpos.get()[i * 19 + 3];
AntDirectX = antpos.get()[i * 19 + 13];// zero doppler
AntDirectY = antpos.get()[i * 19 + 14];
AntDirectZ = antpos.get()[i * 19 + 15];
NormAnt = std::sqrt(AntDirectX * AntDirectX + AntDirectY * AntDirectY + AntDirectZ * AntDirectZ);
AntDirectX = AntDirectX / NormAnt;
AntDirectY = AntDirectY / NormAnt;
AntDirectZ = AntDirectZ / NormAnt;// 归一化
for (long j = 0; j < tempechocol; j++) {
R = (j + startcolidx)*dx + Rnear;
demx(i,j) = Px + AntDirectX * R;
demy(i,j) = Py + AntDirectY * R;
demz(i,j) = Pz + AntDirectZ * R;
}
}
xyzRaster.saveImage(demx, 0, startcolidx, 1);
xyzRaster.saveImage(demy, 0, startcolidx, 2);
xyzRaster.saveImage(demz, 0, startcolidx, 3);
}
}
void TBPImageProcess(QString echofile, QString outImageFolder, QString imagePlanePath,long num_thread)
{
@ -32,6 +104,8 @@ void TBPImageProcess(QString echofile, QString outImageFolder, QString imagePlan
TBPimag.Process(num_thread);
}
void TBPImageAlgCls::setImagePlanePath(QString INimagePlanePath)
{
this->imagePlanePath = INimagePlanePath;
@ -64,14 +138,333 @@ std::shared_ptr<SARSimulationImageL1Dataset> TBPImageAlgCls::getImageL0()
ErrorCode TBPImageAlgCls::Process(long num_thread)
{
qDebug() << u8"创建成像平面的XYZ";
QString outRasterXYZ = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_xyz.tif");
CreatePixelXYZ(this->L0ds, outRasterXYZ);
this->outRasterXYZPath = outRasterXYZ;
// 初始化Raster
qDebug() << u8"初始化影像";
long imageheight = this->L1ds->getrowCount();
long imagewidth = this->L1ds->getcolCount();
long blokline = Memory1GB / 8 / 4 / imageheight * 8;
blokline = blokline < 1000 ? 1000 : blokline;
long startline = 0;
for (startline = 0; startline < imageheight; startline = startline + blokline) {
long templine = blokline;
if (startline + templine >= imageheight) {
templine = imageheight - startline;
}
std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "] imgxyz :\t" << startline << "\t-\t" << startline + templine << " / " << imageheight << std::endl;
std::shared_ptr<std::complex<double>> imageRaster = this->L1ds->getImageRaster(startline, templine);
for (long i = 0; i < templine; i++) {
for (long j = 0; j < imagewidth; j++) {
imageRaster.get()[i * imagewidth + j] = std::complex<double>(0,0);
}
}
this->L1ds->saveImageRaster(imageRaster, startline,templine);
}
qDebug() << u8"开始成像";
if (GPURUN) {
return this->ProcessGPU();
}
else {
return this->ProcessCPU(num_thread);
QMessageBox::information(nullptr,u8"提示",u8"目前只支持显卡");
return ErrorCode::FAIL;
}
}
ErrorCode TBPImageAlgCls::ProcessGPU()
{
// 常用参数
long rowCount = this->L1ds->getrowCount();
long colCount = this->L1ds->getcolCount();
long pixelCount = rowCount * colCount;
long PRFCount = this->L0ds->getPluseCount();
long PlusePoints = this->L0ds->getPlusePoints();
float Rnear = this->L1ds->getNearRange();
float Rfar = this->L1ds->getFarRange();
float fs = this->L1ds->getFs();
double dx = LIGHTSPEED / 2 / fs;
float freq = this->L1ds->getCenterFreq()*1.0*1e9;
double factorj = freq * 4 * M_PI / LIGHTSPEED ;
std::shared_ptr<float> Pxs (new float[this->L0ds->getPluseCount()]);
std::shared_ptr<float> Pys (new float[this->L0ds->getPluseCount()]);
std::shared_ptr<float> Pzs (new float[this->L0ds->getPluseCount()]);
{
std::shared_ptr<double> antpos = this->L0ds->getAntPos();
double time = 0;
double Px = 0;
double Py = 0;
double Pz = 0;
for (long i = 0; i < rowCount; i++) {
time = antpos.get()[i *19 + 0];
Px = antpos.get()[i *19 + 1];
Py = antpos.get()[i *19 + 2];
Pz = antpos.get()[i *19 + 3];
Pxs.get()[i] = Px;
Pys.get()[i] = Py;
Pzs.get()[i] = Pz;
}
antpos.reset();
}
// 按照回波分块，图像分块
long echoBlockline = Memory1GB / 8 / 2 / PlusePoints * 6;
echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;
long imageBlockline = Memory1GB / 8 / 2 / colCount * 2;
imageBlockline = imageBlockline < 1 ? 1 : imageBlockline;
gdalImage imageXYZ(this->outRasterXYZPath);
long startimgrowid = 0;
for (startimgrowid = 0; startimgrowid < rowCount; startimgrowid = startimgrowid + imageBlockline) {
long tempimgBlockline = imageBlockline;
if (startimgrowid + imageBlockline >= rowCount) {
tempimgBlockline = rowCount - startimgrowid;
}
std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "] dem:\t" << startimgrowid << "\t-\t" << startimgrowid + tempimgBlockline << std::endl;
// 提取局部pixel x,y,z
std::shared_ptr<float> img_x = readDataArr<float>(imageXYZ,startimgrowid,0,tempimgBlockline,colCount,1,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<float> img_y = readDataArr<float>(imageXYZ,startimgrowid,0,tempimgBlockline,colCount,2,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<float> img_z = readDataArr<float>(imageXYZ,startimgrowid,0,tempimgBlockline,colCount,3,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<std::complex<double>> imgArr = this->L1ds->getImageRaster(startimgrowid, tempimgBlockline);
// 获取回波
long startechoid = 0;
for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) {
long tempechoBlockline = echoBlockline;
if (startechoid + tempechoBlockline >= PRFCount) {
tempechoBlockline = PRFCount - startechoid;
}
std::shared_ptr<std::complex<double>> echoArr = this->L0ds->getEchoArr(startechoid, tempechoBlockline);
std::shared_ptr<float> antpx(new float[tempechoBlockline*PlusePoints]);
std::shared_ptr<float> antpy(new float[tempechoBlockline* PlusePoints]);
std::shared_ptr<float> antpz(new float[tempechoBlockline* PlusePoints]);
// 复制
for (long anti = 0; anti < tempechoBlockline; anti++) {
antpx.get()[anti] = Pxs.get()[anti + startechoid];
antpy.get()[anti] = Pys.get()[anti + startechoid];
antpz.get()[anti] = Pzs.get()[anti + startechoid];
}
TBPImageGPUAlg(antpx, antpy, antpz,
img_x, img_y, img_z,
echoArr, imgArr,
freq, fs, Rnear, Rfar,
tempimgBlockline, colCount,
tempechoBlockline, PlusePoints );
}
this->L1ds->saveImageRaster(imgArr, startimgrowid, tempimgBlockline);
}
qDebug() << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz") << "] image writing:\t" << this->L1ds->getxmlFilePath();
return ErrorCode::SUCCESS;
}
void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz,
std::shared_ptr<float> imgx, std::shared_ptr<float> imgy, std::shared_ptr<float> imgz,
std::shared_ptr<std::complex<double>> echoArr, std::shared_ptr<std::complex<double>> imgArr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount,
long prfcount, long freqcount
)
{
// 声明GPU变量
float* h_antPx = (float*)mallocCUDAHost(sizeof(float) * prfcount);
float* h_antPy = (float*)mallocCUDAHost(sizeof(float) * prfcount);
float* h_antPz = (float*)mallocCUDAHost(sizeof(float) * prfcount);
float* d_antPx = (float*)mallocCUDADevice(sizeof(float) * prfcount);
float* d_antPy = (float*)mallocCUDADevice(sizeof(float) * prfcount);
float* d_antPz = (float*)mallocCUDADevice(sizeof(float) * prfcount);
float* h_imgx = (float*)mallocCUDAHost(sizeof(float) * rowcount * colcount);
float* h_imgy = (float*)mallocCUDAHost(sizeof(float) * rowcount * colcount);
float* h_imgz = (float*)mallocCUDAHost(sizeof(float) * rowcount * colcount);
float* d_imgx = (float*)mallocCUDADevice(sizeof(float) * rowcount * colcount);
float* d_imgy = (float*)mallocCUDADevice(sizeof(float) * rowcount * colcount);
float* d_imgz = (float*)mallocCUDADevice(sizeof(float) * rowcount * colcount);
cuComplex* h_echoArr = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * prfcount * freqcount);
cuComplex* h_imgArr = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * rowcount * colcount);
cuComplex* d_echoArr = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * prfcount * freqcount);
cuComplex* d_imgArr = (cuComplex*)mallocCUDADevice( sizeof(cuComplex) * rowcount * colcount);
// 初始化
// 天线位置
for (long i = 0; i < prfcount; i++) {
h_antPx[i] = antPx.get()[i];
h_antPy[i] = antPy.get()[i];
h_antPz[i] = antPz.get()[i];
}
// 成像平面
for (long i = 0; i < rowcount; i++) {
for (long j = 0; j < colcount; j++) {
h_imgx[i * colcount + j]=imgx.get()[i * colcount + j];
h_imgy[i * colcount + j]=imgy.get()[i * colcount + j];
h_imgz[i * colcount + j]=imgz.get()[i * colcount + j];
}
}
// 回波
for (long i = 0; i < prfcount; i++) {
for (long j = 0; j < freqcount; j++) {
h_echoArr[i * freqcount + j] = make_cuComplex(echoArr.get()[i * freqcount + j].real(),
echoArr.get()[i * freqcount + j].imag());
}
}
// 图像
for (long i = 0; i < rowcount; i++) {
for (long j = 0; j < colcount; j++) {
h_imgArr[i * colcount + j].x = imgArr.get()[i * colcount + j].real();
h_imgArr[i * colcount + j].y = imgArr.get()[i * colcount + j].imag();
}
}
// Host -> GPU
HostToDevice(h_antPx, d_antPx, sizeof(float) * prfcount);
HostToDevice(h_antPy, d_antPy, sizeof(float) * prfcount);
HostToDevice(h_antPz, d_antPz, sizeof(float) * prfcount);
HostToDevice(h_imgx, d_imgx, sizeof(float) * rowcount * colcount);
HostToDevice(h_imgy, d_imgy, sizeof(float) * rowcount * colcount);
HostToDevice(h_imgz, d_imgz, sizeof(float) * rowcount * colcount);
HostToDevice(h_echoArr, d_echoArr, sizeof(cuComplex) * prfcount * freqcount);
HostToDevice(h_imgArr, d_imgArr, sizeof(cuComplex) * rowcount * colcount);
for (long prfid = 0; prfid < prfcount; prfid++) {
CUDATBPImage(
d_antPx,
d_antPy,
d_antPz,
d_imgx,
d_imgy,
d_imgz,
d_echoArr,
d_imgArr,
freq, fs, Rnear, Rfar,
rowcount, colcount,
prfid, freqcount
);
}
// Device -> Host
DeviceToHost(h_imgArr, d_imgArr, sizeof(cuComplex) * rowcount * colcount);
for (long i = 0; i < rowcount; i++) {
for (long j = 0; j < colcount; j++) {
imgArr.get()[i * colcount + j] = std::complex<double>(h_imgArr[i * colcount + j].x,
h_imgArr[i * colcount + j].y);
}
}
FreeCUDAHost(h_antPx);
FreeCUDAHost(h_antPy);
FreeCUDAHost(h_antPz);
FreeCUDADevice(d_antPx);
FreeCUDADevice(d_antPy);
FreeCUDADevice(d_antPz);
FreeCUDAHost(h_imgx);
FreeCUDAHost(h_imgy);
FreeCUDAHost(h_imgz);
FreeCUDADevice(d_imgx);
FreeCUDADevice(d_imgy);
FreeCUDADevice(d_imgz);
FreeCUDAHost(h_echoArr);
FreeCUDAHost(h_imgArr);
FreeCUDADevice(d_echoArr);
FreeCUDADevice(d_imgArr);
// 释放GPU变量
}
void TBPImageAlgCls::setGPU(bool flag)
{
this->GPURUN = flag;
}
bool TBPImageAlgCls::getGPU( )
{
return this->GPURUN;
}
/// <summary>
/// TBP GPU代码
/// </summary>
/// <param name="antpos_ptr">卫星轨道坐标</param>
/// <param name="echoArr">回波矩阵</param>
/// <param name="img_arr">图像矩阵</param>
void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz, // 天线坐标
std::shared_ptr<float> antVx, std::shared_ptr<float> antVy, std::shared_ptr<float> antVz,
std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z, // 图像坐标
std::shared_ptr<std::complex<double>> echoArr, std::shared_ptr<std::complex<double>> img_arr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount, std::shared_ptr<SARSimulationImageL1Dataset> L1ds) {
float factorj = freq * 4 * PI / LIGHTSPEED;
qDebug() << "factorj:\t" << factorj;
qDebug() << "freq:\t" << freq;
qDebug() << "fs:\t" << fs;
qDebug() << "Rnear:\t" << Rnear;
qDebug() << "Rfar:\t" << Rfar;
qDebug() << "img_x:\t" << img_x.get()[0];
qDebug() << "img_y:\t" << img_y.get()[0];
qDebug() << "img_z:\t" << img_z.get()[0];
long blockline = Memory1MB * 1000 / sizeof(float) / colcount;
blockline = blockline < 10 ? 10 : blockline;
for (long startline = 0; startline < rowcount; startline = startline + blockline) {
long stepline = startline + blockline < rowcount ? blockline : rowcount - startline;
std::cout << startline << " \ " << rowcount << " "<< stepline << " start " << std::endl;
//TBPImageGPUBlock(antPx.get(), antPy.get(), antPz.get(), img_x.get(), img_y.get(), img_z.get(),
// echoArr, rowcount, colcount,
// img_arr,
// freq, fs, Rnear, Rfar, factorj, startline, stepline,
// stepline, colcount);
//std::cout << startline << " \ " << rowcount << " " << stepline << " end " << std::endl;
//L1ds->saveImageRaster(img_arr, 0, rowcount);
}
L1ds->saveImageRaster(img_arr, 0, rowcount);
L1ds->saveToXml();
}
/**
ErrorCode TBPImageAlgCls::ProcessCPU(long num_thread)
{
@ -131,28 +524,28 @@ ErrorCode TBPImageAlgCls::ProcessCPU(long num_thread)
double NormAnt = 0;
for (long i = 0; i < rowCount; i++) {
time = antpos.get()[i *19 + 0];
Px = antpos.get()[i *19 + 1];
Py = antpos.get()[i *19 + 2];
Pz = antpos.get()[i *19 + 3];
Vx = antpos.get()[i *19 + 4];
Vy = antpos.get()[i *19 + 5];
Vz = antpos.get()[i *19 + 6];
AntDirectX = antpos.get()[i *19 + 13]; // Zero doppler
AntDirectY = antpos.get()[i *19 + 14];
AntDirectZ = antpos.get()[i *19 + 15];
AVx = antpos.get()[i *19 + 10];
AVy = antpos.get()[i *19 + 11];
AVz = antpos.get()[i *19 + 12];
time = antpos.get()[i * 19 + 0];
Px = antpos.get()[i * 19 + 1];
Py = antpos.get()[i * 19 + 2];
Pz = antpos.get()[i * 19 + 3];
Vx = antpos.get()[i * 19 + 4];
Vy = antpos.get()[i * 19 + 5];
Vz = antpos.get()[i * 19 + 6];
AntDirectX = antpos.get()[i * 19 + 13]; // Zero doppler
AntDirectY = antpos.get()[i * 19 + 14];
AntDirectZ = antpos.get()[i * 19 + 15];
AVx = antpos.get()[i * 19 + 10];
AVy = antpos.get()[i * 19 + 11];
AVz = antpos.get()[i * 19 + 12];
NormAnt = std::sqrt(AntDirectX * AntDirectX + AntDirectY * AntDirectY + AntDirectZ * AntDirectZ);
AntDirectX = AntDirectX / NormAnt;
AntDirectY = AntDirectY / NormAnt;
AntDirectZ = AntDirectZ / NormAnt;// ¹éÒ»»¯
antpos.get()[i *19 + 13] = AntDirectX;
antpos.get()[i *19 + 14] = AntDirectY;
antpos.get()[i *19 + 15] = AntDirectZ;
antpos.get()[i * 19 + 13] = AntDirectX;
antpos.get()[i * 19 + 14] = AntDirectY;
antpos.get()[i * 19 + 15] = AntDirectZ;
Pxs(i, 0) = Px;
Pys(i, 0) = Py;
Pzs(i, 0) = Pz;
@ -252,175 +645,5 @@ ErrorCode TBPImageAlgCls::ProcessCPU(long num_thread)
return ErrorCode::SUCCESS;
}
ErrorCode TBPImageAlgCls::ProcessGPU()
{
// 常用参数
long rowCount = this->L1ds->getrowCount();
long colCount = this->L1ds->getcolCount();
long pixelCount = rowCount * colCount;
long PRFCount = this->L0ds->getPluseCount();
long PlusePoints = this->L0ds->getPlusePoints();
float Rnear = this->L1ds->getNearRange();
float Rfar = this->L1ds->getFarRange();
float fs = this->L1ds->getFs();
double dx = LIGHTSPEED / 2 / fs;
float freq = this->L1ds->getCenterFreq()*1.0*1e9;
double factorj = freq * 4 * M_PI / LIGHTSPEED ;
std::shared_ptr<float> pixelX(new float[rowCount*colCount],delArrPtr); // 图像成像网格
std::shared_ptr<float> pixelY(new float[rowCount*colCount],delArrPtr);
std::shared_ptr<float> pixelZ(new float[rowCount*colCount],delArrPtr);
std::shared_ptr<float> Pxs (new float[this->L0ds->getPluseCount()]);
std::shared_ptr<float> Pys (new float[this->L0ds->getPluseCount()]);
std::shared_ptr<float> Pzs (new float[this->L0ds->getPluseCount()]);
std::shared_ptr<float> Vxs (new float[this->L0ds->getPluseCount()]);
std::shared_ptr<float> Vys (new float[this->L0ds->getPluseCount()]);
std::shared_ptr<float> Vzs (new float[this->L0ds->getPluseCount()]);
// 图像网格坐标
{
std::shared_ptr<double> antpos = this->L0ds->getAntPos();
double time = 0;
double Px = 0;
double Py = 0;
double Pz = 0;
double Vx = 0;
double Vy = 0;
double Vz = 0;
double AntDirectX = 0;
double AntDirectY = 0;
double AntDirectZ = 0;
double AVx = 0;
double AVy = 0;
double AVz = 0;
double R = 0;
double NormAnt = 0;
for (long i = 0; i < rowCount; i++) {
time = antpos.get()[i *19 + 0];
Px = antpos.get()[i *19 + 1];
Py = antpos.get()[i *19 + 2];
Pz = antpos.get()[i *19 + 3];
Vx = antpos.get()[i *19 + 4];
Vy = antpos.get()[i *19 + 5];
Vz = antpos.get()[i *19 + 6];
AntDirectX = antpos.get()[i *19 + 13];// zero doppler
AntDirectY = antpos.get()[i *19 + 14];
AntDirectZ = antpos.get()[i *19 + 15];
AVx = antpos.get()[i *19 + 10];
AVy = antpos.get()[i *19 + 11];
AVz = antpos.get()[i *19 + 12];
NormAnt = std::sqrt(AntDirectX * AntDirectX + AntDirectY * AntDirectY + AntDirectZ * AntDirectZ);
AntDirectX = AntDirectX / NormAnt;
AntDirectY = AntDirectY / NormAnt;
AntDirectZ = AntDirectZ / NormAnt;// 归一化
antpos.get()[i *19 + 13] = AntDirectX;
antpos.get()[i *19 + 14] = AntDirectY;
antpos.get()[i *19 + 15] = AntDirectZ;
Pxs.get()[i] = Px;
Pys.get()[i] = Py;
Pzs.get()[i] = Pz;
Vxs.get()[i] = Vx;
Vys.get()[i] = Vy;
Vzs.get()[i] = Vz;
for (long j = 0; j < colCount; j++) {
R = j * dx + Rnear;
pixelX.get()[i*colCount+ j] = Px + AntDirectX * R;
pixelY.get()[i*colCount+ j] = Py + AntDirectY * R;
pixelZ.get()[i*colCount+ j] = Pz + AntDirectZ * R;
}
}
qDebug()<<"R: " << R;
this->L1ds->saveAntPos(antpos);
antpos.reset();
}
std::shared_ptr<std::complex<double>> Rasterarr = this->L1ds->getImageRaster();
std::shared_ptr<std::complex<double>> echodataPtr = this->L0ds->getEchoArr();
for (long i = 0; i < rowCount; i++) {
for (long j = 0; j < colCount; j++) {
Rasterarr.get()[i * colCount + j] = Rasterarr.get()[i * colCount + j] * 0.0;
}
}
TBPImageGPUAlg(Pxs, Pys, Pzs, // 天线坐标
Vxs, Vys, Vzs,
pixelX, pixelY, pixelZ, // 图像坐标
echodataPtr, Rasterarr,
freq, fs, Rnear, Rfar,
rowCount, colCount, this->L1ds);
qDebug() << "image writing:\t" << this->L1ds->getxmlFilePath();
return ErrorCode::SUCCESS;
}
void TBPImageAlgCls::setGPU(bool flag)
{
this->GPURUN = flag;
}
bool TBPImageAlgCls::getGPU( )
{
return this->GPURUN;
}
/// <summary>
/// TBP GPU代码
/// </summary>
/// <param name="antpos_ptr">卫星轨道坐标</param>
/// <param name="echoArr">回波矩阵</param>
/// <param name="img_arr">图像矩阵</param>
void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz, // 天线坐标
std::shared_ptr<float> antVx, std::shared_ptr<float> antVy, std::shared_ptr<float> antVz,
std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z, // 图像坐标
std::shared_ptr<std::complex<double>> echoArr, std::shared_ptr<std::complex<double>> img_arr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount, std::shared_ptr<SARSimulationImageL1Dataset> L1ds) {
float factorj = freq * 4 * PI / LIGHTSPEED;
qDebug() << "factorj:\t" << factorj;
qDebug() << "freq:\t" << freq;
qDebug() << "fs:\t" << fs;
qDebug() << "Rnear:\t" << Rnear;
qDebug() << "Rfar:\t" << Rfar;
qDebug() << "img_x:\t" << img_x.get()[0];
qDebug() << "img_y:\t" << img_y.get()[0];
qDebug() << "img_z:\t" << img_z.get()[0];
long blockline = Memory1MB * 1000 / sizeof(float) / colcount;
blockline = blockline < 10 ? 10 : blockline;
for (long startline = 0; startline < rowcount; startline = startline + blockline) {
long stepline = startline + blockline < rowcount ? blockline : rowcount - startline;
std::cout << startline << " \ " << rowcount << " "<< stepline << " start " << std::endl;
//TBPImageGPUBlock(antPx.get(), antPy.get(), antPz.get(), img_x.get(), img_y.get(), img_z.get(),
// echoArr, rowcount, colcount,
// img_arr,
// freq, fs, Rnear, Rfar, factorj, startline, stepline,
// stepline, colcount);
//std::cout << startline << " \ " << rowcount << " " << stepline << " end " << std::endl;
//L1ds->saveImageRaster(img_arr, 0, rowcount);
}
L1ds->saveImageRaster(img_arr, 0, rowcount);
L1ds->saveToXml();
}
*/

17
SimulationSAR/TBPImageAlgCls.h
View File

@ -37,6 +37,7 @@ private:
std::shared_ptr < EchoL0Dataset> L0ds;
QString imagePlanePath;
bool GPURUN;
QString outRasterXYZPath;
public:
void setImagePlanePath(QString imagePlanePath);
@ -53,18 +54,20 @@ public:
void setGPU(bool flag);
bool getGPU( );
private:
ErrorCode ProcessCPU(long num_thread);
//ErrorCode ProcessCPU(long num_thread);
ErrorCode ProcessGPU();
};
void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds,QString outPixelXYZPath);
void TBPImageProcess(QString echofile,QString outImageFolder,QString imagePlanePath,long num_thread);
void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz, // 天线坐标
std::shared_ptr<float> antVx, std::shared_ptr<float> antVy, std::shared_ptr<float> antVz,
std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z, // 图像坐标
void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz,
std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z,
std::shared_ptr<std::complex<double>> echoArr, std::shared_ptr<std::complex<double>> img_arr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount, std::shared_ptr<SARSimulationImageL1Dataset> L1ds);
long rowcount, long colcount,
long prfcount,long freqcount );

2
bugreflex.md Normal file
View File

@ -0,0 +1,2 @@
1. 2024.12.19 bug发现
理论上amp应该都是0，但是从转换的结果来看都是1，这是错误的

2
script/code分析/mainAnayls.m
View File

@ -12,7 +12,7 @@ clear,clc,close all
% 打开二进制文件
fileID = fopen('D:\\Programme\\vs2022\\RasterMergeTest\\TestData\\outData\\GF3_Simulation.gpspos.data', 'rb'); % 假设二进制文件名为 'data.bin'
fileID = fopen('E:\\LAMPCAE_SCANE\\outTestEcho\\GF3_Simulation.gpspos.data', 'rb'); % 假设二进制文件名为 'data.bin'
% 定义每个数据字段的数据类型
% 假设每个数据是双精度浮动数（8字节）