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

同步BP成像算法

pull/6/head
陈增辉 2025-03-03 11:18:50 +08:00
parent 982922506c
commit 94fd13bac8
4 changed files with 304 additions and 380 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

245
Toolbox/SimulationSARTool/SimulationSAR/GPUTBPImage.cu
View File

@ -24,7 +24,6 @@
__global__ void kernel_TimeBPImageGridNet(double* antPx, double* antPy, double* antPz, __global__ void kernel_TimeBPImageGridNet(double* antPx, double* antPy, double* antPz,
double* antDirx, double* antDiry, double* antDirz, double* antDirx, double* antDiry, double* antDirz,
double* imgx, double* imgy, double* imgz, double* imgx, double* imgy, double* imgz,
@ -170,77 +169,193 @@ __global__ void kernel_TimeBPImageGridNet(double* antPx, double* antPy, double*
} }
__device__ double computerR(double& Px, double& Py, double& Pz, double& Tx, double& Ty, double& Tz) {
return sqrt((Px - Tx) * (Px - Tx) + (Py - Ty) * (Py - Ty) + (Pz - Tz) * (Pz - Tz));
}
__device__ void updateBPImage(
long prfid, long pixelidx,double R,double LRrange,
cuComplex* TimeEchoArr, long prfcount, long pointCount,
cuComplex* imgArr,
double startLamda, double Rnear, double dx, double RefRange, double Rfar
) {
double ridx = (R - LRrange) / dx; // 获取范围
if (ridx < 0 || ridx >= pointCount) {
return;
}
else {}
long Ridx0 = floor(ridx);
long Ridx1 = ceil(ridx);
long pid0 = prfid * pointCount + Ridx0;
long pid1 = prfid * pointCount + Ridx1;
cuComplex s0 = TimeEchoArr[pid0];
cuComplex s1 = TimeEchoArr[pid1];
cuComplex s = make_cuComplex(
s0.x + (s1.x - s0.x) * (ridx-Ridx0), // real
s0.y + (s1.y - s0.y) * (ridx-Ridx0) // imag
);
double phi = 4 * LIGHTSPEED / startLamda * (R - RefRange);
// exp(ix)=cos(x)+isin(x)
cuComplex phiCorr = make_cuComplex(cos(phi), sin(phi));
s = cuCmulf(s, phiCorr); // 校正后
imgArr[pixelidx] = cuCaddf(imgArr[pixelidx], s);
return;
}
// 分块计算
__device__ void segmentBPImage(
double* antPx, double* antPy, double* antPz,
double Tx, double Ty, double Tz,
cuComplex* TimeEchoArr, long prfcount, long pointCount,
cuComplex* imgArr,
double startLamda, double Rnear, double dx, double RefRange, double Rfar,
long startSegmentPrfId,long pixelID // 分段起始prfid
) {
// 计算单条脉冲范围
double Rrange = pointCount * dx;// 成像范围
double LRrange = RefRange - Rrange / 2;//左范围
double RRrange = RefRange + Rrange / 2;
long currentprfid = 0;
// 0
currentprfid = currentprfid + 0;
double Px = antPx[currentprfid];
double Py = antPy[currentprfid];
double Pz = antPz[currentprfid];
double R0 = computerR(Px, Py, Pz, Tx, Ty, Tz);
if (Rnear <= R0 && R0 <= RRrange) {
updateBPImage(
currentprfid, pixelID, R0, LRrange,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar
);
}
// 10
currentprfid = currentprfid + 10;
Px = antPx[currentprfid];
Py = antPy[currentprfid];
Pz = antPz[currentprfid];
double R10 = computerR(Px, Py, Pz, Tx, Ty, Tz);
if (Rnear <= R10 && R10 <= RRrange) {
updateBPImage(
currentprfid, pixelID, R10, LRrange,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar
);
}
//19
currentprfid = currentprfid + 19;
Px = antPx[currentprfid];
Py = antPy[currentprfid];
Pz = antPz[currentprfid];
double R19 = computerR(Px, Py, Pz, Tx, Ty, Tz);
if (Rnear <= R19 && R19 <= RRrange) {
updateBPImage(
currentprfid, pixelID, R19, LRrange,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar
);
}
// 判断是否需要处理
if (R0 < LRrange && R10 < LRrange && R19 < LRrange ) { // 越界、不处理
return;
}
else if (R0 > RRrange && R10 > RRrange && R19 > RRrange) {// 越界、不处理
return;
}
else {}
double R = 0;
#pragma unroll
for (long i = 1; i < 10; i++) {
currentprfid = startSegmentPrfId + i;
if (currentprfid < prfcount) {
Px = antPx[currentprfid];
Py = antPy[currentprfid];
Pz = antPz[currentprfid];
R = computerR(Px, Py, Pz, Tx, Ty, Tz);
updateBPImage(
currentprfid, pixelID, R, LRrange,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar
);
}
}
#pragma unroll
for (long i = 11; i < 19; i++) {
currentprfid = startSegmentPrfId + i;
if (currentprfid < prfcount) {
Px = antPx[currentprfid];
Py = antPy[currentprfid];
Pz = antPz[currentprfid];
R = computerR(Px, Py, Pz, Tx, Ty, Tz);
updateBPImage(
currentprfid, pixelID, R, LRrange,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar
);
}
}
}
__global__ void kernel_pixelTimeBP( __global__ void kernel_pixelTimeBP(
double* antPx, double* antPy, double* antPz, double* antPx, double* antPy, double* antPz,
double* imgx, double* imgy, double* imgz, double* imgx, double* imgy, double* imgz,
cuComplex* TimeEchoArr, long prfcount, long pointCount, cuComplex* TimeEchoArr, long prfcount, long pointCount,
cuComplex* imgArr, long imH, long imW, cuComplex* imgArr, long imH, long imW,
double startLamda, double Rnear, double dx,double RefRange double startLamda, double Rnear, double dx,double RefRange,double Rfar
) { ) {
long idx = blockIdx.x * blockDim.x + threadIdx.x; long idx = blockIdx.x * blockDim.x + threadIdx.x;
long pixelcount = imH * imW; long pixelcount = imH * imW;
if (idx < pixelcount) { if (idx < pixelcount) {
double Tx = imgx[idx], Ty = imgy[idx], Tz = imgz[idx], PR = 0; double Tx = imgx[idx]; // 地面坐标点
double Px = 0, Py = 0, Pz = 0; double Ty = imgy[idx];
double Rid = -1; double Tz = imgz[idx];
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]; for (long segid = 0; segid < prfcount; segid = segid + 20) {
Py = antPy[pid]; // 判断范围
Pz = antPz[pid]; segmentBPImage(
antPx, antPy, antPz,
Tx, Ty, Tz,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar,
segid, idx
);
PR = sqrt((Px - Tx) * (Px - Tx) + (Py - Ty) * (Py - Ty) + (Pz - Tz) * (Pz - Tz));
Rid = (PR - Rnear) / dx; // 行数
if (Rid<0 || Rid>maxPointIdx) {
continue;
}
else {}
pid0 = floor(Rid);
pid1 = ceil(Rid);
if (pid1<0 || pid0<0 || pid0>maxPointIdx || pid1>maxPointIdx) {
continue;
}
else {}
// 线性插值
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);
// 相位校正
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); // 校正后
imgArr[idx] = cuCaddf(imgArr[idx], s);
}
} }
} }
} }
@ -254,7 +369,8 @@ extern "C" {
double* antDirx, double* antDiry, double* antDirz, double* antDirx, double* antDiry, double* antDirz,
double* imgx, double* imgy, double* imgz, double* imgx, double* imgy, double* imgz,
long prfcount, long freqpoints, double meanH, long prfcount, long freqpoints, double meanH,
double Rnear, double dx, double RefRange) double Rnear, double dx, double RefRange
)
{ {
long pixelcount = prfcount * freqpoints; long pixelcount = prfcount * freqpoints;
int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE; int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
@ -272,22 +388,23 @@ extern "C" {
void TimeBPImage(double* antPx, double* antPy, double* antPz, void TimeBPImage(
double* antPx, double* antPy, double* antPz,
double* imgx, double* imgy, double* imgz, double* imgx, double* imgy, double* imgz,
cuComplex* TimeEchoArr, long prfcount, long pointCount, cuComplex* TimeEchoArr, long prfcount, long pointCount,
cuComplex* imgArr, long imH, long imW, cuComplex* imgArr, long imH, long imW,
double startLamda, double Rnear, double dx, double RefRange double startLamda, double Rnear, double dx, double RefRange,double Rfar
) )
{ {
long pixelcount = imH * imW; long pixelcount = imH * imW;
int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE; int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
kernel_pixelTimeBP << <1, 1 >> > ( kernel_pixelTimeBP << <grid_size, BLOCK_SIZE >> > (
antPx, antPy, antPz, antPx, antPy, antPz,
imgx, imgy, imgz, imgx, imgy, imgz,
TimeEchoArr, prfcount, pointCount, TimeEchoArr, prfcount, pointCount,
imgArr, imH, imW, imgArr, imH, imW,
startLamda, Rnear, dx, RefRange startLamda, Rnear, dx, RefRange, Rfar
); );
PrintLasterError("TimeBPImage"); PrintLasterError("TimeBPImage");

2
Toolbox/SimulationSARTool/SimulationSAR/GPUTBPImage.cuh
View File

@ -26,7 +26,7 @@ extern "C" void TimeBPImage(
double* imgx, double* imgy, double* imgz, double* imgx, double* imgy, double* imgz,
cuComplex* TimeEchoArr, long prfcount, long pointCount, cuComplex* TimeEchoArr, long prfcount, long pointCount,
cuComplex* imgArr, long imH, long imW, cuComplex* imgArr, long imH, long imW,
double startLamda, double Rnear, double dx, double RefRange double startLamda, double Rnear, double dx, double RefRange,double Rfar
); );

411
Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.cpp
View File

@ -126,7 +126,7 @@ void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZ
d_AntDirectX.get(), d_AntDirectY.get(), d_AntDirectZ.get(), d_AntDirectX.get(), d_AntDirectY.get(), d_AntDirectZ.get(),
d_demx.get(), d_demy.get(), d_demz.get(), d_demx.get(), d_demy.get(), d_demz.get(),
prfcount, tempechocol, 1000, prfcount, tempechocol, 1000,
Rnear+dx* startcolidx, dx, refRange Rnear+dx* startcolidx, dx, refRange // 更新最近修读
); );
DeviceToHost(h_demx.get(), d_demx.get(), sizeof(double) * prfcount * tempechocol); DeviceToHost(h_demx.get(), d_demx.get(), sizeof(double) * prfcount * tempechocol);
@ -252,13 +252,13 @@ ErrorCode TBPImageAlgCls::Process(long num_thread)
qDebug() << u8"频域回波-> 时域回波 结束"; qDebug() << u8"频域回波-> 时域回波 结束";
//if (GPURUN) { if (GPURUN) {
// return this->ProcessGPU(); return this->ProcessGPU();
//} }
//else { else {
// QMessageBox::information(nullptr,u8"提示",u8"目前只支持显卡"); QMessageBox::information(nullptr,u8"提示",u8"目前只支持显卡");
// return ErrorCode::FAIL; return ErrorCode::FAIL;
//} }
return ErrorCode::SUCCESS; return ErrorCode::SUCCESS;
} }
@ -348,8 +348,6 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
qDebug() << "refRange:\t" << refRange; qDebug() << "refRange:\t" << refRange;
// 方位向分辨率 // 方位向分辨率
// 按照回波分块，图像分块
long echoBlockline = Memory1GB / 8 / 2 / PlusePoints * 2; //2GB long echoBlockline = Memory1GB / 8 / 2 / PlusePoints * 2; //2GB
echoBlockline = echoBlockline < 1 ? 1 : echoBlockline; echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;
@ -365,43 +363,41 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
if (startimgrowid + imageBlockline >= rowCount) { if (startimgrowid + imageBlockline >= rowCount) {
tempimgBlockline = rowCount - startimgrowid; tempimgBlockline = rowCount - startimgrowid;
} }
qDebug() << "\r image create Row Range :\t"<<QString("[%1]").arg(startimgrowid*100.0/ rowCount) qDebug() << "\r image create Row Range :\t"<<QString("[%1]").arg(startimgrowid*100.0/ rowCount)<< startimgrowid << "\t-\t" << startimgrowid + tempimgBlockline << "\t/\t" << rowCount;
<< startimgrowid << "\t-\t" << startimgrowid + tempimgBlockline << "\t/\t" << rowCount;
// 提取局部pixel x,y,z // 提取局部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_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_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<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); std::shared_ptr<std::complex<double>> imgArr = this->L1ds->getImageRaster(startimgrowid, tempimgBlockline);
// 获取回波 // 获取回波
long startechoid = 0; long startechoid = 0;
long iffeechoLen = PlusePoints; long iffeechoLen = PlusePoints;
for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) { for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) {
long tempechoBlockline = echoBlockline; long tempechoBlockline = echoBlockline;
if (startechoid + tempechoBlockline >= PRFCount) { if (startechoid + tempechoBlockline >= PRFCount) {
tempechoBlockline = PRFCount - startechoid; tempechoBlockline = PRFCount - startechoid;
} }
std::shared_ptr<std::complex<double>> echoArr = std::shared_ptr<std::complex<double>> echoArr = readDataArrComplex < std::complex<double>>(imagetimeimg,startechoid,long(0), tempechoBlockline, iffeechoLen, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);//; this->L0ds->getEchoArr(startechoid, tempechoBlockline);
readDataArrComplex < std::complex<double>>(imagetimeimg,startechoid,long(0), tempechoBlockline, iffeechoLen, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);//; this->L0ds->getEchoArr(startechoid, tempechoBlockline);
std::shared_ptr<double> antpx(new double[tempechoBlockline],delArrPtr); std::shared_ptr<double> antpx(new double[tempechoBlockline],delArrPtr);
std::shared_ptr<double> antpy(new double[tempechoBlockline], delArrPtr); std::shared_ptr<double> antpy(new double[tempechoBlockline], delArrPtr);
std::shared_ptr<double> antpz(new double[tempechoBlockline], delArrPtr); std::shared_ptr<double> antpz(new double[tempechoBlockline], delArrPtr);
// 复制
for (long anti = 0; anti < tempechoBlockline; anti++) { for (long anti = 0; anti < tempechoBlockline; anti++) { // 天线坐标
antpx.get()[anti] = Pxs.get()[anti + startechoid]; antpx.get()[anti] = Pxs.get()[anti + startechoid];
antpy.get()[anti] = Pys.get()[anti + startechoid]; antpy.get()[anti] = Pys.get()[anti + startechoid];
antpz.get()[anti] = Pzs.get()[anti + startechoid]; antpz.get()[anti] = Pzs.get()[anti + startechoid];
} }
TBPImageGPUAlg2( TBPImageGPUAlg2(
antpx, antpy, antpz, antpx, antpy, antpz, // 坐标
img_x, img_y, img_z, img_x, img_y, img_z, // 图像坐标
echoArr, imgArr, echoArr, // 回波
imgArr, // 图像
startfreq, dx, startfreq, dx,
Rnear, Rfar, refRange, Rnear, Rfar, refRange,
tempimgBlockline, colCount, tempimgBlockline, colCount,
tempechoBlockline, PlusePoints, tempechoBlockline, PlusePoints
startechoid,startimgrowid
); );
qDebug() << QString(" image block PRF:[%1] \t").arg((startechoid + tempechoBlockline) * 100.0 / PRFCount) qDebug() << QString(" image block PRF:[%1] \t").arg((startechoid + tempechoBlockline) * 100.0 / PRFCount)
<< startechoid << "\t-\t" << startechoid + tempechoBlockline; << startechoid << "\t-\t" << startechoid + tempechoBlockline;
@ -414,82 +410,6 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
return ErrorCode::SUCCESS; return ErrorCode::SUCCESS;
} }
void TBPImageAlgCls::EchoFreqToTime( )
{
// 读取数据
long PRFCount = this->L0ds->getPluseCount();
long inColCount = this->L0ds->getPlusePoints();
long outColCount = inColCount;// 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 * 3; //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 * outColCount), FreeCUDAHost);
std::shared_ptr<cuComplex> host_IFFTechoArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex)* tempechoBlockline * outColCount), FreeCUDAHost);
memset(host_echoArr.get(), 0, sizeof(cuComplex) * tempechoBlockline * outColCount);
#pragma omp parallel for
for (long ii = 0; ii < tempechoBlockline ; ii++) {
for (long jj = 0; jj < inColCount; jj++) {
host_echoArr.get()[ii* outColCount +jj] = make_cuComplex(echoArr.get()[ii * inColCount + jj].real(), echoArr.get()[ii * inColCount + jj].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, outColCount, outColCount);
FFTShift1D(device_IFFTechoArr.get(), tempechoBlockline, outColCount);
DeviceToHost(host_IFFTechoArr.get(), device_IFFTechoArr.get(), sizeof(cuComplex) * tempechoBlockline * outColCount);
#pragma omp parallel for
for (long ii = 0; ii < tempechoBlockline ; ii++) {
for (long jj = 0; jj < outColCount; jj++) {
IFFTArr.get()[ii * outColCount + jj] = std::complex<double>(host_IFFTechoArr.get()[ii * outColCount + jj].x, host_IFFTechoArr.get()[ii * outColCount + jj].y);
//IFFTArr.get()[ii * outColCount + jj] = std::complex<double>(host_echoArr.get()[ii * outColCount + jj].x, host_echoArr.get()[ii * outColCount + jj].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, void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antPy, std::shared_ptr<double> antPz,
@ -498,30 +418,11 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
std::shared_ptr<std::complex<double>> img_arr, std::shared_ptr<std::complex<double>> img_arr,
double freq, double dx, double Rnear, double Rfar, double refRange, double freq, double dx, double Rnear, double Rfar, double refRange,
long rowcount, long colcount, long rowcount, long colcount,
long prfcount, long freqcount, 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 * freqcount), 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_echoArrIFFT.get(), sizeof(cuComplex) * prfcount * freqcount);
// 结束傅里叶变换
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_antPx ((double*)mallocCUDAHost(sizeof(double) * prfcount),FreeCUDAHost);
std::shared_ptr<double> h_antPy ((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> h_antPz ((double*)mallocCUDAHost(sizeof(double) * prfcount),FreeCUDAHost);
@ -529,7 +430,7 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
std::shared_ptr<double> d_antPx ((double*)mallocCUDADevice(sizeof(double) * prfcount),FreeCUDADevice); 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_antPy ((double*)mallocCUDADevice(sizeof(double) * prfcount),FreeCUDADevice);
std::shared_ptr<double> d_antPz ((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_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_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> h_imgz((double*)mallocCUDAHost(sizeof(double) * rowcount * colcount),FreeCUDAHost);
@ -538,6 +439,21 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
std::shared_ptr<double> d_imgy ((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); std::shared_ptr<double> d_imgz ((double*)mallocCUDADevice(sizeof(double) * rowcount * colcount),FreeCUDADevice);
// 回波范围
std::shared_ptr<cuComplex> h_echoArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex) * prfcount * freqcount), FreeCUDAHost);
std::shared_ptr<cuComplex> d_echoArr((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * prfcount * freqcount), FreeCUDADevice);
// 成像结果
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 < prfcount; i++) {
for (long j = 0; j < freqcount; j++) {
h_echoArr.get()[i * freqcount + j] = make_cuComplex(echoArr.get()[i * freqcount + j].real(),
echoArr.get()[i * freqcount + j].imag());
}
}
// 天线位置 // 天线位置
for (long i = 0; i < prfcount; i++) { for (long i = 0; i < prfcount; i++) {
@ -555,6 +471,15 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
} }
} }
#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);
HostToDevice(h_echoArr.get(), d_echoArr.get(), sizeof(cuComplex) * prfcount * freqcount);
HostToDevice(h_antPx.get(), d_antPx.get(), sizeof(double) * prfcount); HostToDevice(h_antPx.get(), d_antPx.get(), sizeof(double) * prfcount);
HostToDevice(h_antPy.get(), d_antPy.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_antPz.get(), d_antPz.get(), sizeof(double) * prfcount);
@ -564,28 +489,15 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
HostToDevice(h_imgz.get(), d_imgz.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; double startlamda = LIGHTSPEED / freq;
TimeBPImage( TimeBPImage(
d_antPx.get(), d_antPy.get(), d_antPz.get(), d_antPx.get(), d_antPy.get(), d_antPz.get(),
d_imgx.get(), d_imgy.get(), d_imgz.get(), d_imgx.get(), d_imgy.get(), d_imgz.get(),
d_echoArrIFFT.get(), prfcount, IFFTPadNum, d_echoArr.get(), prfcount, freqcount,
d_imgArr.get(), rowcount, colcount, d_imgArr.get(), rowcount, colcount,
startlamda, Rnear, dx, refRange startlamda, Rnear, dx, refRange,Rfar
); );
// Device -> Host // Device -> Host
@ -605,16 +517,6 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
} }
void TBPImageAlgCls::setGPU(bool flag) void TBPImageAlgCls::setGPU(bool flag)
{ {
this->GPURUN = flag; this->GPURUN = flag;
@ -625,186 +527,79 @@ bool TBPImageAlgCls::getGPU( )
return this->GPURUN; return this->GPURUN;
} }
void TBPImageAlgCls::EchoFreqToTime()
/**
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 PRFCount = this->L0ds->getPluseCount();
long PlusePoints = this->L0ds->getPlusePoints(); long inColCount = this->L0ds->getPlusePoints();
long outColCount = inColCount;// 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 * 3; //1GB
echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;
double Rnear = this->L1ds->getNearRange(); long startechoid = 0;
double Rfar = this->L1ds->getFarRange(); for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) {
double fs = this->L1ds->getFs();
double dx = LIGHTSPEED / 2 / fs;
double factorj = this->L1ds->getCenterFreq() * 4 * M_PI / LIGHTSPEED * 1e9;
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);
Eigen::MatrixXcd echo = Eigen::MatrixXcd::Zero(PRFCount, PlusePoints); std::shared_ptr<cuComplex> host_echoArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex) * tempechoBlockline * outColCount), FreeCUDAHost);
{ std::shared_ptr<cuComplex> host_IFFTechoArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex) * tempechoBlockline * outColCount), FreeCUDAHost);
std::shared_ptr<std::complex<double>> echodata = this->L0ds->getEchoArr();
for (long i = 0; i < PRFCount; i++) { memset(host_echoArr.get(), 0, sizeof(cuComplex) * tempechoBlockline * outColCount);
for (long j = 0; j < PlusePoints; j++) { #pragma omp parallel for
echo(i, j) = echodata.get()[i * PlusePoints + j]; for (long ii = 0; ii < tempechoBlockline; ii++) {
for (long jj = 0; jj < inColCount; jj++) {
host_echoArr.get()[ii * outColCount + jj] = make_cuComplex(echoArr.get()[ii * inColCount + jj].real(), echoArr.get()[ii * inColCount + jj].imag());
} }
} }
echodata.reset(); #pragma omp parallel for
for (long ii = 0; ii < tempechoBlockline * outColCount; ii++) {
host_IFFTechoArr.get()[ii] = make_cuComplex(0, 0);
} }
Eigen::MatrixXd pixelX = Eigen::MatrixXd::Zero(rowCount, colCount); std::shared_ptr<cuComplex> device_echoArr((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * tempechoBlockline * inColCount), FreeCUDADevice);
Eigen::MatrixXd pixelY = Eigen::MatrixXd::Zero(rowCount, colCount); std::shared_ptr<cuComplex> device_IFFTechoArr((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * tempechoBlockline * outColCount), FreeCUDADevice);
Eigen::MatrixXd pixelZ = Eigen::MatrixXd::Zero(rowCount, colCount);
Eigen::MatrixXd Pxs = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1); HostToDevice(host_echoArr.get(), device_echoArr.get(), sizeof(cuComplex) * tempechoBlockline * inColCount);
Eigen::MatrixXd Pys = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1); HostToDevice(host_IFFTechoArr.get(), device_IFFTechoArr.get(), sizeof(cuComplex) * tempechoBlockline * outColCount);
Eigen::MatrixXd Pzs = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1); CUDAIFFT(device_echoArr.get(), device_IFFTechoArr.get(), tempechoBlockline, outColCount, outColCount);
// 图像网格坐标 FFTShift1D(device_IFFTechoArr.get(), tempechoBlockline, outColCount);
{
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; DeviceToHost(host_IFFTechoArr.get(), device_IFFTechoArr.get(), sizeof(cuComplex) * tempechoBlockline * outColCount);
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 #pragma omp parallel for
for (startLine = 0; startLine < rowCount; startLine = startLine + BlockLine) { // 图像大小 for (long ii = 0; ii < tempechoBlockline; ii++) {
long stepLine = startLine + BlockLine < rowCount ? BlockLine : rowCount - startLine; for (long jj = 0; jj < outColCount; jj++) {
long imageRowID = startLine; // IFFTArr.get()[ii * outColCount + jj] = std::complex<double>(host_IFFTechoArr.get()[ii * outColCount + jj].x, host_IFFTechoArr.get()[ii * outColCount + jj].y);
//IFFTArr.get()[ii * outColCount + jj] = std::complex<double>(host_echoArr.get()[ii * outColCount + jj].x, host_echoArr.get()[ii * outColCount + jj].y);
//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);// 信号* 相位校正
} }
} }
} outTimeEchoImg.saveImage(IFFTArr, startechoid, 0, tempechoBlockline, outColCount, 1);
omp_set_lock(&lock); // 保存文件
processValue = processValue + BlockLine; qDebug() << QString(" image block PRF:[%1] \t").arg((startechoid + tempechoBlockline) * 100.0 / PRFCount)
this->L1ds->saveImageRaster(imagarr, 0, rowCount); << startechoid << "\t-\t" << startechoid + tempechoBlockline;
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(); return;
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;
} }
*/

16
Toolbox/SimulationSARTool/SimulationSAR/TBPImageAlgCls.h
View File

@ -79,6 +79,18 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
std::shared_ptr<std::complex<double>> img_arr, std::shared_ptr<std::complex<double>> img_arr,
double freq, double dx, double Rnear, double Rfar,double refRange, double freq, double dx, double Rnear, double Rfar,double refRange,
long rowcount, long colcount, long rowcount, long colcount,
long prfcount, long freqcount, long prfcount, long freqcount
long startPRFId, long startRowID
); );
//
//void TBPImageGridNet(
// 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
//);