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

Merge pull request 'RTPC-BPImage_dev' (#5) from RTPC-BPImage_dev into RTPC-dev 

Reviewed-on: http://123.153.4.249:22779/LAMPSARToolSoftware/RasterProcessTool/pulls/5
pull/6/head^2
chenzenghui 2025-03-02 16:13:32 +08:00
parent b34971e590 982922506c
commit b15a7d0e35
16 changed files with 898 additions and 132 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ALLRelease/ALLRelease.vcxproj
View File

@ -159,6 +159,7 @@
<ConformanceMode>true</ConformanceMode>
<LanguageStandard>stdcpp14</LanguageStandard>
<LanguageStandard_C>stdc11</LanguageStandard_C>
<EnableParallelCodeGeneration>true</EnableParallelCodeGeneration>
</ClCompile>
<Link>
<SubSystem>Console</SubSystem>

7
BaseCommonLibrary/BaseTool/BaseConstVariable.h
View File

@ -66,6 +66,9 @@ inline char* get_cur_time() {
#define EARTHWE 0.000072292115
#define PI 3.141592653589793238462643383279
#define WGS84_A 6378137.0 // 长半轴 (m)
#define WGS84_F (1.0/298.257223563)
#define WGS84_B (WGS84_A*(1-WGS84_F)) // 短半轴 (m)
#define earthRoute 0.000072292115
@ -140,6 +143,10 @@ struct DemBox {
double max_lat;
};
struct Vector3 {
double x, y, z;
};
/*********************************************** FEKO仿真参数 ********************************************************************/
struct SatellitePos {

124
BaseCommonLibrary/BaseTool/ImageOperatorBase.cpp
View File

@ -1207,7 +1207,7 @@ void gdalImage::saveImage(std::shared_ptr<double> data, int start_row, int start
}
else {
poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
this->band_num, GDT_Float32, NULL); // 斤拷锟斤拷
this->band_num, GDT_Float64, NULL); // 斤拷锟斤拷
poDstDS->SetProjection(this->projection.toUtf8().constData());
double gt_ptr[6];
@ -1690,7 +1690,7 @@ gdalImage CreategdalImageDouble(const QString& img_path, int height, int width,
}
gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num,
Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite, bool isEnvi)
Eigen::MatrixXd gt, QString projection,bool need_gt, bool overwrite, bool isEnvi, GDALDataType datetype)
{
if(exists_test(img_path.toUtf8().constData())) {
if(overwrite) {
@ -1713,7 +1713,7 @@ gdalImage CreategdalImage(const QString& img_path, int height, int width, int ba
GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
GDT_Float32, NULL); // 锟斤拷锟斤拷锟斤拷
datetype, NULL); // 锟斤拷锟斤拷锟斤拷
if(need_gt) {
if (!projection.isEmpty()) {
poDstDS->SetProjection(projection.toUtf8().constData());
@ -1822,6 +1822,28 @@ gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int
return result_img;
}
gdalImageComplex BASECONSTVARIABLEAPI CreategdalImageComplexNoProj(const QString& img_path, int height, int width, int band_num, bool overwrite)
{
// 创建新文件
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
GDALDataset* poDstDS = (poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_CFloat64, NULL));
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
//poDstDS.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
gdalImageComplex result_img(img_path);
return result_img;
}
gdalImageComplex CreateEchoComplex(const QString& img_path, int height, int width, int band_num)
{
// 创建图像
@ -2598,6 +2620,60 @@ void gdalImageComplex::saveImage(Eigen::MatrixXcd data, int start_row, int start
omp_destroy_lock(&lock); //
}
void gdalImageComplex::saveImage(std::shared_ptr<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.get()[i*colCount+j].real();
databuffer[i * colCount * 2 + j * 2 + 1] = data.get()[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)
{
@ -2617,8 +2693,8 @@ Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col,
poBand = poDataset->GetRasterBand(1);
// 读取波段信息，假设是复数类型
int nXSize = poBand->GetXSize();
int nYSize = poBand->GetYSize();
int nXSize = cols_count; poBand->GetXSize();
int nYSize = rows_count; poBand->GetYSize();
double* databuffer = new double[nXSize * nYSize * 2];
poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count,
@ -2635,6 +2711,44 @@ Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col,
delete[] databuffer;
return rasterData;
}
std::shared_ptr<std::complex<double>> gdalImageComplex::getDataComplexSharePtr(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
{
GDALDataset* poDataset;
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// 打开TIFF文件
poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
if(poDataset == nullptr) {
QMessageBox::warning(nullptr, u8"错误", u8"无法打开文件：" + this->img_path);
qDebug() << u8"无法打开文件：" + this->img_path;
}
// 获取数据集的第一个波段
GDALRasterBand* poBand;
poBand = poDataset->GetRasterBand(1);
// 读取波段信息，假设是复数类型
double* databuffer = new double[rows_count * cols_count * 2];
poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count,
rows_count, GDT_CFloat64, 0, 0);
GDALClose((GDALDatasetH)poDataset);
std::shared_ptr<std::complex<double>> rasterData(new std::complex<double>[rows_count * cols_count], delArrPtr);
for(size_t i = 0; i < rows_count; i++) {
for(size_t j = 0; j < cols_count; j++) {
rasterData.get()[i*cols_count+ j] = std::complex<double>(databuffer[i * cols_count * 2 + j * 2],
databuffer[i * cols_count * 2 + j * 2 + 1]);
}
}
delete[] databuffer;
return rasterData;
}
void gdalImageComplex::saveImage()
{

5
BaseCommonLibrary/BaseTool/ImageOperatorBase.h
View File

@ -228,7 +228,9 @@ public: // 方法
~gdalImageComplex();
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);
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);
void saveImage() override;
void savePreViewImage();
@ -238,11 +240,12 @@ public:
// 创建影像
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);
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);
gdalImage BASECONSTVARIABLEAPI CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long espgcode, GDALDataType eType = GDT_Float32, bool need_gt = true, bool overwrite = false, bool isENVI = false);
gdalImageComplex BASECONSTVARIABLEAPI CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false);
gdalImageComplex BASECONSTVARIABLEAPI CreategdalImageComplexNoProj(const QString& img_path, int height, int width, int band_num, bool overwrite = true);
gdalImageComplex BASECONSTVARIABLEAPI CreateEchoComplex(const QString& img_path, int height, int width, int band_num);

45
GPUBaseLib/GPUTool/GPUTool.cu
View File

@ -213,6 +213,45 @@ __global__ void CUDA_GridPoint_Linear_Interp1(float* v, float* q, float* qv, lon
}
}
// 一维FFTShift核函数
__global__ void fftshift_1d_kernel(cuComplex* data, int batch_size, int signal_length) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx >= batch_size * signal_length) return;
int batch_id = idx / signal_length;
int signal_id = idx % signal_length;
int half = (signal_length + 1) / 2; // 兼容奇偶长度
if (signal_id >= half) return;
int new_pos = (signal_id + half) % signal_length;
int src_idx = batch_id * signal_length + new_pos;
// 数据交换
cuComplex temp = data[idx];
data[idx] = data[src_idx];
data[src_idx] = temp;
}
// 批量一维FFTShift函数
extern "C" void FFTShift1D(cuComplex* d_data, int batch_size, int signal_length) {
if (signal_length <= 1) return; // 无需处理
// 启动核函数
int total_elements = batch_size * signal_length;
int threads_per_block = 256;
int blocks_per_grid = (total_elements + threads_per_block - 1) / threads_per_block;
fftshift_1d_kernel << <blocks_per_grid, threads_per_block >> > (d_data, batch_size, signal_length);
// 错误检查
PrintLasterError("FFTShift1D");
cudaDeviceSynchronize();
}
extern __global__ void CUDA_D_sin(double* y, double* X, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {
@ -516,11 +555,11 @@ extern "C" void CUDAIFFT(cuComplex* inArr, cuComplex* outArr, long InRowCount,
int rank = 1;
int n[] = { InColCount }; // 每个IFFT处理freqcount点
int inembed[] = { InColCount };
int onembed[] = { InColCount };
int onembed[] = { outColCount };
int istride = 1;
int ostride = 1;
int idist = InColCount; // 输入批次间距
int odist = InColCount; // 输出批次间距
int odist = outColCount; // 输出批次间距
int batch = InRowCount; // 批处理数量
result = cufftPlanMany(&plan, rank, n,
@ -536,7 +575,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;

2
GPUBaseLib/GPUTool/GPUTool.cuh
View File

@ -108,6 +108,6 @@ extern "C" GPUBASELIBAPI void CUDAIFFTScale(cuComplex* inArr, cuComplex* outArr,
extern "C" GPUBASELIBAPI void CUDAIFFT(cuComplex* inArr, cuComplex* outArr, long InRowCount, long InColCount, long outColCount);
extern "C" GPUBASELIBAPI void FFTShift1D(cuComplex* d_data, int batch_size, int signal_length);
#endif
#endif

8
RasterMainWidgetGUI/RasterWidgetMessageShow.cpp
View File

@ -1,5 +1,6 @@
#include "RasterWidgetMessageShow.h"
#include <stdio.h>
#include <iostream>
namespace RasterMessageShow {
RasterWidgetMessageShow* RasterWidgetMessageShow::_instance = nullptr;
@ -27,9 +28,10 @@ namespace RasterMessageShow {
{
if (nullptr != this->textBrowserMessage) {
this->textBrowserMessage->append(Message);
this->textBrowserMessage->moveCursor(QTextCursor::MoveOperation::End);
this->textBrowserMessage->repaint();
std::cout << Message.toLocal8Bit().constData() << std::endl;
}
else {}
}

255
Toolbox/SimulationSARTool/SimulationSAR/GPUBPTool.cu Normal file
View File

@ -0,0 +1,255 @@
#include <iostream>
#include <memory>
#include <cmath>
#include <complex>
#include <device_launch_parameters.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cuComplex.h>
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#include "GPUBPTool.cuh"
#include <cmath>
#include <stdio.h>
#include <cassert>
#define EPSILON 1e-12
#define MAX_ITER 50
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
__device__ __host__ double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees ) {
// 计算点积
double dotProduct = a.x * b.x + a.y * b.y + a.z * b.z;
// 计算模长
double magA = std::sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
double magB = std::sqrt(b.x * b.x + b.y * b.y + b.z * b.z);
// 处理零向量异常
if (magA == 0.0 || magB == 0.0) {
return NAN;
}
double cosTheta = dotProduct / (magA * magB);
cosTheta = cosTheta < -1 ? -1 : cosTheta>1 ? 1 : cosTheta; // 截断到[-1, 1]
double angleRad = std::acos(cosTheta);
return returnDegrees ? angleRad * 180.0 / M_PI : angleRad;
}
// 向量运算
__device__ __host__ Vector3 vec_sub(Vector3 a, Vector3 b) {
return { a.x - b.x, a.y - b.y, a.z - b.z };
}
__device__ __host__ double vec_dot(Vector3 a, Vector3 b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
__device__ __host__ Vector3 vec_cross(Vector3 a, Vector3 b) {
return { a.y * b.z - a.z * b.y,
a.z * b.x - a.x * b.z,
a.x * b.y - a.y * b.x };
}
__device__ __host__ Vector3 vec_normalize(Vector3 v) {
double len = sqrt(vec_dot(v, v));
return (len > 1e-12) ? Vector3{ v.x / len, v.y / len, v.z / len } : v;
}
// 标准正交基底坐标计算
__device__ __host__ Vector3 coordinates_orthonormal_basis(Vector3& A,
Vector3& e1,
Vector3& e2,
Vector3& e3) {
//// 验证基底正交性和单位长度（容差1e-10）
//const double tolerance = 1e-10;
//assert(fabs(dot(e1, e2)) < tolerance);
//assert(fabs(dot(e1, e3)) < tolerance);
//assert(fabs(dot(e2, e3)) < tolerance);
//assert(fabs(norm(e1) - 1.0) < tolerance);
//assert(fabs(norm(e2) - 1.0) < tolerance);
//assert(fabs(norm(e3) - 1.0) < tolerance);
// 计算投影坐标
return Vector3{
vec_dot(A, e1),
vec_dot(A, e2),
vec_dot(A, e3)
};
}
// 计算视线交点T
extern __device__ __host__ Vector3 compute_T(Vector3 S, Vector3 ray, double H) {
Vector3 dir = vec_normalize(ray);
double a_h = WGS84_A + H;
double A = (dir.x * dir.x + dir.y * dir.y) / (a_h * a_h) + dir.z * dir.z / (WGS84_B * WGS84_B); // A > 0
double B = 2.0 * (S.x * dir.x / (a_h * a_h) + S.y * dir.y / (a_h * a_h) + S.z * dir.z / (WGS84_B * WGS84_B));
double C = (S.x * S.x + S.y * S.y) / (a_h * a_h) + S.z * S.z / (WGS84_B * WGS84_B) - 1.0;
double disc = B * B - 4 * A * C;
if (disc < 0) return Vector3{ NAN, NAN, NAN };
double sqrt_d = sqrt(disc);
double t = fmin((-B - sqrt_d) / (2 * A), (-B + sqrt_d) / (2 * A));// 取最小值
return (t > 1e-6) ? Vector3{ S.x + dir.x * t, S.y + dir.y * t, S.z + dir.z * t }
: Vector3{ NAN, NAN, NAN };
}
// 主计算函数A
extern __device__ __host__ Vector3 compute_P(Vector3 S, Vector3 T, double R, double H) {
Vector3 ex, ey, ez; // 平面基函数
Vector3 ST = vec_normalize(vec_sub(T, S));// S->T
Vector3 SO = vec_normalize(vec_sub(Vector3{ 0, 0, 0 }, S)); // S->O
Vector3 st1 = vec_sub(T, S);
double R0 = sqrt(st1.x * st1.x + st1.y * st1.y + st1.z * st1.z);
// S (Z .) --------Y
// |\
// | \
// | \
// X \
// | -> T
// | /
// | /
// | /
// |/
// O
ez = vec_normalize(vec_cross(SO, ST)); // Z 轴
ey = vec_normalize(vec_cross(ez, SO)); // Y 轴 与 ST 同向 --这个结论在星地几何约束，便是显然的；
ex = vec_normalize(SO); //X轴
// 大致的理论推导
// 这里考虑 成像几何，所以点 P 一定在 ex-0-ey 平面上，所以t3=0;
// 定义 SP 的向量与 ex的夹角为 theta , 目标长度为 t
// t1=t*cos(Q);
// t2=t*sin(Q);
// h=(a+H)
// Xp=Sx+t1*ex.x+t2*ey.x +t3*ez.x; //因为 t3=0；
// Yp=Sy+t1*ex.y+t2*ey.y +t3*ez.y;
// Zp=Sz+t1*ex.z+t2*ey.z +t3*ez.z;
// Xp^2+Yp^2 Zp^2 Xp^2+Yp^2 Zp^2
// ---------- + ------- = 1 ==> ---------- + ------- = 1
// (a+H)^2 b^2 h^2 b^2
double h2 = (WGS84_A + H) * (WGS84_A + H);
double b2 = WGS84_B * WGS84_B;
double R2 = R * R;
double A = R2 * ((ex.x * ex.x + ex.y * ex.y) / h2 + (ex.z * ex.z) / b2);
double B = R2 * ((ex.x * ey.x + ex.y * ey.y) / h2 + (ex.z * ey.z) / b2) * 2;
double C = R2 * ((ey.x * ey.x + ey.y * ey.y) / h2 + (ey.z * ey.z) / b2);
double D = 1 - ((S.x * S.x + S.y * S.y) / h2 + (S.z * S.z) / b2);
double E = 2 * R * ((S.x * ex.x + S.y * ex.y) / h2 + (S.z * ex.z) / b2);
double F = 2 * R * ((S.x * ey.x + S.y * ey.y) / h2 + (S.z * ey.z) / b2);
// f(Q)=Acos^2(Q)+Bsin(Q)cos(Q)+Csin^2(Q)+E*cos(Q)+F*sin(Q)-D
// f'(Q)=(C−A)sin(2Q)+2Bcos(2Q)-Esin(Q)+Fcos(Q)
// 牛顿迭代
// f(Q)
// Q(t+1)=Q - -----------
// f'(Q)
// 求解初始值
double Q0 = angleBetweenVectors(SO, ST, false);
double dQ = 0;
double fQ = 0;
double dfQ = 0;
double Q = R < R0 ? Q0 - 1e-3 : Q0 + 1e-3;
// 牛顿迭代法
for (int iter = 0; iter < MAX_ITER * 10; ++iter) {
fQ = A * cos(Q) * cos(Q) + B * sin(Q) * cos(Q) + C * sin(Q) * sin(Q) + E * cos(Q) + F * sin(Q) - D;
dfQ = (C - A) * sin(2 * Q) + B * cos(2 * Q) - E * sin(Q) + F * cos(Q);
dQ = fQ / dfQ;
if (abs(dQ) < 1e-8) {
break;
}
else {
dQ = abs(dQ) < d2r * 2 ? dQ : abs(dQ) / dQ * d2r;
Q = Q - dQ;
}
}
double t1 = R * cos(Q);
double t2 = R * sin(Q);
Vector3 P = {
S.x + t1 * ex.x + t2 * ey.x, //因为 t3=0；
S.y + t1 * ex.y + t2 * ey.y,
S.z + t1 * ex.z + t2 * ey.z,
};
// 椭球验证
double check = (P.x * P.x + P.y * P.y) / ((WGS84_A + H) * (WGS84_A + H))
+ P.z * P.z / (WGS84_B * WGS84_B);
if (isnan(Q) || isinf(Q) || fabs(check - 1.0) > 1e-6) {
return Vector3{ NAN,NAN,NAN };
printf("check value =%f\n", fabs(check - 1.0));
}
return P;
}
//// 参数校验与主函数
//int main() {
// Vector3 S = { -2.8e6, -4.2e6, 3.5e6 }; // 卫星位置 (m)
// Vector3 ray = { 0.6, 0.4, -0.7 }; // 视线方向
// double H = 500.0; // 平均高程
// double R = 1000.0; // 目标距离
//
// // 参数校验
// if (R <= 0 || H < -WGS84_A * 0.1 || H > WGS84_A * 0.1) {
// printf("参数错误：\n H范围：±%.1f km\n R必须>0\n", WGS84_A * 0.1 / 1000);
// return 1;
// }
//
// // Step 1: 计算交点T
// Vector3 T = compute_T(S, ray, H);
// if (isnan(T.x)) {
// printf("错误：视线未与椭球相交\n");
// return 1;
// }
//
// // Step 2: 计算目标点P
// Vector3 P = compute_P(S, T, R, H);
//
// if (!isnan(P.x)) {
// printf("计算结果：\n");
// printf("P = (%.3f, %.3f, %.3f) m\n", P.x, P.y, P.z);
//
// // 验证距离
// Vector3 SP = vec_sub(P, S);
// double dist = sqrt(vec_dot(SP, SP));
// printf("实际距离：%.3f m (期望：%.1f m)\n", dist, R);
//
// // 验证椭球
// double check = (P.x * P.x + P.y * P.y) / ((WGS84_A + H) * (WGS84_A + H))
// + P.z * P.z / (WGS84_B * WGS84_B);
// printf("椭球验证：%.6f (期望：1.0)\n", check);
//
// // 验证最近距离
// Vector3 PT = vec_sub(P, T);
// printf("到T点距离：%.3f m\n", sqrt(vec_dot(PT, PT)));
// }
// else {
// printf("未找到有效解\n");
// }
//
// return 0;
//}
//

22
Toolbox/SimulationSARTool/SimulationSAR/GPUBPTool.cuh Normal file
View File

@ -0,0 +1,22 @@
#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"
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 );
//
//

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

@ -1,4 +1,4 @@

#include <iostream>
#include <memory>
@ -12,10 +12,162 @@
#include "BaseConstVariable.h"
#include "GPUTool.cuh"
#include "GPUTBPImage.cuh"
#include "GPUBPTool.cuh"
#ifdef __CUDANVCC___
#define EPSILON 1e-12
#define MAX_ITER 50
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
__global__ void kernel_TimeBPImageGridNet(double* antPx, double* antPy, double* antPz,
double* antDirx, double* antDiry, double* antDirz,
double* imgx, double* imgy, double* imgz,
long prfcount, long freqpoints, double meanH,
double Rnear, double dx, double RefRange) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
long pixelcount = prfcount * freqpoints;
long prfid = idx / freqpoints;
long Rid = idx % freqpoints;
if (idx < pixelcount) {
// 计算坐标
Vector3 S = { antPx[prfid], antPy[prfid], antPz[prfid] }; // 卫星位置 (m)
Vector3 ray = { antDirx[prfid], antDiry[prfid], antDirz[prfid] }; // 视线方向
double H = meanH; // 平均高程
double R = Rnear + dx * Rid; // 目标距离
// 参数校验
if (R <= 0 || H < -WGS84_A * 0.1 || H > WGS84_A * 0.1) {
//printf("参数错误：\n H范围：±%.1f km\n R必须>0\n", WGS84_A * 0.1 / 1000);
imgx[idx] = NAN;
imgy[idx] = NAN;
imgz[idx] = NAN;
//printf("idx=%d;prfid=%d;Rid=%d;S=[%f , %f ,%f ];ray=[%f ,%f ,%f ];H=%f;R=%f,imgP=[%f ,%f , %f ];Rextend\n",
// idx, prfid, Rid, S.x, S.y, S.z, ray.x, ray.y, ray.z, H, R,imgx[idx],imgy[idx],imgz[idx]);
// 参数校验
return;
}
// Step 1: 计算交点T
Vector3 T = compute_T(S, ray, H);
if (isnan(T.x)) {
imgx[idx] = NAN;
imgy[idx] = NAN;
imgz[idx] = NAN;
//printf("idx=%d;prfid=%d;Rid=%d;Tnan\n",
// idx, prfid, Rid, S.x, S.y, S.z, ray.x, ray.y, ray.z, H, R,T.x,T.y,T.z, imgx[idx], imgy[idx], imgz[idx]);
return;
}
// Step 2: 计算目标点P
Vector3 P;// = compute_P(S, T, R, H);
{ // 计算P
Vector3 ex, ey, ez; // 平面基函数
Vector3 ST = vec_normalize(vec_sub(T, S));// S->T
Vector3 SO = vec_normalize(vec_sub(Vector3{ 0, 0, 0 }, S)); // S->O
Vector3 st1 = vec_sub(T, S);
double R0 = sqrt(st1.x * st1.x + st1.y * st1.y + st1.z * st1.z);
ez = vec_normalize(vec_cross(SO, ST)); // Z 轴
ey = vec_normalize(vec_cross(ez, SO)); // Y 轴 与 ST 同向 --这个结论在星地几何约束，便是显然的；
ex = vec_normalize(SO); //X轴
double h2 = (WGS84_A + H) * (WGS84_A + H);
double b2 = WGS84_B * WGS84_B;
double R2 = R * R;
double A = R2 * ((ex.x * ex.x + ex.y * ex.y) / h2 + (ex.z * ex.z) / b2);
double B = R2 * ((ex.x * ey.x + ex.y * ey.y) / h2 + (ex.z * ey.z) / b2) * 2;
double C = R2 * ((ey.x * ey.x + ey.y * ey.y) / h2 + (ey.z * ey.z) / b2);
double D = 1 - ((S.x * S.x + S.y * S.y) / h2 + (S.z * S.z) / b2);
double E = 2*R * ((S.x * ex.x + S.y * ex.y) / h2 + (S.z * ex.z) / b2);
double F = 2*R * ((S.x * ey.x + S.y * ey.y) / h2 + (S.z * ey.z) / b2);
double Q0 = angleBetweenVectors(SO, ST, false);
double dQ = 0;
double fQ = 0;
double dfQ = 0;
double Q = R < R0 ? Q0 - 1e-3 : Q0 + 1e-3;
//printf("A=%f;B=%f;C=%f;D=%f;E=%f;F=%f;Q=%f;\
// S=[%f , %f ,%f ];\
// T=[%f , %f ,%f ];\
// ex=[%f , %f ,%f ];\
// ey=[%f , %f ,%f ];\
// ez=[%f , %f ,%f ];\
//ray=[%f ,%f ,%f ];\
//H=%f;R=%f;;\n",A,B,C,D,E,F,Q,
// S.x,S.y,S.z,
// T.x,T.y,T.z ,
// ex.x,ex.y,ex.z,
// ey.x,ey.y,ey.z,
// ez.x,ez.y,ez.z,
// ray.x, ray.y, ray.z,
// H, R);
// return;
// 牛顿迭代法
for (int iter = 0; iter < MAX_ITER * 10; ++iter) {
fQ = A * cos(Q) * cos(Q) + B * sin(Q) * cos(Q) + C * sin(Q) * sin(Q) + E * cos(Q) + F * sin(Q) - D;
dfQ = (C - A) * sin(2 * Q) + B * cos(2 * Q) - E * sin(Q) + F * cos(Q);
dQ = fQ / dfQ;
if (abs(dQ) < 1e-8) {
//printf("iter=%d;check Q0=%f;Q=%f;dQ=%f;fQ=%f;dfQ=%f;break\n", iter, Q0, Q, dQ, fQ, dfQ);
break;
}
else {
dQ = (abs(dQ) < d2r * 3) ? dQ :( abs(dQ) / dQ * d2r* 3);
Q = Q - dQ;
//printf("iter=%d;check Q0=%f;Q=%f;dQ=%f;fQ=%f;dfQ=%f;\n", iter, Q0, Q, dQ, fQ, dfQ);
}
}
//printf("check Q0=%f;Q=%f;\n", Q0, Q);
double t1 = R * cos(Q);
double t2 = R * sin(Q);
P = Vector3{
S.x + t1 * ex.x + t2 * ey.x, //因为 t3=0；
S.y + t1 * ex.y + t2 * ey.y,
S.z + t1 * ex.z + t2 * ey.z,
};
double check = (P.x * P.x + P.y * P.y) / ((WGS84_A + H) * (WGS84_A + H))
+ P.z * P.z / (WGS84_B * WGS84_B);
if (isnan(Q) || isinf(Q) || fabs(check - 1.0) > 1e-6) {
P = Vector3{ NAN,NAN,NAN };
}
}
double Rt = sqrt(pow(S.x - T.x, 2) + pow(S.y - T.y, 2) + pow(S.z - T.z, 2));
double Rp = sqrt(pow(S.x - P.x, 2) + pow(S.y - P.y, 2) + pow(S.z - P.z, 2));
double Rop = sqrt(pow( P.x, 2) + pow( P.y, 2) + pow( P.z, 2));
if (!isnan(P.x)&&( Rop>WGS84_A*0.3)&&(Rop<WGS84_A*3)) {
imgx[idx] = P.x;
imgy[idx] = P.y;
imgz[idx] = P.z;
//printf("idx=%d; S=[%f , %f ,%f ]; H=%f;R=%f;RP=%f;Rr=%f;imgT=[%f ,%f ,%f ];imgP=[%f ,%f , %f ]; \n",
// idx, S.x, S.y, S.z, H, R, Rp, Rt,T.x, T.y, T.z, P.x, P.y, P.z);
}
else {
imgx[idx] = NAN;
imgy[idx] = NAN;
imgz[idx] = NAN;
printf("idx=%d; S=[%f , %f ,%f ]; H=%f;R=%f;RP=%f;Rr=%f;imgT=[%f ,%f ,%f ];imgP=[%f ,%f , %f ]; ERROR\n",
idx, S.x, S.y, S.z, H, R, Rp, Rt, T.x, T.y, T.z, P.x, P.y, P.z);
}
}
}
__global__ void kernel_pixelTimeBP(
@ -56,7 +208,7 @@ __global__ void kernel_pixelTimeBP(
Pz = antPz[pid];
PR = sqrt((Px - Tx) * (Px - Tx) + (Py - Ty) * (Py - Ty) + (Pz - Tz) * (Pz - Tz));
Rid = (PR - Rnear) / dx; // 行数
Rid = (PR - Rnear) / dx; // 行数
if (Rid<0 || Rid>maxPointIdx) {
continue;
@ -70,14 +222,14 @@ __global__ void kernel_pixelTimeBP(
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
@ -85,7 +237,7 @@ __global__ void kernel_pixelTimeBP(
phiCorr.x = cos(phi);
phiCorr.y = sin(phi);
s = cuCmulf(s, phiCorr); // 校正后
s = cuCmulf(s, phiCorr); // 校正后
imgArr[idx] = cuCaddf(imgArr[idx], s);
}
@ -97,31 +249,54 @@ __global__ void kernel_pixelTimeBP(
extern "C" 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
)
{
long pixelcount = imH * imW;
int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
extern "C" {
void TIMEBPCreateImageGrid(double* antPx, double* antPy, double* antPz,
double* antDirx, double* antDiry, double* antDirz,
double* imgx, double* imgy, double* imgz,
long prfcount, long freqpoints, double meanH,
double Rnear, double dx, double RefRange)
{
long pixelcount = prfcount * freqpoints;
int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
kernel_pixelTimeBP << <grid_size, BLOCK_SIZE >> > (
antPx, antPy, antPz,
imgx, imgy, imgz,
TimeEchoArr, prfcount, pointCount,
imgArr, imH, imW,
startLamda, Rnear, dx, RefRange
);
kernel_TimeBPImageGridNet << <grid_size, BLOCK_SIZE >> > (
antPx, antPy, antPz,
antDirx, antDiry, antDirz,
imgx, imgy, imgz,
prfcount, freqpoints, meanH,
Rnear, dx, RefRange);
PrintLasterError("TIMEBPCreateImageGrid");
cudaDeviceSynchronize();
}
PrintLasterError("TimeBPImage");
cudaDeviceSynchronize();
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
)
{
long pixelcount = imH * imW;
int grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
kernel_pixelTimeBP << <1, 1 >> > (
antPx, antPy, antPz,
imgx, imgy, imgz,
TimeEchoArr, prfcount, pointCount,
imgArr, imH, imW,
startLamda, Rnear, dx, RefRange
);
PrintLasterError("TimeBPImage");
cudaDeviceSynchronize();
}
}
#endif

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

@ -11,6 +11,15 @@
extern "C" void TIMEBPCreateImageGrid(
double* antPx,double* antPy,double* antPz, // ÎÀÐÇ×ø±ê S
double* antDirx,double* antDiry,double* antDirz, //
double* imgx,double* imgy,double* imgz,
long prfcount,long freqpoints,double meanH,
double Rnear,double dx,double RefRange
);
extern "C" void TimeBPImage(
double* antPx, double* antPy, double* antPz,

4
Toolbox/SimulationSARTool/SimulationSAR/QSimulationBPImage.ui
View File

@ -45,7 +45,7 @@
</size>
</property>
<property name="text">
<string>D:\Programme\vs2022\RasterMergeTest\LAMPCAE_SCANE-all-scane\GF3_Simulation.xml</string>
<string>D:\Programme\vs2022\RasterMergeTest\LAMPCAE_SCANE\GF3_Simulation.xml</string>
</property>
</widget>
</item>
@ -96,7 +96,7 @@
</size>
</property>
<property name="text">
<string>D:\Programme\vs2022\RasterMergeTest\LAMPCAE_SCANE-all-scane\BPImage\GF3BPImage</string>
<string>D:\Programme\vs2022\RasterMergeTest\LAMPCAE_SCANE\BPImage\GF3BPImage</string>
</property>
</widget>
</item>

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

@ -12,6 +12,13 @@
void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZPath)
{
long bandwidth = echoL0ds->getBandwidth();
double Rnear = echoL0ds->getNearRange();
double Rfar = echoL0ds->getFarRange();
double refRange = echoL0ds->getRefPhaseRange();
double dx = LIGHTSPEED / 2.0 / bandwidth; // c/2b
// 创建坐标系统
long prfcount = echoL0ds->getPluseCount();
long freqcount = echoL0ds->getPlusePoints();
@ -22,93 +29,122 @@ void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds, QString outPixelXYZ
gt(1, 0) = 0;
gt(1, 1) = 0;
gt(1, 2) = 1;
gdalImage xyzRaster = CreategdalImage(outPixelXYZPath, prfcount, freqcount, 3, gt, QString(""), false, true,true);
gdalImage xyzRaster = CreategdalImage(outPixelXYZPath, prfcount, freqcount, 3, gt, QString(""), false, true,true,GDT_Float64);
std::shared_ptr<double> antpos = echoL0ds->getAntPos();
double dx = (echoL0ds->getFarRange()-echoL0ds->getNearRange())/(echoL0ds->getPlusePoints()-1);
double Rnear = echoL0ds->getNearRange();
dx = (echoL0ds->getFarRange()-echoL0ds->getNearRange())/(echoL0ds->getPlusePoints()-1);
Rnear = echoL0ds->getNearRange();
double Rref = echoL0ds->getRefPhaseRange();
double centerInc = echoL0ds->getCenterAngle()*d2r;
long echocol = 1073741824 / 8 / 4 / prfcount*8;
long echocol = Memory1GB * 1.0 / 8 / 4 / prfcount * 6;
qDebug() << "echocol:\t " << echocol ;
echocol = echocol < 3000 ? 3000 : echocol;
long startcolidx = 0;
for (startcolidx = 0; startcolidx < freqcount; startcolidx = startcolidx + echocol) {
echocol = echocol < 1 ? 1: echocol;
std::shared_ptr<double> Pxs((double*)mallocCUDAHost(sizeof(double)*prfcount), FreeCUDAHost);
std::shared_ptr<double> Pys((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
std::shared_ptr<double> Pzs((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
std::shared_ptr<double> AntDirectX((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
std::shared_ptr<double> AntDirectY((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
std::shared_ptr<double> AntDirectZ((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
{
std::shared_ptr<double> antpos = echoL0ds->getAntPos();
double time = 0;
double Px = 0;
double Py = 0;
double Pz = 0;
for (long i = 0; i < prfcount; i++) {
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
AntDirectY.get()[i] = antpos.get()[i * 19 + 14];
AntDirectZ.get()[i] = antpos.get()[i * 19 + 15];
double NormAnt = std::sqrt(AntDirectX.get()[i] * AntDirectX.get()[i] +
AntDirectY.get()[i] * AntDirectY.get()[i] +
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;// 归一化
}
antpos.reset();
}
std::shared_ptr<double> d_Pxs((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
std::shared_ptr<double> d_Pys((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
std::shared_ptr<double> d_Pzs((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
std::shared_ptr<double> d_AntDirectX((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
std::shared_ptr<double> d_AntDirectY((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
std::shared_ptr<double> d_AntDirectZ((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
HostToDevice(Pxs.get(), d_Pxs.get(), sizeof(double) * prfcount);
HostToDevice(Pys.get(), d_Pys.get(), sizeof(double) * prfcount);
HostToDevice(Pzs.get(), d_Pzs.get(), sizeof(double) * prfcount);
HostToDevice(AntDirectX.get(), d_AntDirectX.get(), sizeof(double) * prfcount);
HostToDevice(AntDirectY.get(), d_AntDirectY.get(), sizeof(double) * prfcount);
HostToDevice(AntDirectZ.get(), d_AntDirectZ.get(), sizeof(double) * prfcount);
for (long startcolidx = 0; startcolidx < freqcount; startcolidx = startcolidx + echocol) {
long tempechocol = echocol;
if (startcolidx + tempechocol >= freqcount) {
tempechocol = freqcount - startcolidx;
}
qDebug() << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz") << "] imgxyz :\t" << startcolidx << "\t-\t" << startcolidx + tempechocol << " / " << freqcount ;
std::shared_ptr<double> demx = readDataArr<double>(xyzRaster, 0, startcolidx, prfcount, tempechocol, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> demy = readDataArr<double>(xyzRaster, 0, startcolidx, prfcount, tempechocol, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> demz = readDataArr<double>(xyzRaster, 0, startcolidx, prfcount, tempechocol, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
Eigen::MatrixXd demx = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 1);
Eigen::MatrixXd demy = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 2);
Eigen::MatrixXd demz = xyzRaster.getData(0, startcolidx, prfcount, tempechocol, 3);
std::shared_ptr<double> h_demx((double*)mallocCUDAHost(sizeof(double) * prfcount*tempechocol), FreeCUDAHost);
std::shared_ptr<double> h_demy((double*)mallocCUDAHost(sizeof(double) * prfcount*tempechocol), FreeCUDAHost);
std::shared_ptr<double> h_demz((double*)mallocCUDAHost(sizeof(double) * prfcount*tempechocol), FreeCUDAHost);
#pragma omp parallel for
for (long i = 0; i < prfcount; i++) {
double Px = 0;
double Py = 0;
double Pz = 0;
double AntDirectX = 0;
double AntDirectY = 0;
double AntDirectZ = 0;
double R = 0;
double NormAnt = 0;
Px = antpos.get()[i * 19 + 1]; // 卫星坐标
Py = antpos.get()[i * 19 + 2];
Pz = antpos.get()[i * 19 + 3];
AntDirectX = antpos.get()[i * 19 + 13];// zero doppler
AntDirectY = antpos.get()[i * 19 + 14];
AntDirectZ = antpos.get()[i * 19 + 15];
NormAnt = std::sqrt(AntDirectX * AntDirectX + AntDirectY * AntDirectY + AntDirectZ * AntDirectZ);
AntDirectX = AntDirectX / NormAnt;
AntDirectY = AntDirectY / NormAnt;
AntDirectZ = AntDirectZ / NormAnt;// 归一化
// 计算中心参考点
double centerX = Px + Rref * AntDirectX; // T1
double centerY = Py + Rref * AntDirectY;
double centerZ = Pz + Rref * AntDirectZ;
double satH = Rref * std::cos(centerInc); // 卫星高
double PR = sqrt(Px * Px + Py * Py + Pz * Pz);
double satR = PR - satH;
double sPx = satR / PR * Px;
double sPy = satR / PR * Py;
double sPz = satR / PR * Pz;
double dTSx = sPx - centerX;
double dTSy = sPy - centerY;
double dTSz = sPz - centerZ;
double dTSR = sqrt(dTSx * dTSx + dTSy * dTSy + dTSz * dTSz);
dTSx=dTSx/dTSR;
dTSy=dTSy/dTSR;
dTSz=dTSz/dTSR;
for (long j = 0; j < tempechocol; j++) {
R = (j + startcolidx)*dx + Rnear;
double dRp = (Rref - R) / sin(centerInc); // -- 0 +++
demx(i,j) = centerX + dTSx * dRp;
demy(i,j) = centerY + dTSy * dRp;
demz(i,j) = centerZ + dTSz * dRp;
#pragma omp parallel for
for (long ii = 0; ii < prfcount; ii++) {
for (long jj = 0; jj < tempechocol; jj++) {
h_demx.get()[ii*tempechocol+jj]=demx.get()[ii*tempechocol+jj];
h_demy.get()[ii*tempechocol+jj]=demy.get()[ii*tempechocol+jj];
h_demz.get()[ii*tempechocol+jj]=demz.get()[ii*tempechocol+jj];
}
}
xyzRaster.saveImage(demx, 0, startcolidx, 1);
xyzRaster.saveImage(demy, 0, startcolidx, 2);
xyzRaster.saveImage(demz, 0, startcolidx, 3);
std::shared_ptr<double> d_demx((double*)mallocCUDADevice(sizeof(double) * prfcount * tempechocol), FreeCUDADevice);
std::shared_ptr<double> d_demy((double*)mallocCUDADevice(sizeof(double) * prfcount * tempechocol), FreeCUDADevice);
std::shared_ptr<double> d_demz((double*)mallocCUDADevice(sizeof(double) * prfcount * tempechocol), FreeCUDADevice);
HostToDevice(h_demx.get(), d_demx.get(), sizeof(double) * prfcount * tempechocol);
HostToDevice(h_demy.get(), d_demy.get(), sizeof(double) * prfcount * tempechocol);
HostToDevice(h_demz.get(), d_demz.get(), sizeof(double) * prfcount * tempechocol);
TIMEBPCreateImageGrid(
d_Pxs.get(), d_Pys.get(), d_Pzs.get(),
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
);
DeviceToHost(h_demx.get(), d_demx.get(), sizeof(double) * prfcount * tempechocol);
DeviceToHost(h_demy.get(), d_demy.get(), sizeof(double) * prfcount * tempechocol);
DeviceToHost(h_demz.get(), d_demz.get(), sizeof(double) * prfcount * tempechocol);
#pragma omp parallel for
for (long ii = 0; ii < prfcount; ii++) {
for (long jj = 0; jj < tempechocol; jj++) {
demx.get()[ii * tempechocol + jj]=h_demx.get()[ii * tempechocol + jj] ;
demy.get()[ii * tempechocol + jj]=h_demy.get()[ii * tempechocol + jj] ;
demz.get()[ii * tempechocol + jj]=h_demz.get()[ii * tempechocol + jj] ;
}
}
xyzRaster.saveImage(demx, 0, startcolidx,prfcount,tempechocol, 1);
xyzRaster.saveImage(demy, 0, startcolidx,prfcount,tempechocol, 2);
xyzRaster.saveImage(demz, 0, startcolidx,prfcount,tempechocol, 3);
}
}
@ -126,6 +162,7 @@ void TBPImageProcess(QString echofile, QString outImageFolder, QString imagePlan
imagL1->setFarRange(echoL0ds->getFarRange());
imagL1->setFs(echoL0ds->getFs());
imagL1->setLookSide(echoL0ds->getLookSide());
imagL1->OpenOrNew(outImageFolder, echoL0ds->getSimulationTaskName(), echoL0ds->getPluseCount(), echoL0ds->getPlusePoints());
@ -171,11 +208,20 @@ std::shared_ptr<SARSimulationImageL1Dataset> TBPImageAlgCls::getImageL0()
ErrorCode TBPImageAlgCls::Process(long num_thread)
{
qDebug() << u8"开始成像";
qDebug() << u8"创建成像平面的XYZ";
QString outRasterXYZ = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_xyz.bin");
CreatePixelXYZ(this->L0ds, outRasterXYZ);
this->outRasterXYZPath = outRasterXYZ;
qDebug() << u8"频域回波-> 时域回波";
this->TimeEchoDataPath = JoinPath(this->L1ds->getoutFolderPath(), this->L0ds->getSimulationTaskName() + "_Timeecho.bin");
this->EchoFreqToTime();
// 初始化Raster
qDebug() << u8"初始化影像";
long imageheight = this->L1ds->getrowCount();
@ -202,14 +248,18 @@ ErrorCode TBPImageAlgCls::Process(long num_thread)
}
this->L1ds->saveImageRaster(imageRaster, startline,templine);
}
qDebug() << u8"开始成像";
if (GPURUN) {
return this->ProcessGPU();
}
else {
QMessageBox::information(nullptr,u8"提示",u8"目前只支持显卡");
return ErrorCode::FAIL;
}
qDebug() << u8"频域回波-> 时域回波 结束";
//if (GPURUN) {
// return this->ProcessGPU();
//}
//else {
// QMessageBox::information(nullptr,u8"提示",u8"目前只支持显卡");
// return ErrorCode::FAIL;
//}
return ErrorCode::SUCCESS;
}
ErrorCode TBPImageAlgCls::ProcessGPU()
@ -254,9 +304,9 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
qDebug() << "deltaF:\t" << deltaF;
std::shared_ptr<double> Pxs (new double[this->L0ds->getPluseCount()]);
std::shared_ptr<double> Pys (new double[this->L0ds->getPluseCount()]);
std::shared_ptr<double> Pzs (new double[this->L0ds->getPluseCount()]);
std::shared_ptr<double> Pxs (new double[this->L0ds->getPluseCount()],delArrPtr);
std::shared_ptr<double> Pys (new double[this->L0ds->getPluseCount()],delArrPtr);
std::shared_ptr<double> Pzs (new double[this->L0ds->getPluseCount()],delArrPtr);
{
std::shared_ptr<double> antpos = this->L0ds->getAntPos();
@ -307,6 +357,7 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
imageBlockline = imageBlockline < 1 ? 1 : imageBlockline;
gdalImage imageXYZ(this->outRasterXYZPath);
gdalImageComplex imagetimeimg(this->TimeEchoDataPath);
long startimgrowid = 0;
for (startimgrowid = 0; startimgrowid < rowCount; startimgrowid = startimgrowid + imageBlockline) {
@ -324,16 +375,18 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
std::shared_ptr<std::complex<double>> imgArr = this->L1ds->getImageRaster(startimgrowid, tempimgBlockline);
// 获取回波
long startechoid = 0;
long iffeechoLen = PlusePoints;
for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) {
long tempechoBlockline = echoBlockline;
if (startechoid + tempechoBlockline >= PRFCount) {
tempechoBlockline = PRFCount - startechoid;
}
std::shared_ptr<std::complex<double>> echoArr = this->L0ds->getEchoArr(startechoid, tempechoBlockline);
std::shared_ptr<double> antpx(new double[tempechoBlockline*PlusePoints],delArrPtr);
std::shared_ptr<double> antpy(new double[tempechoBlockline* PlusePoints], delArrPtr);
std::shared_ptr<double> antpz(new double[tempechoBlockline* PlusePoints], delArrPtr);
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];
@ -361,9 +414,81 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
return ErrorCode::SUCCESS;
}
void TBPImageAlgCls::EchoFreqToTime( )
{
// 读取数据
long PRFCount = this->L0ds->getPluseCount();
long inColCount = this->L0ds->getPlusePoints();
long outColCount = inColCount;// nextpow2(inColCount);
this->TimeEchoRowCount = PRFCount;
this->TimeEchoColCount = outColCount;
qDebug() << "IFFT : " << this->TimeEchoDataPath;
qDebug() << "PRF Count:\t" << PRFCount;
qDebug() << "inColCount:\t" << inColCount;
qDebug() << "outColCount:\t" << outColCount;
// 创建二进制文件
gdalImageComplex outTimeEchoImg = CreategdalImageComplexNoProj(this->TimeEchoDataPath,this->TimeEchoRowCount,this->TimeEchoColCount,1);
// 分块
long echoBlockline = Memory1GB / 8 / 2 / outColCount * 3; //1GB
echoBlockline = echoBlockline < 1 ? 1 : echoBlockline;
long startechoid = 0;
for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) {
long tempechoBlockline = echoBlockline;
if (startechoid + tempechoBlockline >= PRFCount) {
tempechoBlockline = PRFCount - startechoid;
}
std::shared_ptr<std::complex<double>> echoArr = this->L0ds->getEchoArr(startechoid, tempechoBlockline);
std::shared_ptr<std::complex<double>> IFFTArr = outTimeEchoImg.getDataComplexSharePtr(startechoid, 0, tempechoBlockline, outColCount, 1);
std::shared_ptr<cuComplex> host_echoArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex)* tempechoBlockline * outColCount), FreeCUDAHost);
std::shared_ptr<cuComplex> host_IFFTechoArr((cuComplex*)mallocCUDAHost(sizeof(cuComplex)* tempechoBlockline * outColCount), FreeCUDAHost);
memset(host_echoArr.get(), 0, sizeof(cuComplex) * tempechoBlockline * outColCount);
#pragma omp parallel for
for (long ii = 0; ii < tempechoBlockline ; ii++) {
for (long jj = 0; jj < inColCount; jj++) {
host_echoArr.get()[ii* outColCount +jj] = make_cuComplex(echoArr.get()[ii * inColCount + jj].real(), echoArr.get()[ii * inColCount + jj].imag());
}
}
#pragma omp parallel for
for (long ii = 0; ii < tempechoBlockline * outColCount; ii++) {
host_IFFTechoArr.get()[ii] = make_cuComplex(0, 0);
}
std::shared_ptr<cuComplex> device_echoArr((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * tempechoBlockline * inColCount), FreeCUDADevice);
std::shared_ptr<cuComplex> device_IFFTechoArr((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * tempechoBlockline * outColCount), FreeCUDADevice);
HostToDevice(host_echoArr.get(), device_echoArr.get(), sizeof(cuComplex) * tempechoBlockline * inColCount);
HostToDevice(host_IFFTechoArr.get(), device_IFFTechoArr.get(), sizeof(cuComplex) * tempechoBlockline * outColCount);
CUDAIFFT(device_echoArr.get(), device_IFFTechoArr.get(), tempechoBlockline, outColCount, outColCount);
FFTShift1D(device_IFFTechoArr.get(), tempechoBlockline, outColCount);
DeviceToHost(host_IFFTechoArr.get(), device_IFFTechoArr.get(), sizeof(cuComplex) * tempechoBlockline * outColCount);
#pragma omp parallel for
for (long ii = 0; ii < tempechoBlockline ; ii++) {
for (long jj = 0; jj < outColCount; jj++) {
IFFTArr.get()[ii * outColCount + jj] = std::complex<double>(host_IFFTechoArr.get()[ii * outColCount + jj].x, host_IFFTechoArr.get()[ii * outColCount + jj].y);
//IFFTArr.get()[ii * outColCount + jj] = std::complex<double>(host_echoArr.get()[ii * outColCount + jj].x, host_echoArr.get()[ii * outColCount + jj].y);
}
}
outTimeEchoImg.saveImage(IFFTArr, startechoid, 0, tempechoBlockline, outColCount, 1);
qDebug() << QString(" image block PRF:[%1] \t").arg((startechoid + tempechoBlockline) * 100.0 / PRFCount)
<< startechoid << "\t-\t" << startechoid + tempechoBlockline;
}
return;
}
@ -380,8 +505,8 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
long IFFTPadNum = nextpow2(freqcount);
// 先处理脉冲傅里叶变换
cuComplex* h_echoArr = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * prfcount * freqcount);
cuComplex* d_echoArr = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * prfcount * freqcount);
std::shared_ptr<cuComplex> d_echoArrIFFT((cuComplex*)mallocCUDADevice(sizeof(cuComplex) * prfcount * IFFTPadNum), FreeCUDADevice);
//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++) {
@ -390,11 +515,10 @@ void TBPImageGPUAlg2(std::shared_ptr<double> antPx, std::shared_ptr<double> antP
echoArr.get()[i * freqcount + j].imag());
}
}
HostToDevice(h_echoArr, d_echoArr, sizeof(cuComplex) * prfcount * freqcount);
CUDAIFFT(d_echoArr, d_echoArrIFFT.get(), prfcount, freqcount, IFFTPadNum);
HostToDevice(h_echoArr, d_echoArrIFFT.get(), sizeof(cuComplex) * prfcount * freqcount);
// 结束傅里叶变换
FreeCUDAHost(h_echoArr);
FreeCUDADevice(d_echoArr);
//FreeCUDADevice(d_echoArr);
qDebug() << "IFFT finished!!!";
// 初始化

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

@ -57,6 +57,13 @@ private:
//ErrorCode ProcessCPU(long num_thread);
ErrorCode ProcessGPU();
private://临时成员变量
QString TimeEchoDataPath;
long TimeEchoRowCount;
long TimeEchoColCount;
void EchoFreqToTime( );
};

2
Toolbox/SimulationSARTool/SimulationSARTool.vcxproj
View File

@ -228,6 +228,7 @@
<QtMoc Include="SimulationSAR\QSARLookTableSimualtionGUI.h" />
<QtMoc Include="SimulationSAR\QSimulationBPImage.h" />
<QtMoc Include="SimulationSAR\QSimulationRFPCGUI.h" />
<CudaCompile Include="SimulationSAR\GPUBPTool.cuh" />
<ClInclude Include="SimulationSAR\RFPCProcessCls.h" />
<ClInclude Include="SimulationSAR\SARSatelliteSimulationAbstractCls.h" />
<ClInclude Include="SimulationSAR\SARSimulationTaskSetting.h" />
@ -239,6 +240,7 @@
</ItemGroup>
<ItemGroup>
<CudaCompile Include="PowerSimulationIncoherent\LookTableSimulationComputer.cu" />
<CudaCompile Include="SimulationSAR\GPUBPTool.cu" />
<CudaCompile Include="SimulationSAR\GPURFPC.cu">
<FileType>Document</FileType>
</CudaCompile>

6
Toolbox/SimulationSARTool/SimulationSARTool.vcxproj.filters
View File

@ -172,5 +172,11 @@
<CudaCompile Include="PowerSimulationIncoherent\LookTableSimulationComputer.cuh">
<Filter>PowerSimulationIncoherent</Filter>
</CudaCompile>
<CudaCompile Include="SimulationSAR\GPUBPTool.cu">
<Filter>SimulationSAR</Filter>
</CudaCompile>
<CudaCompile Include="SimulationSAR\GPUBPTool.cuh">
<Filter>SimulationSAR</Filter>
</CudaCompile>
</ItemGroup>
</Project>