matlab代码转换为C++代码

pull/6/head
陈增辉 2025-03-03 17:50:28 +08:00
parent 16b8abb6c5
commit c9cd0c3ee1
8 changed files with 270 additions and 305 deletions

View File

@ -297,12 +297,18 @@ struct RadiationPatternGainPoint {
inline void delArrPtr(void* p)
{
if (nullptr == p || NULL == p) {
return;
}
delete[] p;
p = nullptr;
}
inline void delPointer(void* p)
{
if (nullptr == p || NULL == p) {
return;
}
delete p;
p = nullptr;
}

View File

@ -252,6 +252,14 @@ extern "C" void FFTShift1D(cuComplex* d_data, int batch_size, int signal_length)
cudaDeviceSynchronize();
}
extern "C" GPUBASELIBAPI void shared_complexPtrToHostCuComplex(std::complex<double>* src, cuComplex* dst, long len)
{
for (long i = 0; i < len; i++) {
dst[i] = make_cuComplex(src[i].real(), src[i].imag());
}
return ;
}
extern __global__ void CUDA_D_sin(double* y, double* X, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {

View File

@ -109,5 +109,6 @@ extern "C" GPUBASELIBAPI void CUDAIFFTScale(cuComplex* inArr, cuComplex* outArr,
extern "C" GPUBASELIBAPI void CUDAIFFT(cuComplex* inArr, cuComplex* outArr, long InRowCount, long InColCount, long outColCount);
extern "C" GPUBASELIBAPI void FFTShift1D(cuComplex* d_data, int batch_size, int signal_length);
extern "C" GPUBASELIBAPI void shared_complexPtrToHostCuComplex(std::complex<double>* src, cuComplex* dst, long len);
#endif
#endif

View File

@ -4,13 +4,18 @@
#include <cuda_runtime.h>
#include <iostream>
#include <memory>
#include <vector>
#include <cmath>
#include <complex>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include <cufft.h>
#include <cufftw.h>
#include <cufftXt.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#include "GPUBPTool.cuh"
@ -18,20 +23,20 @@
__global__ void phaseCompensationKernel(cufftComplex* phdata, const float* Freq, float r, int K, int Na) {
__global__ void phaseCompensationKernel(cufftComplex* phdata, const double* Freq, double r, int K, int Na) {
int freqIdx = blockIdx.x * blockDim.x + threadIdx.x;
int pulseIdx = blockIdx.y * blockDim.y + threadIdx.y;
if (freqIdx >= K || pulseIdx >= Na) return;
int idx = pulseIdx * K + freqIdx;
float phase = 4 * PI * Freq[freqIdx] * r / c;
float cos_phase = cosf(phase);
float sin_phase = sinf(phase);
double phase = 4 * PI * Freq[freqIdx] * r / c;
double cos_phase = cosf(phase);
double sin_phase = sinf(phase);
cufftComplex ph = phdata[idx];
float new_real = ph.x * cos_phase - ph.y * sin_phase;
float new_imag = ph.x * sin_phase + ph.y * cos_phase;
double new_real = ph.x * cos_phase - ph.y * sin_phase;
double new_imag = ph.x * sin_phase + ph.y * cos_phase;
phdata[idx] = make_cuComplex(new_real, new_imag);
}
@ -47,9 +52,9 @@ __global__ void fftshiftKernel(cufftComplex* data, int Nfft, int Np) {
}
}
__global__ void processPulseKernel(int nx, int ny, const float* x_mat, const float* y_mat, const float* z_mat,
float AntX, float AntY, float AntZ, float R0, float minF,
const cufftComplex* rc_pulse, float r_start, float dr, int nR,
__global__ void processPulseKernel(int nx, int ny, const double* x_mat, const double* y_mat, const double* z_mat,
double AntX, double AntY, double AntZ, double R0, double minF,
const cufftComplex* rc_pulse, double r_start, double dr, int nR,
cufftComplex* im_final) {
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
@ -57,18 +62,18 @@ __global__ void processPulseKernel(int nx, int ny, const float* x_mat, const flo
if (x >= nx || y >= ny) return;
int idx = x * ny + y;
float dx = AntX - x_mat[idx];
float dy = AntY - y_mat[idx];
float dz = AntZ - z_mat[idx];
float dR = sqrtf(dx * dx + dy * dy + dz * dz) - R0;
double dx = AntX - x_mat[idx];
double dy = AntY - y_mat[idx];
double dz = AntZ - z_mat[idx];
double dR = sqrtf(dx * dx + dy * dy + dz * dz) - R0;
// Range check
if (dR < r_start || dR >= (r_start + dr * (nR - 1))) return;
// Linear interpolation
float pos = (dR - r_start) / dr;
double pos = (dR - r_start) / dr;
int index = (int)floorf(pos);
float weight = pos - index;
double weight = pos - index;
if (index < 0 || index >= nR - 1) return;
@ -79,9 +84,9 @@ __global__ void processPulseKernel(int nx, int ny, const float* x_mat, const flo
rc_interp.y = rc_low.y * (1 - weight) + rc_high.y * weight;
// Phase correction
float phase = 4 * PI * minF * dR / c;
float cos_phase = cosf(phase);
float sin_phase = sinf(phase);
double phase = 4 * PI * minF * dR / c;
double cos_phase = cosf(phase);
double sin_phase = sinf(phase);
cufftComplex phCorr;
phCorr.x = rc_interp.x * cos_phase - rc_interp.y * sin_phase;
@ -99,7 +104,7 @@ void bpBasic0CUDA(GPUDATA& data, int flag) {
dim3 grid((data.K + 15) / 16, (data.Np + 15) / 16);
phaseCompensationKernel << <grid, block >> > (data.phdata, data.Freq, data.R0, data.K, data.Np);
cudaCheckError(cudaDeviceSynchronize());
data.R0 = data.r; // ¼ÙÉèdata.rÒÑÕýÈ·ÉèÖÃ
//data.R0 = data.r; // 假设data.r已正确设置
}
// FFT´¦Àí
@ -115,8 +120,8 @@ void bpBasic0CUDA(GPUDATA& data, int flag) {
cudaCheckError(cudaDeviceSynchronize());
// ͼÏñÖØ½¨
float r_start = data.r_vec[0];
float dr = (data.r_vec[data.Nfft - 1] - r_start) / (data.Nfft - 1);
double r_start = data.r_vec[0];
double dr = (data.r_vec[data.Nfft - 1] - r_start) / (data.Nfft - 1);
dim3 block(16, 16);
dim3 grid((data.nx + 15) / 16, (data.ny + 15) / 16);
@ -136,17 +141,104 @@ void bpBasic0CUDA(GPUDATA& data, int flag) {
cudaCheckError(cudaDeviceSynchronize());
}
void computeRvec(GPUDATA& data) {
// 计算maxWr需要先计算deltaF
double deltaF = data.deltaF; // 从输入参数获取
double maxWr = 299792458.0f / (2.0f * deltaF);
void initGPUData(GPUDATA& h_data, GPUDATA& d_data) {
}
// 生成r_vec主机端
double* r_vec_host=new double[data.Nfft];
const double step = maxWr / data.Nfft;
const double start = -data.Nfft / 2.0f * step;
void freeGPUData(GPUDATA& d_data) {
for (int i = 0; i < data.Nfft; ++i) {
r_vec_host[i] = start + i * step;
}
data.r_vec = r_vec_host;
}
void initGPUData(GPUDATA& h_data, GPUDATA& d_data) {
d_data.AntX = (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
d_data.AntY = (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
d_data.AntZ = (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
d_data.minF = (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
d_data.x_mat = (double*)mallocCUDADevice(sizeof(double) * h_data.nx * h_data.ny);
d_data.y_mat = (double*)mallocCUDADevice(sizeof(double) * h_data.nx * h_data.ny);
d_data.z_mat = (double*)mallocCUDADevice(sizeof(double) * h_data.nx * h_data.ny);
d_data.r_vec = (double*)mallocCUDADevice(sizeof(double) * h_data.Nfft);
d_data.Freq = (double*)mallocCUDADevice(sizeof(double) * h_data.Nfft);
d_data.phdata = (cufftComplex*)mallocCUDADevice(sizeof(cufftComplex) * h_data.K * h_data.Np);
d_data.im_final = (cufftComplex*)mallocCUDADevice(sizeof(cufftComplex) * h_data.nx * h_data.ny);
HostToDevice(h_data.AntX, d_data.AntX,sizeof(double) * h_data.Np);
HostToDevice(h_data.AntY, d_data.AntY,sizeof(double) * h_data.Np);
HostToDevice(h_data.AntZ, d_data.AntZ,sizeof(double) * h_data.Np);
HostToDevice(h_data.minF, d_data.minF,sizeof(double) * h_data.Np);
HostToDevice(h_data.x_mat, d_data.x_mat,sizeof(double) * h_data.nx * h_data.ny);
HostToDevice(h_data.y_mat, d_data.y_mat,sizeof(double) * h_data.nx * h_data.ny);
HostToDevice(h_data.z_mat, d_data.z_mat,sizeof(double) * h_data.nx * h_data.ny);
HostToDevice(h_data.Freq, d_data.Freq, sizeof(double) * h_data.Nfft);
HostToDevice(h_data.r_vec, d_data.r_vec, sizeof(double) * h_data.Nfft);
HostToDevice(h_data.phdata, d_data.phdata, sizeof(cufftComplex) * h_data.K * h_data.Np);
HostToDevice(h_data.im_final, d_data.im_final, sizeof(cufftComplex) * h_data.nx * h_data.ny);
// 拷贝标量参数
d_data.Nfft = h_data.Nfft;
d_data.K = h_data.K;
d_data.Np = h_data.Np;
d_data.nx = h_data.nx;
d_data.ny = h_data.ny;
d_data.R0 = h_data.R0;
d_data.deltaF = h_data.deltaF;
}
void freeGPUData(GPUDATA& d_data) {
FreeCUDADevice((d_data.AntX));
FreeCUDADevice((d_data.AntY));
FreeCUDADevice((d_data.AntZ));
FreeCUDADevice((d_data.minF));
FreeCUDADevice((d_data.x_mat));
FreeCUDADevice((d_data.y_mat));
FreeCUDADevice((d_data.z_mat));
FreeCUDADevice((d_data.r_vec));
FreeCUDADevice((d_data.Freq));
FreeCUDADevice((d_data.phdata));
FreeCUDADevice((d_data.im_final));
}
void freeHostData(GPUDATA& h_data) {
FreeCUDAHost((h_data.AntX));
FreeCUDAHost((h_data.AntY));
FreeCUDAHost((h_data.AntZ));
FreeCUDAHost((h_data.minF));
FreeCUDAHost((h_data.x_mat));
FreeCUDAHost((h_data.y_mat));
FreeCUDAHost((h_data.z_mat));
FreeCUDAHost((h_data.r_vec));
FreeCUDAHost((h_data.Freq));
FreeCUDAHost((h_data.phdata));
FreeCUDAHost((h_data.im_final));
}
void BPBasic0(GPUDATA& h_data)
{
GPUDATA d_data;
initGPUData(h_data, d_data);
bpBasic0CUDA(d_data, 0);
DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * h_data.nx * h_data.ny);
freeGPUData(d_data);
}
//int main() {
// GPUDATA h_data, d_data;
//

View File

@ -17,16 +17,21 @@ inline void gpuAssert(cudaError_t code, const char* file, int line) {
struct GPUDATA {
int Nfft, K, Np, nx, ny; // 傅里叶点数、频点数、脉冲数、图像列、图像行
float* AntX, * AntY, * AntZ, * minF; // 天线坐标、起始频率
float* x_mat, * y_mat, * z_mat;// 地面坐标
float* r_vec; // 坐标范围
double* AntX, * AntY, * AntZ, * minF; // 天线坐标、起始频率
double* x_mat, * y_mat, * z_mat;// 地面坐标
double* r_vec; // 坐标范围
double* Freq;// 频率
cufftComplex* phdata;// 回波
cufftComplex* im_final;// 图像
float R0; // 参考斜距
double R0; // 参考斜距
double deltaF; // 频点范围
};
extern "C" {
void computeRvec(GPUDATA& data);
void initGPUData(GPUDATA& h_data, GPUDATA& d_data);
void freeGPUData(GPUDATA& d_data);
void freeHostData(GPUDATA& d_data);
void BPBasic0(GPUDATA& h_data);
};

View File

@ -1,4 +1,4 @@
#include "stdafx.h"
#include "stdafx.h"
#include "TBPImageAlgCls.h"
#include <QDateTime>
#include <QDebug>
@ -9,6 +9,8 @@
#include "GPUTool.cuh"
#include "GPUTBPImage.cuh"
#include "ImageOperatorBase.h"
#include "BPBasic0_CUDA.cuh"
void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZPath)
{
@ -19,7 +21,7 @@ void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZ
double dx = LIGHTSPEED / 2.0 / bandwidth; // c/2b
// 创建坐标系统
// 创建坐标系统
long prfcount = echoL0ds->getPluseCount();
long freqcount = echoL0ds->getPlusePoints();
Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
@ -55,7 +57,7 @@ void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZ
double Pz = 0;
for (long i = 0; i < prfcount; i++) {
Pxs.get()[i] = antpos.get()[i * 19 + 1]; // 卫星坐标
Pxs.get()[i] = antpos.get()[i * 19 + 1]; // 卫星坐标
Pys.get()[i] = antpos.get()[i * 19 + 2];
Pzs.get()[i] = antpos.get()[i * 19 + 3];
AntDirectX.get()[i] = antpos.get()[i * 19 + 13];// zero doppler
@ -67,7 +69,7 @@ void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZ
AntDirectZ.get()[i] * AntDirectZ.get()[i]);
AntDirectX.get()[i] = AntDirectX.get()[i] / NormAnt;
AntDirectY.get()[i] = AntDirectY.get()[i] / NormAnt;
AntDirectZ.get()[i] = AntDirectZ.get()[i] / NormAnt;// 归一化
AntDirectZ.get()[i] = AntDirectZ.get()[i] / NormAnt;// 归一化
}
antpos.reset();
}
@ -126,7 +128,7 @@ void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZ
d_AntDirectX.get(), d_AntDirectY.get(), d_AntDirectZ.get(),
d_demx.get(), d_demy.get(), d_demz.get(),
prfcount, tempechocol, 1000,
Rnear+dx* startcolidx, dx, refRange // 更新最近修读
Rnear+dx* startcolidx, dx, refRange // 更新最近修读
);
DeviceToHost(h_demx.get(), d_demx.get(), sizeof(double) * prfcount * tempechocol);
@ -209,10 +211,10 @@ std::shared_ptr<SARSimulationImageL1Dataset> TBPImageAlgCls::getImageL0()
ErrorCode TBPImageAlgCls::Process(long num_thread)
{
qDebug() << u8"开始成像";
qDebug() << u8"开始成像";
qDebug() << u8"创建成像平面的XYZ";
qDebug() << u8"创建成像平面的XYZ";
QString outRasterXYZ = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_xyz.bin");
CreatePixelXYZ(this->L0ds, outRasterXYZ);
@ -223,12 +225,12 @@ ErrorCode TBPImageAlgCls::Process(long num_thread)
ErrorCode TBPImageAlgCls::ProcessWithGridNet(long num_thread,QString xyzRasterPath)
{
this->outRasterXYZPath = xyzRasterPath;
qDebug() << u8"频域回波-> 时域回波";
this->TimeEchoDataPath = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_Timeecho.bin");
this->EchoFreqToTime();
//qDebug() << u8"频域回波-> 时域回波";
//this->TimeEchoDataPath = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_Timeecho.bin");
//this->EchoFreqToTime();
// 初始化Raster
qDebug() << u8"初始化影像";
// 初始化Raster
qDebug() << u8"初始化影像";
long imageheight = this->L1ds->getrowCount();
long imagewidth = this->L1ds->getcolCount();
@ -255,13 +257,13 @@ ErrorCode TBPImageAlgCls::ProcessWithGridNet(long num_thread,QString xyzRasterPa
}
qDebug() << u8"频域回波-> 时域回波 结束";
qDebug() << u8"频域回波-> 时域回波 结束";
if (GPURUN) {
return this->ProcessGPU();
}
else {
QMessageBox::information(nullptr, u8"提示", u8"目前只支持显卡");
QMessageBox::information(nullptr, u8"提示", u8"目前只支持显卡");
return ErrorCode::FAIL;
}
return ErrorCode::SUCCESS;
@ -273,49 +275,28 @@ ErrorCode TBPImageAlgCls::ProcessWithGridNet(long num_thread,QString xyzRasterPa
ErrorCode TBPImageAlgCls::ProcessGPU()
{
// 常用参数
long rowCount = this->L1ds->getrowCount();
long colCount = this->L1ds->getcolCount();
long pixelCount = rowCount * colCount;
// 回波大小
long freqpoint = this->L0ds->getPlusePoints();
long PRFCount = this->L0ds->getPluseCount();
long PlusePoints = this->L0ds->getPlusePoints();
long bandwidth = this->L0ds->getBandwidth();
// 图像范围
long imWidth = this->L1ds->getcolCount();
long imHeight = this->L1ds->getrowCount();
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); // 更新斜距距离
// 内存分配大小
long echoBlockline = Memory1GB / 8 / 2 / freqpoint * 4; //2GB
echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;
long imageBlockline = Memory1GB / 8 / 2 / imHeight * 4; //2GB
imageBlockline = imageBlockline < 1 ? 1 : imageBlockline;
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()],delArrPtr);
std::shared_ptr<double> Pys (new double[this->L0ds->getPluseCount()],delArrPtr);
std::shared_ptr<double> Pzs (new double[this->L0ds->getPluseCount()],delArrPtr);
// 天线坐标
std::shared_ptr<double> Pxs(new double[this->L0ds->getPluseCount()], delArrPtr);
std::shared_ptr<double> Pys(new double[this->L0ds->getPluseCount()], delArrPtr);
std::shared_ptr<double> Pzs(new double[this->L0ds->getPluseCount()], delArrPtr);
{
std::shared_ptr<double> antpos = this->L0ds->getAntPos();
@ -324,10 +305,10 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
double Py = 0;
double Pz = 0;
for (long i = 0; i < PRFCount; 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];
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;
@ -336,229 +317,107 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
}
double startFreq=this->L0ds->getCenterFreq() - this->L0ds->getBandwidth() / 2;
gdalImage imageXYZ(this->outRasterXYZPath); // 图像坐标
for (long img_rid = 0; img_rid < imHeight; img_rid = img_rid + imageBlockline) {
// 获取坐标范围
long imrowcount = imageBlockline;
long imcolcount = imWidth;
std::shared_ptr<double> img_x = readDataArr<double>(imageXYZ, img_rid, 0, imrowcount, imcolcount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> img_y = readDataArr<double>(imageXYZ, img_rid, 0, imrowcount, imcolcount, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> img_z = readDataArr<double>(imageXYZ, img_rid, 0, imrowcount, imcolcount, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<std::complex<double>> imgArr = this->L1ds->getImageRaster(img_rid, imrowcount); // 回波值
//shared_complexPtrToHostCuComplex(std::complex<double>*src, cuComplex * dst, long len)
// 处理
GPUDATA h_data;
// 计算成像的基本参数
// 距离向分辨率
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;
// 方位向分辨率
h_data.Nfft = freqpoint;
h_data.K = freqpoint;
h_data.deltaF = this->L0ds->getBandwidth() / (freqpoint - 1);
computeRvec(h_data);
long echoBlockline = Memory1GB / 8 / 2 / PlusePoints * 2; //2GB
echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;
h_data.x_mat = img_x.get();//地面
h_data.y_mat = img_y.get();
h_data.z_mat = img_z.get();
h_data.nx = imcolcount;
h_data.ny = imrowcount;
h_data.R0 = refRange;// 参考斜距
h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * imrowcount * imcolcount);
shared_complexPtrToHostCuComplex(imgArr.get(), h_data.im_final, imrowcount * imcolcount);
for (long prfid = 0; prfid < PRFCount; prfid = prfid + echoBlockline) {
long ehcoprfcount = echoBlockline;
long echofreqpoint = freqpoint;
std::shared_ptr<std::complex<double>> echoArr = this->L0ds->getEchoArr(prfid, ehcoprfcount);
long imageBlockline = Memory1GB / 8 / 2 / colCount * 2; //2GB
imageBlockline = imageBlockline < 1 ? 1 : imageBlockline;
gdalImage imageXYZ(this->outRasterXYZPath);
gdalImageComplex imagetimeimg(this->TimeEchoDataPath);
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;
long iffeechoLen = PlusePoints;
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 = 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> antpy(new double[tempechoBlockline], delArrPtr);
std::shared_ptr<double> antpz(new double[tempechoBlockline], 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];
// 复制天线方向图
std::shared_ptr<double> antpx(new double[ehcoprfcount], delArrPtr);
std::shared_ptr<double> antpy(new double[ehcoprfcount], delArrPtr);
std::shared_ptr<double> antpz(new double[ehcoprfcount], delArrPtr);
for (long anti = 0; anti < ehcoprfcount; anti++) { // ÌìÏß×ø±ê
antpx.get()[anti] = Pxs.get()[anti + prfid];
antpy.get()[anti] = Pys.get()[anti + prfid];
antpz.get()[anti] = Pzs.get()[anti + prfid];
if (abs(antpx.get()[anti]) < 10 || abs(antpy.get()[anti]) < 10 || abs(antpz.get()[anti]) < 10) {
qDebug()<<anti<<":" << antpx.get()[anti] << "," << antpy.get()[anti] << ","<< antpz.get()[anti] ;
qDebug() << anti << ":" << antpx.get()[anti] << "," << antpy.get()[anti] << "," << antpz.get()[anti];
}
}
double maxdouble = 0;
for (long i = 0; i < tempechoBlockline; i++) {
for (long j = 0; j < PlusePoints; j++) {
if (echoArr.get()[i * PlusePoints + j].real() > maxdouble) {
maxdouble = echoArr.get()[i * PlusePoints + j].real();
}
if (echoArr.get()[i * PlusePoints + j].imag() > maxdouble) {
maxdouble = echoArr.get()[i * PlusePoints + j].imag();
}
}
// 起始频率
h_data.minF = (double*)mallocCUDAHost(sizeof(double) * ehcoprfcount);
for (long anti = 0; anti < ehcoprfcount; anti++) {
h_data.minF[anti] = startFreq;
}
h_data.Freq = (double*)mallocCUDAHost(sizeof(double) * freqpoint);
for (long fid = 0; fid < freqpoint; fid++) {
h_data.Freq[fid] = startFreq + fid * h_data.deltaF;
}
h_data.AntX = antpx.get(); // 天线
h_data.AntY = antpy.get();
h_data.AntZ = antpz.get();
h_data.Np = ehcoprfcount;
h_data.phdata= (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * ehcoprfcount * echofreqpoint);
shared_complexPtrToHostCuComplex(echoArr.get(), h_data.phdata, ehcoprfcount * echofreqpoint);
printf("echo %f\n", maxdouble);
BPBasic0(h_data); //BP
TBPImageGPUAlg2(
antpx, antpy, antpz, // 坐标
img_x, img_y, img_z, // 图像坐标
echoArr, // 回波
imgArr, // 图像
startfreq, dx,
Rnear, Rfar, refRange,
tempimgBlockline, colCount,
tempechoBlockline, PlusePoints
);
qDebug() << QString(" image block PRF:[%1] \t").arg((startechoid + tempechoBlockline) * 100.0 / PRFCount)<< startechoid << "\t-\t" << startechoid + tempechoBlockline;
qDebug() << QString(" image block PRF:[%1] \t").arg((img_rid + imrowcount) * 100.0 / imHeight) << QString(" echo [%1] prfidx:\t ").arg((prfid + ehcoprfcount)*100.0/ PRFCount)<< prfid << "\t-\t" << prfid + ehcoprfcount;
this->L1ds->saveImageRaster(imgArr, img_rid, imrowcount);
}
this->L1ds->saveImageRaster(imgArr, startimgrowid, tempimgBlockline);
freeHostData(h_data);
}
qDebug() << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz") << "] image writing:\t" << this->L1ds->getxmlFilePath();
this->L1ds->saveToXml();
this->L1ds->saveToXml();
return ErrorCode::SUCCESS;
}
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
)
ErrorCode TBPImageAlgCls::BPProcessBlockGPU()
{
// 天线
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);
// 回波范围
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
double maxdouble = 0;
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());
if (h_echoArr.get()[i * freqcount + j].x > maxdouble) {
maxdouble = h_echoArr.get()[i * freqcount + j].x;
}
if (h_echoArr.get()[i * freqcount + j].y > maxdouble) {
maxdouble = h_echoArr.get()[i * freqcount + j].y;
}
}
}
printf("h_echoArr max: %e\n", maxdouble);
// 天线位置
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];
}
}
#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_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);
// 直接使用
double startlamda = LIGHTSPEED / freq;
TimeBPImage(
d_antPx.get(), d_antPy.get(), d_antPz.get(),
d_imgx.get(), d_imgy.get(), d_imgz.get(),
d_echoArr.get(), prfcount, freqcount,
d_imgArr.get(), rowcount, colcount,
startlamda, Rnear, dx, refRange,Rfar
);
// Device -> Host
DeviceToHost(h_imgArr.get(), d_imgArr.get(), sizeof(cuComplex)* rowcount* colcount);
maxdouble = 0;
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);
if (img_arr.get()[i * colcount + j].real() > maxdouble) {
maxdouble = img_arr.get()[i * colcount + j].real();
}
if (img_arr.get()[i * colcount + j].imag() > maxdouble) {
maxdouble = img_arr.get()[i * colcount + j].imag();
}
}
}
printf("echoImage max: %f\n", maxdouble);
return ErrorCode::SUCCESS;
}
void TBPImageAlgCls::setGPU(bool flag)
{
@ -572,7 +431,7 @@ bool TBPImageAlgCls::getGPU( )
void TBPImageAlgCls::EchoFreqToTime()
{
// 读取数据
// 读取数据
long PRFCount = this->L0ds->getPluseCount();
long inColCount = this->L0ds->getPlusePoints();
@ -583,10 +442,10 @@ void TBPImageAlgCls::EchoFreqToTime()
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;

View File

@ -51,39 +51,30 @@ public:
public:
ErrorCode Process(long num_thread);
ErrorCode ProcessWithGridNet(long num_thread,QString xyzRasterPath);
ErrorCode ProcessWithGridNet(long num_thread, QString xyzRasterPath);
void setGPU(bool flag);
bool getGPU( );
bool getGPU();
private:
//ErrorCode ProcessCPU(long num_thread);
ErrorCode ProcessGPU();
ErrorCode BPProcessBlockGPU();
private://ÁÙʱ³ÉÔ±±äÁ¿
QString TimeEchoDataPath;
long TimeEchoRowCount;
long TimeEchoColCount;
void EchoFreqToTime( );
void EchoFreqToTime();
};
void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds,QString outPixelXYZPath);
void TBPImageProcess(QString echofile,QString outImageFolder,QString imagePlanePath,long num_thread);
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
);
//
//void TBPImageGridNet(
// std::shared_ptr<double> antPx, std::shared_ptr<double> antPy, std::shared_ptr<double> antPz,

View File

@ -129,6 +129,9 @@
<GenerateRelocatableDeviceCode>true</GenerateRelocatableDeviceCode>
<CodeGeneration>compute_86,sm_86</CodeGeneration>
</CudaCompile>
<Link>
<AdditionalDependencies>cufft.lib;%(AdditionalDependencies);cudart.lib;cudadevrt.lib</AdditionalDependencies>
</Link>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|ARM'">
<ClCompile>