Merge pull request 'RTPC-BPImage_dev' (#8) from RTPC-BPImage_dev into RTPC-dev

Reviewed-on: http://123.153.4.249:22779/LAMPSARToolSoftware/RasterProcessTool/pulls/8
pull/9/head
chenzenghui 2025-03-07 09:16:32 +08:00
commit 338fbb869d
16 changed files with 1206 additions and 55 deletions

View File

@ -148,6 +148,11 @@ struct Vector3 {
};
struct Vector3_single {
float x, y, z;
};
/*********************************************** FEKO仿真参数 ********************************************************************/
struct SatellitePos {
double time;

View File

@ -21,9 +21,9 @@
#include "GPUBPTool.cuh"
#include "BPBasic0_CUDA.cuh"
#include "GPUBpSimulation.cuh"
#include "GPURFPC_single.cuh"
#include "GPURFPC.cuh"
double* getFreqPoints_mallocHost(double startFreq, double endFreq, long freqpoints)
{
@ -68,7 +68,7 @@ __global__ void kernel_RFPCProcess(
// ÈëÉä½Ç
Vector3 TS = vec_sub(S, T);//T-->S
double localIncAngle = angleBetweenVectors(TS, slp, false);// ÈëÉä½Ç
double sigma0 = GPU_getSigma0dB(p1, p2, p3, p4, p5, p6, localIncAngle);// 后向散射系数
double sigma0 = GPU_getSigma0dB_params(p1, p2, p3, p4, p5, p6, localIncAngle);// 后向散射系数
sigma0 = powf(10.0, sigma0 / 10.0);
// ¾àÀë
@ -95,8 +95,6 @@ void RFPCProcess(double Tx, double Ty, double Tz,
double p1, double p2, double p3, double p4, double p5, double p6,
GPUDATA& d_data)
{
double* AntX = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
double* AntY = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
double* AntZ = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
@ -128,3 +126,125 @@ void RFPCProcess(double Tx, double Ty, double Tz,
/** 单精度成像测试 *****************************************************************************************************************************************/
__global__ void kernel_RFPCProcess_single(
float* Sx, float* Sy, float* Sz,
float Tx, float Ty, float Tz,
float Tslx, float Tsly, float Tslz, // 目标的坡面向量
float p1, float p2, float p3, float p4, float p5, float p6,
long Np, long Nf,
float minF, float dFreq, float RefRange,
cuComplex* phdata
) {
long prfid = blockIdx.x * blockDim.x + threadIdx.x; // 获取当前的线程编码
if (prfid >= Np || prfid < 0) { return; }
else {}
// 距离
Vector3_single S{ Sx[prfid],Sy[prfid],Sz[prfid] };
Vector3_single T{ Tx,Ty,Tz };
Vector3_single slp{ Tslx,Tsly,Tslz };
//printf("S=(%e,%e,%e)\n", S.x, S.y, S.z);
// 入射角
Vector3_single TS = vec_sub_single(S, T);//T-->S
float localIncAngle = angleBetweenVectors_single(TS, slp, false);// 入射角
float sigma0 = GPU_getSigma0dB_single(p1, p2, p3, p4, p5, p6, localIncAngle);// 后向散射系数
sigma0 = powf(10.0, sigma0 / 10.0);
// 距离
float R = sqrtf(vec_dot_single(TS, TS));
// 计算增益
float amp_echo = sigma0 / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4)); // 反射强度
float phi = 0;
for (long fid = 0; fid < Nf; fid++) {
phi = -4.0 * PI / LIGHTSPEED * (minF + dFreq * fid) * (R - RefRange);
phdata[prfid * Nf + fid].x += amp_echo * cosf(phi);
phdata[prfid * Nf + fid].y += amp_echo * sinf(phi);
//printf("phdata(%d,%d)=complex(%e,%e)\n", prfid, fid, phdata[prfid * Nf + fid].x, phdata[prfid * Nf + fid].y);
}
//printf("Nf:%d\n", Nf);
//printf("amp_echo:%f\n", amp_echo);
//printf("sigma0:%f\n", sigma0);
//printf("R:%e\n", R);
}
float* getFreqPoints_mallocHost_single(float startFreq, float endFreq, long freqpoints)
{
float dfreq = (endFreq - startFreq) / (freqpoints - 1);
float* freqlist = (float*)mallocCUDAHost(sizeof(float) * freqpoints);
for (long i = 0; i < freqpoints; i++) {
freqlist[i] = startFreq + dfreq * i;
}
return freqlist;
}
void RFPCProcess_single(
float Tx, float Ty, float Tz,
float Tslx, float Tsly, float Tslz,
float p1, float p2, float p3, float p4, float p5, float p6,
GPUDATA_single& d_data)
{
float* AntX = (float*)mallocCUDADevice(sizeof(float) * d_data.Np);
float* AntY = (float*)mallocCUDADevice(sizeof(float) * d_data.Np);
float* AntZ = (float*)mallocCUDADevice(sizeof(float) * d_data.Np);
HostToDevice(d_data.AntX, AntX, sizeof(float) * d_data.Np); printf("antX host to device finished!!\n");
HostToDevice(d_data.AntY, AntY, sizeof(float) * d_data.Np); printf("antY host to device finished!!\n");
HostToDevice(d_data.AntZ, AntZ, sizeof(float) * d_data.Np); printf("antZ host to device finished!!\n");
float minF = d_data.minF[0];
long grid_size = (d_data.Np + BLOCK_SIZE - 1) / BLOCK_SIZE;
float dfreq = d_data.deltaF;
float R0 = d_data.R0;
kernel_RFPCProcess_single << <grid_size, BLOCK_SIZE >> > (
AntX, AntY, AntZ,
Tx, Ty, Tz,
Tslx, Tsly, Tslz, // 目标的坡面向量
p1, p2, p3, p4, p5, p6,
d_data.Np, d_data.Nfft,
minF, dfreq, R0,
d_data.phdata
);
PrintLasterError("RFPCProcess");
FreeCUDADevice(AntX);
FreeCUDADevice(AntY);
FreeCUDADevice(AntZ);
return;
}

View File

@ -30,6 +30,17 @@ extern "C" void RFPCProcess(
GPUDATA& d_data
);
extern "C" float* getFreqPoints_mallocHost_single(float startFreq, float endFreq, long freqpoints);
extern "C" void RFPCProcess_single(
float Tx, float Ty, float Tz, // 目标点坐标
float Tslx, float Tsly, float Tslz, // 目标的坡面向量
float p1, float p2, float p3, float p4, float p5, float p6,// 地面目标后向散射系数与入射角关系 系数
GPUDATA_single& d_data
);
#endif

View File

@ -13,19 +13,19 @@
#include "LookTableSimulationComputer.cuh"
extern __device__ __host__ double getNumberDopplerCenterRate(double R, double r0, double r1, double r2, double r3, double r4, double reftime)
extern __device__ double getNumberDopplerCenterRate(double R, double r0, double r1, double r2, double r3, double r4, double reftime)
{
double t=R / LIGHTSPEED - reftime;
double dopplerCenterRate = r0 + r1*t + r2*t*t + r3*t*t*t + r4*t*t*t*t;
return dopplerCenterRate;
}
__device__ __host__ double getDopplerCenterRate(double Rx, double Ry, double Rz, double Vx, double Vy, double Vz, double fact_lamda)
__device__ double getDopplerCenterRate(double Rx, double Ry, double Rz, double Vx, double Vy, double Vz, double fact_lamda)
{
return -2*(Rx*Vx+Ry*Vy+Rz*Vz)* fact_lamda;
}
__device__ __host__ double getPolyfitNumber(double x, double a0, double a1, double a2, double a3, double a4, double a5)
__device__ double getPolyfitNumber(double x, double a0, double a1, double a2, double a3, double a4, double a5)
{
return a0 + a1 * x + a2 * x * x + a3 * x * x * x + a4 * x * x * x * x + a5 * x * x * x * x * x;
}

View File

@ -8,11 +8,11 @@
extern __device__ __host__ double getNumberDopplerCenterRate(double R, double r0, double r1, double r2, double r3, double r4,double reftime);
extern __device__ double getNumberDopplerCenterRate(double R, double r0, double r1, double r2, double r3, double r4,double reftime);
extern __device__ __host__ double getDopplerCenterRate(double Rx,double Ry,double Rz,double Vx,double Vy,double Vz,double fact_lamda); //fact_lamda lamda ľÄľšĘý
extern __device__ double getDopplerCenterRate(double Rx,double Ry,double Rz,double Vx,double Vy,double Vz,double fact_lamda); //fact_lamda lamda ľÄľšĘý
extern __device__ __host__ double getPolyfitNumber(double x, double a0, double a1, double a2, double a3, double a4,double a5);
extern __device__ double getPolyfitNumber(double x, double a0, double a1, double a2, double a3, double a4,double a5);

View File

@ -185,9 +185,8 @@ void bpBasic0CUDA(GPUDATA& data, int flag,double* h_R) {
data.im_final
//,d_R
);
PrintLasterError("processPulseKernel");
if (ii % 1000==0) {
printfinfo("\rPRF(%f %) %d / %d\t\t\t\t",(ii*100.0/data.Np), ii,data.Np);
if (ii % 10000 == 0) {
printfinfo("PRF[%f]: %d / %d", ii * 100.0 / data.Np, ii, data.Np);
}
}
//FreeCUDADevice(d_R);
@ -276,6 +275,10 @@ void BPBasic0(GPUDATA& h_data)
freeGPUData(d_data);
}
//int main() {
// GPUDATA h_data, d_data;
//

View File

@ -30,6 +30,21 @@ struct GPUDATA {
double deltaF; // 频点范围
};
struct GPUDATA_single {
long Nfft, K, Np, nx, ny; // 傅里叶点数、频点数、脉冲数、图像列、图像行
float* AntX, * AntY, * AntZ, * minF; // 天线坐标、起始频率
float* x_mat, * y_mat, * z_mat;// 地面坐标
float* r_vec; // 坐标范围
float* Freq;// 频率
cuComplex* phdata;// 回波
cuComplex* im_final;// 图像
float R0; // 参考斜距
float deltaF; // 频点范围
};
extern "C" {
void bpBasic0CUDA(GPUDATA& data, int flag,double* h_R=nullptr);
void initGPUData(GPUDATA& h_data, GPUDATA& d_data);
@ -38,4 +53,15 @@ extern "C" {
void BPBasic0(GPUDATA& h_data);
};
extern "C" {
// 函数废弃,因为这个成像精度不够
//void bpBasic0CUDA_single(GPUDATA_single& data, int flag, float* h_R = nullptr);
//void initGPUData_single(GPUDATA_single& h_data, GPUDATA_single& d_data);
//void freeGPUData_single(GPUDATA_single& d_data);
//void freeHostData_single(GPUDATA_single& d_data);
//void BPBasic0_single(GPUDATA_single& h_data);
};
#endif

View File

@ -22,13 +22,13 @@
#define M_PI 3.14159265358979323846
#endif
__device__ __host__ double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees ) {
__device__ double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees ) {
// 计算点积
double dotProduct = a.x * b.x + a.y * b.y + a.z * b.z;
// 计算模长
double magA = std::sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
double magB = std::sqrt(b.x * b.x + b.y * b.y + b.z * b.z);
double magA = sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
double magB = sqrt(b.x * b.x + b.y * b.y + b.z * b.z);
// 处理零向量异常
if (magA == 0.0 || magB == 0.0) {
@ -37,33 +37,33 @@ __device__ __host__ double angleBetweenVectors(Vector3 a, Vector3 b, bool retur
double cosTheta = dotProduct / (magA * magB);
cosTheta = cosTheta < -1 ? -1 : cosTheta>1 ? 1 : cosTheta; // 截断到[-1, 1]
double angleRad = std::acos(cosTheta);
double angleRad = acos(cosTheta);
return returnDegrees ? angleRad * 180.0 / M_PI : angleRad;
}
// 向量运算
__device__ __host__ Vector3 vec_sub(Vector3 a, Vector3 b) {
__device__ Vector3 vec_sub(Vector3 a, Vector3 b) {
return { a.x - b.x, a.y - b.y, a.z - b.z };
}
__device__ __host__ double vec_dot(Vector3 a, Vector3 b) {
__device__ double vec_dot(Vector3 a, Vector3 b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
__device__ __host__ Vector3 vec_cross(Vector3 a, Vector3 b) {
__device__ Vector3 vec_cross(Vector3 a, Vector3 b) {
return { a.y * b.z - a.z * b.y,
a.z * b.x - a.x * b.z,
a.x * b.y - a.y * b.x };
}
__device__ __host__ Vector3 vec_normalize(Vector3 v) {
__device__ Vector3 vec_normalize(Vector3 v) {
double len = sqrt(vec_dot(v, v));
return (len > 1e-12) ? Vector3{ v.x / len, v.y / len, v.z / len } : v;
}
// 标准正交基底坐标计算
__device__ __host__ Vector3 coordinates_orthonormal_basis(Vector3& A,
__device__ Vector3 coordinates_orthonormal_basis(Vector3& A,
Vector3& e1,
Vector3& e2,
Vector3& e3) {
@ -87,7 +87,7 @@ __device__ __host__ Vector3 coordinates_orthonormal_basis(Vector3& A,
// 计算视线交点T
extern __device__ __host__ Vector3 compute_T(Vector3 S, Vector3 ray, double H) {
extern __device__ Vector3 compute_T(Vector3 S, Vector3 ray, double H) {
Vector3 dir = vec_normalize(ray);
double a_h = WGS84_A + H;
@ -105,7 +105,7 @@ extern __device__ __host__ Vector3 compute_T(Vector3 S, Vector3 ray, double H) {
}
// 主计算函数A
extern __device__ __host__ Vector3 compute_P(Vector3 S, Vector3 T, double R, double H) {
extern __device__ Vector3 compute_P(Vector3 S, Vector3 T, double R, double H) {
Vector3 ex, ey, ez; // 平面基函数
Vector3 ST = vec_normalize(vec_sub(T, S));// S->T
Vector3 SO = vec_normalize(vec_sub(Vector3{ 0, 0, 0 }, S)); // S->O
@ -203,7 +203,44 @@ extern __device__ __host__ Vector3 compute_P(Vector3 S, Vector3 T, double R, dou
return P;
}
__device__ double angleBetweenVectors_single(Vector3_single a, Vector3_single b, bool returnDegrees)
{
// 计算点积
float dotProduct = a.x * b.x + a.y * b.y + a.z * b.z;
// 计算模长
float magA = sqrtf(a.x * a.x + a.y * a.y + a.z * a.z);
float magB = sqrtf(b.x * b.x + b.y * b.y + b.z * b.z);
// 处理零向量异常
if (magA == 0.0 || magB == 0.0) {
return NAN;
}
float cosTheta = dotProduct / (magA * magB);
cosTheta = cosTheta < -1 ? -1 : cosTheta>1 ? 1 : cosTheta; // 截断到[-1, 1]
float angleRad = acosf(cosTheta);
return returnDegrees ? angleRad * 180.0 / M_PI : angleRad;
}
__device__ Vector3_single vec_sub_single(Vector3_single a, Vector3_single b)
{
return { a.x - b.x, a.y - b.y, a.z - b.z };
}
__device__ float vec_dot_single(Vector3_single a, Vector3_single b)
{
return a.x * b.x + a.y * b.y + a.z * b.z;
}
__device__ Vector3_single vec_cross_single(Vector3_single a, Vector3_single b)
{
return Vector3_single{ a.y * b.z - a.z * b.y,
a.z * b.x - a.x * b.z,
a.x * b.y - a.y * b.x };
}
//// 参数校验与主函数
//int main() {

View File

@ -8,12 +8,16 @@
extern __device__ __host__ double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees = false);
extern __device__ __host__ Vector3 vec_sub(Vector3 a, Vector3 b);
extern __device__ __host__ double vec_dot(Vector3 a, Vector3 b);
extern __device__ __host__ Vector3 vec_cross(Vector3 a, Vector3 b);
extern __device__ __host__ Vector3 vec_normalize(Vector3 v);
extern __device__ __host__ Vector3 compute_T(Vector3 S, Vector3 ray_dir, double H);
//
extern __device__ __host__ Vector3 compute_P(Vector3 S, Vector3 T, double R, double H );
extern __device__ double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees = false);
extern __device__ Vector3 vec_sub(Vector3 a, Vector3 b);
extern __device__ double vec_dot(Vector3 a, Vector3 b);
extern __device__ Vector3 vec_cross(Vector3 a, Vector3 b);
extern __device__ Vector3 vec_normalize(Vector3 v);
extern __device__ Vector3 compute_T(Vector3 S, Vector3 ray_dir, double H);
extern __device__ Vector3 compute_P(Vector3 S, Vector3 T, double R, double H );
extern __device__ double angleBetweenVectors_single(Vector3_single a, Vector3_single b, bool returnDegrees = false);
extern __device__ Vector3_single vec_sub_single(Vector3_single a, Vector3_single b);
extern __device__ float vec_dot_single(Vector3_single a, Vector3_single b);
extern __device__ Vector3_single vec_cross_single(Vector3_single a, Vector3_single b);

View File

@ -1,9 +1,4 @@
#include <time.h>
#include <iostream>
#include <memory>
#include <cmath>
#include <complex>
#include <cuda.h>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
@ -19,12 +14,13 @@
/* »úÆ÷º¯Êý ****************************************************************************************************************************/
extern __host__ __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta) {//ÏßÐÔÖµ
extern __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta) {//ÏßÐÔÖµ
double sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
return sigma;
}
extern __host__ __device__ double GPU_getSigma0dB(
__device__ double GPU_getSigma0dB_params(
const double p1, const double p2, const double p3, const double p4, const double p5, const double p6,
double theta) {//ÏßÐÔÖµ
return p1 + p2 * exp(-p3 * theta) + p4 * cos(p5 * theta + p6);
@ -32,7 +28,7 @@ extern __host__ __device__ double GPU_getSigma0dB(
extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
extern __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
double RstX, double RstY, double RstZ,
double AntXaxisX, double AntXaxisY, double AntXaxisZ,
double AntYaxisX, double AntYaxisY, double AntYaxisZ,
@ -111,7 +107,7 @@ extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
return result;
}
extern __host__ __device__ double GPU_BillerInterpAntPattern(double* antpattern,
extern __device__ double GPU_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints,
double searththeta, double searchphi) {

View File

@ -20,15 +20,15 @@ extern "C" struct CUDASigmaParam {
double p6;
};
extern __host__ __device__ double GPU_getSigma0dB(
extern __device__ double GPU_getSigma0dB_params(
const double p1, const double p2, const double p3, const double p4, const double p5, const double p6,
double theta);
extern __host__ __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta);
extern __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta);
extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
extern __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
double RstX, double RstY, double RstZ,
double AntXaxisX, double AntXaxisY, double AntXaxisZ,
double AntYaxisX, double AntYaxisY, double AntYaxisZ,
@ -36,7 +36,7 @@ extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
double AntDirectX, double AntDirectY, double AntDirectZ
);
extern __host__ __device__ double GPU_BillerInterpAntPattern(double* antpattern,
extern __device__ double GPU_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints,
double searththeta, double searchphi);

View File

@ -0,0 +1,468 @@
#include <cuda.h>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "BaseConstVariable.h"
#include "GPURFPC_single.cuh"
#ifdef __CUDANVCC___
/* 机器函数 ****************************************************************************************************************************/
extern __device__ float GPU_getSigma0dB_single(CUDASigmaParam_single param, float theta) {//线性值
float sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
return sigma;
}
extern __device__ float GPU_getSigma0dB_single(const float p1, const float p2, const float p3, const float p4, const float p5, const float p6, float theta)
{
return p1 + p2 * expf(-p3 * theta) + p4 * cosf(p5 * theta + p6);
}
extern __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal_single(
float RstX, float RstY, float RstZ,
float AntXaxisX, float AntXaxisY, float AntXaxisZ,
float AntYaxisX, float AntYaxisY, float AntYaxisZ,
float AntZaxisX, float AntZaxisY, float AntZaxisZ,
float AntDirectX, float AntDirectY, float AntDirectZ
) {
CUDAVectorEllipsoidal result{ 0,0,-1 };
// 求解天线增益
float Xst = -1 * RstX; // 卫星 --> 地面
float Yst = -1 * RstY;
float Zst = -1 * RstZ;
// 归一化
float RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
float AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
float AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
float AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
float Rx = Xst / RstNorm;
float Ry = Yst / RstNorm;
float Rz = Zst / RstNorm;
float Xx = AntXaxisX / AntXaxisNorm;
float Xy = AntXaxisY / AntXaxisNorm;
float Xz = AntXaxisZ / AntXaxisNorm;
float Yx = AntYaxisX / AntYaxisNorm;
float Yy = AntYaxisY / AntYaxisNorm;
float Yz = AntYaxisZ / AntYaxisNorm;
float Zx = AntZaxisX / AntZaxisNorm;
float Zy = AntZaxisY / AntZaxisNorm;
float Zz = AntZaxisZ / AntZaxisNorm;
float Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
// 计算theta 与 phi
float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
float Zn = Zant / Norm;
float ThetaAnt = ( - 1 > Zn) ? PI : (Zn > 1 ? 0 : acos(Zn));// acosf(Zant / Norm); // theta 与 Z轴的夹角
float PhiAnt = abs(Xant)<PRECISIONTOLERANCE ?0: atanf(Yant / Xant); // -pi/2 ~pi/2
if (abs(Yant) < PRECISIONTOLERANCE) { // X轴上
PhiAnt = 0;
}
else if (abs(Xant) < PRECISIONTOLERANCE) { // Y轴上原点
if (Yant > 0) {
PhiAnt = PI / 2;
}
else {
PhiAnt = -PI / 2;
}
}
else if (Xant < 0) {
if (Yant > 0) {
PhiAnt = PI + PhiAnt;
}
else {
PhiAnt = -PI + PhiAnt;
}
}
else { // Xant>0 X 正轴
}
if (isnan(PhiAnt)) {
printf("V=[%f,%f,%f];norm=%f;thetaAnt=%f;phiAnt=%f;\n", Xant, Yant, Zant, Norm, ThetaAnt, PhiAnt);
}
result.theta = ThetaAnt;
result.phi = PhiAnt;
result.Rho = Norm;
return result;
}
extern __device__ float GPU_BillerInterpAntPattern_single(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
long thetapoints, long phipoints,
float searththeta, float searchphi) {
float stheta = searththeta;
float sphi = searchphi;
if (stheta > 90) {
return 0;
}
else {}
float pthetaid = (stheta - starttheta) / dtheta;//
float pphiid = (sphi - startphi) / dphi;
long lasttheta = floorf(pthetaid);
long nextTheta = lasttheta + 1;
long lastphi = floorf(pphiid);
long nextPhi = lastphi + 1;
if (lasttheta < 0 || nextTheta < 0 || lastphi < 0 || nextPhi < 0 ||
lasttheta >= thetapoints || nextTheta >= thetapoints || lastphi >= phipoints || nextPhi >= phipoints)
{
return 0;
}
else {
float x = stheta;
float y = sphi;
float x1 = lasttheta * dtheta + starttheta;
float x2 = nextTheta * dtheta + starttheta;
float y1 = lastphi * dphi + startphi;
float y2 = nextPhi * dphi + startphi;
float z11 = antpattern[lasttheta * phipoints + lastphi];
float z12 = antpattern[lasttheta * phipoints + nextPhi];
float z21 = antpattern[nextTheta * phipoints + lastphi];
float z22 = antpattern[nextTheta * phipoints + nextPhi];
//z11 = powf(10, z11 / 10); // dB-> 线性
//z12 = powf(10, z12 / 10);
//z21 = powf(10, z21 / 10);
//z22 = powf(10, z22 / 10);
float GainValue = (z11 * (x2 - x) * (y2 - y)
+ z21 * (x - x1) * (y2 - y)
+ z12 * (x2 - x) * (y - y1)
+ z22 * (x - x1) * (y - y1));
GainValue = GainValue / ((x2 - x1) * (y2 - y1));
return GainValue;
}
}
/* 核函数 ****************************************************************************************************************************/
// 计算每块
__global__ void CUDA_Kernel_Computer_R_amp_single(
float* antX, float* antY, float* antZ,
float* antXaxisX, float* antXaxisY, float* antXaxisZ,
float* antYaxisX, float* antYaxisY, float* antYaxisZ,
float* antZaxisX, float* antZaxisY, float* antZaxisZ,
float* antDirectX, float* antDirectY, float* antDirectZ,
long PRFCount, // 整体的脉冲数,
float* targetX, float* targetY, float* targetZ, long* demCls,
float* demSlopeX, float* demSlopeY, float* demSlopeZ ,
long startPosId, long pixelcount,
CUDASigmaParam_single* sigma0Paramslist, long sigmaparamslistlen,
float Pt,
float refPhaseRange,
float* TransAntpattern,
float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints,
float* ReceiveAntpattern,
float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,
float maxTransAntPatternValue, float maxReceiveAntPatternValue,
float NearR, float FarR,
float* d_temp_R, float* d_temp_amps// 计算输出
) {
long idx = blockIdx.x * blockDim.x + threadIdx.x; // 获取当前的线程编码
long prfId = idx / SHAREMEMORY_FLOAT_HALF;
long posId = idx % SHAREMEMORY_FLOAT_HALF+ startPosId; // 当前线程对应的影像点
if (prfId < PRFCount && posId < pixelcount) {
float RstX = antX[prfId] - targetX[posId]; // 计算坐标矢量
float RstY = antY[prfId] - targetY[posId];
float RstZ = antZ[prfId] - targetZ[posId];
float RstR = sqrt(RstX * RstX + RstY * RstY + RstZ * RstZ); // 矢量距离
if (RstR<NearR || RstR>FarR) {
d_temp_R[idx] = 0;
d_temp_amps[idx] = 0;
return;
}
else {
float slopeX = demSlopeX[posId];
float slopeY = demSlopeY[posId];
float slopeZ = demSlopeZ[posId];
float slopR = sqrtf(slopeX * slopeX + slopeY * slopeY + slopeZ * slopeZ); //
if (abs(slopR - 0) > 1e-3) {
float dotAB = RstX * slopeX + RstY * slopeY + RstZ * slopeZ;
float localangle = acos(dotAB / (RstR * slopR));
if (localangle < 0 || localangle >= LAMP_CUDA_PI / 2|| isnan(localangle)) {
d_temp_R[idx] = 0;
d_temp_amps[idx] = 0;
return;
}
else {}
float ampGain = 0;
// 求解天线方向图指向
CUDAVectorEllipsoidal antVector = GPU_SatelliteAntDirectNormal_single(
RstX, RstY, RstZ,
antXaxisX[prfId], antXaxisY[prfId], antXaxisZ[prfId],
antYaxisX[prfId], antYaxisY[prfId], antYaxisZ[prfId],
antZaxisX[prfId], antZaxisY[prfId], antZaxisZ[prfId],
antDirectX[prfId], antDirectY[prfId], antDirectZ[prfId]
);
antVector.theta = antVector.theta * r2d;
antVector.phi = antVector.phi * r2d;
//printf("theta: %f , phi: %f \n", antVector.theta, antVector.phi);
if (antVector.Rho > 0) {
//float TansantPatternGain = GPU_BillerInterpAntPattern(
// TransAntpattern,
// Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
// antVector.theta, antVector.phi);
//float antPatternGain = GPU_BillerInterpAntPattern(
// ReceiveAntpattern,
// Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
// antVector.theta, antVector.phi);
float sigma0 = 0;
{
long clsid = demCls[posId];
//printf("clsid=%d\n", clsid);
CUDASigmaParam_single tempsigma = sigma0Paramslist[clsid];
if (abs(tempsigma.p1) < PRECISIONTOLERANCE &&
abs(tempsigma.p2) < PRECISIONTOLERANCE &&
abs(tempsigma.p3) < PRECISIONTOLERANCE &&
abs(tempsigma.p4) < PRECISIONTOLERANCE &&
abs(tempsigma.p5) < PRECISIONTOLERANCE &&
abs(tempsigma.p6) < PRECISIONTOLERANCE
) {
sigma0 = 0;
}
else {
float sigma = GPU_getSigma0dB_single(tempsigma, localangle);
sigma0 = powf(10.0, sigma / 10.0);
}
}
//ampGain = TansantPatternGain * antPatternGain;
ampGain = 1;
//if (10 * log10(ampGain / maxReceiveAntPatternValue / maxTransAntPatternValue) < -3) { // 小于-3dB
// d_temp_R[idx] = 0;
// d_temp_amps[idx] = 0;
// return;
//}
//else {}
ampGain = ampGain / (powf(4 * LAMP_CUDA_PI, 2) * powf(RstR, 4)); // 反射强度
float temp_amp = float(ampGain * Pt * sigma0);
float temp_R = float(RstR - refPhaseRange);
if (isnan(temp_amp) || isnan(temp_R)|| isinf(temp_amp) || isinf(temp_R)) {
printf("amp is nan or R is nan,amp=%f;R=%f; \n", temp_amp, temp_R);
d_temp_R[idx] = 0;
d_temp_amps[idx] = 0;
return;
}
else {}
d_temp_amps[idx] = temp_amp;
d_temp_R[idx] = temp_R;
return;
}
else {
d_temp_R[idx] = 0;
d_temp_amps[idx] = 0;
return;
}
}
else {
d_temp_R[idx] = 0;
d_temp_amps[idx] = 0;
return;
}
}
}
}
__global__ void CUDA_Kernel_Computer_echo_single(
float* d_temp_R, float* d_temp_amps, long posNum,
float f0, float dfreq,
long FreqPoints, // 当前频率的分块
long maxfreqnum, // 最大脉冲值
float* d_temp_echo_real, float* d_temp_echo_imag,
long temp_PRF_Count
) {
__shared__ float s_R[SHAREMEMORY_FLOAT_HALF]; // 注意一个完整的block_size 共享相同内存
__shared__ float s_amp[SHAREMEMORY_FLOAT_HALF];
long tid = threadIdx.x;
long bid = blockIdx.x;
long idx = bid * blockDim.x + tid;
long prfId = idx / FreqPoints; // 脉冲ID
long fId = idx % FreqPoints;//频率ID
long psid = 0;
long pixelId = 0;
for (long ii = 0; ii < SHAREMEMORY_FLOAT_HALF_STEP; ii++) { // SHAREMEMORY_FLOAT_HALF_STEP * BLOCK_SIZE=SHAREMEMORY_FLOAT_HALF
psid = tid * SHAREMEMORY_FLOAT_HALF_STEP + ii;
pixelId = prfId * posNum + psid; //
if (psid < posNum) {
s_R[psid] = d_temp_R[pixelId];
s_amp[psid] = d_temp_amps[pixelId];
}
else {
s_R[psid] = 0;
s_amp[psid] = 0;
}
}
__syncthreads(); // 确定所有待处理数据都已经进入程序中
if (fId < maxfreqnum && prfId < temp_PRF_Count) {
long echo_ID = prfId * maxfreqnum + fId; // 计算对应的回波位置
float factorjTemp = RFPCPIDIVLIGHT * (f0 + fId * dfreq);
float temp_real = 0;
float temp_imag = 0;
float temp_phi = 0;
float temp_amp = 0;
for (long dataid = 0; dataid < SHAREMEMORY_FLOAT_HALF; dataid++) {
temp_phi = s_R[dataid] * factorjTemp;
temp_amp = s_amp[dataid];
temp_real += (temp_amp * cosf(temp_phi));
temp_imag += (temp_amp * sinf(temp_phi));
//if (dataid > 5000) {
// printf("echo_ID=%d; dataid=%d;ehodata=(%f,%f);R=%f;amp=%f;\n", echo_ID, dataid, temp_real, temp_imag, s_R[0], s_amp[0]);
//}
if (isnan(temp_phi) || isnan(temp_amp) || isnan(temp_real) || isnan(temp_imag)
|| isinf(temp_phi) || isinf(temp_amp) || isinf(temp_real) || isinf(temp_imag)
) {
printf("[amp,phi,real,imag]=[%f,%f,%f,%f];\n",temp_amp,temp_phi,temp_real,temp_imag);
}
}
//printf("echo_ID=%d; ehodata=(%f,%f)\n", echo_ID, temp_real, temp_imag);
//printf("(%f %f %f) ", factorjTemp, s_amp[0], s_R[0]);
d_temp_echo_real[echo_ID] += /*d_temp_echo_real[echo_ID] + */temp_real;
d_temp_echo_imag[echo_ID] += /*d_temp_echo_imag[echo_ID] +*/ temp_imag;
}
}
/**
* 分块计算主流程
*/
void CUDA_RFPC_MainProcess_single(
float* antX, float* antY, float* antZ,
float* antXaxisX, float* antXaxisY, float* antXaxisZ,
float* antYaxisX, float* antYaxisY, float* antYaxisZ,
float* antZaxisX, float* antZaxisY, float* antZaxisZ,
float* antDirectX, float* antDirectY, float* antDirectZ,
long PRFCount, long FreqNum,
float f0, float dfreq,
float Pt,
float refPhaseRange,
float* TransAntpattern,
float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints,
float* ReceiveAntpattern,
float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,
float maxTransAntPatternValue, float maxReceiveAntPatternValue,
float NearR, float FarR,
float* targetX, float* targetY, float* targetZ, long* demCls, long TargetNumber,
float* demSlopeX, float* demSlopeY, float* demSlopeZ,
CUDASigmaParam_single* sigma0Paramslist, long sigmaparamslistlen,
float* out_echoReal, float* out_echoImag,
float* d_temp_R, float* d_temp_amp
)
{
long BLOCK_FREQNUM = NextBlockPad(FreqNum, BLOCK_SIZE); // 256*freqBlockID
long cudaBlocknum = 0;
long freqpoints = BLOCK_FREQNUM;
printf("freqpoints:%d\n", freqpoints);
long process = 0;
for (long sTi = 0; sTi < TargetNumber; sTi = sTi + SHAREMEMORY_FLOAT_HALF) {
cudaBlocknum = (PRFCount * SHAREMEMORY_FLOAT_HALF + BLOCK_SIZE - 1) / BLOCK_SIZE;
CUDA_Kernel_Computer_R_amp_single << <cudaBlocknum, BLOCK_SIZE >> > (
antX, antY, antZ,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ,
PRFCount,
targetX, targetY, targetZ, demCls,
demSlopeX, demSlopeY, demSlopeZ,
sTi, TargetNumber,
sigma0Paramslist, sigmaparamslistlen,
Pt,
refPhaseRange,
TransAntpattern,
Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
ReceiveAntpattern,
Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
maxTransAntPatternValue, maxReceiveAntPatternValue,
NearR, FarR,
d_temp_R, d_temp_amp// 计算输出
);
PrintLasterError("CUDA_Kernel_Computer_R_amp");
cudaBlocknum = (PRFCount * BLOCK_FREQNUM + BLOCK_SIZE - 1) / BLOCK_SIZE;
CUDA_Kernel_Computer_echo_single << <cudaBlocknum, BLOCK_SIZE >> > (
d_temp_R, d_temp_amp, SHAREMEMORY_FLOAT_HALF,
f0, dfreq,
freqpoints, FreqNum,
out_echoReal, out_echoImag,
PRFCount
);
PrintLasterError("CUDA_Kernel_Computer_echo");
if ((sTi * 100.0 / TargetNumber ) - process >= 1) {
process = sTi * 100.0 / TargetNumber;
PRINT("TargetID [%f]: %d / %d finished\n", sTi*100.0/ TargetNumber,sTi, TargetNumber);
}
}
cudaDeviceSynchronize();
}
#endif

View File

@ -0,0 +1,85 @@
#ifndef _GPURFPC_SINGLE_H_
#define _GPURFPC_SINGLE_H_
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#define RFPCPIDIVLIGHT -4*PI/(LIGHTSPEED/1e9)
extern "C" struct CUDASigmaParam_single {
float p1;
float p2;
float p3;
float p4;
float p5;
float p6;
};
extern __device__ float GPU_getSigma0dB_single(const float p1, const float p2, const float p3, const float p4, const float p5, const float p6, float theta)
;
extern __device__ float GPU_getSigma0dB_single(CUDASigmaParam_single param, float theta);
extern __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal_single(
float RstX, float RstY, float RstZ,
float AntXaxisX, float AntXaxisY, float AntXaxisZ,
float AntYaxisX, float AntYaxisY, float AntYaxisZ,
float AntZaxisX, float AntZaxisY, float AntZaxisZ,
float AntDirectX, float AntDirectY, float AntDirectZ
);
extern __device__ float GPU_BillerInterpAntPattern_single(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
long thetapoints, long phipoints,
float searththeta, float searchphi);
extern "C" void CUDA_RFPC_MainProcess_single(
// 天线
float* antX, float* antY, float* antZ, // 天线坐标
float* antXaxisX, float* antXaxisY, float* antXaxisZ, // 天线坐标系的X轴
float* antYaxisX, float* antYaxisY, float* antYaxisZ,// 天线坐标系的Y轴
float* antZaxisX, float* antZaxisY, float* antZaxisZ,// 天线坐标系的Z轴
float* antDirectX, float* antDirectY, float* antDirectZ,// 天线的指向
long PRFCount, long FreqNum, // 脉冲数量,频率数量
float f0, float dfreq,// 起始频率,终止频率
float Pt,// 发射能量
float refPhaseRange,
// 天线方向图
float* TransAntpattern, float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
float* ReceiveAntpattern, float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
float maxTransAntPatternValue,float maxReceiveAntPatternValue,
float NearR, float FarR, // 距离范围
// 地面
float* targetX, float* targetY, float* targetZ, long* demCls, long TargetPixelNumber, // 地面坐标、地表覆盖类型,像素数
float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
CUDASigmaParam_single* sigma0Paramslist, long sigmaparamslistlen,// 插值图像
float* out_echoReal, float* out_echoImag,// 输出回波
float* d_temp_R, float* d_temp_amp
);
#endif

View File

@ -222,6 +222,7 @@
<ClCompile Include="SimulationSAR\TBPImageAlgCls.cpp" />
<ClCompile Include="UnitTestMain.cpp" />
<CudaCompile Include="GPUBpSimulation.cu" />
<CudaCompile Include="SimulationSAR\GPURFPC_single.cu" />
<CudaCompile Include="SimulationSAR\GPUTBPImage.cu" />
<QtMoc Include="PowerSimulationIncoherent\QSimulationSARPolynomialOrbitModel.h" />
<CudaCompile Include="PowerSimulationIncoherent\LookTableSimulationComputer.cuh" />
@ -238,6 +239,7 @@
<CudaCompile Include="SimulationSAR\BPBasic0_CUDA.cu" />
<CudaCompile Include="SimulationSAR\GPUBPTool.cuh" />
<ClInclude Include="SimulationSAR\BPBasic0_CUDA.cuh" />
<ClInclude Include="SimulationSAR\GPURFPC_single.cuh" />
<ClInclude Include="SimulationSAR\GPUTBPImage.cuh" />
<ClInclude Include="SimulationSAR\RFPCProcessCls.h" />
<ClInclude Include="SimulationSAR\SARSatelliteSimulationAbstractCls.h" />

View File

@ -68,6 +68,9 @@
<ClInclude Include="GPUBpSimulation.cuh">
<Filter>SimulationSAR</Filter>
</ClInclude>
<ClInclude Include="SimulationSAR\GPURFPC_single.cuh">
<Filter>SimulationSAR</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ClCompile Include="SimulationSAR\QImageSARRFPC.cpp">
@ -193,5 +196,8 @@
<CudaCompile Include="GPUBpSimulation.cu">
<Filter>SimulationSAR</Filter>
</CudaCompile>
<CudaCompile Include="SimulationSAR\GPURFPC_single.cu">
<Filter>SimulationSAR</Filter>
</CudaCompile>
</ItemGroup>
</Project>

View File

@ -44,7 +44,7 @@ void testSimualtionEchoPoint() {
.arg(h_data.AntZ[gpspoints - 1]);
/** 2. 成像网格 **************************************************************************************************/
qDebug() << u8"轨道文件读取结束\n2.成像网格读取。。。";
QString demxyzPath = u8"D:\\Programme\\vs2022\\RasterMergeTest\\simulationData\\demdataset\\demxyz.bin";
QString demxyzPath = u8"C:\\Users\\30453\\Desktop\\script\\已修改GF3场景\\data\\demxyz.bin";
gdalImage demgridimg(demxyzPath);
long dem_rowCount = demgridimg.height;
long dem_ColCount = demgridimg.width;
@ -115,7 +115,7 @@ void testSimualtionEchoPoint() {
/** 7. 目标 **************************************************************************************************/
double Tx = -2029086.618142, Ty = 4139594.934504, Tz = 4392846.782027;
double Tslx = -2029086.618142, Tsly = 4139594.934504, Tslz = 4392846.782027;
double p1 = 33.3, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
double p1 = 1, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
/** 7. 构建回波 **************************************************************************************************/
GPUDATA d_data;
@ -404,17 +404,405 @@ void testBpImage() {
}
//
//
//void testSimualtionEchoPoint_single() {
// GPUDATA_single h_data{};
// /** 1. 轨道 **************************************************************************************************/
// qDebug() << u8"1.轨道文件读取中。。。";
// QString inGPSPath = u8"C:\\Users\\30453\\Desktop\\script\\data\\GF3_Simulation.gpspos.data";
// long gpspoints = gdalImage(inGPSPath).height;
// std::shared_ptr<SatelliteAntPos> antpos = SatelliteAntPosOperator::readAntPosFile(inGPSPath, gpspoints);
// h_data.AntX = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
// h_data.AntY = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
// h_data.AntZ = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
//
// for (long i = 0; i < gpspoints; i = i + 1) {
// h_data.AntX[i] = antpos.get()[i].Px;
// h_data.AntY[i] = antpos.get()[i].Py;
// h_data.AntZ[i] = antpos.get()[i].Pz;
// }
// //gpspoints = gpspoints / 2;
// qDebug() << "GPS points Count:\t" << gpspoints;
// qDebug() << "PRF 0:\t " << QString("%1,%2,%3").arg(h_data.AntX[0]).arg(h_data.AntY[0]).arg(h_data.AntZ[0]);
// qDebug() << "PRF " << gpspoints - 1 << ":\t " << QString("%1,%2,%3")
// .arg(h_data.AntX[gpspoints - 1])
// .arg(h_data.AntY[gpspoints - 1])
// .arg(h_data.AntZ[gpspoints - 1]);
// /** 2. 成像网格 **************************************************************************************************/
// qDebug() << u8"轨道文件读取结束\n2.成像网格读取。。。";
// QString demxyzPath = u8"C:\\Users\\30453\\Desktop\\script\\已修改GF3场景\\data\\demxyz.bin";
// gdalImage demgridimg(demxyzPath);
// long dem_rowCount = demgridimg.height;
// long dem_ColCount = demgridimg.width;
// std::shared_ptr<float> demX = readDataArr<float>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
// std::shared_ptr<float> demY = readDataArr<float>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
// std::shared_ptr<float> demZ = readDataArr<float>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
//
// h_data.x_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
// h_data.y_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
// h_data.z_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
//
// for (long i = 0; i < dem_rowCount; i++) {
// for (long j = 0; j < dem_ColCount; j++) {
// h_data.x_mat[i * dem_ColCount + j] = demX.get()[i * dem_ColCount + j];
// h_data.y_mat[i * dem_ColCount + j] = demY.get()[i * dem_ColCount + j];
// h_data.z_mat[i * dem_ColCount + j] = demZ.get()[i * dem_ColCount + j];
// }
// }
// qDebug() << "dem row Count:\t" << dem_rowCount;
// qDebug() << "dem col Count:\t" << dem_ColCount;
//
// qDebug() << u8"成像网格读取结束";
// /** 3. 频率网格 **************************************************************************************************/
// qDebug() << u8"3.处理频率";
// float centerFreq = 5.3e9;
// float bandwidth = 40e6;
// long freqpoints = 2048;
// std::shared_ptr<float> freqlist(getFreqPoints_mallocHost_single(centerFreq - bandwidth / 2, centerFreq + bandwidth / 2, freqpoints), FreeCUDAHost);
// std::shared_ptr<float> minF(new float[gpspoints], delArrPtr);
// for (long i = 0; i < gpspoints; i++) {
// minF.get()[i] = freqlist.get()[0];
// }
// h_data.deltaF = bandwidth / (freqpoints - 1);
// qDebug() << "start Freq:\t" << centerFreq - bandwidth / 2;
// qDebug() << "end Freq:\t" << centerFreq + bandwidth / 2;
// qDebug() << "freq points:\t" << freqpoints;
// qDebug() << "delta freq:\t" << freqlist.get()[1] - freqlist.get()[0];
// qDebug() << u8"频率结束";
// /** 4. 初始化回波 **************************************************************************************************/
// qDebug() << u8"4.初始化回波";
// std::shared_ptr<cuComplex> phdata(createEchoPhase_mallocHost(gpspoints, freqpoints), FreeCUDAHost);
// qDebug() << "Azimuth Points:\t" << gpspoints;
// qDebug() << "Range Points:\t" << freqpoints;
// qDebug() << u8"初始化回波结束";
// /** 5. 初始化图像 **************************************************************************************************/
// qDebug() << u8"5.初始化图像";
// h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * dem_rowCount * dem_ColCount);
// qDebug() << "Azimuth Points:\t" << gpspoints;
// qDebug() << "Range Points:\t" << freqpoints;
// qDebug() << u8"初始化图像结束";
//
// /** 6. 模型参数初始化 **************************************************************************************************/
// qDebug() << u8"6.模型参数初始化";
//
// h_data.nx = dem_ColCount;
// h_data.ny = dem_rowCount;
//
// h_data.Np = gpspoints;
// h_data.Freq = freqlist.get();
// h_data.minF = minF.get();
// h_data.Nfft = freqpoints;
// h_data.K = h_data.Nfft;
// h_data.phdata = phdata.get();
// h_data.R0 = 900000;
// qDebug() << u8"模型参数结束";
//
//
// /** 7. 目标 **************************************************************************************************/
// float Tx = -2029086.618142, Ty = 4139594.934504, Tz = 4392846.782027;
// float Tslx = -2029086.618142, Tsly = 4139594.934504, Tslz = 4392846.782027;
// float p1 = 1, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
//
// /** 7. 构建回波 **************************************************************************************************/
// GPUDATA_single d_data;
// initGPUData_single(h_data, d_data);
//
// RFPCProcess_single(Tx, Ty, Tz,
// Tslx, Tsly, Tslz, // 目标的坡面向量
// p1, p2, p3, p4, p5, p6,
// d_data);
//
//
// /** 8. 展示回波 **************************************************************************************************/
// {
// DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// cuComplex* h_ifftphdata = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// cuComplex* d_ifftphdata = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// CUDAIFFT(d_data.phdata, d_ifftphdata, d_data.Np, d_data.Nfft, d_data.Nfft);
// FFTShift1D(d_ifftphdata, d_data.Np, d_data.Nfft);
// DeviceToHost(h_ifftphdata, d_ifftphdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
// {
// for (long i = 0; i < d_data.Np; i++) {
// for (long j = 0; j < d_data.Nfft; j++) {
// ifftdata.get()[i * d_data.Nfft + j] =
// std::complex<double>(
// h_ifftphdata[i * d_data.Nfft + j].x,
// h_ifftphdata[i * d_data.Nfft + j].y);
// }
// }
// }
// testOutComplexDoubleArr(QString("echo_ifft_single.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
//
// FreeCUDADevice(d_ifftphdata);
// FreeCUDAHost(h_ifftphdata);
//
// }
// /** 9. 成像 **************************************************************************************************/
//
// // 计算maxWr需要先计算deltaF
// float deltaF = h_data.deltaF; // 从输入参数获取
// float maxWr = 299792458.0f / (2.0f * deltaF);
// qDebug() << "maxWr :\t" << maxWr;
// float* r_vec_host = (float*)mallocCUDAHost(sizeof(float) * h_data.Nfft);// new float[data.Nfft];
// const float step = maxWr / h_data.Nfft;
// const float start = -1 * h_data.Nfft / 2.0f * step;
// printf("nfft=%d\n", h_data.Nfft);
//
// for (int i = 0; i < h_data.Nfft; ++i) {
// r_vec_host[i] = start + i * step;
// }
//
// h_data.r_vec = r_vec_host;
// d_data.r_vec = h_data.r_vec;
//
// qDebug() << "r_vec [0]:\t" << h_data.r_vec[0];
// qDebug() << "r_vec step:\t" << step;
//
// bpBasic0CUDA_single(d_data, 0);
//
// DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * d_data.nx * d_data.ny);
//
// {
// DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
// {
// for (long i = 0; i < d_data.Np; i++) {
// for (long j = 0; j < d_data.Nfft; j++) {
// ifftdata.get()[i * d_data.Nfft + j] =
// std::complex<double>(
// h_data.phdata[i * d_data.Nfft + j].x,
// h_data.phdata[i * d_data.Nfft + j].y);
// }
// }
// }
// testOutComplexDoubleArr(QString("echo_ifft_BPBasic_single.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
//
//
//
// std::shared_ptr<std::complex<double>> im_finals(new std::complex<double>[d_data.nx * d_data.ny], delArrPtr);
// {
// for (long i = 0; i < d_data.ny; i++) {
// for (long j = 0; j < d_data.nx; j++) {
// im_finals.get()[i * d_data.nx + j] = std::complex<double>(
// h_data.im_final[i * d_data.nx + j].x,
// h_data.im_final[i * d_data.nx + j].y);
// }
// }
// }
// testOutComplexDoubleArr(QString("im_finals_single.bin"), im_finals.get(), d_data.ny, d_data.nx);
// }
//}
//
//
//void testBpImage_single() {
//
// GPUDATA_single h_data{};
// /** 0. 轨道 **************************************************************************************************/
// QString echoPath = "D:\\Programme\\vs2022\\RasterMergeTest\\LAMPCAE_SCANE-all-scane\\GF3_Simulation.xml";
// std::shared_ptr<EchoL0Dataset> echoL0ds(new EchoL0Dataset);
// echoL0ds->Open(echoPath);
// qDebug() << u8"加载回拨文件:\t" << echoPath;
//
// /** 1. 轨道 **************************************************************************************************/
// qDebug() << u8"1.轨道文件读取中。。。";
// QString inGPSPath = echoL0ds->getGPSPointFilePath();
// long gpspoints = gdalImage(inGPSPath).height;
// std::shared_ptr<SatelliteAntPos> antpos = SatelliteAntPosOperator::readAntPosFile(inGPSPath, gpspoints);
// h_data.AntX = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
// h_data.AntY = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
// h_data.AntZ = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
//
// for (long i = 0; i < gpspoints; i = i + 1) {
// h_data.AntX[i] = antpos.get()[i].Px;
// h_data.AntY[i] = antpos.get()[i].Py;
// h_data.AntZ[i] = antpos.get()[i].Pz;
// }
// //gpspoints = gpspoints / 2;
// qDebug() << "GPS points Count:\t" << gpspoints;
// qDebug() << "PRF 0:\t " << QString("%1,%2,%3").arg(h_data.AntX[0]).arg(h_data.AntY[0]).arg(h_data.AntZ[0]);
// qDebug() << "PRF " << gpspoints - 1 << ":\t " << QString("%1,%2,%3")
// .arg(h_data.AntX[gpspoints - 1])
// .arg(h_data.AntY[gpspoints - 1])
// .arg(h_data.AntZ[gpspoints - 1]);
// /** 2. 成像网格 **************************************************************************************************/
// qDebug() << u8"轨道文件读取结束\n2.成像网格读取。。。";
// QString demxyzPath = u8"D:\\Programme\\vs2022\\RasterMergeTest\\simulationData\\demdataset\\demxyz.bin";
// gdalImage demgridimg(demxyzPath);
// long dem_rowCount = demgridimg.height;
// long dem_ColCount = demgridimg.width;
// std::shared_ptr<double> demX = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
// std::shared_ptr<double> demY = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
// std::shared_ptr<double> demZ = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
//
// h_data.x_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
// h_data.y_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
// h_data.z_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
//
// for (long i = 0; i < dem_rowCount; i++) {
// for (long j = 0; j < dem_ColCount; j++) {
// h_data.x_mat[i * dem_ColCount + j] = demX.get()[i * dem_ColCount + j];
// h_data.y_mat[i * dem_ColCount + j] = demY.get()[i * dem_ColCount + j];
// h_data.z_mat[i * dem_ColCount + j] = demZ.get()[i * dem_ColCount + j];
// }
// }
// qDebug() << "dem row Count:\t" << dem_rowCount;
// qDebug() << "dem col Count:\t" << dem_ColCount;
//
// qDebug() << u8"成像网格读取结束";
// /** 3. 频率网格 **************************************************************************************************/
// qDebug() << u8"3.处理频率";
// float centerFreq = echoL0ds->getCenterFreq();
// float bandwidth = echoL0ds->getBandwidth();
// size_t freqpoints = echoL0ds->getPlusePoints();
// std::shared_ptr<float> freqlist(getFreqPoints_mallocHost_single(centerFreq - bandwidth / 2, centerFreq + bandwidth / 2, freqpoints), FreeCUDAHost);
// std::shared_ptr<float> minF(new float[gpspoints], delArrPtr);
// for (long i = 0; i < gpspoints; i++) {
// minF.get()[i] = freqlist.get()[0];
// }
// h_data.deltaF = bandwidth / (freqpoints - 1);
// qDebug() << "start Freq:\t" << centerFreq - bandwidth / 2;
// qDebug() << "end Freq:\t" << centerFreq + bandwidth / 2;
// qDebug() << "freq points:\t" << freqpoints;
// qDebug() << "delta freq:\t" << freqlist.get()[1] - freqlist.get()[0];
// qDebug() << u8"频率结束";
// /** 4. 初始化回波 **************************************************************************************************/
// qDebug() << u8"4.初始化回波";
// std::shared_ptr<std::complex<double>> echoData = echoL0ds->getEchoArr();
// size_t echosize = sizeof(cuComplex) * gpspoints * freqpoints;
// qDebug() << "echo data size (byte) :\t" << echosize;
// h_data.phdata = (cuComplex*)mallocCUDAHost(echosize);
// shared_complexPtrToHostCuComplex(echoData.get(), h_data.phdata, gpspoints * freqpoints);
//
// qDebug() << "Azimuth Points:\t" << gpspoints;
// qDebug() << "Range Points:\t" << freqpoints;
// qDebug() << u8"初始化回波结束";
// /** 5. 初始化图像 **************************************************************************************************/
// qDebug() << u8"5.初始化图像";
// h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * dem_rowCount * dem_ColCount);
// qDebug() << "Azimuth Points:\t" << gpspoints;
// qDebug() << "Range Points:\t" << freqpoints;
// qDebug() << u8"初始化图像结束";
//
// /** 6. 模型参数初始化 **************************************************************************************************/
// qDebug() << u8"6.模型参数初始化";
//
// h_data.nx = dem_ColCount;
// h_data.ny = dem_rowCount;
//
// h_data.Np = gpspoints;
// h_data.Freq = freqlist.get();
// h_data.minF = minF.get();
// h_data.Nfft = freqpoints;
// h_data.K = h_data.Nfft;
//
// h_data.R0 = echoL0ds->getRefPhaseRange();
// qDebug() << u8"模型参数结束";
//
//
// /** 7. 目标 **************************************************************************************************/
// double Tx = -2029086.618142, Ty = 4139594.934504, Tz = 4392846.782027;
// double Tslx = -2029086.618142, Tsly = 4139594.934504, Tslz = 4392846.782027;
// double p1 = 33.3, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
//
// /** 7. 构建回波 **************************************************************************************************/
// GPUDATA_single d_data;
// initGPUData_single(h_data, d_data);
//
// /** 8. 展示回波 **************************************************************************************************/
// {
// DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// cuComplex* h_ifftphdata = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// cuComplex* d_ifftphdata = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// CUDAIFFT(d_data.phdata, d_ifftphdata, d_data.Np, d_data.Nfft, d_data.Nfft);
// FFTShift1D(d_ifftphdata, d_data.Np, d_data.Nfft);
// DeviceToHost(h_ifftphdata, d_ifftphdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
// {
// for (long i = 0; i < d_data.Np; i++) {
// for (long j = 0; j < d_data.Nfft; j++) {
// ifftdata.get()[i * d_data.Nfft + j] =
// std::complex<double>(
// h_ifftphdata[i * d_data.Nfft + j].x,
// h_ifftphdata[i * d_data.Nfft + j].y);
// }
// }
// }
// testOutComplexDoubleArr(QString("echo_ifft_single.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
//
// FreeCUDADevice(d_ifftphdata);
// FreeCUDAHost(h_ifftphdata);
//
// }
// /** 9. 成像 **************************************************************************************************/
//
// // 计算maxWr需要先计算deltaF
// double deltaF = h_data.deltaF; // 从输入参数获取
// double maxWr = 299792458.0f / (2.0f * deltaF);
// qDebug() << "maxWr :\t" << maxWr;
// float* r_vec_host = (float*)mallocCUDAHost(sizeof(float) * h_data.Nfft);// new double[data.Nfft];
// const double step = maxWr / h_data.Nfft;
// const double start = -1 * h_data.Nfft / 2.0f * step;
// printf("nfft=%d\n", h_data.Nfft);
//
// for (int i = 0; i < h_data.Nfft; ++i) {
// r_vec_host[i] = start + i * step;
// }
//
// h_data.r_vec = r_vec_host;
// d_data.r_vec = h_data.r_vec;
//
// qDebug() << "r_vec [0]:\t" << h_data.r_vec[0];
// qDebug() << "r_vec step:\t" << step;
//
// bpBasic0CUDA_single(d_data, 0);
//
// DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * d_data.nx * d_data.ny);
// {
// DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
// std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
// {
// for (long i = 0; i < d_data.Np; i++) {
// for (long j = 0; j < d_data.Nfft; j++) {
// ifftdata.get()[i * d_data.Nfft + j] =
// std::complex<double>(
// h_data.phdata[i * d_data.Nfft + j].x,
// h_data.phdata[i * d_data.Nfft + j].y);
// }
// }
// }
// testOutComplexDoubleArr(QString("echo_ifft_BPBasic_single.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
//
// std::shared_ptr<std::complex<double>> im_finals(new std::complex<double>[d_data.nx * d_data.ny], delArrPtr);
// {
// for (long i = 0; i < d_data.ny; i++) {
// for (long j = 0; j < d_data.nx; j++) {
// im_finals.get()[i * d_data.nx + j] = std::complex<double>(
// h_data.im_final[i * d_data.nx + j].x,
// h_data.im_final[i * d_data.nx + j].y);
// }
// }
// }
// testOutComplexDoubleArr(QString("im_finals_single.bin"), im_finals.get(), d_data.ny, d_data.nx);
// }
//
//
//
//}
//
//
//int main(int argc, char* argv[]) {
//
// QApplication a(argc, argv);
//
// testBpImage();
//
//
// return 0;
//}
//
//