#include "stdafx.h"
#include "TBPImageAlgCls.h"
#include <QDateTime>
#include <QDebug>
#include <QString>
#include <cmath>
#include <QProgressDialog>
#include <QMessageBox>
#include "GPUTool.cuh"
#include "GPUTBPImage.cuh"
#include "ImageOperatorBase.h"

void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZPath)
{
	// 创建坐标系统
	long prfcount = echoL0ds->getPluseCount();
	long freqcount = echoL0ds->getPlusePoints();
	Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
	gt(0, 0) = 0;
	gt(0, 1) = 1;
	gt(0, 2) = 0;
	gt(1, 0) = 0;
	gt(1, 1) = 0;
	gt(1, 2) = 1;
	gdalImage xyzRaster = CreategdalImage(outPixelXYZPath, prfcount, freqcount, 3, gt, QString(""), false, true,true);
	std::shared_ptr<double> antpos = echoL0ds->getAntPos();
	double dx = (echoL0ds->getFarRange()-echoL0ds->getNearRange())/(echoL0ds->getPlusePoints()-1);
	double Rnear = echoL0ds->getNearRange();
	double Rref = echoL0ds->getRefPhaseRange();
	double centerInc = echoL0ds->getCenterAngle()*d2r;
	long echocol = 1073741824  / 8 / 4 / prfcount*8;
	qDebug() << "echocol:\t " << echocol ;
	echocol = echocol < 3000 ? 3000 : echocol;
	long startcolidx = 0;
	for (startcolidx = 0; startcolidx < freqcount; startcolidx = startcolidx + echocol) {

		long tempechocol = echocol;
		if (startcolidx + tempechocol >= freqcount) {
			tempechocol = freqcount - startcolidx;
		}
		qDebug() << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz")  << "]  imgxyz :\t" << startcolidx << "\t-\t" << startcolidx + tempechocol << " / " << freqcount  ;
		
		Eigen::MatrixXd demx = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 1);
		Eigen::MatrixXd demy = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 2);
		Eigen::MatrixXd demz = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 3);

		#pragma omp parallel for 
		for (long i = 0; i < prfcount; i++) {

			double Px = 0;
			double Py = 0;
			double Pz = 0;

			double AntDirectX = 0;
			double AntDirectY = 0;
			double AntDirectZ = 0;
			double R = 0;
			double NormAnt = 0;

			Px = antpos.get()[i * 19 + 1]; // 卫星坐标
			Py = antpos.get()[i * 19 + 2];
			Pz = antpos.get()[i * 19 + 3];
			AntDirectX = antpos.get()[i * 19 + 13];// zero doppler
			AntDirectY = antpos.get()[i * 19 + 14];
			AntDirectZ = antpos.get()[i * 19 + 15];

			NormAnt = std::sqrt(AntDirectX * AntDirectX + AntDirectY * AntDirectY + AntDirectZ * AntDirectZ);
			AntDirectX = AntDirectX / NormAnt;
			AntDirectY = AntDirectY / NormAnt;
			AntDirectZ = AntDirectZ / NormAnt;// 归一化

			// 计算中心参考点
			double centerX =  Px + Rref * AntDirectX; // T1
			double centerY =  Py + Rref * AntDirectY;
			double centerZ =  Pz + Rref * AntDirectZ;

			
			double satH = Rref * std::cos(centerInc); // 卫星高

			double PR = sqrt(Px * Px + Py * Py + Pz * Pz);

			double satR = PR - satH;
			double sPx = satR / PR * Px;
			double sPy = satR / PR * Py;
			double sPz = satR / PR * Pz;


			double dTSx  = sPx  - centerX;
			double dTSy  = sPy  - centerY;
			double dTSz  = sPz  - centerZ;

			double dTSR = sqrt(dTSx * dTSx + dTSy * dTSy + dTSz * dTSz);
			dTSx=dTSx/dTSR;
			dTSy=dTSy/dTSR;
			dTSz=dTSz/dTSR;

			for (long j = 0; j < tempechocol; j++) {
				R = (j + startcolidx)*dx + Rnear;

				double dRp = (Rref - R) / sin(centerInc); //  -- 0 +++


				demx(i,j) = centerX + dTSx * dRp;
				demy(i,j) = centerY + dTSy * dRp;
				demz(i,j) = centerZ + dTSz * dRp;
			}
		}

		xyzRaster.saveImage(demx, 0, startcolidx, 1);
		xyzRaster.saveImage(demy, 0, startcolidx, 2);
		xyzRaster.saveImage(demz, 0, startcolidx, 3);
	}
	
}

void TBPImageProcess(QString echofile, QString outImageFolder, QString imagePlanePath,long num_thread)
{
	std::shared_ptr<EchoL0Dataset> echoL0ds(new EchoL0Dataset);
	echoL0ds->Open(echofile);

	std::shared_ptr< SARSimulationImageL1Dataset> imagL1(new SARSimulationImageL1Dataset);
	imagL1->setCenterAngle(echoL0ds->getCenterAngle());
	imagL1->setCenterFreq(echoL0ds->getCenterFreq());
	imagL1->setNearRange(echoL0ds->getNearRange());
	imagL1->setRefRange((echoL0ds->getNearRange() + echoL0ds->getFarRange()) / 2);
	imagL1->setFarRange(echoL0ds->getFarRange());
	imagL1->setFs(echoL0ds->getFs());
	imagL1->setLookSide(echoL0ds->getLookSide());
	imagL1->OpenOrNew(outImageFolder, echoL0ds->getSimulationTaskName(), echoL0ds->getPluseCount(), echoL0ds->getPlusePoints());

	
	TBPImageAlgCls TBPimag;
	TBPimag.setEchoL0(echoL0ds);
	TBPimag.setImageL1(imagL1);
	TBPimag.setImagePlanePath(imagePlanePath);

	TBPimag.Process(num_thread);
}



void TBPImageAlgCls::setImagePlanePath(QString INimagePlanePath)
{
	this->imagePlanePath = INimagePlanePath;
}

QString TBPImageAlgCls::getImagePlanePath()
{
	return this->imagePlanePath;
}

void TBPImageAlgCls::setEchoL0(std::shared_ptr<EchoL0Dataset> inL0ds)
{
	this->L0ds = inL0ds;
}

void TBPImageAlgCls::setImageL1(std::shared_ptr<SARSimulationImageL1Dataset> inL1ds)
{
	this->L1ds = inL1ds;
}

std::shared_ptr<EchoL0Dataset> TBPImageAlgCls::getEchoL1()
{
	return  this->L0ds;
}

std::shared_ptr<SARSimulationImageL1Dataset> TBPImageAlgCls::getImageL0()
{
	return this->L1ds;
}

ErrorCode TBPImageAlgCls::Process(long num_thread)
{
	qDebug() << u8"创建成像平面的XYZ";
	QString outRasterXYZ = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_xyz.bin");
	CreatePixelXYZ(this->L0ds, outRasterXYZ);
	this->outRasterXYZPath = outRasterXYZ;

	// 初始化Raster
	qDebug() << u8"初始化影像";
	long imageheight = this->L1ds->getrowCount();
	long imagewidth = this->L1ds->getcolCount();

	 
	long blokline = Memory1GB  / 8 / 4 / imageheight * 32;
	blokline = blokline < 1000 ? 1000 : blokline;

	long startline = 0;
	for (startline = 0; startline < imageheight; startline = startline + blokline) {
		long templine = blokline;
		if (startline + templine >= imageheight) {
			templine = imageheight - startline;
		}
		qDebug() << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz")  << "]  imgxyz :\t" << startline << "\t-\t" << startline + templine << " / " << imageheight  ;

		std::shared_ptr<std::complex<double>> imageRaster = this->L1ds->getImageRaster(startline, templine);

		for (long i = 0; i < templine; i++) {
			for (long j = 0; j < imagewidth; j++) {
				imageRaster.get()[i * imagewidth + j] = std::complex<double>(0,0);
			}
		}
		this->L1ds->saveImageRaster(imageRaster, startline,templine);
	}
	qDebug() << u8"开始成像";
	qDebug() << u8"频域回波-> 时域回波";
	this->TimeEchoDataPath = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_Timeecho.bin");
	this->EchoFreqToTime();

	qDebug() << u8"频域回波-> 时域回波 结束";

	if (GPURUN) {
		return this->ProcessGPU();
	}
	else {
		QMessageBox::information(nullptr,u8"提示",u8"目前只支持显卡");
		return	ErrorCode::FAIL;
	}
}

ErrorCode TBPImageAlgCls::ProcessGPU()
{
	// 常用参数
	long rowCount = this->L1ds->getrowCount();
	long colCount = this->L1ds->getcolCount();
	long pixelCount = rowCount * colCount;
	long PRFCount = this->L0ds->getPluseCount();
	long PlusePoints = this->L0ds->getPlusePoints();
	long bandwidth = this->L0ds->getBandwidth();

	double Rnear = this->L1ds->getNearRange();
	double Rfar = this->L1ds->getFarRange();
	double refRange = this->L0ds->getRefPhaseRange();
	double dx = LIGHTSPEED / 2.0 / bandwidth; // c/2b
	Rfar = Rnear + dx * (PlusePoints - 1); // 更新斜距距离


	float freq =  this->L1ds->getCenterFreq();
	double factorj = freq * 4 * M_PI / LIGHTSPEED  ;

	qDebug() << "------------------------------------------------";
	qDebug() << "TBP params:";
	qDebug() << "Rnear:\t" << Rnear;
	qDebug() << "Rfar:\t" << Rfar;
	qDebug() << "refRange:\t" << this->getEchoL1()->getRefPhaseRange();
	qDebug() << "dx:\t" << dx;
	qDebug() << "freq:\t" << freq;
	qDebug() << "rowCount:\t" << rowCount;
	qDebug() << "colCount:\t" << colCount;
	qDebug() << "PRFCount:\t" << PRFCount;
	qDebug() << "PlusePoints:\t" << PlusePoints;
	qDebug() << "bandwidth:\t" << bandwidth;
	
	// 反方向计算起始相位
	
	double deltaF = bandwidth / (PlusePoints - 1);
	double startfreq = freq - bandwidth / 2;
	
	double startlamda = LIGHTSPEED / startfreq;
	qDebug() << "deltaF:\t" << deltaF;


	std::shared_ptr<double>  Pxs (new double[this->L0ds->getPluseCount()]);
	std::shared_ptr<double>  Pys (new double[this->L0ds->getPluseCount()]);
	std::shared_ptr<double>  Pzs (new double[this->L0ds->getPluseCount()]);

	{
		std::shared_ptr<double> antpos = this->L0ds->getAntPos();
		double time = 0;
		double Px = 0;
		double Py = 0;
		double Pz = 0;
		for (long i = 0; i < rowCount; i++) {
			time = antpos.get()[i *19 + 0];
			Px = antpos.get()[i *19 + 1];
			Py = antpos.get()[i *19 + 2];
			Pz = antpos.get()[i *19 + 3];
			Pxs.get()[i] = Px;
			Pys.get()[i] = Py;
			Pzs.get()[i] = Pz;
		}
		antpos.reset();
	}




	// 计算成像的基本参数
	// 距离向分辨率
	double dr = sqrt(pow(Pxs.get()[2]- Pxs.get()[1],2)+pow(Pys.get()[2] - Pys.get()[1],2)+pow(Pzs.get()[2] - Pzs.get()[1],2));
	qDebug() << "-------  resolution ----------------------------------";
	qDebug() << "Range Resolution (m):\t" << dx ;
	qDebug() << "Cross Resolution (m):\t" << dr;
	qDebug() << "Range Range (m):\t" << dx*PlusePoints;
	qDebug() << "Cross Range (m):\t" << dr*PRFCount;
	qDebug() << "start Freq (Hz):\t" << startfreq;
	qDebug() << "start lamda (m):\t" << startlamda;
	qDebug() << "rowCount:\t" << rowCount;
	qDebug() << "colCount:\t" << colCount;
	qDebug() << "PRFCount:\t" << PRFCount;
	qDebug() << "PlusePoints:\t" << PlusePoints;
	qDebug() << "Rnear:\t" << Rnear;
	qDebug() << "Rfar:\t" << Rfar;
	qDebug() << "refRange:\t" << refRange;
	// 方位向分辨率


	// 按照回波分块，图像分块
	long echoBlockline = Memory1GB  / 8 / 2 / PlusePoints * 2; //2GB
	echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;

	long imageBlockline = Memory1GB  / 8 / 2 / colCount * 2; //2GB
	imageBlockline = imageBlockline < 1 ? 1 : imageBlockline;

	gdalImage imageXYZ(this->outRasterXYZPath);

	long startimgrowid = 0;
	for (startimgrowid = 0; startimgrowid < rowCount; startimgrowid = startimgrowid + imageBlockline) {
		long tempimgBlockline = imageBlockline;
		if (startimgrowid + imageBlockline >= rowCount) {
			tempimgBlockline = rowCount - startimgrowid;
		}
		qDebug() << "\r image create Row Range  :\t"<<QString("[%1]").arg(startimgrowid*100.0/ rowCount)
			<< startimgrowid << "\t-\t" << startimgrowid + tempimgBlockline << "\t/\t" << rowCount;
		// 提取局部pixel x,y,z
		std::shared_ptr<double> img_x = readDataArr<double>(imageXYZ,startimgrowid,0,tempimgBlockline,colCount,1,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
		std::shared_ptr<double> img_y = readDataArr<double>(imageXYZ,startimgrowid,0,tempimgBlockline,colCount,2,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
		std::shared_ptr<double> img_z = readDataArr<double>(imageXYZ,startimgrowid,0,tempimgBlockline,colCount,3,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);

		std::shared_ptr<std::complex<double>> imgArr = this->L1ds->getImageRaster(startimgrowid, tempimgBlockline);
		// 获取回波
		long startechoid = 0;
		for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) {
			
			long tempechoBlockline = echoBlockline;
			if (startechoid + tempechoBlockline >= PRFCount) {
				tempechoBlockline = PRFCount - startechoid;
			}
			std::shared_ptr<std::complex<double>> echoArr = this->L0ds->getEchoArr(startechoid, tempechoBlockline);
			std::shared_ptr<double> antpx(new double[tempechoBlockline*PlusePoints],delArrPtr);
			std::shared_ptr<double> antpy(new double[tempechoBlockline* PlusePoints], delArrPtr);
			std::shared_ptr<double> antpz(new double[tempechoBlockline* PlusePoints], delArrPtr);
			// 复制
			for (long anti = 0; anti < tempechoBlockline; anti++) {
				antpx.get()[anti] = Pxs.get()[anti + startechoid];
				antpy.get()[anti] = Pys.get()[anti + startechoid];
				antpz.get()[anti] = Pzs.get()[anti + startechoid];
			}
			TBPImageGPUAlg2(
				antpx, antpy, antpz,
				img_x, img_y, img_z,
				echoArr, imgArr,
				startfreq, dx,
				Rnear, Rfar, refRange,
				tempimgBlockline, colCount, 
				tempechoBlockline, PlusePoints,
				startechoid,startimgrowid
				);
			qDebug() << QString(" image block PRF:[%1] \t").arg((startechoid + tempechoBlockline) * 100.0 / PRFCount)
				<< startechoid << "\t-\t" << startechoid + tempechoBlockline;
		}
		this->L1ds->saveImageRaster(imgArr, startimgrowid, tempimgBlockline);
	}
 
	qDebug() << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz")  << "] image writing:\t" << this->L1ds->getxmlFilePath();
	this->L1ds->saveToXml();
	return ErrorCode::SUCCESS;
}

void TBPImageAlgCls::EchoFreqToTime( )
{
	// 读取数据

	long PRFCount = this->L0ds->getPluseCount();
	long inColCount = this->L0ds->getPlusePoints();
	long outColCount = nextpow2(inColCount);
	this->TimeEchoRowCount = PRFCount;
	this->TimeEchoColCount = outColCount;
	qDebug() << "IFFT : " << this->TimeEchoDataPath;
	qDebug() << "PRF Count:\t" << PRFCount;
	qDebug() << "inColCount:\t" << inColCount;
	qDebug() << "outColCount:\t" << outColCount;
	// 创建二进制文件
	gdalImageComplex outTimeEchoImg = CreategdalImageComplexNoProj(this->TimeEchoDataPath,this->TimeEchoRowCount,this->TimeEchoColCount,1);

	// 分块
	long echoBlockline = Memory1GB / 8 / 2 / outColCount * 1; //1GB
	echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;


	long startechoid = 0;
	for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) {

		long tempechoBlockline = echoBlockline;
		if (startechoid + tempechoBlockline >= PRFCount) {
			tempechoBlockline = PRFCount - startechoid;
		}
		std::shared_ptr<std::complex<double>> echoArr = this->L0ds->getEchoArr(startechoid, tempechoBlockline);
		std::shared_ptr<std::complex<double>> IFFTArr = outTimeEchoImg.getDataComplexSharePtr(startechoid, 0, tempechoBlockline, outColCount, 1);

		std::shared_ptr<cuComplex> host_echoArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex)* tempechoBlockline * inColCount), FreeCUDAHost);
		std::shared_ptr<cuComplex> host_IFFTechoArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex)* tempechoBlockline * outColCount), FreeCUDAHost);
	
#pragma omp parallel for
		for (long ii = 0; ii < tempechoBlockline * inColCount; ii++) {
			host_echoArr.get()[ii] = make_cuComplex(echoArr.get()[ii].real(), echoArr.get()[ii].imag());
		}
#pragma omp parallel for
		for (long ii = 0; ii < tempechoBlockline * outColCount; ii++) {
			host_IFFTechoArr.get()[ii] = make_cuComplex(0,0);
		}

		std::shared_ptr<cuComplex> device_echoArr((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * tempechoBlockline * inColCount), FreeCUDADevice);
		std::shared_ptr<cuComplex> device_IFFTechoArr((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * tempechoBlockline * outColCount), FreeCUDADevice);

		HostToDevice(host_echoArr.get(), device_echoArr.get(), sizeof(cuComplex) * tempechoBlockline * inColCount);
		HostToDevice(host_IFFTechoArr.get(), device_IFFTechoArr.get(), sizeof(cuComplex) * tempechoBlockline * outColCount);
		CUDAIFFT(device_echoArr.get(), device_IFFTechoArr.get(), tempechoBlockline, inColCount, outColCount);

		DeviceToHost(host_IFFTechoArr.get(), device_IFFTechoArr.get(), sizeof(cuComplex) * tempechoBlockline * outColCount);

#pragma omp parallel for
		for (long ii = 0; ii < tempechoBlockline * outColCount; ii++) {
			IFFTArr.get()[ii] = std::complex<double>(host_IFFTechoArr.get()[ii].x, host_IFFTechoArr.get()[ii].y);
		}

		outTimeEchoImg.saveImage(IFFTArr, startechoid, 0, tempechoBlockline, outColCount, 1);

		qDebug() << QString(" image block PRF:[%1] \t").arg((startechoid + tempechoBlockline) * 100.0 / PRFCount)
			<< startechoid << "\t-\t" << startechoid + tempechoBlockline;
	}

	return;


}







void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antPy, std::shared_ptr<double> antPz,
	std::shared_ptr<double> img_x, std::shared_ptr<double> img_y, std::shared_ptr<double> img_z,
	std::shared_ptr<std::complex<double>> echoArr,
	std::shared_ptr<std::complex<double>> img_arr,
	double freq, double dx, double Rnear, double Rfar, double refRange,
	long rowcount, long colcount,
	long prfcount, long freqcount,
	long startPRFId, long startRowID
)
{
	long IFFTPadNum = nextpow2(freqcount);
	// 先处理脉冲傅里叶变换
	cuComplex* h_echoArr = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * prfcount * freqcount);
	cuComplex* d_echoArr = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * prfcount * freqcount);
	std::shared_ptr<cuComplex> d_echoArrIFFT((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * prfcount * IFFTPadNum), FreeCUDADevice);

	// 回波赋值
	for (long i = 0; i < prfcount; i++) {
		for (long j = 0; j < freqcount; j++) {
			h_echoArr[i * freqcount + j] = make_cuComplex(echoArr.get()[i * freqcount + j].real(),
				echoArr.get()[i * freqcount + j].imag());
		}
	}
	HostToDevice(h_echoArr, d_echoArr, sizeof(cuComplex) * prfcount * freqcount);
	CUDAIFFT(d_echoArr, d_echoArrIFFT.get(), prfcount, freqcount, IFFTPadNum);
	// 结束傅里叶变换
	FreeCUDAHost(h_echoArr);
	FreeCUDADevice(d_echoArr);

	qDebug() << "IFFT finished!!!";
	// 初始化
	std::shared_ptr<double> h_antPx ((double*)mallocCUDAHost(sizeof(double) * prfcount),FreeCUDAHost);
	std::shared_ptr<double> h_antPy ((double*)mallocCUDAHost(sizeof(double) * prfcount),FreeCUDAHost);
	std::shared_ptr<double> h_antPz ((double*)mallocCUDAHost(sizeof(double) * prfcount),FreeCUDAHost);

	std::shared_ptr<double> d_antPx ((double*)mallocCUDADevice(sizeof(double) * prfcount),FreeCUDADevice);
	std::shared_ptr<double> d_antPy ((double*)mallocCUDADevice(sizeof(double) * prfcount),FreeCUDADevice);
	std::shared_ptr<double> d_antPz ((double*)mallocCUDADevice(sizeof(double) * prfcount),FreeCUDADevice);

	std::shared_ptr<double> h_imgx((double*)mallocCUDAHost(sizeof(double) * rowcount * colcount),FreeCUDAHost);
	std::shared_ptr<double> h_imgy((double*)mallocCUDAHost(sizeof(double) * rowcount * colcount),FreeCUDAHost);
	std::shared_ptr<double> h_imgz((double*)mallocCUDAHost(sizeof(double) * rowcount * colcount),FreeCUDAHost);

	std::shared_ptr<double> d_imgx ((double*)mallocCUDADevice(sizeof(double) * rowcount * colcount),FreeCUDADevice);
	std::shared_ptr<double> d_imgy ((double*)mallocCUDADevice(sizeof(double) * rowcount * colcount),FreeCUDADevice);
	std::shared_ptr<double> d_imgz ((double*)mallocCUDADevice(sizeof(double) * rowcount * colcount),FreeCUDADevice);
 
	
	// 天线位置
	for (long i = 0; i < prfcount; i++) {
		h_antPx.get()[i] = antPx.get()[i];
		h_antPy.get()[i] = antPy.get()[i];
		h_antPz.get()[i] = antPz.get()[i];
	}

#pragma omp parallel for
	for (long i = 0; i < rowcount; i++) {
		for (long j = 0; j < colcount; j++) {
			h_imgx.get()[i * colcount + j] = img_x.get()[i * colcount + j];
			h_imgy.get()[i * colcount + j] = img_y.get()[i * colcount + j];
			h_imgz.get()[i * colcount + j] = img_z.get()[i * colcount + j];
		}
	}

	HostToDevice(h_antPx.get(), d_antPx.get(), sizeof(double) * prfcount);
	HostToDevice(h_antPy.get(), d_antPy.get(), sizeof(double) * prfcount);
	HostToDevice(h_antPz.get(), d_antPz.get(), sizeof(double) * prfcount);

	HostToDevice(h_imgx.get(), d_imgx.get(), sizeof(double) * rowcount * colcount);
	HostToDevice(h_imgy.get(), d_imgy.get(), sizeof(double) * rowcount * colcount);
	HostToDevice(h_imgz.get(), d_imgz.get(), sizeof(double) * rowcount * colcount);


	std::shared_ptr<cuComplex> h_imgArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex) * rowcount * colcount), FreeCUDAHost);
	std::shared_ptr<cuComplex> d_imgArr((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * rowcount * colcount), FreeCUDADevice);


#pragma omp parallel for
	for (long i = 0; i < rowcount; i++) {
		for (long j = 0; j < colcount; j++) {
			h_imgArr.get()[i * colcount + j].x = img_arr.get()[i * colcount + j].real();
			h_imgArr.get()[i * colcount + j].y = img_arr.get()[i * colcount + j].imag();
		}
	}
	HostToDevice(h_imgArr.get(), d_imgArr.get(), sizeof(cuComplex) * rowcount * colcount);

	// 直接使用
	double startlamda = LIGHTSPEED / freq;
 
	TimeBPImage(
		d_antPx.get(), d_antPy.get(), d_antPz.get(),
		d_imgx.get(), d_imgy.get(), d_imgz.get(),
		d_echoArrIFFT.get(), prfcount, IFFTPadNum,
		d_imgArr.get(), rowcount, colcount,
		startlamda, Rnear, dx, refRange
	);

	// Device -> Host
	DeviceToHost(h_imgArr.get(), d_imgArr.get(), sizeof(cuComplex)* rowcount* colcount);
 


#pragma omp parallel for
	for (long i = 0; i < rowcount; i++) {
		for (long j = 0; j < colcount; j++) {
			img_arr.get()[i * colcount + j] = std::complex<double>(h_imgArr.get()[i * colcount + j].x, h_imgArr.get()[i * colcount + j].y);
		}
	}

 

}












void TBPImageAlgCls::setGPU(bool flag)
{
	this->GPURUN = flag;
}

bool TBPImageAlgCls::getGPU( )
{
	return this->GPURUN;
}


/**
ErrorCode TBPImageAlgCls::ProcessCPU(long num_thread)
{
	omp_set_num_threads(num_thread);
	// 常用参数
	long rowCount = this->L1ds->getrowCount();
	long colCount = this->L1ds->getcolCount();
	long pixelCount = rowCount * colCount;
	long PRFCount = this->L0ds->getPluseCount();
	long PlusePoints = this->L0ds->getPlusePoints();


	double Rnear = this->L1ds->getNearRange();
	double Rfar = this->L1ds->getFarRange();
	double fs = this->L1ds->getFs();
	double dx = LIGHTSPEED / 2 / fs;
	double factorj = this->L1ds->getCenterFreq() * 4 * M_PI / LIGHTSPEED * 1e9;


	Eigen::MatrixXcd echo = Eigen::MatrixXcd::Zero(PRFCount, PlusePoints);
	{
		std::shared_ptr<std::complex<double>> echodata = this->L0ds->getEchoArr();
		for (long i = 0; i < PRFCount; i++) {
			for (long j = 0; j < PlusePoints; j++) {
				echo(i, j) = echodata.get()[i * PlusePoints + j];
			}
		}
		echodata.reset();
	}

	Eigen::MatrixXd pixelX = Eigen::MatrixXd::Zero(rowCount, colCount);
	Eigen::MatrixXd pixelY = Eigen::MatrixXd::Zero(rowCount, colCount);
	Eigen::MatrixXd pixelZ = Eigen::MatrixXd::Zero(rowCount, colCount);

	Eigen::MatrixXd Pxs = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1);
	Eigen::MatrixXd Pys = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1);
	Eigen::MatrixXd Pzs = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1);

	// 图像网格坐标 
	{
		std::shared_ptr<double> antpos = this->L0ds->getAntPos();
		double time = 0;
		double Px = 0;
		double Py = 0;
		double Pz = 0;
		double Vx = 0;
		double Vy = 0;
		double Vz = 0;
		double AntDirectX = 0;
		double AntDirectY = 0;
		double AntDirectZ = 0;
		double AVx = 0;
		double AVy = 0;
		double AVz = 0;

		double R = 0;
		double NormAnt = 0;

		for (long i = 0; i < rowCount; i++) {
			time = antpos.get()[i * 19 + 0];
			Px = antpos.get()[i * 19 + 1];
			Py = antpos.get()[i * 19 + 2];
			Pz = antpos.get()[i * 19 + 3];
			Vx = antpos.get()[i * 19 + 4];
			Vy = antpos.get()[i * 19 + 5];
			Vz = antpos.get()[i * 19 + 6];
			AntDirectX = antpos.get()[i * 19 + 13]; // Zero doppler 
			AntDirectY = antpos.get()[i * 19 + 14];
			AntDirectZ = antpos.get()[i * 19 + 15];
			AVx = antpos.get()[i * 19 + 10];
			AVy = antpos.get()[i * 19 + 11];
			AVz = antpos.get()[i * 19 + 12];

			NormAnt = std::sqrt(AntDirectX * AntDirectX + AntDirectY * AntDirectY + AntDirectZ * AntDirectZ);
			AntDirectX = AntDirectX / NormAnt;
			AntDirectY = AntDirectY / NormAnt;
			AntDirectZ = AntDirectZ / NormAnt;// 归一化

			antpos.get()[i * 19 + 13] = AntDirectX;
			antpos.get()[i * 19 + 14] = AntDirectY;
			antpos.get()[i * 19 + 15] = AntDirectZ;
			Pxs(i, 0) = Px;
			Pys(i, 0) = Py;
			Pzs(i, 0) = Pz;

			for (long j = 0; j < colCount; j++) {
				R = j * dx + Rnear;
				pixelX(i, j) = Px + AntDirectX * R;
				pixelY(i, j) = Py + AntDirectY * R;
				pixelZ(i, j) = Pz + AntDirectZ * R;
			}
		}

		this->L1ds->saveAntPos(antpos);
		antpos.reset();
	}

	// BP成像
	long BlockLine = Memory1MB * 10 / 16 / rowCount;
	if (rowCount / BlockLine / num_thread < 3) {
		BlockLine = rowCount / num_thread / 3;
	}
	BlockLine = BlockLine > 10 ? BlockLine : 10;

	std::shared_ptr<std::complex<double>> imagarr = this->L1ds->getImageRaster();
	{
		for (long i = 0; i < pixelCount; i++) {
			imagarr.get()[i] = imagarr.get()[i];
		}
	}
	omp_lock_t lock; // 定义锁
	omp_init_lock(&lock); // 初始化锁

	long writeImageCount = 0;
	qDebug() << "block line:\t" << BlockLine;


	long startLine = 0;
	long processValue = 0;
	long processNumber = 0;
	QProgressDialog progressDialog(u8"RFPC回波生成中", u8"终止", 0, rowCount);
	progressDialog.setWindowTitle(u8"RFPC回波生成中");
	progressDialog.setWindowModality(Qt::WindowModal);
	progressDialog.setAutoClose(true);
	progressDialog.setValue(0);
	progressDialog.setMaximum(rowCount);
	progressDialog.setMinimum(0);
	progressDialog.show();

#pragma omp parallel for
	for (startLine = 0; startLine < rowCount; startLine = startLine + BlockLine) { // 图像大小
		long stepLine = startLine + BlockLine < rowCount ? BlockLine : rowCount - startLine;
		long imageRowID = startLine; // 

		//Eigen::MatrixXd R = Eigen::MatrixXd::Zero(rowCount, 1);
		long pluseId = 0;
		std::complex<double> factPhas(0, 0);
		std::complex < double> sign(0, 0);
		Eigen::MatrixXd R = Eigen::MatrixXd::Zero(rowCount, 1);
		Eigen::MatrixXcd Rphi = Eigen::MatrixXd::Zero(rowCount, 1);

		long PluseIDs = 0;
		double mask = 0;

		for (long i = 0; i < stepLine; i++) { // 图像行
			imageRowID = startLine + i;
			for (long j = 0; j < colCount; j++) { //图像列
				R = ((pixelX(i, j) - Pxs.array()).array().pow(2) + (pixelY(i, j) - Pys.array()).array().pow(2) + (pixelZ(i, j) - Pzs.array()).array().pow(2)).array().sqrt();
				Rphi = Rphi.array() * 0;
				Rphi.imag() = R.array() * factorj;
				Rphi = Rphi.array().exp();

				for (long prfid = 0; prfid < rowCount; prfid++) { // 脉冲行
					PluseIDs = std::floor((R(prfid, 0) - Rnear) / dx);
					mask = (PluseIDs < 0 || PluseIDs >= PlusePoints) ? 0 : 1;
					PluseIDs = (PluseIDs < 0 || PluseIDs >= PlusePoints) ? 0 : PluseIDs;
					imagarr.get()[imageRowID * colCount + j] =
						imagarr.get()[imageRowID * colCount + j] +
						mask * echo(prfid, PluseIDs) * Rphi(prfid, 0);// 信号* 相位校正
				}
			}

		}
		omp_set_lock(&lock); // 保存文件
		processValue = processValue + BlockLine;
		this->L1ds->saveImageRaster(imagarr, 0, rowCount);
		qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz").toUtf8().constData() << "\t" << processValue * 100.0 / rowCount << "%\t" << startLine << "\t-\t" << startLine + BlockLine << "\tend\t\t";
		processNumber = processNumber + BlockLine;
		processNumber = processNumber < progressDialog.maximum() ? processNumber : progressDialog.maximum();
		progressDialog.setValue(processNumber);
		omp_unset_lock(&lock); // 解锁
	}
	omp_destroy_lock(&lock); // 销毁锁

	this->L1ds->saveImageRaster(imagarr, 0, rowCount);
	this->L1ds->saveToXml();
	progressDialog.close();
	return ErrorCode::SUCCESS;
}


*/