增加GPU的RTPC代码

GPUTool.cu

@@ -19,7 +19,6 @@
 #define LAMP_CUDA_PI 3.141592653589793238462643383279
 
 
-
 // 定义参数
  __device__  cuComplex cuCexpf(cuComplex x)
 {
@@ -27,6 +26,134 @@
 	return make_cuComplex(factor * cos(x.y), factor * sin(x.y));
 }
 
+// 定义仿真所需参数
+
+ __device__ float GPU_getSigma0dB(CUDASigmaParam param,float theta) {
+	 return  param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
+}
+
+ __device__ CUDAVector GPU_VectorAB(CUDAVector A, CUDAVector B) {
+	 CUDAVector C;
+	 C.x = B.x - A.x;
+	 C.y = B.y - A.y;
+	 C.z = B.z - A.z;
+	 return C;
+}
+
+ __device__ float GPU_VectorNorm2(CUDAVector A) {
+	 return sqrtf(A.x * A.x + A.y * A.y + A.z * A.z);
+ }
+
+ __device__ float GPU_dotVector(CUDAVector A, CUDAVector B) {
+	 return A.x * B.x + A.y * B.y + A.z * B.z;
+}
+
+ __device__ float GPU_CosAngle_VectorA_VectorB(CUDAVector A, CUDAVector B) {
+	 return GPU_dotVector(A, B) / (GPU_VectorNorm2(A)*GPU_VectorNorm2(B));
+}
+
+ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(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 AntXaxisX = antXaxisX;
+	 float AntXaxisY = antXaxisY;
+	 float AntXaxisZ = antXaxisZ;
+	 float AntYaxisX = antYaxisX;
+	 float AntYaxisY = antYaxisY;
+	 float AntYaxisZ = antYaxisZ;
+	 float AntZaxisX = antZaxisX;
+	 float AntZaxisY = antZaxisY;
+	 float AntZaxisZ = antZaxisZ;
+	 // 天线指向在天线坐标系下的值
+	 float Xant = (Xst * (AntYaxisY * AntZaxisZ - AntYaxisZ * AntZaxisY) + Xst * (AntXaxisZ * AntZaxisY - AntXaxisY * AntZaxisZ) + Xst * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
+	 float Yant = (Yst * (AntYaxisZ * AntZaxisX - AntYaxisX * AntZaxisZ) + Yst * (AntXaxisX * AntZaxisZ - AntXaxisZ * AntZaxisX) + Yst * (AntYaxisX * AntXaxisZ - AntXaxisX * AntYaxisZ)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
+	 float Zant = (Zst * (AntYaxisX * AntZaxisY - AntYaxisY * AntZaxisX) + Zst * (AntXaxisY * AntZaxisX - AntXaxisX * AntZaxisY) + Zst * (AntXaxisX * AntYaxisY - AntYaxisX * AntXaxisY)) / (AntXaxisX * (AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * (AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * (AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
+	 // 计算theta 与 phi
+	 float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
+	 float ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
+	 float YsinTheta = Yant / sinf(ThetaAnt);
+	 float PhiAnt = (YsinTheta / abs(YsinTheta)) * acosf(Xant / (Norm * sinf(ThetaAnt)));
+	 result.theta = ThetaAnt;
+	 result.phi = PhiAnt;
+	 result.pho = Norm;
+	 return result;
+}
+
+ /**
+天线方向图插值方法，以双线性插值算法为基础，由theta与phi组合得到的矩阵图为基础数据，通过插值计算的方法获取目标点的数据。
+其中行是theta、列是phi
+*/
+ __device__ float GPU_BillerInterpAntPattern(float* antpattern,
+	 float starttheta, float startphi, float dtheta, float dphi,
+	 long thetapoints, long phipoints,
+	 float searththeta, float searchphi) {
+	 float stheta = searththeta;
+	 float sphi = searchphi;
+	 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);
+		 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;
+	 }
+ }
+
+
+ __device__ cuComplex  GPU_calculationEcho(float sigma0, float TransAnt, float ReciveAnt,
+	 float localangle, float R, float slopeangle,float Pt, float lamda ) {
+	 float r = R;
+	 float amp = Pt * TransAnt * ReciveAnt;
+	 amp = amp * sigma0;
+	 amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // 反射强度
+	 float phi = (-4 * LAMP_CUDA_PI / lamda) * r;
+	 cuComplex echophi = make_cuComplex(0, phi);
+	 cuComplex echophiexp = cuCexpf(echophi);
+	 cuComplex echo;
+	 echo.x = echophiexp.x * amp;
+	 echo.y = echophiexp.y * amp;
+	 return echo;
+ }
+
 __global__ void CUDA_DistanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz,float *R, long len) {
 	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
 	if (idx < len) {
@@ -144,6 +271,127 @@ __global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* Reci
 	}
 }
 
+__global__ void CUDA_BillerInterpAntPattern(float* antpattern,
+	float starttheta, float startphi, float dtheta, float dphi,
+	long thetapoints, long phipoints,
+	float* searththeta, float* searchphi,float* searchantpattern,
+	long len) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < len) {
+		float stheta = searththeta[idx];
+		float sphi = searchphi[idx];
+		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)
+		{
+			searchantpattern[idx] = 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);
+			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));
+			searchantpattern[idx] = GainValue;
+		}
+	}
+}
+
+
+
+__global__ void CUDA_RTPC_Kernel(
+	float antPx, float antPy, float antPz,// 天线坐标
+	float  antXaxisX, float  antXaxisY, float  antXaxisZ, // 天线坐标系
+	float  antYaxisX, float  antYaxisY, float  antYaxisZ, //
+	float  antZaxisX, float  antZaxisY, float  antZaxisZ,
+	float  antDirectX, float  antDirectY, float  antDirectZ,// 天线指向
+	float* demx, float* demy, float* demz, long* demcls, // 地面坐标
+	float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,// 地面坡度
+	float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, long Tthetapoints, long Tphipoints,// 天线方向图相关
+	float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, long Rthetapoints, long Rphipoints,// 天线方向图相关
+	float lamda, float fs, float nearrange, float Pt, long Freqnumbers, // 参数
+	CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby
+	cuComplex* outecho, long* d_echoAmpFID,
+	long len
+ ) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < len) {
+		long clsid = demcls[idx];
+		CUDAVector Rs{ antPx,antPy,antPz };
+		CUDAVector Rt{ demx[idx],demy[idx],demz[idx] };
+		CUDAVector Rst{ Rs.x - Rt.x,Rs.y - Rt.y,Rs.z - Rt.z };
+		CUDAVector Vslope{ demslopex[idx],demslopey[idx],demslopez[idx] };
+		float R = GPU_VectorNorm2(Rst); // 斜距
+
+		CUDAVectorEllipsoidal Rtanttheta = GPU_SatelliteAntDirectNormal( // 地面目标在天线的位置
+			Rst.x, Rst.y, Rst.z,
+			antXaxisX, antXaxisY, antXaxisZ,
+			antYaxisX, antYaxisY, antYaxisZ,
+			antZaxisX, antZaxisY, antZaxisZ,
+			antDirectX, antDirectY, antDirectZ);
+
+		float localangle=GPU_CosAngle_VectorA_VectorB(Rst, Vslope); // 距地入射角
+		float sigma = GPU_getSigma0dB(sigma0Paramslist[clsid], localangle); 
+		sigma = powf(10.0, sigma / 10.0);// 后向散射系数
+
+		// 发射方向图
+		float transPattern = GPU_BillerInterpAntPattern(Tantpattern,
+			Tstarttheta, Tstartphi, Tdtheta, Tdphi,Tthetapoints, Tphipoints,
+			Rtanttheta.theta, Rtanttheta.phi);
+		
+		// 接收方向图
+		float receivePattern = GPU_BillerInterpAntPattern(Rantpattern,
+			Rstarttheta, Rstartphi, Rdtheta, Rdphi, Rthetapoints, Rphipoints,
+			Rtanttheta.theta, Rtanttheta.phi);
+		// 计算振幅、相位
+		float amp = Pt * transPattern * receivePattern * sigma / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4));
+		float phi = (-4 * LAMP_CUDA_PI / lamda) * R;
+
+		// 构建回波
+		cuComplex echophi = make_cuComplex(0, phi);
+		cuComplex echophiexp = cuCexpf(echophi);
+		float timeR = 2 * (R - nearrange) / LIGHTSPEED * fs;
+		long timeID = floorf(timeR);
+		if (timeID < 0 || timeID >= Freqnumbers) {
+			timeID = 0;
+			amp = 0;
+			
+		}
+		else {}
+
+		cuComplex echo;
+		echo.x = echophiexp.x * amp;
+		echo.y = echophiexp.y * amp;
+		outecho[idx] = echo;
+		d_echoAmpFID[idx] = timeID;
+	}
+}
 
 
 //错误提示
@@ -156,8 +404,10 @@ void checkCudaError(cudaError_t err, const char* msg) {
 }
  
 // 主机参数内存声明
-extern "C"  void mallocCUDAHost(void* ptr, long memsize) {
+extern "C"  void* mallocCUDAHost( long memsize) {
+	void* ptr;
 	cudaMallocHost(&ptr, memsize);
+	return ptr;
 }
 
 // 主机参数内存释放
@@ -166,8 +416,10 @@ extern "C"  void FreeCUDAHost(void* ptr) {
 }
 
 // GPU参数内存声明
-extern "C" void mallocCUDADevice(void* ptr, long memsize) {
+extern "C" void* mallocCUDADevice( long memsize) {
+	void* ptr;
 	cudaMalloc(&ptr, memsize);
+	return ptr;
 }
 
 // GPU参数内存释放
@@ -268,7 +520,6 @@ extern "C" void calculationEcho(float* sigma0,float* TransAnt,float* ReciveAnt,
 		echoAmp, FreqID,
 		len);
 	cudaDeviceSynchronize();
- 
 }
 
 
@@ -276,4 +527,44 @@ extern "C" void calculationEcho(float* sigma0,float* TransAnt,float* ReciveAnt,
  
 
 
+extern "C" void CUDA_RTPC_SiglePRF(
+	float antPx, float antPy, float antPZ,// 天线坐标
+	float  antXaxisX, float  antXaxisY, float  antXaxisZ, // 天线坐标系
+	float  antYaxisX, float  antYaxisY, float  antYaxisZ, //
+	float  antZaxisX, float  antZaxisY, float  antZaxisZ,
+	float  antDirectX, float  antDirectY, float  antDirectZ,// 天线指向
+	float* demx, float* demy, float* demz, long* demcls, // 地面坐标
+	float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,// 地面坡度
+	float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, long Tthetapoints, long Tphipoints,// 天线方向图相关
+	float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, long Rthetapoints, long Rphipoints,// 天线方向图相关
+	float lamda, float fs, float nearrange, float Pt, long Freqnumbers, // 参数
+	CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby
+	cuComplex* outecho, long* d_echoAmpFID,
+	long len
+) {
+	int blockSize = 256; // 每个块的线程数
+	int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
+	// 调用 CUDA 核函数
+	CUDA_RTPC_Kernel<<<blockSize,numBlocks>>>(
+		antPx, antPy, antPZ,// 天线坐标
+		antXaxisX, antXaxisY, antXaxisZ, // 天线坐标系
+		antYaxisX, antYaxisY, antYaxisZ, //
+		antZaxisX, antZaxisY, antZaxisZ,
+		antDirectX, antDirectY, antDirectZ,// 天线指向
+		demx, demy, demz, demcls, // 地面坐标
+		demslopex, demslopey, demslopez, demslopeangle,// 地面坡度
+		Tantpattern, Tstarttheta, Tstartphi, Tdtheta, Tdphi, Tthetapoints, Tphipoints,// 天线方向图相关
+		Rantpattern, Rstarttheta, Rstartphi, Rdtheta, Rdphi, Rthetapoints, Rphipoints,// 天线方向图相关
+		lamda, fs, nearrange, Pt, Freqnumbers, // 参数
+		sigma0Paramslist, sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby
+		outecho, d_echoAmpFID,
+		len
+		);
+		cudaDeviceSynchronize();
+}
+
+
+
+
+
 #endif

GPUTool.cuh

@@ -17,14 +17,33 @@ enum LAMPGPUDATETYPE {
 };
 
 
+extern "C" struct  CUDASigmaParam {
+	float p1;
+	float p2;
+	float p3;
+	float p4;
+	float p5;
+	float p6;
+};
 
+extern "C" struct CUDAVector {
+	float x;
+	float y;
+	float z;
+};
+
+extern "C" struct CUDAVectorEllipsoidal {
+	float theta;
+	float phi;
+	float pho;
+};
 
 
 
 // GPU 内存函数
-extern "C" void mallocCUDAHost(void* ptr, long memsize); // 主机内存声明
+extern "C" void* mallocCUDAHost(  long memsize); // 主机内存声明
 extern "C" void FreeCUDAHost(void* ptr);
-extern "C" void mallocCUDADevice(void* ptr, long memsize); // GPU内存声明
+extern "C" void* mallocCUDADevice(  long memsize); // GPU内存声明
 extern "C" void FreeCUDADevice(void* ptr);
 extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移  设备 -> GPU
 extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移  GPU -> 设备
@@ -38,11 +57,35 @@ extern "C" void Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long memb
 extern "C" void cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos, long len);
 extern "C" void SatelliteAntDirectNormal(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, float* thetaAnt, float* phiAnt, long len);
 extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,float* localangle, float* R, float* slopeangle,float nearRange, float Fs, float pt, float lamda, long FreqIDmax,cuComplex* echoAmp, long* FreqID, long len);
+extern "C" void CUDA_RTPC_SiglePRF(
+	float antPx, float antPy, float antPZ,// 天线坐标
+	float  antXaxisX, float  antXaxisY, float  antXaxisZ, // 天线坐标系
+	float  antYaxisX, float  antYaxisY, float  antYaxisZ, //
+	float  antZaxisX, float  antZaxisY, float  antZaxisZ,
+	float  antDirectX, float  antDirectY, float  antDirectZ,// 天线指向
+	float* demx, float* demy, float* demz,long* demcls, // 地面坐标
+	float* demslopex, float* demslopey, float* demslopez, float* demslopeangle,// 地面坡度
+	float* Tantpattern, float Tstarttheta, float Tstartphi, float Tdtheta, float Tdphi, long Tthetapoints, long Tphipoints,// 天线方向图相关
+	float* Rantpattern, float Rstarttheta, float Rstartphi, float Rdtheta, float Rdphi, long Rthetapoints, long Rphipoints,// 天线方向图相关
+	float lamda, float fs, float nearrange, float Pt, long Freqnumbers, // 参数
+	CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 地表覆盖类型-sigma插值对应函数-ulaby
+	cuComplex* outecho,long* d_echoAmpFID,
+	long len
+);
+
+
 
 #endif
 #endif
 
 
+
+
+
+
+
+
+
 /**
 *