Merge pull request '合并修复BP成像bug' (#6) from RTPC-BPImage_dev into RTPC-dev

Reviewed-on: http://123.153.4.249:22779/LAMPSARToolSoftware/RasterProcessTool/pulls/6
pull/7/head
chenzenghui 2025-03-06 16:45:12 +08:00
commit a83c39999c
33 changed files with 2166 additions and 639 deletions

1
.gitignore vendored
View File

@ -377,3 +377,4 @@ FodyWeavers.xsd
/script/landArr.bin.enp
/script/landArr.hdr
/script/.ipynb_checkpoints
/LAMPDataProcessEXE

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

@ -130,23 +130,68 @@ void initializeMatrixWithSSE2(Eigen::MatrixXd& mat, const double* data, long row
void initializeMatrixWithSSE2(Eigen::MatrixXf& mat, const float* data, long rowcount, long colcount);
/** 模板函数类 ***********************************************************************************************************/
template<typename T>
inline void BASECONSTVARIABLEAPI memsetInitArray(std::shared_ptr<T> ptr, long arrcount,T ti) {
inline void memsetInitArray(std::shared_ptr<T> ptr, long arrcount,T ti) {
for (long i = 0; i < arrcount; i++) {
ptr.get()[i] = ti;
}
}
template<typename T>
inline void BASECONSTVARIABLEAPI memcpyArray(std::shared_ptr<T> srct, std::shared_ptr<T> dest, long arrcount) {
inline void memcpyArray(std::shared_ptr<T> srct, std::shared_ptr<T> dest, long arrcount) {
for (long i = 0; i < arrcount; i++) {
dest.get()[i] = srct.get()[i];
}
}
template<typename T>
inline void minValueInArr(T* ptr, long arrcount, T& minvalue) {
if (arrcount == 0)return;
minvalue = ptr[0];
for (long i = 0; i < arrcount; i++) {
if (minvalue > ptr[i]) {
minvalue = ptr[i];
}
}
}
template<typename T>
inline void maxValueInArr(T* ptr, long arrcount, T& maxvalue) {
if (arrcount == 0)return;
maxvalue = ptr[0];
for (long i = 0; i < arrcount; i++) {
if (maxvalue < ptr[i]) {
maxvalue = ptr[i];
}
}
}
/** 常用SAR工具 ***********************************************************************************************************/
template<typename T>
inline T complexAbs(std::complex<T> ccdata) {
return T(sqrt(pow(ccdata.real(), 2) + pow(ccdata.imag(), 2)));
}
template<typename T>
inline void complex2dB(std::complex<T>* ccdata, T* outdata, long long count) {
for (long long i = 0; i < count; i++) {
outdata[i] = 20 * log10(complexAbs(ccdata[i]));
}
}

View File

@ -251,6 +251,16 @@ QString EchoL0Dataset::getEchoDataFilename()
return GPSPointFilePath;
}
QString EchoL0Dataset::getGPSPointFilePath()
{
return this->GPSPointFilePath;
}
QString EchoL0Dataset::getEchoDataFilePath()
{
return this->echoDataFilePath;
}
void EchoL0Dataset::initEchoArr(std::complex<double> init0)
{
long blockline = Memory1MB / 8 / 2 / this->PlusePoints * 8000;
@ -531,7 +541,7 @@ std::shared_ptr<double> EchoL0Dataset::getAntPos()
return temp;
}
std::shared_ptr<std::complex<double>> EchoL0Dataset::getEchoArr(long startPRF, long PRFLen)
std::shared_ptr<std::complex<double>> EchoL0Dataset::getEchoArr(long startPRF, long& PRFLen)
{
if (!(startPRF < this->PluseCount)) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE))<<startPRF<<" "<<this->PluseCount;
@ -558,6 +568,9 @@ std::shared_ptr<std::complex<double>> EchoL0Dataset::getEchoArr(long startPRF, l
long width = rasterDataset->GetRasterXSize();
long height = rasterDataset->GetRasterYSize();
long band_num = rasterDataset->GetRasterCount();
PRFLen = (PRFLen + startPRF) < height ? PRFLen : height - startPRF;
std::shared_ptr<std::complex<double>> temp = nullptr;
if (height != this->PluseCount || width != this->PlusePoints) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
@ -629,8 +642,6 @@ ErrorCode EchoL0Dataset::saveEchoArr(std::shared_ptr<std::complex<double>> echoP
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->echoDataFilePath.toUtf8().constData(), GDALAccess::GA_Update));
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* poBand = rasterDataset->GetRasterBand(1);
@ -667,4 +678,69 @@ ErrorCode EchoL0Dataset::saveEchoArr(std::shared_ptr<std::complex<double>> echoP
omp_destroy_lock(&lock); //
return ErrorCode::SUCCESS;
}
std::shared_ptr<SatelliteAntPos> SatelliteAntPosOperator::readAntPosFile(QString filepath, long& count)
{
gdalImage antimg(filepath);
long rowcount = count;
long colcount = 19;
std::shared_ptr<double> antlist = readDataArr<double>(antimg, 0, 0, rowcount, colcount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<SatelliteAntPos> antpos(new SatelliteAntPos[rowcount], delArrPtr);
for (long i = 0; i < rowcount; i++) {
antpos.get()[i].time = antlist.get()[i * 19 + 0];
antpos.get()[i].Px = antlist.get()[i * 19 + 1];
antpos.get()[i].Py = antlist.get()[i * 19 + 2];
antpos.get()[i].Pz = antlist.get()[i * 19 + 3];
antpos.get()[i].Vx = antlist.get()[i * 19 + 4];
antpos.get()[i].Vy = antlist.get()[i * 19 + 5];
antpos.get()[i].Vz = antlist.get()[i * 19 + 6]; //7
antpos.get()[i].AntDirectX = antlist.get()[i * 19 + 7];
antpos.get()[i].AntDirectY = antlist.get()[i * 19 + 8];
antpos.get()[i].AntDirectZ = antlist.get()[i * 19 + 9];
antpos.get()[i].AVx = antlist.get()[i * 19 + 10];
antpos.get()[i].AVy = antlist.get()[i * 19 + 11];
antpos.get()[i].AVz = antlist.get()[i * 19 + 12];
antpos.get()[i].ZeroAntDiectX = antlist.get()[i * 19 + 13];
antpos.get()[i].ZeroAntDiectY = antlist.get()[i * 19 + 14];
antpos.get()[i].ZeroAntDiectZ = antlist.get()[i * 19 + 15];
antpos.get()[i].lon = antlist.get()[i * 19 + 16];
antpos.get()[i].lat = antlist.get()[i * 19 + 17];
antpos.get()[i].ati = antlist.get()[i * 19 + 18]; // 19
}
return antpos;
}
void SatelliteAntPosOperator::writeAntPosFile(QString filepath, std::shared_ptr<SatelliteAntPos> data, const long count)
{
gdalImage antimg=CreategdalImageDouble(filepath,count,19,1,true,true);
long rowcount = count;
long colcount = 19;
std::shared_ptr<double> antpos(new double[rowcount*19], delArrPtr);
for (long i = 0; i < colcount; i++) {
antpos.get()[i * 19 + 1] = data.get()[i].time;
antpos.get()[i * 19 + 2] = data.get()[i].Px;
antpos.get()[i * 19 + 3] = data.get()[i].Py;
antpos.get()[i * 19 + 4] = data.get()[i].Pz;
antpos.get()[i * 19 + 5] = data.get()[i].Vx;
antpos.get()[i * 19 + 6] = data.get()[i].Vy;
antpos.get()[i * 19 + 7] = data.get()[i].Vz;
antpos.get()[i * 19 + 8] = data.get()[i].AntDirectX;
antpos.get()[i * 19 + 9] = data.get()[i].AntDirectY;
antpos.get()[i * 19 + 10] = data.get()[i].AntDirectZ;
antpos.get()[i * 19 + 11] = data.get()[i].AVx;
antpos.get()[i * 19 + 12] = data.get()[i].AVy;
antpos.get()[i * 19 + 13] = data.get()[i].AVz;
antpos.get()[i * 19 + 14] = data.get()[i].ZeroAntDiectX;
antpos.get()[i * 19 + 15] = data.get()[i].ZeroAntDiectY;
antpos.get()[i * 19 + 16] = data.get()[i].ZeroAntDiectZ;
antpos.get()[i * 19 + 17] = data.get()[i].lon;
antpos.get()[i * 19 + 18] = data.get()[i].lat;
antpos.get()[i * 19 + 19] = data.get()[i].ati;
}
antimg.saveImage(antpos, 0,0,rowcount, colcount, 1);
return ;
}

View File

@ -92,6 +92,31 @@ struct PluseAntPos {
std::shared_ptr<PluseAntPos> BASECONSTVARIABLEAPI CreatePluseAntPosArr(long pluseCount);
class BASECONSTVARIABLEAPI SatelliteAntPosOperator {
public:
static std::shared_ptr<SatelliteAntPos> readAntPosFile(QString filepath,long& count);
static void writeAntPosFile(QString filepath, std::shared_ptr< SatelliteAntPos> data,const long count);
};
// 定义L0级数据
class BASECONSTVARIABLEAPI EchoL0Dataset {
@ -107,7 +132,8 @@ public:
QString getxmlName();
QString getGPSPointFilename();
QString getEchoDataFilename();
QString getGPSPointFilePath();
QString getEchoDataFilePath();
void initEchoArr(std::complex<double> init0);
@ -184,7 +210,7 @@ public: //
public:
// 读取文件
std::shared_ptr<double> getAntPos();
std::shared_ptr<std::complex<double>> getEchoArr(long startPRF, long PRFLen);
std::shared_ptr<std::complex<double>> getEchoArr(long startPRF, long& PRFLen);
std::shared_ptr<std::complex<double>> getEchoArr();
//保存文件
ErrorCode saveAntPos(std::shared_ptr<double> ptr); // 注意这个方法很危险,请写入前检查数据是否正确

View File

@ -477,7 +477,7 @@ gdalImage::gdalImage(const QString& raster_path)
rasterDataset = NULL; // 指矫匡拷
this->InitInv_gt();
delete[] gt;
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
}
@ -624,7 +624,7 @@ Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count,
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return datamatrix;
}
@ -753,7 +753,7 @@ Eigen::MatrixXf gdalImage::getDataf(int start_row, int start_col, int rows_count
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return datamatrix;
}
@ -882,7 +882,7 @@ Eigen::MatrixXi gdalImage::getDatai(int start_row, int start_col, int rows_count
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return datamatrix;
}
@ -1039,7 +1039,7 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
}
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
delete[] databuffer;
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
@ -1120,7 +1120,7 @@ void gdalImage::saveImage(Eigen::MatrixXf data, int start_row = 0, int start_col
}
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
delete[] databuffer;
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
@ -1179,7 +1179,7 @@ void gdalImage::saveImage(Eigen::MatrixXi data, int start_row, int start_col, in
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
delete[] databuffer;
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
@ -1240,7 +1240,7 @@ void gdalImage::saveImage(std::shared_ptr<double> data, int start_row, int start
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
delete[] databuffer;
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
@ -1303,7 +1303,7 @@ void gdalImage::saveImage(std::shared_ptr<float> data, int start_row, int start_
GDALClose((GDALDatasetH)poDstDS);
//delete poDstDS;
//poDstDS = nullptr;
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
delete[] databuffer;
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
@ -1365,7 +1365,7 @@ void gdalImage::saveImage(std::shared_ptr<int> data, int start_row, int start_co
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
delete[] databuffer;
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
@ -1642,6 +1642,39 @@ RasterExtend gdalImage::getExtend()
return extend;
}
gdalImage BASECONSTVARIABLEAPI CreategdalImageDouble(QString& img_path, int height, int width, int band_num, bool overwrite, bool isEnvi)
{
if (exists_test(img_path.toUtf8().constData())) {
if (overwrite) {
gdalImage result_img(img_path);
return result_img;
}
else {
throw "file has exist!!!";
exit(1);
}
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
GDALDriver* poDriver = nullptr;
if (isEnvi) {
poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
}
else {
poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
}
GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,GDT_Float64, NULL); // 锟斤拷锟斤拷锟斤拷
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImage result_img(img_path);
return result_img;
}
gdalImage CreategdalImageDouble(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite, bool isEnvi)
{
if (exists_test(img_path.toUtf8().constData())) {
@ -1684,7 +1717,7 @@ gdalImage CreategdalImageDouble(const QString& img_path, int height, int width,
}
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImage result_img(img_path);
return result_img;
}
@ -1731,7 +1764,7 @@ gdalImage CreategdalImage(const QString& img_path, int height, int width, int ba
}
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImage result_img(img_path);
return result_img;
}
@ -1777,7 +1810,7 @@ gdalImage CreategdalImage(const QString& img_path, int height, int width, int ba
}
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImage result_img(img_path);
return result_img;
}
@ -1817,7 +1850,7 @@ gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int
//}
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImageComplex result_img(img_path);
return result_img;
}
@ -1962,7 +1995,7 @@ int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int
// GDALDestroyWarpOptions(psWo);
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return 0;
}
@ -2551,7 +2584,7 @@ gdalImageComplex::gdalImageComplex(const QString& raster_path)
rasterDataset = NULL; // 指矫匡拷
this->InitInv_gt();
delete[] gt;
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
}
@ -2615,7 +2648,7 @@ void gdalImageComplex::saveImage(Eigen::MatrixXcd data, int start_row, int start
GDALFlushCache(poDstDS);
delete databuffer;
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
@ -2668,12 +2701,65 @@ void gdalImageComplex::saveImage(std::shared_ptr<std::complex<double>> data, lon
GDALFlushCache(poDstDS);
delete databuffer;
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
void gdalImageComplex::saveImage(std::complex<double>* data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
if (start_row + rowCount > this->height || start_col + colCount > this->width) {
QString tip = u8"file path: " + this->img_path;
qDebug() << tip;
throw std::exception(tip.toUtf8().constData());
return;
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
GDALDataset* poDstDS = nullptr;
if (exists_test(this->img_path)) {
poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
}
else {
poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
this->band_num, GDT_CFloat64, NULL); // 斤拷锟斤拷
poDstDS->SetProjection(this->projection.toUtf8().constData());
double gt_ptr[6];
for (int i = 0; i < this->gt.rows(); i++) {
for (int j = 0; j < this->gt.cols(); j++) {
gt_ptr[i * 3 + j] = this->gt(i, j);
}
}
poDstDS->SetGeoTransform(gt_ptr);
//delete[] gt_ptr;
}
double* databuffer = new double[rowCount * colCount * 2];
for (long i = 0; i < rowCount; i++) {
for (long j = 0; j < colCount; j++) {
databuffer[i * colCount * 2 + j * 2] = data[i * colCount + j].real();
databuffer[i * colCount * 2 + j * 2 + 1] = data[i * colCount + j].imag();
}
}
// poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
// datarows, GDT_Float32,band_ids, num,0,0,0);
poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, colCount, rowCount,
databuffer, colCount, rowCount, GDT_CFloat64, 0, 0);
GDALFlushCache(poDstDS);
delete databuffer;
GDALClose((GDALDatasetH)poDstDS);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col, int rows_count,
int cols_count, int band_ids)
{
@ -3525,6 +3611,29 @@ void testOutClsArr(QString filename, long* amp, long rowcount, long colcount) {
}
void BASECONSTVARIABLEAPI testOutComplexDoubleArr(QString filename, std::complex<double>* data, long rowcount, long colcount)
{
gdalImageComplex compleximg= CreateEchoComplex(filename, rowcount, colcount, 1);
compleximg.saveImage( data, 0, 0, rowcount, colcount, 1);
return void BASECONSTVARIABLEAPI();
}
void BASECONSTVARIABLEAPI testOutDataArr(QString filename, double* data, long rowcount, long colcount)
{
return testOutAmpArr(filename, data, rowcount, colcount);
}
void BASECONSTVARIABLEAPI testOutDataArr(QString filename, float* data, long rowcount, long colcount)
{
return testOutAmpArr(filename, data, rowcount, colcount);
}
void BASECONSTVARIABLEAPI testOutDataArr(QString filename, long* data, long rowcount, long colcount)
{
return testOutClsArr(filename,data,rowcount,colcount);
}
void testOutAntPatternTrans(QString antpatternfilename, double* antPatternArr,
double starttheta, double deltetheta,
double startphi, double deltaphi,
@ -3905,7 +4014,7 @@ void ResampleByReferenceRasterB(QString pszSrcFile, QString RefrasterBPath, QStr
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return ;
}

View File

@ -229,6 +229,8 @@ public: // 方法
void setData(Eigen::MatrixXcd);
void saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids);
void saveImage(std::shared_ptr<std::complex<double>> data, long start_row, long start_col, long rowCount, long colCount, int band_ids);
void saveImage(std::complex<double>* data, long start_row, long start_col, long rowcount, long colcount, int banids);
Eigen::MatrixXcd getDataComplex(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
std::shared_ptr<std::complex<double>> getDataComplexSharePtr(int start_row, int start_col, int rows_count, int cols_count, int band_ids);
@ -239,6 +241,7 @@ public:
};
// 创建影像
gdalImage BASECONSTVARIABLEAPI CreategdalImageDouble(QString& img_path, int height, int width, int band_num, bool overwrite = false, bool isEnvi = false);
gdalImage BASECONSTVARIABLEAPI CreategdalImageDouble(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false, bool isEnvi = false);
gdalImage BASECONSTVARIABLEAPI CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection,bool need_gt = true, bool overwrite = false, bool isEnvi = false, GDALDataType datetype = GDT_Float32);
@ -301,8 +304,12 @@ bool BASECONSTVARIABLEAPI saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString
void BASECONSTVARIABLEAPI testOutAntPatternTrans(QString antpatternfilename, double* antPatternArr, double starttheta, double deltetheta, double startphi, double deltaphi, long thetanum, long phinum);
void BASECONSTVARIABLEAPI testOutAmpArr(QString filename, float* amp, long rowcount, long colcount);
void BASECONSTVARIABLEAPI testOutAmpArr(QString filename, double* amp, long rowcount, long colcount);
void BASECONSTVARIABLEAPI testOutClsArr(QString filename, long* amp, long rowcount, long colcount);
void BASECONSTVARIABLEAPI testOutClsArr(QString filename, long* amp, long rowcount, long colcount);
void BASECONSTVARIABLEAPI testOutComplexDoubleArr(QString filename, std::complex<double>* data, long rowcount, long colcount);
void BASECONSTVARIABLEAPI testOutDataArr(QString filename, double* data, long rowcount, long colcount);
void BASECONSTVARIABLEAPI testOutDataArr(QString filename, float* data, long rowcount, long colcount);
void BASECONSTVARIABLEAPI testOutDataArr(QString filename, long* data, long rowcount, long colcount);
void BASECONSTVARIABLEAPI CreateSARIntensityByLookTable(QString IntensityRasterPath, QString LookTableRasterPath, QString SARIntensityPath, long min_rid, long max_rid, long min_cid, long max_cid, std::function<void(long, long)> processBarShow = {});
@ -504,7 +511,7 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return result;
}
@ -577,7 +584,7 @@ inline std::shared_ptr<T> readDataArrComplex(gdalImageComplex& imgds, long start
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return result;
}

View File

@ -6,8 +6,31 @@ BASECONSTVARIABLEAPI void PrintMsgToQDebug(char* msg)
return ;
}
BASECONSTVARIABLEAPI void PrintfToQDebug(const char* msg)
{
qDebug() << QString(msg);
return;
}
BASECONSTVARIABLEAPI void PrintTipMsgToQDebug(const char* tip, const char* msg)
{
qDebug() <<QString(tip)<<"\t:\t" << QString(msg);
return;
}
// 自定义的 printf 风格函数
BASECONSTVARIABLEAPI void printfinfo(const char* format, ...) {
char buffer[256]; // 假设最大输出长度为 256 字节
va_list args;
// 使用 va_start 获取可变参数列表
va_start(args, format);
// 使用 vsnprintf 将格式化字符串写入 buffer
vsnprintf(buffer, sizeof(buffer), format, args);
// 结束可变参数列表的使用
va_end(args);
// 将格式化后的字符串转发给 PrintfToQDebug
PrintfToQDebug(buffer);
}

View File

@ -2,9 +2,11 @@
#ifndef PRINTMSGTOQDEBUG_H_
#define PRINTMSGTOQDEBUG_H_
#include "BaseConstVariable.h"
#include <format>
extern "C" BASECONSTVARIABLEAPI void PrintMsgToQDebug(char* msg);
extern "C" BASECONSTVARIABLEAPI void PrintfToQDebug(const char* msg);
extern "C" BASECONSTVARIABLEAPI void PrintTipMsgToQDebug(const char* tip, const char* msg);
extern "C" BASECONSTVARIABLEAPI void printfinfo(const char* format, ...);
#endif // !PRINTMSGTOQDEBUG_H_

View File

@ -0,0 +1,104 @@
<?xml version="1.0" encoding="utf-8"?>
<Project DefaultTargets="Build" ToolsVersion="15.0" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
<ItemGroup Label="ProjectConfigurations">
<ProjectConfiguration Include="Debug|x64">
<Configuration>Debug</Configuration>
<Platform>x64</Platform>
</ProjectConfiguration>
<ProjectConfiguration Include="Release|x64">
<Configuration>Release</Configuration>
<Platform>x64</Platform>
</ProjectConfiguration>
</ItemGroup>
<PropertyGroup Label="Globals">
<ProjectGuid>9233788c-bd43-41aa-b157-d77e92616d00</ProjectGuid>
<RootNamespace>GPUBPSimulation</RootNamespace>
</PropertyGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
<ConfigurationType>Application</ConfigurationType>
<UseDebugLibraries>true</UseDebugLibraries>
<CharacterSet>MultiByte</CharacterSet>
<PlatformToolset>v143</PlatformToolset>
</PropertyGroup>
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
<ConfigurationType>Application</ConfigurationType>
<UseDebugLibraries>false</UseDebugLibraries>
<WholeProgramOptimization>true</WholeProgramOptimization>
<CharacterSet>MultiByte</CharacterSet>
<PlatformToolset>v143</PlatformToolset>
</PropertyGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
<ImportGroup Label="ExtensionSettings">
<Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 12.6.props" />
</ImportGroup>
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
</ImportGroup>
<ImportGroup Label="PropertySheets" Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
<Import Project="$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props" Condition="exists('$(UserRootDir)\Microsoft.Cpp.$(Platform).user.props')" Label="LocalAppDataPlatform" />
</ImportGroup>
<PropertyGroup Label="UserMacros" />
<PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
<LinkIncremental>true</LinkIncremental>
</PropertyGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">
<ClCompile>
<WarningLevel>Level3</WarningLevel>
<Optimization>Disabled</Optimization>
<PreprocessorDefinitions>WIN32;WIN64;_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
</ClCompile>
<Link>
<GenerateDebugInformation>true</GenerateDebugInformation>
<SubSystem>Console</SubSystem>
<AdditionalDependencies>cudart_static.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)</AdditionalDependencies>
</Link>
<CudaCompile>
<TargetMachinePlatform>64</TargetMachinePlatform>
</CudaCompile>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
<ClCompile>
<WarningLevel>Level3</WarningLevel>
<Optimization>MaxSpeed</Optimization>
<FunctionLevelLinking>true</FunctionLevelLinking>
<IntrinsicFunctions>true</IntrinsicFunctions>
<PreprocessorDefinitions>WIN32;WIN64;NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
</ClCompile>
<Link>
<GenerateDebugInformation>true</GenerateDebugInformation>
<EnableCOMDATFolding>true</EnableCOMDATFolding>
<OptimizeReferences>true</OptimizeReferences>
<SubSystem>Console</SubSystem>
<AdditionalDependencies>cudart_static.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)</AdditionalDependencies>
</Link>
<CudaCompile>
<TargetMachinePlatform>64</TargetMachinePlatform>
</CudaCompile>
</ItemDefinitionGroup>
<ItemGroup>
<CudaCompile Include="kernel.cu" />
</ItemGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
<ImportGroup Label="ExtensionTargets">
<Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 12.6.targets" />
</ImportGroup>
</Project>

121
GPUBPSimulation/kernel.cu Normal file
View File

@ -0,0 +1,121 @@

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
cudaError_t addWithCuda(int *c, const int *a, const int *b, unsigned int size);
__global__ void addKernel(int *c, const int *a, const int *b)
{
int i = threadIdx.x;
c[i] = a[i] + b[i];
}
int main()
{
const int arraySize = 5;
const int a[arraySize] = { 1, 2, 3, 4, 5 };
const int b[arraySize] = { 10, 20, 30, 40, 50 };
int c[arraySize] = { 0 };
// Add vectors in parallel.
cudaError_t cudaStatus = addWithCuda(c, a, b, arraySize);
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "addWithCuda failed!");
return 1;
}
printf("{1,2,3,4,5} + {10,20,30,40,50} = {%d,%d,%d,%d,%d}\n",
c[0], c[1], c[2], c[3], c[4]);
// cudaDeviceReset must be called before exiting in order for profiling and
// tracing tools such as Nsight and Visual Profiler to show complete traces.
cudaStatus = cudaDeviceReset();
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaDeviceReset failed!");
return 1;
}
return 0;
}
// Helper function for using CUDA to add vectors in parallel.
cudaError_t addWithCuda(int *c, const int *a, const int *b, unsigned int size)
{
int *dev_a = 0;
int *dev_b = 0;
int *dev_c = 0;
cudaError_t cudaStatus;
// Choose which GPU to run on, change this on a multi-GPU system.
cudaStatus = cudaSetDevice(0);
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaSetDevice failed! Do you have a CUDA-capable GPU installed?");
goto Error;
}
// Allocate GPU buffers for three vectors (two input, one output) .
cudaStatus = cudaMalloc((void**)&dev_c, size * sizeof(int));
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaMalloc failed!");
goto Error;
}
cudaStatus = cudaMalloc((void**)&dev_a, size * sizeof(int));
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaMalloc failed!");
goto Error;
}
cudaStatus = cudaMalloc((void**)&dev_b, size * sizeof(int));
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaMalloc failed!");
goto Error;
}
// Copy input vectors from host memory to GPU buffers.
cudaStatus = cudaMemcpy(dev_a, a, size * sizeof(int), cudaMemcpyHostToDevice);
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaMemcpy failed!");
goto Error;
}
cudaStatus = cudaMemcpy(dev_b, b, size * sizeof(int), cudaMemcpyHostToDevice);
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaMemcpy failed!");
goto Error;
}
// Launch a kernel on the GPU with one thread for each element.
addKernel<<<1, size>>>(dev_c, dev_a, dev_b);
// Check for any errors launching the kernel
cudaStatus = cudaGetLastError();
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "addKernel launch failed: %s\n", cudaGetErrorString(cudaStatus));
goto Error;
}
// cudaDeviceSynchronize waits for the kernel to finish, and returns
// any errors encountered during the launch.
cudaStatus = cudaDeviceSynchronize();
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching addKernel!\n", cudaStatus);
goto Error;
}
// Copy output vector from GPU buffer to host memory.
cudaStatus = cudaMemcpy(c, dev_c, size * sizeof(int), cudaMemcpyDeviceToHost);
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "cudaMemcpy failed!");
goto Error;
}
Error:
cudaFree(dev_c);
cudaFree(dev_a);
cudaFree(dev_b);
return cudaStatus;
}

View File

@ -252,6 +252,31 @@ extern "C" void FFTShift1D(cuComplex* d_data, int batch_size, int signal_length)
cudaDeviceSynchronize();
}
extern "C" void shared_complexPtrToHostCuComplex(std::complex<double>* src, cuComplex* dst, size_t len)
{
for (long i = 0; i < len; i++) {
dst[i] = make_cuComplex(src[i].real(), src[i].imag());
}
return ;
}
extern "C" void HostCuComplexToshared_complexPtr( cuComplex* src, std::complex<double>* dst, size_t len)
{
double maxvalue = src[0].x;
for (long i = 0; i < len; i++) {
dst[i] = std::complex<double>(src[i].x, src[i].y);
if (maxvalue < src[i].x) {
maxvalue = src[i].x;
}
if (maxvalue < src[i].y) {
maxvalue = src[i].y;
}
}
printf("max value %e\n", maxvalue);
return;
}
extern __global__ void CUDA_D_sin(double* y, double* X, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {
@ -293,14 +318,14 @@ extern "C" void checkCudaError(cudaError_t err, const char* msg) {
}
// 主机参数内存声明
extern "C" void* mallocCUDAHost(long memsize) {
extern "C" void* mallocCUDAHost(size_t memsize) {
void* ptr;
cudaMallocHost(&ptr, memsize);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("mallocCUDAHost CUDA Error: %s\n", cudaGetErrorString(err));
printf("mallocCUDAHost CUDA Error: %s, malloc memory : %d byte\n", cudaGetErrorString(err),memsize);
exit(2);
}
#endif // __CUDADEBUG__
@ -310,25 +335,32 @@ extern "C" void* mallocCUDAHost(long memsize) {
// 主机参数内存释放
extern "C" void FreeCUDAHost(void* ptr) {
if (nullptr == ptr||NULL==ptr) {
return;
}
cudaFreeHost(ptr);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("FreeCUDAHost CUDA Error: %s\n", cudaGetErrorString(err));
printf("FreeCUDAHost CUDA Error: %s,\n", cudaGetErrorString(err));
exit(2);
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
ptr = nullptr;
}
// GPU参数内存声明
extern "C" void* mallocCUDADevice(long memsize) {
extern "C" void* mallocCUDADevice(size_t memsize) {
void* ptr;
cudaMalloc(&ptr, memsize);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("mallocCUDADevice CUDA Error: %s\n", cudaGetErrorString(err));
printf("mallocCUDADevice CUDA Error: %s, malloc memory : %d byte\n", cudaGetErrorString(err), memsize);
exit(2);
}
#endif // __CUDADEBUG__
@ -338,6 +370,9 @@ extern "C" void* mallocCUDADevice(long memsize) {
// GPU参数内存释放
extern "C" void FreeCUDADevice(void* ptr) {
if (nullptr == ptr || NULL == ptr) {
return;
}
cudaFree(ptr);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
@ -347,10 +382,11 @@ extern "C" void FreeCUDADevice(void* ptr) {
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
ptr = nullptr;
}
// GPU 内存数据转移
extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize) {
extern "C" void HostToDevice(void* hostptr, void* deviceptr, size_t memsize) {
cudaMemcpy(deviceptr, hostptr, memsize, cudaMemcpyHostToDevice);
#ifdef __CUDADEBUG__
@ -364,7 +400,7 @@ extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize) {
cudaDeviceSynchronize();
}
extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize) {
extern "C" void DeviceToHost(void* hostptr, void* deviceptr, size_t memsize) {
cudaMemcpy(hostptr, deviceptr, memsize, cudaMemcpyDeviceToHost);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
@ -376,7 +412,7 @@ extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize) {
cudaDeviceSynchronize();
}
void DeviceToDevice(void* s_deviceptr, void* t_deviceptr, long memsize)
void DeviceToDevice(void* s_deviceptr, void* t_deviceptr, size_t memsize)
{
cudaMemcpy(t_deviceptr, s_deviceptr, memsize, cudaMemcpyDeviceToDevice);
#ifdef __CUDADEBUG__
@ -535,6 +571,7 @@ long NextBlockPad(long num, long blocksize)
void PrintLasterError(const char* s)
{
cudaDeviceSynchronize();
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
//printf("%s: %s\n", s, cudaGetErrorString(err));
@ -548,6 +585,7 @@ void PrintLasterError(const char* s)
extern "C" void CUDAIFFT(cuComplex* inArr, cuComplex* outArr, long InRowCount, long InColCount, long outColCount) {
cufftHandle plan;
cufftResult result;
@ -575,7 +613,7 @@ extern "C" void CUDAIFFT(cuComplex* inArr, cuComplex* outArr, long InRowCount,
cuComplex* in_ptr = inArr;
cuComplex* out_ptr = outArr;
result = cufftExecC2C(plan, (cufftComplex*)in_ptr, (cufftComplex*)out_ptr, CUFFT_INVERSE);
if (result != CUFFT_SUCCESS) {
cufftDestroy(plan);
return;
@ -585,7 +623,8 @@ extern "C" void CUDAIFFT(cuComplex* inArr, cuComplex* outArr, long InRowCount,
cudaDeviceSynchronize();
cufftDestroy(plan);
}
}

View File

@ -75,13 +75,13 @@ extern "C" GPUBASELIBAPI void printDeviceInfo(int deviceId);
extern "C" GPUBASELIBAPI void checkCudaError(cudaError_t err, const char* msg);
// GPU 内存函数
extern "C" GPUBASELIBAPI void* mallocCUDAHost(long memsize); // 主机内存声明
extern "C" GPUBASELIBAPI void* mallocCUDAHost(size_t memsize); // 主机内存声明
extern "C" GPUBASELIBAPI void FreeCUDAHost(void* ptr);
extern "C" GPUBASELIBAPI void* mallocCUDADevice(long memsize); // GPU内存声明
extern "C" GPUBASELIBAPI void* mallocCUDADevice(size_t memsize); // GPU内存声明
extern "C" GPUBASELIBAPI void FreeCUDADevice(void* ptr);
extern "C" GPUBASELIBAPI void HostToDevice(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移 设备 -> GPU
extern "C" GPUBASELIBAPI void DeviceToHost(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移 GPU -> 设备
extern "C" GPUBASELIBAPI void DeviceToDevice(void* s_deviceptr, void* t_deviceptr, long memsize);//GPU 内存数据转移 GPU -> 设备
extern "C" GPUBASELIBAPI void HostToDevice(void* hostptr, void* deviceptr, size_t memsize);//GPU 内存数据转移 设备 -> GPU
extern "C" GPUBASELIBAPI void DeviceToHost(void* hostptr, void* deviceptr, size_t memsize);//GPU 内存数据转移 GPU -> 设备
extern "C" GPUBASELIBAPI void DeviceToDevice(void* s_deviceptr, void* t_deviceptr, size_t memsize);//GPU 内存数据转移 GPU -> 设备
extern "C" GPUBASELIBAPI void CUDA_MemsetBlock(cuComplex* data, cuComplex init0, long len);
// 矢量基础运算函数
@ -109,5 +109,7 @@ 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, size_t len);
extern "C" GPUBASELIBAPI void HostCuComplexToshared_complexPtr(cuComplex* src, std::complex<double>* dst, size_t len);
#endif
#endif

View File

@ -36,12 +36,16 @@ ImageShowDialogClass ::ImageShowDialogClass(QWidget *parent)
ImageShowDialogClass::~ImageShowDialogClass()
{}
{
if (nullptr != this->desCursor) {
this->desCursor->close();
}
}
void ImageShowDialogClass::on_action_cursor_enable_trigged()
{
this->desCursor = new ImageShowCursorDesClass(this);
this->desCursor = new ImageShowCursorDesClass;
this->desCursorflag = true;
for (size_t i = 0; i < this->getGraphCount(); i++) {
@ -313,11 +317,14 @@ void ImageShowDialogClass::updateCursor(QMouseEvent *event)
QPoint pos = event->pos();
double x=this->m_plot->xAxis->pixelToCoord(pos.x()); // 将鼠标位置映射到图表坐标系中
double y = this->m_plot->yAxis->pixelToCoord(pos.y());
this->statusbar->showMessage(u8"X: "+QString::number(x,'f', 6)+" y: "+QString::number(y, 'f', 6));
QCPColorMap* colorMap = dynamic_cast<QCPColorMap*>(this->ui->m_plot->plottable(0));
double dataValue = colorMap->data()->data(x, y);
this->statusbar->showMessage(u8"X: "+QString::number(x,'f', 6)+" y: "+QString::number(y, 'f', 6)+
" value: "+QString::number(dataValue,'e',6));
if (this->desCursorflag) {
if (this->graphclass == LAMPDATASHOWCLASS::LAMPColorMap) {
QCPColorMap* colorMap = dynamic_cast<QCPColorMap*>(this->ui->m_plot->plottable(0));
double dataValue = colorMap->data()->data(x, y);
this->desCursor->updateCursorContent(u8"X: " + QString::number(x, 'f', 6) + " y: " + QString::number(y, 'f', 6) +u8"\n" +u8"DataValue: "+QString::number(dataValue));

View File

@ -85,7 +85,7 @@
<x>0</x>
<y>0</y>
<width>906</width>
<height>22</height>
<height>23</height>
</rect>
</property>
<widget class="QMenu" name="projectMenu">
@ -367,7 +367,7 @@ p, li { white-space: pre-wrap; }
<string>信息窗口</string>
</property>
<attribute name="dockWidgetArea">
<number>1</number>
<number>2</number>
</attribute>
<widget class="QWidget" name="dockWidgetContents_4">
<layout class="QVBoxLayout" name="verticalLayout_5">

View File

@ -0,0 +1,130 @@
#include <cstdio>
#include <cufft.h>
#include <cmath>
#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"
#include "BPBasic0_CUDA.cuh"
#include "GPUBpSimulation.cuh"
#include "GPURFPC.cuh"
double* getFreqPoints_mallocHost(double startFreq, double endFreq, long freqpoints)
{
long double dfreq = (endFreq - startFreq) / (freqpoints - 1);
double* freqlist = (double*)mallocCUDAHost(sizeof(double) * freqpoints);
for (long i = 0; i < freqpoints; i++) {
freqlist[i] = startFreq + dfreq * i;
}
return freqlist;
}
cuComplex* createEchoPhase_mallocHost(long Np, long Nf)
{
cuComplex* phdata = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * Np * Nf);
for (long i = 0; i < Np; i++) {
for (long j = 0; j < Nf; j++) {
phdata[i * Nf + j] = make_cuComplex(0, 0);
}
}
return phdata;
}
__global__ void kernel_RFPCProcess(
double* Sx,double* Sy,double* Sz,
double Tx, double Ty,double Tz,
double Tslx,double Tsly,double Tslz, // 目标的坡面向量
double p1,double p2, double p3, double p4, double p5, double p6,
long Np,long Nf,
double minF,double dFreq,double RefRange,
cuComplex* phdata
) {
long prfid = blockIdx.x * blockDim.x + threadIdx.x; // 获取当前的线程编码
if (prfid >= Np || prfid < 0) { return; }
else {}
// 距离
Vector3 S{ Sx[prfid],Sy[prfid],Sz[prfid] };
Vector3 T{ Tx,Ty,Tz };
Vector3 slp{ Tslx,Tsly,Tslz };
//printf("S=(%e,%e,%e)\n", S.x, S.y, S.z);
// 入射角
Vector3 TS = vec_sub(S, T);//T-->S
double localIncAngle = angleBetweenVectors(TS, slp, false);// 入射角
double sigma0 = GPU_getSigma0dB(p1, p2, p3, p4, p5, p6, localIncAngle);// 后向散射系数
sigma0 = powf(10.0, sigma0 / 10.0);
// 距离
double R = sqrt(vec_dot(TS, TS));
// 计算增益
double amp_echo = sigma0 / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4)); // 反射强度
double phi = 0;
for (long fid = 0; fid < Nf; fid++) {
phi = -4.0 * PI / LIGHTSPEED * (minF + dFreq * fid) * (R - RefRange);
phdata[prfid * Nf + fid].x += amp_echo * cos(phi);
phdata[prfid * Nf + fid].y += amp_echo * sin(phi);
//printf("phdata(%d,%d)=complex(%e,%e)\n", prfid, fid, phdata[prfid * Nf + fid].x, phdata[prfid * Nf + fid].y);
}
//printf("Nf:%d\n", Nf);
//printf("amp_echo:%f\n", amp_echo);
//printf("sigma0:%f\n", sigma0);
//printf("R:%e\n", R);
}
void RFPCProcess(double Tx, double Ty, double Tz,
double Tslx, double Tsly, double Tslz, // 目标的坡面向量
double p1, double p2, double p3, double p4, double p5, double p6,
GPUDATA& d_data)
{
double* AntX = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
double* AntY = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
double* AntZ = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
HostToDevice(d_data.AntX, AntX, sizeof(double) * d_data.Np); printf("antX host to device finished!!\n");
HostToDevice(d_data.AntY, AntY, sizeof(double) * d_data.Np); printf("antY host to device finished!!\n");
HostToDevice(d_data.AntZ, AntZ, sizeof(double) * d_data.Np); printf("antZ host to device finished!!\n");
double minF = d_data.minF[0];
long grid_size = (d_data.Np + BLOCK_SIZE - 1) / BLOCK_SIZE;
double dfreq = d_data.deltaF;
double R0 = d_data.R0;
kernel_RFPCProcess<<<grid_size , BLOCK_SIZE >>>(
AntX, AntY, AntZ,
Tx, Ty, Tz,
Tslx, Tsly, Tslz, // 目标的坡面向量
p1, p2, p3, p4, p5, p6,
d_data.Np, d_data.Nfft,
minF, dfreq,R0,
d_data.phdata
);
PrintLasterError("RFPCProcess");
FreeCUDADevice(AntX);
FreeCUDADevice(AntY);
FreeCUDADevice(AntZ);
return ;
}

View File

@ -0,0 +1,35 @@
/*****************************************************************//**
* \file GPUBpSimulation.cuh
* \brief GPU的局部仿真代码
*
* \author 30453
* \date March 2025
*********************************************************************/
#ifndef _GPUBPSIMUALTION_CUDA_H_
#define _GPUBPSIMUALTION_CUDA_H_
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "GPUTool.cuh"
#include "BPBasic0_CUDA.cuh"
extern "C" double* getFreqPoints_mallocHost(double startFreq, double endFreq, long freqpoints);
extern "C" cuComplex* createEchoPhase_mallocHost(long Np, long Nf);
extern "C" void RFPCProcess(
double Tx, double Ty, double Tz, // 目标点坐标
double Tslx,double Tsly,double Tslz, // 目标的坡面向量
double p1, double p2, double p3, double p4, double p5, double p6,// 地面目标后向散射系数与入射角关系 系数
GPUDATA& d_data
);
#endif

View File

@ -0,0 +1,301 @@
#include <cstdio>
#include <cufft.h>
#include <cmath>
#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"
#include "BPBasic0_CUDA.cuh"
#include <PrintMsgToQDebug.h>
#define c LIGHTSPEED
__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;
double phase = 4 * PI * Freq[freqIdx] * r / c;
double cos_phase = cos(phase);
double sin_phase = sin(phase);
cufftComplex ph = phdata[idx];
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);
}
__global__ void fftshiftKernel(cufftComplex* data, int Nfft, int Np) {
int pulse = blockIdx.x * blockDim.x + threadIdx.x;
if (pulse >= Np) return;
int half = Nfft / 2;
for (int i = 0; i < half; ++i) {
cufftComplex temp = data[pulse * Nfft + i];
data[pulse * Nfft + i] = data[pulse * Nfft + i + half];
data[pulse * Nfft + i + half] = temp;
}
}
__global__ void processPulseKernel(
long prfid,
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,
const double r_start, const double dr, const int nR,
cufftComplex* im_final
) {
//
long long idx = blockIdx.x * blockDim.x + threadIdx.x;
long long pixelcount = nx * ny;
if (idx >= pixelcount) return;
//printf("processPulseKernel start!!\n");
//if (x >= nx || y >= ny) return;
//int idx = x * ny + y;
double dx = AntX - x_mat[idx];
double dy = AntY - y_mat[idx];
double dz = AntZ - z_mat[idx];
//printf("processPulseKernel xmat !!\n");
double R = sqrt(dx * dx + dy * dy + dz * dz);
double dR = R - R0;
if (dR < r_start || dR >= (r_start + dr * (nR - 1))) return;
// Linear interpolation
double pos = (dR - r_start) / dr;
int index = (int)floor(pos);
double weight = pos - index;
if (index < 0 || index >= nR - 1) return;
cufftComplex rc_low = rc_pulse[prfid * nR +index];
cufftComplex rc_high = rc_pulse[prfid * nR+index + 1];
cufftComplex rc_interp;
rc_interp.x = rc_low.x * (1 - weight) + rc_high.x * weight;
rc_interp.y = rc_low.y * (1 - weight) + rc_high.y * weight;
// Phase correction
double phase = 4 * PI * minF / c * dR;
double cos_phase = cos(phase);
double sin_phase = sin(phase);
cufftComplex phCorr;
phCorr.x = rc_interp.x * cos_phase - rc_interp.y * sin_phase;
phCorr.y = rc_interp.x * sin_phase + rc_interp.y * cos_phase;
// Accumulate
im_final[idx].x += phCorr.x;
im_final[idx].y += phCorr.y;
//printf("r_start=%e;dr=%e;nR=%d\n", r_start, dr, nR);
if (abs(phCorr.x) > 1e-100 || abs(phCorr.y > 1e-100)) {
//printf(
// "[DEBUG] prfid=%-4ld | idx=%-8lld\n"
// " Ant: X=%-18.10e Y=%-18.10e Z=%-18.10e\n"
// " Pix: X=%-18.10e Y=%-18.10e Z=%-18.10e\n"
// " R=%-18.10e|dR=%-18.10e | pos=%-8.4f[%-6d+%-8.6f]\n"
// " RC: low=(%-18.10e,%-18.10e) high=(%-18.10e,%-18.10e)\n"
// " => interp=(%-18.10e,%-18.10e)\n"
// " Phase: val=%-18.10e | corr=(%-18.10e,%-18.10e)\n"
// " Final: im=(%-18.10e,%-18.10e)\n"
// "----------------------------------------\n",
// prfid, idx,
// AntX, AntY, AntZ,
// x_mat[idx], y_mat[idx], z_mat[idx],
// R,dR,
// pos, index, weight,
// rc_low.x, rc_low.y,
// rc_high.x, rc_high.y,
// rc_interp.x, rc_interp.y,
// phase,
// phCorr.x, phCorr.y,
// im_final[idx].x, im_final[idx].y
//);
}
}
void bpBasic0CUDA(GPUDATA& data, int flag,double* h_R) {
// Phase compensation
if (flag == 1) {
dim3 block(16, 16);
dim3 grid((data.K + 15) / 16, (data.Np + 15) / 16);
phaseCompensationKernel << <grid, block >> > (data.phdata, data.Freq, data.R0, data.K, data.Np);
PrintLasterError("bpBasic0CUDA Phase compensation");
//data.R0 = data.r; // 假设data.r已正确设置
}
// FFT处理
cufftHandle plan;
cufftPlan1d(&plan, data.Nfft, CUFFT_C2C, data.Np);
cufftExecC2C(plan, data.phdata, data.phdata, CUFFT_INVERSE);
cufftDestroy(plan);
// FFT移位
dim3 blockShift(256);
dim3 gridShift((data.Np + 255) / 256);
fftshiftKernel << <gridShift, blockShift >> > (data.phdata, data.Nfft, data.Np);
PrintLasterError("bpBasic0CUDA Phase FFT Process");
printfinfo("fft finished!!\n");
// 图像重建
double r_start = data.r_vec[0];
double dr = (data.r_vec[data.Nfft - 1] - r_start) / (data.Nfft - 1);
printfinfo("dr = %f\n",dr);
long pixelcount = data.nx* data.ny;
long grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
printfinfo("grid finished!!\n");
//double* d_R = (double*)mallocCUDADevice(sizeof(double) * data.nx * data.ny);
printfinfo("r_start=%e;dr=%e;nR=%d\n", r_start, dr, data.Nfft);
printfinfo("BPimage .....\n");
for (long ii = 0; ii < data.Np; ++ii) {
processPulseKernel << <grid_size, BLOCK_SIZE >> > (
ii,
data.nx, data.ny,
data.x_mat, data.y_mat, data.z_mat,
data.AntX[ii], data.AntY[ii], data.AntZ[ii],
data.R0, data.minF[ii],
data.phdata,
r_start, dr, data.Nfft,
data.im_final
//,d_R
);
PrintLasterError("processPulseKernel");
if (ii % 1000==0) {
printfinfo("\rPRF(%f %) %d / %d\t\t\t\t",(ii*100.0/data.Np), ii,data.Np);
}
}
//FreeCUDADevice(d_R);
PrintLasterError("bpBasic0CUDA Phase BPimage Process finished!!");
}
void initGPUData(GPUDATA& h_data, GPUDATA& d_data) {
d_data.AntX =h_data.AntX; //(double*)mallocCUDADevice(sizeof(double) * h_data.Np);
d_data.AntY = h_data.AntY;//(double*)mallocCUDADevice(sizeof(double) * h_data.Np);
d_data.AntZ = h_data.AntZ;// (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
d_data.minF = h_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 = h_data.r_vec;// (double*)mallocCUDADevice(sizeof(double) * h_data.Nfft);
d_data.Freq = (double*)mallocCUDADevice(sizeof(double) * h_data.Nfft);
d_data.phdata = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * h_data.Nfft * h_data.Np);
d_data.im_final = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * 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); printf("image X Copy finished!!!\n");
HostToDevice(h_data.y_mat, d_data.y_mat,sizeof(double) * h_data.nx * h_data.ny); printf("image Y Copy finished!!!\n");
HostToDevice(h_data.z_mat, d_data.z_mat, sizeof(double) * h_data.nx * h_data.ny); printf("image Z Copy finished!!!\n");
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(cuComplex) * h_data.Nfft * h_data.Np); printf("image echo Copy finished!!!\n");
HostToDevice(h_data.im_final, d_data.im_final, sizeof(cuComplex) * h_data.nx * h_data.ny); printf("image data Copy finished!!!\n");
// 拷贝标量参数
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;
//
// // 初始化主机数据
// h_data.Nfft = 1024;
// h_data.K = 512;
// // ... 其他参数初始化
//
// // 初始化设备内存
// initGPUData(h_data, d_data);
//
// // 执行算法
// bpBasic0CUDA(d_data, 0);
//
// // 拷贝结果回主机
// cudaCheckError(cudaMemcpy(h_data.im_final, d_data.im_final,
// sizeof(cufftComplex) * h_data.nx * h_data.ny, cudaMemcpyDeviceToHost));
//
// // 释放资源
// freeGPUData(d_data);
//
// return 0;
//}

View File

@ -0,0 +1,41 @@
#ifndef _BPBASIC0_CUDA_H_
#define _BPBASIC0_CUDA_H_
#include <cstdio>
#include <cufft.h>
#include <cmath>
#include <cuda_runtime.h>
#include "BaseConstVariable.h"
//#define cudaCheckError(ans) { gpuAssert((ans), __FILE__, __LINE__); }
//inline void gpuAssert(cudaError_t code, const char* file, int line) {
// if (code != cudaSuccess) {
// fprintf(stderr, "CUDA Error: %s %s %d\n", cudaGetErrorString(code), file, line);
// exit(code);
// }
//}
//
//#define c LIGHTSPEED
struct GPUDATA {
long Nfft, K, Np, nx, ny; // 傅里叶点数、频点数、脉冲数、图像列、图像行
double* AntX, * AntY, * AntZ, * minF; // 天线坐标、起始频率
double* x_mat, * y_mat, * z_mat;// 地面坐标
double* r_vec; // 坐标范围
double* Freq;// 频率
cuComplex* phdata;// 回波
cuComplex* im_final;// 图像
double R0; // 参考斜距
double deltaF; // 频点范围
};
extern "C" {
void bpBasic0CUDA(GPUDATA& data, int flag,double* h_R=nullptr);
void initGPUData(GPUDATA& h_data, GPUDATA& d_data);
void freeGPUData(GPUDATA& d_data);
void freeHostData(GPUDATA& d_data);
void BPBasic0(GPUDATA& h_data);
};
#endif

View File

@ -8,15 +8,12 @@
extern __device__ __host__ double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees = false);
extern __device__ __host__ double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees = false);
extern __device__ __host__ Vector3 vec_sub(Vector3 a, Vector3 b);
extern __device__ __host__ double vec_dot(Vector3 a, Vector3 b);
extern __device__ __host__ Vector3 vec_cross(Vector3 a, Vector3 b);
extern __device__ __host__ Vector3 vec_normalize(Vector3 v);
extern __device__ __host__ Vector3 vec_normalize(Vector3 v);
extern __device__ __host__ Vector3 compute_T(Vector3 S, Vector3 ray_dir, double H);
//
extern __device__ __host__ Vector3 compute_P(Vector3 S, Vector3 T, double R, double H );
//
//

View File

@ -19,13 +19,20 @@
/* »úÆ÷º¯Êý ****************************************************************************************************************************/
__device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta) {//线性值
extern __host__ __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta) {//线性值
double sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
return sigma;
}
extern __host__ __device__ double GPU_getSigma0dB(
const double p1, const double p2, const double p3, const double p4, const double p5, const double p6,
double theta) {//线性值
return p1 + p2 * exp(-p3 * theta) + p4 * cos(p5 * theta + p6);
}
__device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
double RstX, double RstY, double RstZ,
double AntXaxisX, double AntXaxisY, double AntXaxisZ,
double AntYaxisX, double AntYaxisY, double AntYaxisZ,
@ -104,7 +111,7 @@ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
return result;
}
__device__ double GPU_BillerInterpAntPattern(double* antpattern,
extern __host__ __device__ double GPU_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints,
double searththeta, double searchphi) {

View File

@ -22,15 +22,24 @@ extern "C" struct CUDASigmaParam {
extern __host__ __device__ double GPU_getSigma0dB(
const double p1, const double p2, const double p3, const double p4, const double p5, const double p6,
double theta);
extern __host__ __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta);
extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
double RstX, double RstY, double RstZ,
double AntXaxisX, double AntXaxisY, double AntXaxisZ,
double AntYaxisX, double AntYaxisY, double AntYaxisZ,
double AntZaxisX, double AntZaxisY, double AntZaxisZ,
double AntDirectX, double AntDirectY, double AntDirectZ
);
extern __host__ __device__ double GPU_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints,
double searththeta, double searchphi);

View File

@ -24,7 +24,6 @@
__global__ void kernel_TimeBPImageGridNet(double* antPx, double* antPy, double* antPz,
double* antDirx, double* antDiry, double* antDirz,
double* imgx, double* imgy, double* imgz,
@ -170,77 +169,230 @@ __global__ void kernel_TimeBPImageGridNet(double* antPx, double* antPy, double*
}
__device__ double computerR(double& Px, double& Py, double& Pz, double& Tx, double& Ty, double& Tz) {
//double R= sqrt((Px - Tx) * (Px - Tx) + (Py - Ty) * (Py - Ty) + (Pz - Tz) * (Pz - Tz));
//if (R > 900000) {
// printf("R=%f\n", R);
//}
//return R;
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];
if (isinf(s0.x) || isinf(s0.y) || isinf(s1.x) || isinf(s1.y)) {
return;
}
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 * PI / 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);
//imgArr[pixelidx] = cuCaddf(imgArr[pixelidx], make_cuComplex(1,1));
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 = startSegmentPrfId + 0;
double Px = antPx[currentprfid];
double Py = antPy[currentprfid];
double Pz = antPz[currentprfid];
double R0 = computerR(Px, Py, Pz, Tx, Ty, Tz);
if (LRrange <= R0 && R0 <= RRrange) {
updateBPImage(
currentprfid, pixelID, R0, LRrange,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar
);
}
// 10
currentprfid = startSegmentPrfId + 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 = startSegmentPrfId + 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 < 20; 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(
double* antPx, double* antPy, double* antPz,
double* imgx, double* imgy, double* imgz,
cuComplex* TimeEchoArr, long prfcount, long pointCount,
cuComplex* imgArr, long imH, long imW,
double startLamda, double Rnear, double dx,double RefRange
double startLamda, double Rnear, double dx,double RefRange,double Rfar
) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
long pixelcount = imH * imW;
if (idx < pixelcount) {
double Tx = imgx[idx]; // 地面坐标点
double Ty = imgy[idx];
double Tz = imgz[idx];
double Rrange = pointCount * dx;// 成像范围
double LRrange = RefRange - Rrange / 2;//左范围
double RRrange = RefRange + Rrange / 2;
double Tx = imgx[idx], Ty = imgy[idx], Tz = imgz[idx], PR = 0;
double Px = 0, Py = 0, Pz = 0;
double Rid = -1;
long maxPointIdx = pointCount - 1;
long pid = 0;
long pid0 = 0;
long pid1 = 0;
double phi = 0;
cuComplex s0=make_cuComplex(0,0), s1=make_cuComplex(0,0);
cuComplex s = make_cuComplex(0, 0);
cuComplex phiCorr = make_cuComplex(0, 0);
for (long spid = 0; spid < prfcount; spid = spid + 8) {
#pragma unroll
for (long sid = 0; sid < 8; sid++)
{
pid = spid + sid;
if (pid < prfcount) {}
else {
continue;
for (long segid = 0; segid < prfcount; segid = segid + 20) {
long seglen = prfcount - segid;
if (seglen < 20) {
for (long i = 1; i < 10; i++) {
long currentprfid = segid + i;
if (currentprfid < prfcount) {
double Px = antPx[currentprfid];
double Py = antPy[currentprfid];
double Pz = antPz[currentprfid];
double R = computerR(Px, Py, Pz, Tx, Ty, Tz);
updateBPImage(
currentprfid, idx, R, LRrange,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar
);
}
}
Px = antPx[pid];
Py = antPy[pid];
Pz = antPz[pid];
PR = sqrt((Px - Tx) * (Px - Tx) + (Py - Ty) * (Py - Ty) + (Pz - Tz) * (Pz - Tz));
Rid = (PR - Rnear) / dx; // 行数
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);
}
else {
// 判断范围
segmentBPImage(
antPx, antPy, antPz,
Tx, Ty, Tz,
TimeEchoArr, prfcount, pointCount,
imgArr,
startLamda, Rnear, dx, RefRange, Rfar,
segid, idx
);
}
}
}
}
@ -254,7 +406,8 @@ extern "C" {
double* antDirx, double* antDiry, double* antDirz,
double* imgx, double* imgy, double* imgz,
long prfcount, long freqpoints, double meanH,
double Rnear, double dx, double RefRange)
double Rnear, double dx, double RefRange
)
{
long pixelcount = prfcount * freqpoints;
int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
@ -272,22 +425,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,
cuComplex* TimeEchoArr, long prfcount, long pointCount,
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;
int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
kernel_pixelTimeBP << <1, 1 >> > (
kernel_pixelTimeBP << <grid_size, BLOCK_SIZE >> > (
antPx, antPy, antPz,
imgx, imgy, imgz,
TimeEchoArr, prfcount, pointCount,
imgArr, imH, imW,
startLamda, Rnear, dx, RefRange
startLamda, Rnear, dx, RefRange, Rfar
);
PrintLasterError("TimeBPImage");

View File

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

View File

@ -15,8 +15,13 @@ QSimulationBPImage::QSimulationBPImage(QWidget *parent)
QObject::connect(ui->pushButtonEchoSelect, SIGNAL(clicked()), this, SLOT(onpushButtonEchoSelectClicked()));
QObject::connect(ui->pushButtonImageSelect, SIGNAL(clicked()), this, SLOT(onpushButtonImageSelectClicked()));
QObject::connect(ui->GridNetBtn, SIGNAL(clicked()), this, SLOT(onpushButtonGridNetSelectClicked()));
QObject::connect(ui->buttonBox, SIGNAL(accepted()), this, SLOT(onbtnaccepted()));
QObject::connect(ui->buttonBox, SIGNAL(rejected()), this, SLOT(onbtnrejected()));
QObject::connect(ui->checkBox, SIGNAL(stateChanged(int)), this, SLOT(oncheckboxTrigged(int)));
}
@ -62,6 +67,7 @@ void QSimulationBPImage::onpushButtonImageSelectClicked()
void QSimulationBPImage::onbtnaccepted()
{
this->hide();
QString echofile = this->ui->lineEditEchoPath->text().trimmed();
QString outImageFolder = getParantFromPath(this->ui->lineEditImagePath->text().trimmed());
QString imagename = getFileNameFromPath(this->ui->lineEditImagePath->text().trimmed());
@ -76,7 +82,14 @@ void QSimulationBPImage::onbtnaccepted()
imagL1->setFarRange(echoL0ds->getFarRange());
imagL1->setFs(echoL0ds->getFs());
imagL1->setLookSide(echoL0ds->getLookSide());
imagL1->OpenOrNew(outImageFolder, imagename, echoL0ds->getPluseCount(), echoL0ds->getPlusePoints());
if (ui->checkBox->isChecked()) {
gdalImage imgxyzimg(ui->lineEdit->text().trimmed());
imagL1->OpenOrNew(outImageFolder, imagename, imgxyzimg.height, imgxyzimg.width);
}
else {
imagL1->OpenOrNew(outImageFolder, imagename, echoL0ds->getPluseCount(), echoL0ds->getPlusePoints());
}
TBPImageAlgCls TBPimag;
@ -84,7 +97,13 @@ void QSimulationBPImage::onbtnaccepted()
TBPimag.setImageL1(imagL1);
long cpucore_num = std::thread::hardware_concurrency();
TBPimag.setGPU(true);
TBPimag.Process(cpucore_num);
if (ui->checkBox->isChecked()) {
TBPimag.ProcessWithGridNet(cpucore_num,ui->lineEdit->text().trimmed());
}
else {
TBPimag.Process(cpucore_num);
}
this->show();
QMessageBox::information(this,u8"³ÉÏñ",u8"³ÉÏñ½áÊø");
}
@ -93,3 +112,28 @@ void QSimulationBPImage::onbtnrejected()
{
this->close();
}
void QSimulationBPImage::oncheckboxTrigged(int)
{
this->ui->lineEdit->setEnabled(this->ui->checkBox->isChecked());
this->ui->GridNetBtn->setEnabled(this->ui->checkBox->isChecked());
}
void QSimulationBPImage::onpushButtonGridNetSelectClicked( )
{
QString fileNames = QFileDialog::getOpenFileName(
this, // 父窗口
tr(u8"选择影像文件"), // 标题
QString(), // 默认路径
tr(u8"All Files(*);;dat(*.dat);;tif(*.tif);;tiff(*.tiff)") // 文件过滤器
);
// 如果用户选择了文件
if (!fileNames.isEmpty()) {
this->ui->lineEdit->setText(fileNames);
}
else {
QMessageBox::information(this, tr(u8"没有选择文件"), tr(u8"没有选择任何文件。"));
}
}

View File

@ -23,6 +23,9 @@ public slots:
void onbtnaccepted();
void onbtnrejected();
void oncheckboxTrigged(int);
void onpushButtonGridNetSelectClicked( );
private:
Ui::QSimulationBPImageClass* ui;
};

View File

@ -45,7 +45,7 @@
</size>
</property>
<property name="text">
<string>D:\Programme\vs2022\RasterMergeTest\LAMPCAE_SCANE\GF3_Simulation.xml</string>
<string>D:\Programme\vs2022\RasterMergeTest\LAMPCAE_SCANE-all-scane\GF3_Simulation.xml</string>
</property>
</widget>
</item>
@ -65,6 +65,66 @@
</layout>
</widget>
</item>
<item>
<widget class="QFrame" name="frame_3">
<property name="frameShape">
<enum>QFrame::StyledPanel</enum>
</property>
<property name="frameShadow">
<enum>QFrame::Raised</enum>
</property>
<layout class="QHBoxLayout" name="horizontalLayout_3">
<item>
<widget class="QCheckBox" name="checkBox">
<property name="minimumSize">
<size>
<width>0</width>
<height>30</height>
</size>
</property>
<property name="text">
<string>图像网格</string>
</property>
<property name="checked">
<bool>true</bool>
</property>
</widget>
</item>
<item>
<widget class="QLineEdit" name="lineEdit">
<property name="enabled">
<bool>true</bool>
</property>
<property name="minimumSize">
<size>
<width>0</width>
<height>30</height>
</size>
</property>
<property name="text">
<string>D:\Programme\vs2022\RasterMergeTest\simulationData\demdataset\demxyz.bin</string>
</property>
</widget>
</item>
<item>
<widget class="QPushButton" name="GridNetBtn">
<property name="enabled">
<bool>true</bool>
</property>
<property name="minimumSize">
<size>
<width>0</width>
<height>30</height>
</size>
</property>
<property name="text">
<string>选择</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item>
<widget class="QFrame" name="frame_2">
<property name="frameShape">
@ -96,7 +156,7 @@
</size>
</property>
<property name="text">
<string>D:\Programme\vs2022\RasterMergeTest\LAMPCAE_SCANE\BPImage\GF3BPImage</string>
<string>D:\Programme\vs2022\RasterMergeTest\LAMPCAE_SCANE-all-scane\BPImage\GF3BPImage</string>
</property>
</widget>
</item>

View File

@ -1,4 +1,4 @@
#include "stdafx.h"
#include "stdafx.h"
#include "TBPImageAlgCls.h"
#include <QDateTime>
#include <QDebug>
@ -9,6 +9,10 @@
#include "GPUTool.cuh"
#include "GPUTBPImage.cuh"
#include "ImageOperatorBase.h"
#include "BPBasic0_CUDA.cuh"
#include "BaseTool.h"
#include <GPUBpSimulation.cuh>
void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZPath)
{
@ -19,7 +23,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 +59,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 +71,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 +130,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,26 +213,31 @@ 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);
this->outRasterXYZPath = outRasterXYZ;
return this->ProcessWithGridNet(num_thread, outRasterXYZ);
}
qDebug() << u8"频域回波-> 时域回波";
this->TimeEchoDataPath = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_Timeecho.bin");
this->EchoFreqToTime();
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();
// 初始化Raster
qDebug() << u8"初始化影像";
// 初始化Raster
qDebug() << u8"初始化影像";
long imageheight = this->L1ds->getrowCount();
long imagewidth = this->L1ds->getcolCount();
long blokline = Memory1GB / 8 / 4 / imageheight * 32;
long blokline = Memory1GB / 8 / 4 / imageheight * 32;
blokline = blokline < 1000 ? 1000 : blokline;
long startline = 0;
@ -237,186 +246,211 @@ ErrorCode TBPImageAlgCls::Process(long num_thread)
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 ;
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);
imageRaster.get()[i * imagewidth + j] = std::complex<double>(0, 0);
}
}
this->L1ds->saveImageRaster(imageRaster, startline,templine);
this->L1ds->saveImageRaster(imageRaster, startline, templine);
}
qDebug() << u8"频域回波-> 时域回波 结束";
qDebug() << u8"频域回波-> 时域回波 结束";
//if (GPURUN) {
// return this->ProcessGPU();
//}
//else {
// QMessageBox::information(nullptr,u8"提示",u8"目前只支持显卡");
// return ErrorCode::FAIL;
//}
if (GPURUN) {
return this->ProcessGPU();
}
else {
QMessageBox::information(nullptr, u8"提示", u8"目前只支持显卡");
return ErrorCode::FAIL;
}
return ErrorCode::SUCCESS;
}
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();
std::shared_ptr<SARSimulationImageL1Dataset> L1ds=this->L1ds;
std::shared_ptr < EchoL0Dataset> L0ds=this->L0ds;
QString inGPSPath = L0ds->getGPSPointFilePath();
QString echoFilePath = L0ds->getEchoDataFilePath();
QString imgXYZPath = this->outRasterXYZPath;
QString outimgDataPath = L1ds->getImageRasterPath();
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;
size_t prfcount = L0ds->getPluseCount();
size_t freqpoints = L0ds->getPlusePoints();
size_t block_pfrcount = Memory1MB / freqpoints / 8 * 2000;// 4GB -- 可以分配内存
// 反方向计算起始相位
double deltaF = bandwidth / (PlusePoints - 1);
double startfreq = freq - bandwidth / 2;
double startlamda = LIGHTSPEED / startfreq;
qDebug() << "deltaF:\t" << deltaF;
size_t img_rowCont = L1ds->getrowCount();
size_t img_colCont = L1ds->getcolCount();
size_t block_imgRowCount = Memory1MB / img_colCont / 8 / 3 * 2000;// 4GB-- 可以分配内存
gdalImage demgridimg(imgXYZPath);
gdalImageComplex im_finalds(outimgDataPath);
gdalImageComplex echods(echoFilePath);
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);
// 加载 GPS it should be same as prfcount;
long gpspoints = prfcount;
std::shared_ptr<SatelliteAntPos> antpos = SatelliteAntPosOperator::readAntPosFile(inGPSPath, gpspoints);
GPUDATA h_data{};// 任务参数
{
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;
//回波 频率参数
h_data.Nfft = freqpoints;
qDebug() << u8"3.proces echo params:";
double centerFreq = L0ds->getCenterFreq();
double bandwidth = L0ds->getBandwidth();
size_t freqpoints = L0ds->getPlusePoints();
h_data.Freq = getFreqPoints_mallocHost(centerFreq - bandwidth / 2, centerFreq + bandwidth / 2, freqpoints);
h_data.minF = (double*)mallocCUDAHost(sizeof(double) * gpspoints);
for (long i = 0; i < gpspoints; i++) {
h_data.minF[i] = h_data.Freq[0];
}
antpos.reset();
h_data.deltaF = bandwidth / (freqpoints - 1);
h_data.K = h_data.Nfft;
h_data.R0 = L0ds->getRefPhaseRange();
double deltaF = h_data.deltaF; // 从输入参数获取
double maxWr = 299792458.0 / (2.0 * deltaF);
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.0 * step;
for (int i = 0; i < h_data.Nfft; ++i) {
r_vec_host[i] = start + i * step;
}
h_data.r_vec = r_vec_host;
// 处理天线坐标
h_data.AntX = (double*)mallocCUDAHost(sizeof(double) * block_pfrcount);
h_data.AntY = (double*)mallocCUDAHost(sizeof(double) * block_pfrcount);
h_data.AntZ = (double*)mallocCUDAHost(sizeof(double) * block_pfrcount);
qDebug() << "r_vec [0]:\t" << h_data.r_vec[0];
qDebug() << "Range resolution(m):\t" << step;
qDebug() << "range Scence(m):\t" << maxWr;
qDebug() << "start Freq:\t" << centerFreq - bandwidth / 2;
qDebug() << "end Freq:\t" << centerFreq + bandwidth / 2;
qDebug() << "freq points:\t" << freqpoints;
qDebug() << "delta freq:\t" << h_data.deltaF;
qDebug() << "prf count:\t" << gpspoints;
qDebug() << "-----------------------------------------------------------";
}
qDebug() << "BP......";
for (long img_start_rid = 0; img_start_rid < img_rowCont; img_start_rid = img_start_rid + block_imgRowCount) {
long img_blockRowCount = block_imgRowCount;
long img_blockColCount = img_colCont;
std::shared_ptr<double> demX = readDataArr<double>(demgridimg, img_start_rid, 0, img_blockRowCount, img_blockColCount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> demY = readDataArr<double>(demgridimg, img_start_rid, 0, img_blockRowCount, img_blockColCount, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> demZ = readDataArr<double>(demgridimg, img_start_rid, 0, img_blockRowCount, img_blockColCount, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
h_data.x_mat = (double*)mallocCUDAHost(sizeof(double) * img_blockRowCount * img_blockColCount); // 成像网格
h_data.y_mat = (double*)mallocCUDAHost(sizeof(double) * img_blockRowCount * img_blockColCount);
h_data.z_mat = (double*)mallocCUDAHost(sizeof(double) * img_blockRowCount * img_blockColCount);
h_data.nx = img_blockColCount;
h_data.ny = img_blockRowCount;
// 计算成像的基本参数
// 距离向分辨率
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);
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;
for (long i = 0; i < h_data.ny; i++) {
for (long j = 0; j < h_data.nx; j++) {
h_data.x_mat[i * h_data.nx + j] = demX.get()[i * h_data.nx + j];
h_data.y_mat[i * h_data.nx + j] = demY.get()[i * h_data.nx + j];
h_data.z_mat[i * h_data.nx + j] = demZ.get()[i * h_data.nx + j];
}
}
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) {
std::shared_ptr<std::complex<double>> imgArr = im_finalds.getDataComplexSharePtr(img_start_rid, 0, img_blockRowCount, img_blockColCount, 1);
h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex)* img_blockRowCount* img_blockColCount);
for (long echo_start_pid = 0; echo_start_pid < prfcount; echo_start_pid = echo_start_pid + block_pfrcount) {
long echo_blockPRFCount = block_pfrcount;
long echo_blockFreqCount = freqpoints;
// 读取回波
std::shared_ptr<std::complex<double>> echoData = readDataArrComplex<std::complex<double>>(echods, echo_start_pid,0, echo_blockPRFCount, echo_blockFreqCount,1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
size_t echosize = sizeof(cuComplex) * echo_blockPRFCount * echo_blockFreqCount;
h_data.phdata = (cuComplex*)mallocCUDAHost(echosize);
shared_complexPtrToHostCuComplex(echoData.get(), h_data.phdata, size_t(echo_blockPRFCount * echo_blockFreqCount));
// 处理天线坐标
h_data.Np = echo_blockPRFCount;
for (long i = 0; i < h_data.Np; i++) {
h_data.AntX[i] = antpos.get()[i+echo_start_pid].Px;
h_data.AntY[i] = antpos.get()[i+echo_start_pid].Py;
h_data.AntZ[i] = antpos.get()[i+echo_start_pid].Pz;
}
{
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]);
}
GPUDATA d_data;
initGPUData(h_data, d_data);
bpBasic0CUDA(d_data, 0);
DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * d_data.nx * d_data.ny);
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];
}
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;
freeGPUData(d_data);
FreeCUDAHost(h_data.phdata);
qDebug() << QString(u8"\nimg proces [%1 precent ] , echo process [%2 precess]\t:img row [%3%4,echo pfr [%5,%6]")
.arg((img_start_rid+ img_blockRowCount)*100.0/ img_rowCont)
.arg((echo_start_pid + echo_blockPRFCount) * 100.0 / prfcount)
.arg(img_start_rid)
.arg(img_start_rid + img_blockRowCount)
.arg(echo_start_pid)
.arg(echo_start_pid + echo_blockPRFCount)
;
//HostCuComplexToshared_complexPtr(h_data.im_final, imgArr.get(), size_t(h_data.nx)* size_t(h_data.ny));
//testOutComplexDoubleArr(QString("im_final.bin"), imgArr.get(), h_data.ny, h_data.nx);
}
this->L1ds->saveImageRaster(imgArr, startimgrowid, tempimgBlockline);
HostCuComplexToshared_complexPtr(h_data.im_final, imgArr.get(), size_t(h_data.nx)* size_t(h_data.ny));
im_finalds.saveImage(imgArr, img_start_rid, 0,img_blockRowCount, img_blockColCount, 1);
}
qDebug() << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz") << "] image writing:\t" << this->L1ds->getxmlFilePath();
this->L1ds->saveToXml();
qDebug() << QString("img proces [100 precent ] , echo process [100 precess] ");
L1ds->saveToXml();
qDebug() << "bp Image Result write to file :\t" << L1ds->getoutFolderPath();
return ErrorCode::SUCCESS;
}
void TBPImageAlgCls::EchoFreqToTime( )
ErrorCode TBPImageAlgCls::BPProcessBlockGPU()
{
// 读取数据
return ErrorCode::SUCCESS;
}
void TBPImageAlgCls::setGPU(bool flag)
{
this->GPURUN = flag;
}
bool TBPImageAlgCls::getGPU( )
{
return this->GPURUN;
}
void TBPImageAlgCls::EchoFreqToTime()
{
// 读取数据
long PRFCount = this->L0ds->getPluseCount();
long inColCount = this->L0ds->getPlusePoints();
@ -427,10 +461,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);
// 创建二进制文件
gdalImageComplex outTimeEchoImg = CreategdalImageComplexNoProj(this->TimeEchoDataPath, this->TimeEchoRowCount, this->TimeEchoColCount, 1);
// 分块
// 分块
long echoBlockline = Memory1GB / 8 / 2 / outColCount * 3; //1GB
echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;
@ -445,14 +479,14 @@ void TBPImageAlgCls::EchoFreqToTime( )
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);
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 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());
host_echoArr.get()[ii * outColCount + jj] = make_cuComplex(echoArr.get()[ii * inColCount + jj].real(), echoArr.get()[ii * inColCount + jj].imag());
}
}
#pragma omp parallel for
@ -472,7 +506,7 @@ void TBPImageAlgCls::EchoFreqToTime( )
DeviceToHost(host_IFFTechoArr.get(), device_IFFTechoArr.get(), sizeof(cuComplex) * tempechoBlockline * outColCount);
#pragma omp parallel for
for (long ii = 0; ii < tempechoBlockline ; ii++) {
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);
@ -490,321 +524,16 @@ void TBPImageAlgCls::EchoFreqToTime( )
}
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
)
bool TBPImageAlgCls::checkZeros(double* data, long long len)
{
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_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];
}
bool flag = true;
#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];
for (long long i = 0; i < len; i++) {
if (abs(data[i]) < 1e-7) {
flag= false;
break;
}
}
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);
}
}
return flag;
}
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;
}
*/

View File

@ -51,34 +51,41 @@ public:
public:
ErrorCode Process(long num_thread);
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();
private:
bool checkZeros(double* data, long long len);
};
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,
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
//);

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>
@ -179,7 +182,7 @@
<IntrinsicFunctions>true</IntrinsicFunctions>
</ClCompile>
<Link>
<SubSystem>Windows</SubSystem>
<SubSystem>Console</SubSystem>
<GenerateDebugInformation>DebugFull</GenerateDebugInformation>
<EnableCOMDATFolding>true</EnableCOMDATFolding>
<OptimizeReferences>true</OptimizeReferences>
@ -217,8 +220,12 @@
<ClCompile Include="SimulationSAR\SatelliteOribtModel.cpp" />
<ClCompile Include="SimulationSAR\SigmaDatabase.cpp" />
<ClCompile Include="SimulationSAR\TBPImageAlgCls.cpp" />
<ClCompile Include="UnitTestMain.cpp" />
<CudaCompile Include="GPUBpSimulation.cu" />
<CudaCompile Include="SimulationSAR\GPUTBPImage.cu" />
<QtMoc Include="PowerSimulationIncoherent\QSimulationSARPolynomialOrbitModel.h" />
<CudaCompile Include="PowerSimulationIncoherent\LookTableSimulationComputer.cuh" />
<ClInclude Include="GPUBpSimulation.cuh" />
<ClInclude Include="PowerSimulationIncoherent\OribtModelOperator.h" />
<QtMoc Include="PowerSimulationIncoherent\QSimulationLookTableDialog.h" />
<QtMoc Include="PowerSimulationIncoherent\QCreateSARIntensityByLookTableDialog.h" />
@ -228,7 +235,10 @@
<QtMoc Include="SimulationSAR\QSARLookTableSimualtionGUI.h" />
<QtMoc Include="SimulationSAR\QSimulationBPImage.h" />
<QtMoc Include="SimulationSAR\QSimulationRFPCGUI.h" />
<CudaCompile Include="SimulationSAR\BPBasic0_CUDA.cu" />
<CudaCompile Include="SimulationSAR\GPUBPTool.cuh" />
<ClInclude Include="SimulationSAR\BPBasic0_CUDA.cuh" />
<ClInclude Include="SimulationSAR\GPUTBPImage.cuh" />
<ClInclude Include="SimulationSAR\RFPCProcessCls.h" />
<ClInclude Include="SimulationSAR\SARSatelliteSimulationAbstractCls.h" />
<ClInclude Include="SimulationSAR\SARSimulationTaskSetting.h" />
@ -245,10 +255,6 @@
<FileType>Document</FileType>
</CudaCompile>
<CudaCompile Include="SimulationSAR\GPURFPC.cuh" />
<CudaCompile Include="SimulationSAR\GPUTBPImage.cu">
<FileType>Document</FileType>
</CudaCompile>
<CudaCompile Include="SimulationSAR\GPUTBPImage.cuh" />
</ItemGroup>
<ItemGroup>
<QtUic Include="PowerSimulationIncoherent\QCreateSARIntensityByLookTableDialog.ui" />

View File

@ -59,6 +59,15 @@
<ClInclude Include="PowerSimulationIncoherent\OribtModelOperator.h">
<Filter>PowerSimulationIncoherent</Filter>
</ClInclude>
<ClInclude Include="SimulationSAR\BPBasic0_CUDA.cuh">
<Filter>SimulationSAR</Filter>
</ClInclude>
<ClInclude Include="SimulationSAR\GPUTBPImage.cuh">
<Filter>SimulationSAR</Filter>
</ClInclude>
<ClInclude Include="GPUBpSimulation.cuh">
<Filter>SimulationSAR</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<ClCompile Include="SimulationSAR\QImageSARRFPC.cpp">
@ -103,6 +112,9 @@
<ClCompile Include="PowerSimulationIncoherent\QCreateSARIntensityByLookTableDialog.cpp">
<Filter>PowerSimulationIncoherent</Filter>
</ClCompile>
<ClCompile Include="UnitTestMain.cpp">
<Filter>Source Files</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<QtUic Include="SimulationSAR\QImageSARRFPC.ui">
@ -160,12 +172,6 @@
<CudaCompile Include="SimulationSAR\GPURFPC.cuh">
<Filter>SimulationSAR</Filter>
</CudaCompile>
<CudaCompile Include="SimulationSAR\GPUTBPImage.cu">
<Filter>SimulationSAR</Filter>
</CudaCompile>
<CudaCompile Include="SimulationSAR\GPUTBPImage.cuh">
<Filter>SimulationSAR</Filter>
</CudaCompile>
<CudaCompile Include="PowerSimulationIncoherent\LookTableSimulationComputer.cu">
<Filter>PowerSimulationIncoherent</Filter>
</CudaCompile>
@ -178,5 +184,14 @@
<CudaCompile Include="SimulationSAR\GPUBPTool.cuh">
<Filter>SimulationSAR</Filter>
</CudaCompile>
<CudaCompile Include="SimulationSAR\BPBasic0_CUDA.cu">
<Filter>SimulationSAR</Filter>
</CudaCompile>
<CudaCompile Include="SimulationSAR\GPUTBPImage.cu">
<Filter>SimulationSAR</Filter>
</CudaCompile>
<CudaCompile Include="GPUBpSimulation.cu">
<Filter>SimulationSAR</Filter>
</CudaCompile>
</ItemGroup>
</Project>

View File

@ -0,0 +1,421 @@
/*
* 仿
*
*
*/
#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;
//}