RasterProcessTool/Toolbox/SimulationSARTool/UnitTestMain.cpp

427 lines
17 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

/*
* 仿真测试代码
*
*
*/
#include <stdio.h>
#include <QString>
#include "EchoDataFormat.h"
#include "ImageShowDialogClass.h"
#include "GPUBaseLibAPI.h"
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#include "GPUBPTool.cuh"
#include "BPBasic0_CUDA.cuh"
#include "ImageOperatorBase.h"
#include "GPUBpSimulation.cuh"
#include <QApplication>
void testSimualtionEchoPoint() {
GPUDATA h_data{};
/** 1. 轨道 **************************************************************************************************/
qDebug() << u8"1.轨道文件读取中。。。";
QString inGPSPath = u8"C:\\Users\\30453\\Desktop\\script\\data\\GF3_Simulation.gpspos.data";
long gpspoints = gdalImage(inGPSPath).height;
std::shared_ptr<SatelliteAntPos> antpos = SatelliteAntPosOperator::readAntPosFile(inGPSPath, gpspoints);
h_data.AntX = (double*)mallocCUDAHost(sizeof(double) * gpspoints);
h_data.AntY = (double*)mallocCUDAHost(sizeof(double) * gpspoints);
h_data.AntZ = (double*)mallocCUDAHost(sizeof(double) * gpspoints);
for (long i = 0; i < gpspoints; i = i + 1) {
h_data.AntX[i] = antpos.get()[i].Px;
h_data.AntY[i] = antpos.get()[i].Py;
h_data.AntZ[i] = antpos.get()[i].Pz;
}
//gpspoints = gpspoints / 2;
qDebug() << "GPS points Count:\t" << gpspoints;
qDebug() << "PRF 0:\t " << QString("%1,%2,%3").arg(h_data.AntX[0]).arg(h_data.AntY[0]).arg(h_data.AntZ[0]);
qDebug() << "PRF " << gpspoints - 1 << ":\t " << QString("%1,%2,%3")
.arg(h_data.AntX[gpspoints - 1])
.arg(h_data.AntY[gpspoints - 1])
.arg(h_data.AntZ[gpspoints - 1]);
/** 2. 成像网格 **************************************************************************************************/
qDebug() << u8"轨道文件读取结束\n2.成像网格读取。。。";
QString demxyzPath = u8"D:\\Programme\\vs2022\\RasterMergeTest\\simulationData\\demdataset\\demxyz.bin";
gdalImage demgridimg(demxyzPath);
long dem_rowCount = demgridimg.height;
long dem_ColCount = demgridimg.width;
std::shared_ptr<double> demX = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> demY = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> demZ = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
h_data.x_mat = (double*)mallocCUDAHost(sizeof(double) * dem_rowCount * dem_ColCount);
h_data.y_mat = (double*)mallocCUDAHost(sizeof(double) * dem_rowCount * dem_ColCount);
h_data.z_mat = (double*)mallocCUDAHost(sizeof(double) * dem_rowCount * dem_ColCount);
for (long i = 0; i < dem_rowCount; i++) {
for (long j = 0; j < dem_ColCount; j++) {
h_data.x_mat[i * dem_ColCount + j] = demX.get()[i * dem_ColCount + j];
h_data.y_mat[i * dem_ColCount + j] = demY.get()[i * dem_ColCount + j];
h_data.z_mat[i * dem_ColCount + j] = demZ.get()[i * dem_ColCount + j];
}
}
qDebug() << "dem row Count:\t" << dem_rowCount;
qDebug() << "dem col Count:\t" << dem_ColCount;
qDebug() << u8"成像网格读取结束";
/** 3. 频率网格 **************************************************************************************************/
qDebug() << u8"3.处理频率";
double centerFreq = 5.3e9;
double bandwidth = 40e6;
long freqpoints = 2048;
std::shared_ptr<double> freqlist(getFreqPoints_mallocHost(centerFreq - bandwidth / 2, centerFreq + bandwidth / 2, freqpoints), FreeCUDAHost);
std::shared_ptr<double> minF(new double[gpspoints], delArrPtr);
for (long i = 0; i < gpspoints; i++) {
minF.get()[i] = freqlist.get()[0];
}
h_data.deltaF = bandwidth / (freqpoints - 1);
qDebug() << "start Freq:\t" << centerFreq - bandwidth / 2;
qDebug() << "end Freq:\t" << centerFreq + bandwidth / 2;
qDebug() << "freq points:\t" << freqpoints;
qDebug() << "delta freq:\t" << freqlist.get()[1] - freqlist.get()[0];
qDebug() << u8"频率结束";
/** 4. 初始化回波 **************************************************************************************************/
qDebug() << u8"4.初始化回波";
std::shared_ptr<cuComplex> phdata(createEchoPhase_mallocHost(gpspoints, freqpoints), FreeCUDAHost);
qDebug() << "Azimuth Points:\t" << gpspoints;
qDebug() << "Range Points:\t" << freqpoints;
qDebug() << u8"初始化回波结束";
/** 5. 初始化图像 **************************************************************************************************/
qDebug() << u8"5.初始化图像";
h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * dem_rowCount * dem_ColCount);
qDebug() << "Azimuth Points:\t" << gpspoints;
qDebug() << "Range Points:\t" << freqpoints;
qDebug() << u8"初始化图像结束";
/** 6. 模型参数初始化 **************************************************************************************************/
qDebug() << u8"6.模型参数初始化";
h_data.nx = dem_ColCount;
h_data.ny = dem_rowCount;
h_data.Np = gpspoints;
h_data.Freq = freqlist.get();
h_data.minF = minF.get();
h_data.Nfft = freqpoints;
h_data.K = h_data.Nfft;
h_data.phdata = phdata.get();
h_data.R0 = 900000;
qDebug() << u8"模型参数结束";
/** 7. 目标 **************************************************************************************************/
double Tx = -2029086.618142, Ty = 4139594.934504, Tz = 4392846.782027;
double Tslx = -2029086.618142, Tsly = 4139594.934504, Tslz = 4392846.782027;
double p1 = 33.3, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
/** 7. 构建回波 **************************************************************************************************/
GPUDATA d_data;
initGPUData(h_data, d_data);
RFPCProcess(Tx, Ty, Tz,
Tslx, Tsly, Tslz, // 目标的坡面向量
p1, p2, p3, p4, p5, p6,
d_data);
/** 8. 展示回波 **************************************************************************************************/
{
DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
cuComplex* h_ifftphdata = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
cuComplex* d_ifftphdata = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
CUDAIFFT(d_data.phdata, d_ifftphdata, d_data.Np, d_data.Nfft, d_data.Nfft);
FFTShift1D(d_ifftphdata, d_data.Np, d_data.Nfft);
DeviceToHost(h_ifftphdata, d_ifftphdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
{
for (long i = 0; i < d_data.Np; i++) {
for (long j = 0; j < d_data.Nfft; j++) {
ifftdata.get()[i * d_data.Nfft + j] =
std::complex<double>(
h_ifftphdata[i * d_data.Nfft + j].x,
h_ifftphdata[i * d_data.Nfft + j].y);
}
}
}
testOutComplexDoubleArr(QString("echo_ifft.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
FreeCUDADevice(d_ifftphdata);
FreeCUDAHost(h_ifftphdata);
}
/** 9. 成像 **************************************************************************************************/
// 计算maxWr需要先计算deltaF
double deltaF = h_data.deltaF; // 从输入参数获取
double maxWr = 299792458.0f / (2.0f * deltaF);
qDebug() << "maxWr :\t" << maxWr;
double* r_vec_host = (double*)mallocCUDAHost(sizeof(double) * h_data.Nfft);// new double[data.Nfft];
const double step = maxWr / h_data.Nfft;
const double start = -1 * h_data.Nfft / 2.0f * step;
printf("nfft=%d\n", h_data.Nfft);
for (int i = 0; i < h_data.Nfft; ++i) {
r_vec_host[i] = start + i * step;
}
h_data.r_vec = r_vec_host;
d_data.r_vec = h_data.r_vec;
qDebug() << "r_vec [0]:\t" << h_data.r_vec[0];
qDebug() << "r_vec step:\t" << step;
bpBasic0CUDA(d_data, 0);
DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * d_data.nx * d_data.ny);
{
DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
{
for (long i = 0; i < d_data.Np; i++) {
for (long j = 0; j < d_data.Nfft; j++) {
ifftdata.get()[i * d_data.Nfft + j] =
std::complex<double>(
h_data.phdata[i * d_data.Nfft + j].x,
h_data.phdata[i * d_data.Nfft + j].y);
}
}
}
testOutComplexDoubleArr(QString("echo_ifft_BPBasic.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
std::shared_ptr<std::complex<double>> im_finals(new std::complex<double>[d_data.nx * d_data.ny], delArrPtr);
{
for (long i = 0; i < d_data.ny; i++) {
for (long j = 0; j < d_data.nx; j++) {
im_finals.get()[i * d_data.nx + j] = std::complex<double>(
h_data.im_final[i * d_data.nx + j].x,
h_data.im_final[i * d_data.nx + j].y);
}
}
}
testOutComplexDoubleArr(QString("im_finals.bin"), im_finals.get(), d_data.ny, d_data.nx);
}
}
void testBpImage() {
GPUDATA h_data{};
/** 0. 轨道 **************************************************************************************************/
QString echoPath = "D:\\Programme\\vs2022\\RasterMergeTest\\LAMPCAE_SCANE-all-scane\\GF3_Simulation.xml";
std::shared_ptr<EchoL0Dataset> echoL0ds(new EchoL0Dataset);
echoL0ds->Open(echoPath);
qDebug() << u8"加载回拨文件:\t" << echoPath;
/** 1. 轨道 **************************************************************************************************/
qDebug() << u8"1.轨道文件读取中。。。";
QString inGPSPath = echoL0ds->getGPSPointFilePath();
long gpspoints = gdalImage(inGPSPath).height;
std::shared_ptr<SatelliteAntPos> antpos = SatelliteAntPosOperator::readAntPosFile(inGPSPath, gpspoints);
h_data.AntX = (double*)mallocCUDAHost(sizeof(double) * gpspoints);
h_data.AntY = (double*)mallocCUDAHost(sizeof(double) * gpspoints);
h_data.AntZ = (double*)mallocCUDAHost(sizeof(double) * gpspoints);
for (long i = 0; i < gpspoints; i = i + 1) {
h_data.AntX[i] = antpos.get()[i].Px;
h_data.AntY[i] = antpos.get()[i].Py;
h_data.AntZ[i] = antpos.get()[i].Pz;
}
//gpspoints = gpspoints / 2;
qDebug() << "GPS points Count:\t" << gpspoints;
qDebug() << "PRF 0:\t " << QString("%1,%2,%3").arg(h_data.AntX[0]).arg(h_data.AntY[0]).arg(h_data.AntZ[0]);
qDebug() << "PRF " << gpspoints - 1 << ":\t " << QString("%1,%2,%3")
.arg(h_data.AntX[gpspoints - 1])
.arg(h_data.AntY[gpspoints - 1])
.arg(h_data.AntZ[gpspoints - 1]);
/** 2. 成像网格 **************************************************************************************************/
qDebug() << u8"轨道文件读取结束\n2.成像网格读取。。。";
QString demxyzPath = u8"D:\\Programme\\vs2022\\RasterMergeTest\\simulationData\\demdataset\\demxyz.bin";
gdalImage demgridimg(demxyzPath);
long dem_rowCount = demgridimg.height;
long dem_ColCount = demgridimg.width;
std::shared_ptr<double> demX = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> demY = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> demZ = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
h_data.x_mat = (double*)mallocCUDAHost(sizeof(double) * dem_rowCount * dem_ColCount);
h_data.y_mat = (double*)mallocCUDAHost(sizeof(double) * dem_rowCount * dem_ColCount);
h_data.z_mat = (double*)mallocCUDAHost(sizeof(double) * dem_rowCount * dem_ColCount);
for (long i = 0; i < dem_rowCount; i++) {
for (long j = 0; j < dem_ColCount; j++) {
h_data.x_mat[i * dem_ColCount + j] = demX.get()[i * dem_ColCount + j];
h_data.y_mat[i * dem_ColCount + j] = demY.get()[i * dem_ColCount + j];
h_data.z_mat[i * dem_ColCount + j] = demZ.get()[i * dem_ColCount + j];
}
}
qDebug() << "dem row Count:\t" << dem_rowCount;
qDebug() << "dem col Count:\t" << dem_ColCount;
qDebug() << u8"成像网格读取结束";
/** 3. 频率网格 **************************************************************************************************/
qDebug() << u8"3.处理频率";
double centerFreq = echoL0ds->getCenterFreq();
double bandwidth = echoL0ds->getBandwidth();
size_t freqpoints = echoL0ds->getPlusePoints();
std::shared_ptr<double> freqlist(getFreqPoints_mallocHost(centerFreq - bandwidth / 2, centerFreq + bandwidth / 2, freqpoints), FreeCUDAHost);
std::shared_ptr<double> minF(new double[gpspoints], delArrPtr);
for (long i = 0; i < gpspoints; i++) {
minF.get()[i] = freqlist.get()[0];
}
h_data.deltaF = bandwidth / (freqpoints - 1);
qDebug() << "start Freq:\t" << centerFreq - bandwidth / 2;
qDebug() << "end Freq:\t" << centerFreq + bandwidth / 2;
qDebug() << "freq points:\t" << freqpoints;
qDebug() << "delta freq:\t" << freqlist.get()[1] - freqlist.get()[0];
qDebug() << u8"频率结束";
/** 4. 初始化回波 **************************************************************************************************/
qDebug() << u8"4.初始化回波";
std::shared_ptr<std::complex<double>> echoData = echoL0ds->getEchoArr();
size_t echosize = sizeof(cuComplex) * gpspoints * freqpoints;
qDebug() << "echo data size (byte) :\t" << echosize;
h_data.phdata = (cuComplex*)mallocCUDAHost(echosize);
shared_complexPtrToHostCuComplex(echoData.get(), h_data.phdata, gpspoints * freqpoints);
qDebug() << "Azimuth Points:\t" << gpspoints;
qDebug() << "Range Points:\t" << freqpoints;
qDebug() << u8"初始化回波结束";
/** 5. 初始化图像 **************************************************************************************************/
qDebug() << u8"5.初始化图像";
h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * dem_rowCount * dem_ColCount);
qDebug() << "Azimuth Points:\t" << gpspoints;
qDebug() << "Range Points:\t" << freqpoints;
qDebug() << u8"初始化图像结束";
/** 6. 模型参数初始化 **************************************************************************************************/
qDebug() << u8"6.模型参数初始化";
h_data.nx = dem_ColCount;
h_data.ny = dem_rowCount;
h_data.Np = gpspoints;
h_data.Freq = freqlist.get();
h_data.minF = minF.get();
h_data.Nfft = freqpoints;
h_data.K = h_data.Nfft;
h_data.R0 = echoL0ds->getRefPhaseRange();
qDebug() << u8"模型参数结束";
/** 7. 目标 **************************************************************************************************/
double Tx = -2029086.618142, Ty = 4139594.934504, Tz = 4392846.782027;
double Tslx = -2029086.618142, Tsly = 4139594.934504, Tslz = 4392846.782027;
double p1 = 33.3, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
/** 7. 构建回波 **************************************************************************************************/
GPUDATA d_data;
initGPUData(h_data, d_data);
/** 8. 展示回波 **************************************************************************************************/
{
DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
cuComplex* h_ifftphdata = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
cuComplex* d_ifftphdata = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
CUDAIFFT(d_data.phdata, d_ifftphdata, d_data.Np, d_data.Nfft, d_data.Nfft);
FFTShift1D(d_ifftphdata, d_data.Np, d_data.Nfft);
DeviceToHost(h_ifftphdata, d_ifftphdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
{
for (long i = 0; i < d_data.Np; i++) {
for (long j = 0; j < d_data.Nfft; j++) {
ifftdata.get()[i * d_data.Nfft + j] =
std::complex<double>(
h_ifftphdata[i * d_data.Nfft + j].x,
h_ifftphdata[i * d_data.Nfft + j].y);
}
}
}
testOutComplexDoubleArr(QString("echo_ifft.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
FreeCUDADevice(d_ifftphdata);
FreeCUDAHost(h_ifftphdata);
}
/** 9. 成像 **************************************************************************************************/
// 计算maxWr需要先计算deltaF
double deltaF = h_data.deltaF; // 从输入参数获取
double maxWr = 299792458.0f / (2.0f * deltaF);
qDebug() << "maxWr :\t" << maxWr;
double* r_vec_host = (double*)mallocCUDAHost(sizeof(double) * h_data.Nfft);// new double[data.Nfft];
const double step = maxWr / h_data.Nfft;
const double start = -1 * h_data.Nfft / 2.0f * step;
printf("nfft=%d\n", h_data.Nfft);
for (int i = 0; i < h_data.Nfft; ++i) {
r_vec_host[i] = start + i * step;
}
h_data.r_vec = r_vec_host;
d_data.r_vec = h_data.r_vec;
qDebug() << "r_vec [0]:\t" << h_data.r_vec[0];
qDebug() << "r_vec step:\t" << step;
bpBasic0CUDA(d_data, 0);
DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * d_data.nx * d_data.ny);
{
DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
{
for (long i = 0; i < d_data.Np; i++) {
for (long j = 0; j < d_data.Nfft; j++) {
ifftdata.get()[i * d_data.Nfft + j] =
std::complex<double>(
h_data.phdata[i * d_data.Nfft + j].x,
h_data.phdata[i * d_data.Nfft + j].y);
}
}
}
testOutComplexDoubleArr(QString("echo_ifft_BPBasic.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
std::shared_ptr<std::complex<double>> im_finals(new std::complex<double>[d_data.nx * d_data.ny], delArrPtr);
{
for (long i = 0; i < d_data.ny; i++) {
for (long j = 0; j < d_data.nx; j++) {
im_finals.get()[i * d_data.nx + j] = std::complex<double>(
h_data.im_final[i * d_data.nx + j].x,
h_data.im_final[i * d_data.nx + j].y);
}
}
}
testOutComplexDoubleArr(QString("im_finals.bin"), im_finals.get(), d_data.ny, d_data.nx);
}
}
int main(int argc, char* argv[]) {
QApplication a(argc, argv);
testBpImage();
return 0;
}