RasterProcessTool/Toolbox/SimulationSARTool/SimulationSAR/GPUTBPImage.cu

#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 "GPUTool.cuh"
#include "GPUTBPImage.cuh"
#include "GPUBPTool.cuh"

#ifdef __CUDANVCC___





__global__ void kernel_TimeBPImageGridNet(double* antPx, double* antPy, double* antPz,
	double* antDirx, double* antDiry, double* antDirz,
	double* imgx, double* imgy, double* imgz,
	long prfcount, long freqpoints, double meanH,
	double Rnear, double dx, double RefRange) {

	long idx = blockIdx.x * blockDim.x + threadIdx.x;
	long pixelcount = prfcount * freqpoints;
	long prfid = idx / freqpoints;
	long Rid = idx % freqpoints;
	if (idx < pixelcount) {
		// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>

		Vector3 S = { antPx[prfid], antPy[prfid], antPz[prfid] }; // <20><><EFBFBD><EFBFBD>λ<EFBFBD><CEBB> (m)
		Vector3 ray = { antDirx[prfid], antDiry[prfid], antDirz[prfid] };      // <20><><EFBFBD>߷<EFBFBD><DFB7><EFBFBD>
		double H = meanH;                    // ƽ<><C6BD><EFBFBD>߳<EFBFBD>
		double R = Rnear+ dx*Rid;                   // Ŀ<><C4BF><EFBFBD><EFBFBD><EFBFBD><EFBFBD>

		    // <20><><EFBFBD><EFBFBD>У<EFBFBD><D0A3>
		    if (R <= 0 || H < -WGS84_A * 0.1 || H > WGS84_A * 0.1) {
		        //printf("<22><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n  H<><48>Χ<EFBFBD><CEA7><EFBFBD><EFBFBD>%.1f km\n  R<><52><EFBFBD><EFBFBD>>0\n", WGS84_A * 0.1 / 1000);
				imgx[idx] = 1.0/0;
				imgy[idx] = 1.0/0;
				imgz[idx] = 1.0/0;

		        return ;
		    }
		
		    // Step 1: <20><><EFBFBD>㽻<EFBFBD><E3BDBB>T
		    Vector3 T = compute_T(S, ray, H);
		    if (isnan(T.x)) {
				imgx[idx] = 1.0/0;
				imgy[idx] = 1.0/0;
				imgz[idx] = 1.0/0;
		        return ;
		    }
		
		    // Step 2: <20><><EFBFBD><EFBFBD>Ŀ<EFBFBD><C4BF><EFBFBD><EFBFBD>P
		    Vector3 P = compute_P(S, T, R, H);
		
		    if (!isnan(P.x)) {

				imgx[idx] = P.x;
				imgy[idx] = P.y;
				imgz[idx] = P.z;
 
		    }
		    else {
				imgx[idx] = 1.0/0;
				imgy[idx] = 1.0/0;
				imgz[idx] = 1.0/0;
		        //printf("δ<>ҵ<EFBFBD><D2B5><EFBFBD>Ч<EFBFBD><D0A7>\n");
		    }
	}

}


__global__ void kernel_pixelTimeBP(
	double* antPx, double* antPy, double* antPz,
	double* imgx, double* imgy, double* imgz,
	cuComplex* TimeEchoArr, long prfcount, long pointCount,
	cuComplex* imgArr, long imH, long imW,
	double startLamda, double Rnear, double dx,double RefRange
) {
	long idx = blockIdx.x * blockDim.x + threadIdx.x;
	long pixelcount = imH * imW;
	if (idx < pixelcount) {

		double Tx = imgx[idx], Ty = imgy[idx], Tz = imgz[idx], PR = 0;
		double Px = 0, Py = 0, Pz = 0;
		double Rid = -1;
		long maxPointIdx = pointCount - 1;
		long pid = 0;
		long pid0 = 0;
		long pid1 = 0;
		double phi = 0;
		cuComplex s0=make_cuComplex(0,0), s1=make_cuComplex(0,0);
		cuComplex s = make_cuComplex(0, 0);
		cuComplex phiCorr = make_cuComplex(0, 0);
		for (long spid = 0; spid < prfcount; spid = spid + 8) {

#pragma unroll
			for (long sid = 0; sid < 8; sid++)
			{
				pid = spid + sid;
				if (pid < prfcount) {}
				else {
					continue;
				}

				Px = antPx[pid];
				Py = antPy[pid];
				Pz = antPz[pid];

				PR = sqrt((Px - Tx) * (Px - Tx) + (Py - Ty) * (Py - Ty) + (Pz - Tz) * (Pz - Tz));
				Rid = (PR - Rnear) / dx; // <20><><EFBFBD><EFBFBD>

				if (Rid<0 || Rid>maxPointIdx) {
					continue;
				}
				else {}

				pid0 = floor(Rid);
				pid1 = ceil(Rid);

				if (pid1<0 || pid0<0 || pid0>maxPointIdx || pid1>maxPointIdx) {
					continue;
				}
				else {}
				// <20><><EFBFBD>Բ<EFBFBD>ֵ
				s0 = TimeEchoArr[pid0];
				s1 = TimeEchoArr[pid1];

				s.x = s0.x * (Rid - pid0) + s1.x * (pid1 - Rid);
				s.y = s0.y * (Rid - pid0) + s1.y * (pid1 - Rid);

				// <20><>λУ<CEBB><D0A3>

				phi = 4 * LIGHTSPEED/startLamda* (PR - RefRange) ; // 4PI/lamda * R

				// exp(ix)=cos(x)+isin(x)
				phiCorr.x = cos(phi);
				phiCorr.y = sin(phi);
				
				s = cuCmulf(s, phiCorr); // У<><D0A3><EFBFBD><EFBFBD>

				imgArr[idx] = cuCaddf(imgArr[idx], s);
			}
		}
	}
}





extern "C" {
	void TIMEBPCreateImageGrid(double* antPx, double* antPy, double* antPz,
		double* antDirx, double* antDiry, double* antDirz, 
		double* imgx, double* imgy, double* imgz, 
		long prfcount, long freqpoints,   double meanH,
		double Rnear, double dx, double RefRange)
	{
		long pixelcount = prfcount * freqpoints;
		int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;

		kernel_TimeBPImageGridNet << <grid_size, BLOCK_SIZE >> > (
			antPx, antPy, antPz,
			antDirx, antDiry, antDirz,
			imgx, imgy, imgz,
			prfcount, freqpoints, meanH,
			Rnear, dx, RefRange);
		PrintLasterError("TIMEBPCreateImageGrid");
		cudaDeviceSynchronize();
	}



 
	void TimeBPImage(double* antPx, double* antPy, double* antPz,
		double* imgx, double* imgy, double* imgz,
		cuComplex* TimeEchoArr, long prfcount, long pointCount,
		cuComplex* imgArr, long imH, long imW,
		double startLamda, double Rnear, double dx, double RefRange
	)
	{
		long pixelcount = imH * imW;
		int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;

		kernel_pixelTimeBP << <grid_size, BLOCK_SIZE >> > (
			antPx, antPy, antPz,
			imgx, imgy, imgz,
			TimeEchoArr, prfcount, pointCount,
			imgArr, imH, imW,
			startLamda, Rnear, dx, RefRange
			);

		PrintLasterError("TimeBPImage");
		cudaDeviceSynchronize();

	}

}

#endif