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

更新RFPC代码

pull/3/head
陈增辉 2025-01-15 18:48:43 +08:00
parent 46b3c47889
commit 4d9d16bb65
8 changed files with 155 additions and 512 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
GPUTool/GPUGarbage.cu Normal file
View File

@ -0,0 +1,17 @@
#include <time.h>
#include <iostream>
#include <memory>
#include <cmath>
#include <complex>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "BaseConstVariable.h"
#include "GPUGarbage.cuh"

38
GPUTool/GPUGarbage.cuh Normal file
View File

@ -0,0 +1,38 @@
#ifndef _GPUGARBAGE_H_
#define _GPUGARBAGE_H_
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "GPURFPC.cuh"
#endif

469
GPUTool/GPURFPC.cu
View File

@ -331,35 +331,6 @@ __global__ void CUDA_BillerInterpAntPattern(double* antpattern,
} }
} }
__global__ void CUDA_calculationEcho(double* sigma0, double* TransAnt, double* ReciveAnt,
double* localangle, double* R, double* slopeangle,
double nearRange, double Fs, double Pt, double lamda, long FreqIDmax,
cuComplex* echoArr, long* FreqID,
long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
double r = R[idx];
double amp = Pt * TransAnt[idx] * ReciveAnt[idx];
amp = amp * sigma0[idx];
amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // 反射强度
// 处理相位
double phi = (-4 * LAMP_CUDA_PI / lamda) * r;
cuComplex echophi = make_cuComplex(0, phi);
cuComplex echophiexp = cuCexpf(echophi);
double timeR = 2 * (r - nearRange) / LIGHTSPEED * Fs;
long timeID = floorf(timeR);
//if (timeID < 0 || timeID >= FreqIDmax) {
// timeID = 0;
// amp = 0;
//}
cuComplex echo = make_cuComplex(echophiexp.x, echophiexp.y);
echoArr[idx] = echo;
FreqID[idx] = timeID;
}
}
__global__ void CUDA_AntPatternInterpGain(double* anttheta, double* antphi, double* gain, __global__ void CUDA_AntPatternInterpGain(double* anttheta, double* antphi, double* gain,
@ -412,126 +383,6 @@ __global__ void CUDA_InterpSigma(
__global__ void CUDAKernel_RFPC_Caluation_R_Gain(
double antX, double antY, double antZ, // 天线的坐标
double* targetX, double* targetY, double* targetZ, long len, // 地面坐标
long* demCls,
double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量
double antXaxisX, double antXaxisY, double antXaxisZ, // 天线坐标系的X轴
double antYaxisX, double antYaxisY, double antYaxisZ,// 天线坐标系的Y轴
double antZaxisX, double antZaxisY, double antZaxisZ,// 天线坐标系的Z轴
double antDirectX, double antDirectY, double antDirectZ,// 天线的指向
double Pt,// 发射能量
double refPhaseRange,
double* TransAntpattern, double Transtarttheta, double Transstartphi, double Transdtheta, double Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
double* ReceiveAntpattern, double Receivestarttheta, double Receivestartphi, double Receivedtheta, double Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
double NearR, double FarR, // 距离范围
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
double* factorj, long freqnum,
double* outR, // 输出距离
//double* outAmp // 输出增益
double* PRFEcho_real, double* PRFEcho_imag, long prfid
) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
double tx = targetX[idx];
double ty = targetY[idx];
double tz = targetZ[idx];
double RstX = antX - tx; // 计算坐标矢量
double RstY = antY - ty;
double RstZ = antZ - tz;
double slopeX = demSlopeX[idx];
double slopeY = demSlopeY[idx];
double slopeZ = demSlopeZ[idx];
double RstR2 = RstX * RstX + RstY * RstY + RstZ * RstZ;
double RstR = sqrt(RstR2); // 矢量距离
//printf("antX=%f;antY=%f;antZ=%f;targetX=%f;targetY=%f;targetZ=%f;RstR=%.6f;diffR=%.6f;\n",antX,antY,antZ,targetX,targetY,targetZ,RstR, RstR - 9.010858499003178e+05);
if (RstR<NearR || RstR>FarR) {
}
else {
// 求解坡度
double slopR = sqrtf(slopeX * slopeX + slopeY * slopeY + slopeZ * slopeZ); //
double dotAB = RstX * slopeX + RstY * slopeY + RstZ * slopeZ;
double localangle = acosf(dotAB / (RstR * slopR)); // 局地入射角
double ampGain = 0;
// 求解天线方向图指向
CUDAVectorEllipsoidal antVector = GPU_SatelliteAntDirectNormal(
RstX, RstY, RstZ,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ
);
if (antVector.Rho > 0) {
// 发射方向图
double temptheta = antVector.theta * r2d;
double tempphi = antVector.phi * r2d;
double TansantPatternGain =
GPU_BillerInterpAntPattern(
TransAntpattern,
Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
temptheta, tempphi);
// 接收方向图
double antPatternGain = GPU_BillerInterpAntPattern(
ReceiveAntpattern,
Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
temptheta, tempphi);
// 计算
double sigma0 = 0;
{
long clsid = demCls[idx];
//printf("clsid=%d\n", clsid);
CUDASigmaParam tempsigma = sigma0Paramslist[clsid];
if (localangle < 0 || localangle >= LAMP_CUDA_PI / 2) {
sigma0 = 0;
}
else {}
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 {
double sigma = GPU_getSigma0dB(tempsigma, localangle);
sigma0 = powf(10.0, sigma / 10.0);// 后向散射系数
}
}
ampGain = TansantPatternGain * antPatternGain;
ampGain = ampGain / (powf(4 * LAMP_CUDA_PI, 2) * powf(RstR, 4)); // 反射强度
double outAmp_temp = ampGain * Pt * sigma0;
double tempR = RstR- refPhaseRange;
outR[idx] = RstR ;
for (long ii = 0; ii < freqnum; ii++) {
double phi= tempR * factorj[ii]; // 相位
// Eular; exp(ix)=cos(x)+isin(x)
double real = outAmp_temp * cos(phi); // 实部
double imag = outAmp_temp * sin(phi); // 虚部
atomicAdd(&PRFEcho_real[prfid * freqnum+ ii], real);
atomicAdd(&PRFEcho_imag[prfid * freqnum+ ii], imag);
}
}
else {
}
}
}
}
__global__ void CUDAKernel_RFPC_Computer_R_Gain( __global__ void CUDAKernel_RFPC_Computer_R_Gain(
double antX, double antY, double antZ, // 天线的坐标 double antX, double antY, double antZ, // 天线的坐标
double* targetX, double* targetY, double* targetZ, long len, // 地面坐标 double* targetX, double* targetY, double* targetZ, long len, // 地面坐标
@ -569,6 +420,8 @@ __global__ void CUDAKernel_RFPC_Computer_R_Gain(
//printf("antX=%f;antY=%f;antZ=%f;targetX=%f;targetY=%f;targetZ=%f;RstR=%.6f;diffR=%.6f;\n",antX,antY,antZ,targetX,targetY,targetZ,RstR, RstR - 9.010858499003178e+05); //printf("antX=%f;antY=%f;antZ=%f;targetX=%f;targetY=%f;targetZ=%f;RstR=%.6f;diffR=%.6f;\n",antX,antY,antZ,targetX,targetY,targetZ,RstR, RstR - 9.010858499003178e+05);
if (RstR<NearR || RstR>FarR) { if (RstR<NearR || RstR>FarR) {
outAmp[idx] = 0;
outR[idx] = 0;
} }
else { else {
// 求解坡度 // 求解坡度
@ -638,88 +491,51 @@ __global__ void CUDAKernel_RFPC_Computer_R_Gain(
} }
__global__ void CUDAKernel_PRF_FreqEcho(double* temp_R,double factorj,double* temp_real,double* temp_imag,long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
double phi = factorj * temp_R[idx];
temp_real[idx] = cos(phi);
temp_imag[idx] = sin(phi);
}
}
__global__ void CUDAKernel_PRF_GeneratorEcho(float* Rarr, float* ampArr,
long pixelcount,
__global__ void CUDAKernel_PRF_CalFreqEcho( float f0, float dfreq,long freqnum,
double* Rarr, double* ampArr, long pixelcount, double* echo_real,double* echo_imag, long prfid)
double* factorj, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < freqnum) {
double fatorj = factorj[idx];
double phi = 0;
double amptemp = 0;
cuComplex tempfreqEcho = PRFEcho[prfid * freqnum + idx];
for (long i = 0; i < pixelcount; i++) { // 区域积分
//phi = (R = R - (floor(R / lamda) - 1) * lamda)* fatorj; // 相位
double phi = Rarr[i] * factorj[idx]; // 相位
amptemp = ampArr[i];
//printf("amp=%f\n", amptemp);
// Eular; exp(ix)=cos(x)+isin(x)
tempfreqEcho.x = tempfreqEcho.x + amptemp * cos(phi); // 实部
tempfreqEcho.y = tempfreqEcho.y + amptemp * sin(phi); // 虚部
//printf("freqid=%d;fatorj=%.12f;d_R=%.10f;phi=%.10f;echo=complex(%.5f,%.5f)\n", idx, fatorj, Rarr[i], phi, tempfreqEcho.x, tempfreqEcho.y);
}
PRFEcho[prfid*freqnum+idx] = tempfreqEcho;
}
}
__global__ void CUDAKernel_PRF_GeneratorEcho(float* Rarr, float* ampArr, long blocknum, long pixelcount, double* factorj, long freqnum,
double nearR, double farR, double* echo_real, double* echo_imag, long prfid) //11
{ {
//// 假定共享内存大小为49152 byte //// 假定共享内存大小为49152 byte
//// 假定每个Block 线程数大小为 32 //// 假定每个Block 线程数大小为 32
__shared__ float s_R[GPU_SHARE_MEMORY]; // 距离 32*12 * 8= 49.2kb __shared__ float s_R[GPU_SHARE_MEMORY]; // 距离 32*12 * 8= 49.2kb
__shared__ float s_Amp[GPU_SHARE_MEMORY]; // 振幅 3072 * 8= 49.2kb 49.2*2 = 98.4 < 100 KB __shared__ float s_Amp[GPU_SHARE_MEMORY]; // 振幅 3072 * 8= 49.2kb 49.2*2 = 98.4 < 100 KB
const int bid = blockIdx.x; // 获取 grid网格编号ID int idx = blockIdx.x * blockDim.x + threadIdx.x;; // 获取当前的线程编码
const int tid = threadIdx.x;// 获取 单个 block 中的线程ID int tid = threadIdx.x;// 获取 单个 block 中的线程ID
const int startPIX = bid * GPU_SHARE_STEP; const long startPIX = idx * GPU_SHARE_STEP; // 计算偏移
int curthreadidx = 0; int curthreadidx = 0;
for (long i = 0; i < GPU_SHARE_STEP; i++) { for (long i = 0; i < GPU_SHARE_STEP; i++) {
curthreadidx = i * blockDim.x + tid; // 计算分块 curthreadidx = i * BLOCK_SIZE + tid; // 计算分块
s_R[curthreadidx] = (startPIX + i) < pixelcount ? Rarr[startPIX + i] : 0.0; s_R[curthreadidx] = (startPIX + i) < pixelcount ? Rarr[startPIX + i] : 0.0;
s_Amp[curthreadidx] = (startPIX + i) < pixelcount ? ampArr[startPIX + i] : 0.0; s_Amp[curthreadidx] = (startPIX + i) < pixelcount ? ampArr[startPIX + i] : 0.0;
} }
//__syncthreads(); // 确定所有待处理数据都已经进入程序中 //__syncthreads(); // 确定所有待处理数据都已经进入程序中
long freqnumblock = (freqnum + 32 - 1) / 32; //16
if (startPIX < pixelcount) { // 存在可能处理的计算 if (startPIX < pixelcount) { // 存在可能处理的计算
double temp_real = 0; float temp_real = 0;
double temp_imag = 0; float temp_imag = 0;
double factorjTemp = 0; float factorjTemp = 0;
double temp_phi = 0; float temp_phi = 0;
double temp_amp = 0; float temp_amp = 0;
long dataid = 0; long dataid = 0;
curthreadidx = 0; curthreadidx = 0;
for (long i = 0; i < freqnumblock; i++) { for (long fid = 0; fid < freqnum; fid++) {
curthreadidx = tid * freqnumblock + i; // 获取当前频率 factorjTemp = RFPCPIDIVLIGHT *(f0+ fid* dfreq);
if (curthreadidx < freqnum) { // 存在频率 for (long j = 0; j < GPU_SHARE_STEP; j++) {
factorjTemp = factorj[curthreadidx]; dataid = j * BLOCK_SIZE + tid; //
for (long j = 0; j < GPU_SHARE_STEP; j++) { temp_phi = s_R[dataid] * factorjTemp;
dataid = j * blockDim.x + tid; // 数据编辑 temp_amp = s_Amp[dataid];
temp_phi = s_R[dataid] * factorjTemp;
temp_amp = s_Amp[dataid];
temp_real = temp_real + temp_amp * cos(temp_phi); temp_real += temp_amp * cosf(temp_phi);
temp_imag = temp_imag + temp_amp * sin(temp_phi); temp_imag += temp_amp * sinf(temp_phi);
}
atomicAdd(&echo_real[prfid * freqnum + curthreadidx], temp_real); // 更新实部
atomicAdd(&echo_imag[prfid * freqnum + curthreadidx], temp_imag); // 更新虚部
} }
atomicAdd(&echo_real[prfid * freqnum + fid], temp_real); // 更新实部
atomicAdd(&echo_imag[prfid * freqnum + fid], temp_imag); // 更新虚部
} }
} }
} }
@ -731,213 +547,6 @@ __global__ void CUDAKernel_PRF_GeneratorEcho(float* Rarr, float* ampArr, long bl
/** 对外封装接口 *******************************************************************************************************/ /** 对外封装接口 *******************************************************************************************************/
extern "C" void SatelliteAntDirectNormal(double* RstX, double* RstY, double* RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ,
double* thetaAnt, double* phiAnt,
long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_SatelliteAntDirectNormal << <numBlocks, blockSize >> > (RstX, RstY, RstZ,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ,
thetaAnt, phiAnt
, len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RFPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void AntPatternInterpGain(double* anttheta, double* antphi, double* gain,
double* antpattern, double starttheta, double startphi, double dtheta, double dphi, int thetapoints, int phipoints, long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
//printf("\nCUDA_RFPC_SiglePRF blockSize:%d ,numBlock:%d\n", blockSize, numBlocks);
CUDA_AntPatternInterpGain << <numBlocks, blockSize >> > (anttheta, antphi, gain,
antpattern,
starttheta, startphi, dtheta, dphi, thetapoints, phipoints,
len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RFPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void calculationEcho(double* sigma0, double* TransAnt, double* ReciveAnt,
double* localangle, double* R, double* slopeangle,
double nearRange, double Fs, double pt, double lamda, long FreqIDmax,
cuComplex* echoAmp, long* FreqID,
long len)
{
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_calculationEcho << <numBlocks, blockSize >> > (sigma0, TransAnt, ReciveAnt,
localangle, R, slopeangle,
nearRange, Fs, pt, lamda, FreqIDmax,
echoAmp, FreqID,
len);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RFPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void CUDAInterpSigma(
long* demcls, double* sigmaAmp, double* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) {// 地表覆盖类型-sigma插值对应函数-ulaby
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_InterpSigma << <numBlocks, blockSize >> > (
demcls, sigmaAmp, localanglearr, len,
sigma0Paramslist, sigmaparamslistlen
);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_RFPC_SiglePRF CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void CUDARFPC_Caluation_R_Gain(
double antX, double antY, double antZ, // 天线的坐标
double* targetX, double* targetY, double* targetZ, long TargetPixelNumber, // 地面坐标
long* demCls,
double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量
double antXaxisX, double antXaxisY, double antXaxisZ, // 天线坐标系的X轴
double antYaxisX, double antYaxisY, double antYaxisZ,// 天线坐标系的Y轴
double antZaxisX, double antZaxisY, double antZaxisZ,// 天线坐标系的Z轴
double antDirectX, double antDirectY, double antDirectZ,// 天线的指向
double Pt,// 发射能量
double refPhaseRange,
double* TransAntpattern, double Transtarttheta, double Transstartphi, double Transdtheta, double Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
double* ReceiveAntpattern, double Receivestarttheta, double Receivestartphi, double Receivedtheta, double Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
double NearR, double FarR, // 距离范围
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
double* factorj, long freqnum,
double* outR, // 输出距离
//double* outAmp // 输出增益
double* PRFEcho_real, double* PRFEcho_imag, long prfid
)
{
int blockSize = 256; // 每个块的线程数
int numBlocks = (TargetPixelNumber + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDAKernel_RFPC_Caluation_R_Gain << <numBlocks, blockSize >> > (
antX, antY, antZ,
targetX, targetY, targetZ, TargetPixelNumber,
demCls,
demSlopeX, demSlopeY, demSlopeZ,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ,
Pt,
refPhaseRange,
TransAntpattern,
Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
ReceiveAntpattern,
Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
NearR, FarR,
sigma0Paramslist, sigmaparamslistlen,
factorj, freqnum,
outR,
//outAmp
PRFEcho_real, PRFEcho_imag, prfid
);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDARFPC_Caluation_R_Gain CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void CUDA_PRF_CalFreqEcho(
double* Rarr, double* ampArr, long pixelcount,
double* factorj, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid)
{
int blockSize = 256; // 每个块的线程数
int numBlocks = (freqnum + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
CUDAKernel_PRF_CalFreqEcho << <numBlocks, blockSize >> > (
Rarr, ampArr, pixelcount,
factorj, freqnum,
dx,nearR,
PRFEcho, prfid
);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_PRF_CalFreqEcho CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}
extern "C" void CUDA_PRF_GeneratorEcho(cublasHandle_t handle,double* Rarr, double* ampArr, long pixelcount, double* factorj, long freqnum, double nearR, double farR, double* echo_real, double* echo_imag, long prfid)
{
//cublasHandle_t handle;
//cublasStatus_t status = cublasCreate(&handle);
long blocknum = pixelcount / GPU_SHARE_MEMORY + 1;
int blockSize = 256; // 每个块的线程数
int numBlocks = (pixelcount + GPU_SHARE_MEMORY - 1) / GPU_SHARE_MEMORY; // 网格数量
CUDAKernel_PRF_GeneratorEcho << <numBlocks, blockSize >> > (Rarr, ampArr, blocknum, pixelcount,
factorj, freqnum,
nearR, farR,
echo_real, echo_imag, prfid);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_PRF_GeneratorEcho CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
//cublasDestroy(handle);
}
extern "C" void CUDA_RFPC_MainBlock( extern "C" void CUDA_RFPC_MainBlock(
double* antX, double* antY, double* antZ, // 天线的坐标 double* antX, double* antY, double* antZ, // 天线的坐标
double* antXaxisX, double* antXaxisY, double* antXaxisZ, // 天线坐标系的X轴 double* antXaxisX, double* antXaxisY, double* antXaxisZ, // 天线坐标系的X轴
@ -945,7 +554,7 @@ extern "C" void CUDA_RFPC_MainBlock(
double* antZaxisX, double* antZaxisY, double* antZaxisZ,// 天线坐标系的Z轴 double* antZaxisX, double* antZaxisY, double* antZaxisZ,// 天线坐标系的Z轴
double* antDirectX, double* antDirectY, double* antDirectZ,// 天线的指向 double* antDirectX, double* antDirectY, double* antDirectZ,// 天线的指向
long PRFCount, // 脉冲数 long PRFCount, // 脉冲数
double* freqpoints,double* factorj ,long freqnum,// 频率数 float f0, float dfreq, long freqnum, // 频率数
double* targetX, double* targetY, double* targetZ, long TargetPixelNumber, // 地面坐标 double* targetX, double* targetY, double* targetZ, long TargetPixelNumber, // 地面坐标
long* demCls, // 地表类别 long* demCls, // 地表类别
double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量 double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量
@ -962,13 +571,12 @@ extern "C" void CUDA_RFPC_MainBlock(
long blocknum = 0; long blocknum = 0;
long pixelcount=TargetPixelNumber; long pixelcount=TargetPixelNumber;
int numBlocks = 0;
for(long pid=0;pid<PRFCount;pid++){ for(long pid=0;pid<PRFCount;pid++){
numBlocks = (TargetPixelNumber + BLOCK_SIZE - 1) / BLOCK_SIZE; // 根据 pixelcount 计算网格大小
int blockSize = 256; // 每个块的线程数
int numBlocks = (TargetPixelNumber + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
CUDAKernel_RFPC_Computer_R_Gain<<<numBlocks , blockSize >>>( CUDAKernel_RFPC_Computer_R_Gain<<<numBlocks , BLOCK_SIZE >>>(
antX[pid], antY[pid], antZ[pid], antX[pid], antY[pid], antZ[pid],
targetX, targetY, targetZ, TargetPixelNumber, targetX, targetY, targetZ, TargetPixelNumber,
demCls, demCls,
@ -990,19 +598,16 @@ extern "C" void CUDA_RFPC_MainBlock(
temp_amp temp_amp
//out_echoReal, out_echoImag, pid // 输出振幅 //out_echoReal, out_echoImag, pid // 输出振幅
); );
cudaDeviceSynchronize();
blocknum = (pixelcount+ GPU_SHARE_STEP-1) / GPU_SHARE_STEP ; blocknum = (pixelcount + GPU_SHARE_STEP - 1) / GPU_SHARE_STEP;
blockSize = 32; // 每个块的线程数 numBlocks = (blocknum + BLOCK_SIZE - 1) / BLOCK_SIZE; // 网格数量
numBlocks = (pixelcount + GPU_SHARE_STEP - 1) / GPU_SHARE_STEP; // 网格数量 CUDAKernel_PRF_GeneratorEcho << <numBlocks, BLOCK_SIZE >> >
(temp_R, temp_amp, pixelcount,
CUDAKernel_PRF_GeneratorEcho << <numBlocks, blockSize >> > (temp_R, temp_amp, blocknum, pixelcount, f0,dfreq,freqnum,
factorj, freqnum,
NearR, FarR,
out_echoReal, out_echoImag, pid); out_echoReal, out_echoImag, pid);
cudaDeviceSynchronize();
}
}
cudaDeviceSynchronize();
} }

101
GPUTool/GPURFPC.cuh
View File

@ -8,6 +8,12 @@
#include <cublas_v2.h> #include <cublas_v2.h>
#include <cuComplex.h> #include <cuComplex.h>
#define RFPCPIDIVLIGHT -4*PI/(LIGHTSPEED/1e9)
#define GPU_SHARE_MEMORY 5888
#define GPU_SHARE_STEP 23
extern "C" struct CUDASigmaParam { extern "C" struct CUDASigmaParam {
double p1; double p1;
@ -44,23 +50,6 @@ extern __global__ void CUDA_InterpSigma(
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen); CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen);
extern __global__ void CUDAKernel_PRF_CalFreqEcho(
double* Rarr, double* ampArr, long pixelcount,
double* factorj, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid);
extern __device__ cuComplex GPU_calculationEcho(double sigma0, double TransAnt, double ReciveAnt,
double localangle, double R, double slopeangle, double Pt, double lamda);
extern __global__ void CUDA_calculationEcho(double* sigma0, double* TransAnt, double* ReciveAnt,
double* localangle, double* R, double* slopeangle,
double nearRange, double Fs, double Pt, double lamda, long FreqIDmax,
cuComplex* echoArr, long* FreqID,
long len);
extern __global__ void CUDA_BillerInterpAntPattern(double* antpattern, extern __global__ void CUDA_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi, double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints, long thetapoints, long phipoints,
@ -68,51 +57,12 @@ extern __global__ void CUDA_BillerInterpAntPattern(double* antpattern,
long len); long len);
extern __global__ void CUDA_SatelliteAntDirectNormal(double* RstX, double* RstY, double* RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ,
double* thetaAnt, double* phiAnt
, long len);
extern __global__ void CUDAKernel_RFPC_Computer_R_Gain(
extern __global__ void CUDAKernel_PRF_GeneratorEcho(float* Rarr, float* ampArr, long blocknum, long pixelcount, double* factorj, long freqnum,
double nearR, double farR, double* echo_real, double* echo_imag, long prfid);
extern "C" void SatelliteAntDirectNormal(double* RstX, double* RstY, double* RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ,
double* thetaAnt, double* phiAnt, long len);
extern "C" void AntPatternInterpGain(double* anttheta, double* antphi, double* gain,
double* antpattern,
double starttheta, double startphi, double dtheta, double dphi, int thetapoints, int phipoints,
long len);
extern "C" void calculationEcho(double* sigma0, double* TransAnt, double* ReciveAnt,
double* localangle, double* R, double* slopeangle,
double nearRange, double Fs, double pt, double lamda, long FreqIDmax,
cuComplex* echoAmp, long* FreqID,
long len);
extern "C" void CUDAInterpSigma(
long* demcls, double* sigmaAmp, double* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen);
// 计算坐标的 距离、增益
extern "C" void CUDARFPC_Caluation_R_Gain(
double antX, double antY, double antZ, // 天线的坐标 double antX, double antY, double antZ, // 天线的坐标
double* targetX, double* targetY, double* targetZ, long TargetPixelNumber, // 地面坐标 double* targetX, double* targetY, double* targetZ, long len, // 地面坐标
long* demCls, long* demCls,
double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量 double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量
double antXaxisX, double antXaxisY, double antXaxisZ, // 天线坐标系的X轴 double antXaxisX, double antXaxisY, double antXaxisZ, // 天线坐标系的X轴
@ -125,26 +75,31 @@ extern "C" void CUDARFPC_Caluation_R_Gain(
double* ReceiveAntpattern, double Receivestarttheta, double Receivestartphi, double Receivedtheta, double Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图 double* ReceiveAntpattern, double Receivestarttheta, double Receivestartphi, double Receivedtheta, double Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
double NearR, double FarR, // 距离范围 double NearR, double FarR, // 距离范围
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图 CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
double* factorj, long freqnum, float* outR, // 输出距离
double* outR, // 输出距离 float* outAmp
//double* outAmp // 输出增益
double* PRFEcho_real, double* PRFEcho_imag, long prfid
); );
extern "C" void CUDA_PRF_CalFreqEcho( extern __global__ void CUDAKernel_PRF_GeneratorEcho(float* Rarr, float* ampArr,
double* Rarr, double* ampArr, long pixelcount, long pixelcount,
double* factorj, long freqnum, float f0, float dfreq, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid);
extern "C" void CUDA_PRF_GeneratorEcho(cublasHandle_t handle,float* Rarr, float* ampArr, long pixelcount,
double* factorj, long freqnum,
double nearR, double farR,
double* echo_real, double* echo_imag, long prfid); double* echo_real, double* echo_imag, long prfid);
extern "C" void CUDA_RFPC_MainBlock( extern "C" void CUDA_RFPC_MainBlock(
double* antX, double* antY, double* antZ, // 天线的坐标 double* antX, double* antY, double* antZ, // 天线的坐标
double* antXaxisX, double* antXaxisY, double* antXaxisZ, // 天线坐标系的X轴 double* antXaxisX, double* antXaxisY, double* antXaxisZ, // 天线坐标系的X轴
@ -152,7 +107,7 @@ extern "C" void CUDA_RFPC_MainBlock(
double* antZaxisX, double* antZaxisY, double* antZaxisZ,// 天线坐标系的Z轴 double* antZaxisX, double* antZaxisY, double* antZaxisZ,// 天线坐标系的Z轴
double* antDirectX, double* antDirectY, double* antDirectZ,// 天线的指向 double* antDirectX, double* antDirectY, double* antDirectZ,// 天线的指向
long PRFCount, // 脉冲数 long PRFCount, // 脉冲数
double* freqpoints, double* factorj ,long freqnunm,// 频率数 float f0, float dfreq, long freqnum,// 频率数
double* targetX, double* targetY, double* targetZ, long TargetPixelNumber, // 地面坐标 double* targetX, double* targetY, double* targetZ, long TargetPixelNumber, // 地面坐标
long* demCls, // 地表类别 long* demCls, // 地表类别
double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量 double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量

30
GPUTool/GPUTool.cuh
View File

@ -17,9 +17,7 @@
#define BLOCK_SIZE 256 #define BLOCK_SIZE 256
// 默认显存分布
#define GPU_SHARE_MEMORY 5632
#define GPU_SHARE_STEP 176
@ -55,6 +53,18 @@ extern __device__ float GPU_dotVector(CUDAVector A, CUDAVector B);
extern __device__ float GPU_CosAngle_VectorA_VectorB(CUDAVector A, CUDAVector B); extern __device__ float GPU_CosAngle_VectorA_VectorB(CUDAVector A, CUDAVector B);
// 定义全局函数 // 定义全局函数
extern __global__ void CUDA_DistanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* R, long len); extern __global__ void CUDA_DistanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* R, long len);
extern __global__ void CUDA_B_DistanceA(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long len); extern __global__ void CUDA_B_DistanceA(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long len);
@ -66,6 +76,20 @@ extern __global__ void CUDA_D_sin(double* y, double* X, int n);
extern __global__ void CUDA_D_cos(double* y, double* X, int n); extern __global__ void CUDA_D_cos(double* y, double* X, int n);
// 误差处理函数 // 误差处理函数
extern "C" void checkCudaError(cudaError_t err, const char* msg); extern "C" void checkCudaError(cudaError_t err, const char* msg);

2
RasterProcessTool.vcxproj
View File

@ -210,6 +210,7 @@
<ClInclude Include="BaseTool\RasterToolBase.h" /> <ClInclude Include="BaseTool\RasterToolBase.h" />
<ClInclude Include="BaseTool\SARSimulationImageL1.h" /> <ClInclude Include="BaseTool\SARSimulationImageL1.h" />
<ClInclude Include="BaseTool\stdafx.h" /> <ClInclude Include="BaseTool\stdafx.h" />
<ClInclude Include="GPUTool\GPUGarbage.cuh" />
<ClInclude Include="GPUTool\GPURFPC.cuh" /> <ClInclude Include="GPUTool\GPURFPC.cuh" />
<ClInclude Include="GPUTool\GPUTBPImage.cuh" /> <ClInclude Include="GPUTool\GPUTBPImage.cuh" />
<ClInclude Include="GPUTool\GPUTool.cuh" /> <ClInclude Include="GPUTool\GPUTool.cuh" />
@ -229,6 +230,7 @@
<None Include="cpp.hint" /> <None Include="cpp.hint" />
</ItemGroup> </ItemGroup>
<ItemGroup> <ItemGroup>
<CudaCompile Include="GPUTool\GPUGarbage.cu" />
<CudaCompile Include="GPUTool\GPURFPC.cu"> <CudaCompile Include="GPUTool\GPURFPC.cu">
<GenerateRelocatableDeviceCode Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</GenerateRelocatableDeviceCode> <GenerateRelocatableDeviceCode Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</GenerateRelocatableDeviceCode>
<CodeGeneration Condition="'$(Configuration)|$(Platform)'=='Release|x64'">compute_86,sm_86;%(CodeGeneration)</CodeGeneration> <CodeGeneration Condition="'$(Configuration)|$(Platform)'=='Release|x64'">compute_86,sm_86;%(CodeGeneration)</CodeGeneration>

6
SimulationSAR/RFPCProcessCls.cpp
View File

@ -493,6 +493,10 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
} }
testOutAmpArr("freqlist.bin", (double*)(freqPtr.get()), freqnum, 1); testOutAmpArr("freqlist.bin", (double*)(freqPtr.get()), freqnum, 1);
float f0 = float(freqlist[0] / 1e9);
float dfreq = float((freqlist[1] - freqlist[0]) / 1e9);
long double imageStarttime = 0; long double imageStarttime = 0;
@ -917,7 +921,7 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
antdirectx, antdirecty, antdirectz,// 天线的指向 antdirectx, antdirecty, antdirectz,// 天线的指向
templine, // 脉冲数 templine, // 脉冲数
//h_freqlist, h_factorj, PlusePoint,// 频率数 //h_freqlist, h_factorj, PlusePoint,// 频率数
d_freqlist, d_factorj, PlusePoint,// 频率数 f0, dfreq, PlusePoint,// 频率数
d_dem_x, d_dem_y, d_dem_z, pixelcount, // 地面坐标 d_dem_x, d_dem_y, d_dem_z, pixelcount, // 地面坐标
d_demcls, d_demcls,
d_demsloper_x, d_demsloper_y, d_demsloper_z, // 地表坡度矢量 d_demsloper_x, d_demsloper_y, d_demsloper_z, // 地表坡度矢量

2
SimulationSAR/SARSatelliteSimulationAbstractCls.cpp
View File

@ -207,8 +207,6 @@ QVector<double> AbstractSARSatelliteModel::getFreqList()
double farR = this->getFarRange(); double farR = this->getFarRange();
// 计算分辨率 // 计算分辨率
long long centerFreq_long = long long(centerFreq); long long centerFreq_long = long long(centerFreq);
long long bandwidth_long = long long(bandwidth); long long bandwidth_long = long long(bandwidth);
long long bandhalf = ceil(bandwidth_long / 2); long long bandhalf = ceil(bandwidth_long / 2);