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"

#ifdef __CUDANVCC___




__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 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