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

增加RTPC的GPU工具代码

main
陈增辉 2024-11-29 23:29:45 +08:00
parent 458be41d96
commit ef94aa388a
8 changed files with 737 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
BaseConstVariable.h
View File

@ -3,7 +3,7 @@
#ifndef BASECONSTVARIABLE_H #ifndef BASECONSTVARIABLE_H
#define BASECONSTVARIABLE_H #define BASECONSTVARIABLE_H
#define EIGEN_USE_MKL_ALL #define EIGEN_USE_MKL_ALL
//#define EIGEN_NO_DEBUG #define EIGEN_NO_DEBUG
#define EIGEN_USE_BLAS #define EIGEN_USE_BLAS
@ -12,11 +12,17 @@
//#define DEBUGSHOWDIALOG //#define DEBUGSHOWDIALOG
#define __CUDANVCC___ // 定义CUDA函数
//#include <mkl.h> //#include <mkl.h>
#include <complex> #include <complex>
#include <math.h> #include <math.h>
#include <complex> #include <complex>
#define MATPLOTDRAWIMAGE #define MATPLOTDRAWIMAGE
@ -120,16 +126,19 @@ struct SatelliteAntPos {
double Pz; double Pz;
double Vx; double Vx;
double Vy; double Vy;
double Vz; double Vz; //7
double AntDirectX; double AntDirectX;
double AntDirectY; double AntDirectY;
double AntDirectZ; double AntDirectZ;
double AVx; double AVx;
double AVy; double AVy;
double AVz; double AVz;//13
double ZeroAntDiectX;
double ZeroAntDiectY;
double ZeroAntDiectZ;
double lon; double lon;
double lat; double lat;
double ati; // 15 double ati; // 19
}; };
@ -150,6 +159,4 @@ inline void delPointer(void* p)
} }
#endif #endif

BIN
BaseLibraryCPP.rar Normal file
View File

Binary file not shown.

43
EchoDataFormat.cpp
View File

@ -145,7 +145,7 @@ ErrorCode EchoL0Dataset::OpenOrNew(QString folder, QString filename, long PluseC
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 16, PluseCount, 1, GDT_Float64, NULL)); std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 19, PluseCount, 1, GDT_Float64, NULL));
GDALFlushCache((GDALDatasetH)poDstDS.get()); GDALFlushCache((GDALDatasetH)poDstDS.get());
poDstDS.reset(); poDstDS.reset();
omp_unset_lock(&lock); // omp_unset_lock(&lock); //
@ -296,22 +296,25 @@ SatelliteAntPos EchoL0Dataset::getSatelliteAntPos(long prf_id)
{ {
std::shared_ptr<double> antpos = this->getAntPos(); std::shared_ptr<double> antpos = this->getAntPos();
SatelliteAntPos prfpos{}; SatelliteAntPos prfpos{};
prfpos.time = antpos.get()[prf_id * 16 + 0]; prfpos.time = antpos.get()[prf_id *19 + 0];
prfpos.Px = antpos.get()[prf_id * 16 + 1]; prfpos.Px = antpos.get()[prf_id *19 + 1];
prfpos.Py = antpos.get()[prf_id * 16 + 2]; prfpos.Py = antpos.get()[prf_id *19 + 2];
prfpos.Pz = antpos.get()[prf_id * 16 + 3]; prfpos.Pz = antpos.get()[prf_id *19 + 3];
prfpos.Vx = antpos.get()[prf_id * 16 + 4]; prfpos.Vx = antpos.get()[prf_id *19 + 4];
prfpos.Vy = antpos.get()[prf_id * 16 + 5]; prfpos.Vy = antpos.get()[prf_id *19 + 5];
prfpos.Vz = antpos.get()[prf_id * 16 + 6]; prfpos.Vz = antpos.get()[prf_id *19 + 6];
prfpos.AntDirectX = antpos.get()[prf_id * 16 + 7]; prfpos.AntDirectX = antpos.get()[prf_id *19 + 7];
prfpos.AntDirectY = antpos.get()[prf_id * 16 + 8]; prfpos.AntDirectY = antpos.get()[prf_id *19 + 8];
prfpos.AntDirectZ = antpos.get()[prf_id * 16 + 9]; prfpos.AntDirectZ = antpos.get()[prf_id *19 + 9];
prfpos.AVx = antpos.get()[prf_id * 16 + 10]; prfpos.AVx = antpos.get()[prf_id *19 + 10];
prfpos.AVy = antpos.get()[prf_id * 16 + 11]; prfpos.AVy = antpos.get()[prf_id *19 + 11];
prfpos.AVz =antpos.get()[prf_id * 16 + 12]; prfpos.AVz =antpos.get()[prf_id *19 + 12];
prfpos.lon =antpos.get()[prf_id * 16 + 13]; prfpos.ZeroAntDiectX = antpos.get()[prf_id *19 + 13];
prfpos.lat =antpos.get()[prf_id * 16 + 14]; prfpos.ZeroAntDiectY = antpos.get()[prf_id *19 + 14];
prfpos.ati =antpos.get()[prf_id * 16 + 15]; prfpos.ZeroAntDiectZ = antpos.get()[prf_id *19 + 15];
prfpos.lon =antpos.get()[prf_id *19 + 16];
prfpos.lat =antpos.get()[prf_id *19 + 17];
prfpos.ati =antpos.get()[prf_id *19 + 18];
return prfpos; return prfpos;
} }
@ -467,8 +470,8 @@ std::shared_ptr<double> EchoL0Dataset::getAntPos()
std::shared_ptr<double> temp = nullptr; std::shared_ptr<double> temp = nullptr;
if (gdal_datatype == GDT_Float64) { if (gdal_datatype == GDT_Float64) {
temp=std::shared_ptr<double>(new double[this->PluseCount * 16],delArrPtr); temp=std::shared_ptr<double>(new double[this->PluseCount * 19],delArrPtr);
demBand->RasterIO(GF_Read, 0, 0, 16, this->PluseCount, temp.get(), 16, this->PluseCount, gdal_datatype, 0, 0); demBand->RasterIO(GF_Read, 0, 0, 19, this->PluseCount, temp.get(), 19, this->PluseCount, gdal_datatype, 0, 0);
} }
else { else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)) ; qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)) ;
@ -548,7 +551,7 @@ ErrorCode EchoL0Dataset::saveAntPos(std::shared_ptr<double> ptr)
long band_num = rasterDataset->GetRasterCount(); long band_num = rasterDataset->GetRasterCount();
if (gdal_datatype == GDT_Float64) { if (gdal_datatype == GDT_Float64) {
demBand->RasterIO(GF_Write, 0, 0, 16, this->PluseCount, ptr.get(), 16, this->PluseCount, gdal_datatype, 0, 0); demBand->RasterIO(GF_Write, 0, 0, 19, this->PluseCount, ptr.get(), 19, this->PluseCount, gdal_datatype, 0, 0);
} }
else { else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)); qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR));

279
GPUTool.cu Normal file
View File

@ -0,0 +1,279 @@
#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"
#ifdef __CUDANVCC___
#define CUDAMEMORY Memory1MB*100
#define LAMP_CUDA_PI 3.141592653589793238462643383279
// 定义参数
__device__ cuComplex cuCexpf(cuComplex x)
{
float factor = exp(x.x);
return make_cuComplex(factor * cos(x.y), factor * sin(x.y));
}
__global__ void CUDA_DistanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz,float *R, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
R[idx] = sqrtf(powf(Ax[idx]-Bx[idx], 2) + powf(Ay[idx] - By[idx], 2) + powf(Az[idx] - Bz[idx], 2));
}
}
__global__ void CUDA_B_DistanceA(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
R[idx] = sqrtf(powf(Ax[idx] - Bx, 2) + powf(Ay[idx] - By, 2) + powf(Az[idx] - Bz, 2));
}
}
__global__ void CUDA_make_VectorA_B(float sX, float sY, float sZ, float* tX, float* tY, float* tZ, float* RstX, float* RstY, float* RstZ, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
RstX[idx] = sX - tX[idx];
RstY[idx] = sY - tY[idx];
RstZ[idx] = sZ - tZ[idx];
}
}
__global__ void CUDA_Norm_Vector(float* Vx, float* Vy, float* Vz,float *R, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
R[idx] = sqrtf(powf(Vx[idx],2)+powf(Vy[idx],2)+powf(Vz[idx], 2));
}
}
__global__ void CUDA_cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos,long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float tAx = Ax[idx];
float tAy = Ay[idx];
float tAz = Az[idx];
float tBx = Bx[idx];
float tBy = By[idx];
float tBz = Bz[idx];
float AR = sqrtf(powf(tAx,2) + powf(tAy,2) + powf(tAz,2));
float BR = sqrtf(powf(tBx,2) + powf(tBy,2) + powf(tBz,2));
float dotAB = tAx * tBx + tAy * tBy + tAz * tBz;
float result =acosf( dotAB / (AR * BR));
anglecos[idx] = result;
}
}
__global__ void CUDA_SatelliteAntDirectNormal(float* RstX,float* RstY,float* RstZ,
float antXaxisX,float antXaxisY,float antXaxisZ,
float antYaxisX,float antYaxisY,float antYaxisZ,
float antZaxisX,float antZaxisY,float antZaxisZ,
float antDirectX,float antDirectY,float antDirectZ,
float* thetaAnt,float* phiAnt
, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float Xst = -1*RstX[idx]; // 卫星 --> 地面
float Yst = -1*RstY[idx];
float Zst = -1*RstZ[idx];
float AntXaxisX=antXaxisX ;
float AntXaxisY=antXaxisY ;
float AntXaxisZ=antXaxisZ ;
float AntYaxisX=antYaxisX ;
float AntYaxisY=antYaxisY ;
float AntYaxisZ=antYaxisZ ;
float AntZaxisX=antZaxisX ;
float AntZaxisY=antZaxisY ;
float AntZaxisZ=antZaxisZ ;
// 天线指向在天线坐标系下的值
float Xant = (Xst * (AntYaxisY * AntZaxisZ - AntYaxisZ * AntZaxisY) + Xst * ( AntXaxisZ * AntZaxisY - AntXaxisY * AntZaxisZ) + Xst * ( AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY)) / ( AntXaxisX * ( AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * ( AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * ( AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
float Yant = (Yst * (AntYaxisZ * AntZaxisX - AntYaxisX * AntZaxisZ) + Yst * ( AntXaxisX * AntZaxisZ - AntXaxisZ * AntZaxisX) + Yst * ( AntYaxisX * AntXaxisZ - AntXaxisX * AntYaxisZ)) / ( AntXaxisX * ( AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * ( AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * ( AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
float Zant = (Zst * (AntYaxisX * AntZaxisY - AntYaxisY * AntZaxisX) + Zst * ( AntXaxisY * AntZaxisX - AntXaxisX * AntZaxisY) + Zst * ( AntXaxisX * AntYaxisY - AntYaxisX * AntXaxisY)) / ( AntXaxisX * ( AntYaxisY * AntZaxisZ - AntZaxisY * AntYaxisZ) - AntYaxisX * ( AntXaxisY * AntZaxisZ - AntXaxisZ * AntZaxisY) + AntZaxisX * ( AntXaxisY * AntYaxisZ - AntXaxisZ * AntYaxisY));
// 计算theta 与 phi
float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
float ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
float YsinTheta = Yant / sinf(ThetaAnt);
float PhiAnt = (YsinTheta/abs(YsinTheta)) * acosf( Xant / (Norm * sinf(ThetaAnt)));
thetaAnt[idx] = ThetaAnt;
phiAnt[idx] = PhiAnt;
}
}
__global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R,float* slopeangle,
float nearRange, float Fs,float Pt,float lamda,long FreqIDmax,
cuComplex* echoArr , long* FreqID,
long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float r = R[idx];
float amp = Pt * TransAnt[idx] * ReciveAnt[idx];
amp= amp * sigma0[idx];
amp = amp / (powf(4* LAMP_CUDA_PI,2)*powf(r,4)); // 反射强度
// 处理相位
float phi = (-4 * LAMP_CUDA_PI / lamda) * r;
cuComplex echophi = make_cuComplex(0, phi) ;
cuComplex echophiexp = cuCexpf(echophi);
float timeR = 2 * (r - nearRange) / LIGHTSPEED * Fs;
long timeID = floorf(timeR);
if (timeID < 0 || timeID >= FreqIDmax) {
timeID = 0;
amp = 0;
}
cuComplex echo;
echo.x = echophiexp.x * amp;
echo.y = echophiexp.y * amp;
echoArr[idx] = echo;
}
}
//错误提示
void checkCudaError(cudaError_t err, const char* msg) {
if (err != cudaSuccess) {
std::cerr << "CUDA error: " << msg << " (" << cudaGetErrorString(err) << ")" << std::endl;
exit(EXIT_FAILURE);
}
}
// 主机参数内存声明
extern "C" void mallocCUDAHost(void* ptr, long memsize) {
cudaMallocHost(&ptr, memsize);
}
// 主机参数内存释放
extern "C" void FreeCUDAHost(void* ptr) {
cudaFreeHost(ptr);
}
// GPU参数内存声明
extern "C" void mallocCUDADevice(void* ptr, long memsize) {
cudaMalloc(&ptr, memsize);
}
// GPU参数内存释放
extern "C" void FreeCUDADevice(void* ptr) {
cudaFree(ptr);
}
// GPU 内存数据转移
extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize) {
cudaMemcpy(deviceptr, hostptr, memsize, cudaMemcpyHostToDevice);
}
extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize) {
cudaMemcpy(hostptr, deviceptr, memsize, cudaMemcpyDeviceToHost);
}
extern "C" void distanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* R,long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_DistanceAB << <blockSize, numBlocks >> > ( Ax, Ay, Az, Bx, By, Bz, R, len);
}
extern "C" void BdistanceAs(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_B_DistanceA << <blockSize, numBlocks >> > (Ax, Ay, Az, Bx, By, Bz, R, len);
cudaDeviceSynchronize();
}
extern "C" void make_VectorA_B(float sX, float sY, float sZ, float* tX, float* tY, float* tZ, float* RstX, float* RstY, float* RstZ, long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_make_VectorA_B << <blockSize, numBlocks >> > (sX, sY, sZ,tX, tY, tZ, RstX,RstY, RstZ, len);
cudaDeviceSynchronize();
}
extern "C" void Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long len) {
// 设置 CUDA 核函数的网格和块的尺寸
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_Norm_Vector << <blockSize, numBlocks >> > (Vx,Vy,Vz,R, len);
cudaDeviceSynchronize();
}
extern "C" void cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos, long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_cosAngle_VA_AB << <blockSize, numBlocks >> > (Ax, Ay, Az, Bx, By, Bz, anglecos, len);
cudaDeviceSynchronize();
}
extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* thetaAnt, float* phiAnt
, long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_SatelliteAntDirectNormal << <blockSize, numBlocks >> > ( RstX, RstY, RstZ,
antXaxisX, antXaxisY, antXaxisZ,
antYaxisX, antYaxisY, antYaxisZ,
antZaxisX, antZaxisY, antZaxisZ,
antDirectX, antDirectY, antDirectZ,
thetaAnt, phiAnt
, len);
cudaDeviceSynchronize();
}
extern "C" void calculationEcho(float* sigma0,float* TransAnt,float* ReciveAnt,
float* localangle,float* R, float* slopeangle,
float nearRange,float Fs, float pt, float lamda, long FreqIDmax,
cuComplex* echoAmp,long* FreqID,
long len)
{
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
// 调用 CUDA 核函数
CUDA_calculationEcho << <blockSize, numBlocks >> > ( sigma0, TransAnt,ReciveAnt,
localangle, R, slopeangle,
nearRange, Fs, pt, lamda, FreqIDmax,
echoAmp, FreqID,
len);
cudaDeviceSynchronize();
}
#endif

382
GPUTool.cuh Normal file
View File

@ -0,0 +1,382 @@
#ifndef GPUTOOL_H
#define GPUTOOL_H
#ifdef __CUDANVCC___
#include "BaseConstVariable.h"
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cublas_v2.h>
#include <cuComplex.h>
// 默认显存分布
enum LAMPGPUDATETYPE {
LAMP_LONG,
LAMP_FLOAT,
LAMP_COMPLEXFLOAT
};
// GPU 内存函数
extern "C" void mallocCUDAHost(void* ptr, long memsize); // 主机内存声明
extern "C" void FreeCUDAHost(void* ptr);
extern "C" void mallocCUDADevice(void* ptr, long memsize); // GPU内存声明
extern "C" void FreeCUDADevice(void* ptr);
extern "C" void HostToDevice(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移 设备 -> GPU
extern "C" void DeviceToHost(void* hostptr, void* deviceptr, long memsize);//GPU 内存数据转移 GPU -> 设备
// 仿真所需的常用函数
extern "C" void distanceAB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* R, long member);
extern "C" void BdistanceAs(float* Ax, float* Ay, float* Az, float Bx, float By, float Bz, float* R, long member);
extern "C" void make_VectorA_B(float sX, float sY, float sZ, float* tX, float* tY, float* tZ, float* RstX, float* RstY, float* RstZ, long member);
extern "C" void Norm_Vector(float* Vx, float* Vy, float* Vz, float* R, long member);
extern "C" void cosAngle_VA_AB(float* Ax, float* Ay, float* Az, float* Bx, float* By, float* Bz, float* anglecos, long len);
extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ, float antXaxisX, float antXaxisY, float antXaxisZ, float antYaxisX, float antYaxisY, float antYaxisZ, float antZaxisX, float antZaxisY, float antZaxisZ, float antDirectX, float antDirectY, float antDirectZ, float* thetaAnt, float* phiAnt, long len);
extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,float* localangle, float* R, float* slopeangle,float nearRange, float Fs, float pt, float lamda, long FreqIDmax,cuComplex* echoAmp, long* FreqID, long len);
#endif
#endif
/**
*
double* databuffer = new double[nXSize * nYSize * 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);
Eigen::MatrixXcd rasterData(nYSize, nXSize); // 使用Eigen的MatrixXcd
for(size_t i = 0; i < nYSize; i++) {
for(size_t j = 0; j < nXSize; j++) {
rasterData(i, j) = std::complex<double>(databuffer[i * nXSize * 2 + j * 2],
databuffer[i * nXSize * 2 + j * 2 + 1]);
}
}
delete[] databuffer;
gdalImage demxyz(demxyzPath);// 地面点坐标
gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型
gdalImage demsloperxyz(this->demsloperPath);// 地面坡向
omp_lock_t lock; // 定义锁
omp_init_lock(&lock); // 初始化锁
long start_ids = 1250;
for (start_ids = 1; start_ids < demxyz.height; start_ids = start_ids + line_invert) { // 8+ 17 + 0.3 MB
QDateTime current = QDateTime::currentDateTime();
long pluseStep = Memory1MB * 100 / 3 / PlusePoint;
if (pluseStep * num_thread * 3 > this->PluseCount) {
pluseStep = this->PluseCount / num_thread / 3;
}
pluseStep = pluseStep > 50 ? pluseStep : 50;
qDebug() << current.toString("yyyy-MM-dd HH:mm:ss.zzz") << " \tstart \t " << start_ids << " - " << start_ids + line_invert << "\t" << demxyz.height << "\t pluseCount:\t" << pluseStep;
// 文件读取
Eigen::MatrixXd dem_x = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 1); //
Eigen::MatrixXd dem_y = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 2); //
Eigen::MatrixXd dem_z = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 3); //
// 地表覆盖
std::shared_ptr<long> dem_landcls = readDataArr<long>(demlandcls, start_ids - 1, 0, line_invert + 1, demxyz.width, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD); // 地表覆盖类型
long* dem_landcls_ptr = dem_landcls.get();
double localAngle = 30.0;
bool sigmaNoZeroFlag = true;
for (long ii = 0; ii < dem_x.rows(); ii++) {
for (long jj = 0; jj < dem_y.cols(); jj++) {
if (0 != this->SigmaDatabasePtr->getAmp(dem_landcls_ptr[dem_x.cols() * ii + jj], localAngle, polartype)) {
sigmaNoZeroFlag = false;
break;
}
}
if (!sigmaNoZeroFlag) {
break;
}
}
if (sigmaNoZeroFlag) {
continue;
}
//#ifdef DEBUGSHOWDIALOG
// dialog->load_double_MatrixX_data(dem_z, "dem_z");
//#endif
Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 1); //
Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 2); //
Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 3); //
Eigen::MatrixXd sloperAngle = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 4); //
sloperAngle = sloperAngle.array() * T180_PI;
long dem_rows = dem_x.rows();
long dem_cols = dem_x.cols();
long freqidx = 0;//
#ifdef DEBUGSHOWDIALOG
ImageShowDialogClass* dialog = new ImageShowDialogClass(nullptr);
dialog->show();
Eigen::MatrixXd landaArr = Eigen::MatrixXd::Zero(dem_rows, dem_cols);
for (long i = 0; i < dem_rows; i++) {
for (long j = 0; j < dem_cols; j++) {
landaArr(i, j) = dem_landcls.get()[i * dem_cols + j];
}
}
dialog->load_double_MatrixX_data(landaArr, "landCover");
#endif
//qDebug() << " pluse bolck size :\t " << pluseStep << " all size:\t" << this->PluseCount;
long processNumber = 0;
#pragma omp parallel for
for (long startprfidx = 0; startprfidx < this->PluseCount; startprfidx = startprfidx + pluseStep) { // 17 + 0.3 MB
long prfcount_step = startprfidx + pluseStep < this->PluseCount ? pluseStep : this->PluseCount - startprfidx;
Eigen::MatrixXcd echoPluse = Eigen::MatrixXcd::Zero(prfcount_step, PlusePoint); // 当前脉冲的回波积分情况
// 内存预分配
Eigen::MatrixXd Rst_x = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Rst_y = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Rst_z = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd R = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd localangle = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Vst_x = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Vst_y = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Vst_z = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd fde = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd fr = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd Rx = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd sigam = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd echoAmp = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols()).array() + Pt;
Eigen::MatrixXd Rphi = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd TimeRange = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd TransAnt = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd ReciveAnt = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd AntTheta = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
Eigen::MatrixXd AntPhi = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
double minR = 0, maxR = 0;
double minLocalAngle = 0, maxLocalAngle = 0;
Vector3D Rt = { 0,0,0 };
SatelliteOribtNode oRs = SatelliteOribtNode{ 0 };;
Vector3D p0 = {}, slopeVector = {}, sateAntDirect = {};
Vector3D Rs = {}, Vs = {}, Ast = {};
SatelliteAntDirect antdirectNode = {};
std::complex<double> echofreq;
std::complex<double> Imag1(0, 1);
double TAntPattern = 1; // 发射天线方向图
double RAntPanttern = 1;// 接收天线方向图
double maxechoAmp = 1;
double tempAmp = 1;
for (long prfidx = 0; prfidx < prfcount_step; prfidx++)
{
oRs = sateOirbtNodes[prfidx + startprfidx];
// 计算天线方向图
for (long jj = 1; jj < dem_cols - 1; jj++) {
for (long ii = 1; ii < dem_rows - 1; ii++) {
p0.x = dem_x(ii, jj);
p0.y = dem_y(ii, jj);
p0.z = dem_z(ii, jj);
this->TaskSetting->getSatelliteAntDirectNormal(oRs, p0, antdirectNode);
//antdirectNode.ThetaAnt = antdirectNode.ThetaAnt * r2d;
//antdirectNode.PhiAnt = antdirectNode.PhiAnt * r2d;
AntTheta(ii, jj) = antdirectNode.ThetaAnt * r2d;
AntPhi(ii, jj) = antdirectNode.PhiAnt * r2d;
}
}
// 计算发射天线方向图
for (long jj = 1; jj < dem_cols - 1; jj++) {
for (long ii = 1; ii < dem_rows - 1; ii++) {
TransformPattern->getGainLinear(AntTheta(ii, jj), AntPhi(ii, jj), TransAnt(ii, jj));
//TransAnt(ii, jj) = TAntPattern;
}
}
// 计算接收天线方向图
for (long jj = 1; jj < dem_cols - 1; jj++) {
for (long ii = 1; ii < dem_rows - 1; ii++) {
TransformPattern->getGainLinear(AntTheta(ii, jj), AntPhi(ii, jj), ReciveAnt(ii, jj));
//ReciveAnt(ii, jj) = RAntPanttern;
}
}
// 计算经过增益的能量
echoAmp = Pt * TransAnt.array() * ReciveAnt.array();
maxechoAmp = echoAmp.maxCoeff();
if (std::abs(maxechoAmp) < PRECISIONTOLERANCE) { // 这种情况下，不在合成孔径范围中
continue;
}
Rs.x = sateOirbtNodes[prfidx + startprfidx].Px; // 卫星位置
Rs.y = sateOirbtNodes[prfidx + startprfidx].Py;
Rs.z = sateOirbtNodes[prfidx + startprfidx].Pz;
Vs.x = sateOirbtNodes[prfidx + startprfidx].Vx; // 卫星速度
Vs.y = sateOirbtNodes[prfidx + startprfidx].Vy;
Vs.z = sateOirbtNodes[prfidx + startprfidx].Vz;
Ast.x = sateOirbtNodes[prfidx + startprfidx].AVx;// 卫星加速度
Ast.y = sateOirbtNodes[prfidx + startprfidx].AVy;
Ast.z = sateOirbtNodes[prfidx + startprfidx].AVz;
Rst_x = Rs.x - dem_x.array(); // Rst = Rs - Rt;
Rst_y = Rs.y - dem_y.array();
Rst_z = Rs.z - dem_z.array();
R = (Rst_x.array().pow(2) + Rst_y.array().pow(2) + Rst_z.array().pow(2)).array().sqrt(); // R
minR = R.minCoeff();
maxR = R.maxCoeff();
//qDebug() << "minR:\t" << minR << " maxR:\t" << maxR;
if (maxR<NearRange || minR>FarRange) {
continue;
}
else {}
// getCosAngle
// double c = dot(a, b) / (getlength(a) * getlength(b));
// return acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d;
// localangle = getCosAngle(Rst, slopeVector) * T180_PI; // 注意这个只能实时计算，因为非实时计算代价太大
localangle = (Rst_x.array() * demsloper_x.array() + Rst_y.array() * demsloper_y.array() + Rst_z.array() * demsloper_z.array()).array(); // dot(a, b)
localangle = localangle.array() / R.array();
localangle = localangle.array() / (demsloper_x.array().pow(2) + demsloper_y.array().pow(2) + demsloper_z.array().pow(2)).array().sqrt().array();
localangle = localangle.array().acos(); // 弧度值
minLocalAngle = localangle.minCoeff();
maxLocalAngle = localangle.maxCoeff();
if (maxLocalAngle<0 || minLocalAngle>PI / 2) {
continue;
}
else {}
//Vst_x = Vs.x + 1 * earthRoute * dem_y.array(); // Vst = Vs - Vt;
//Vst_y = Vs.y - 1 * earthRoute * dem_x.array();
//Vst_z = Vs.z - Eigen::MatrixXd::Zero(dem_x.rows(), dem_y.cols()).array();
//// 计算多普勒中心频率 Rst, Vst : ( - 2 / lamda) * dot(Rs - Rt, Vs - Vt) / R; // 星载合成孔径雷达原始回波数据模拟研究 3.18
//fde = (-2 / lamda) * (Rst_x.array() * Vst_x.array() + Rst_y.array() * Vst_y.array() + Rst_z.array() * Vst_z.array()).array() / (R.array());
//// 计算多普勒频率斜率 // 星载合成孔径雷达原始回波数据模拟研究 3.19
//// -(2/lamda)*( dot(Vs - Vt, Vs - Vt)/R + dot(Ast, Rs - Rt)/R - std::pow(dot(Vs - Vt, Rs - Rt),2 )/std::pow(R,3));
//fr = (-2 / lamda) *
// (Vst_x.array() * Vst_x.array() + Vst_y.array() * Vst_y.array() + Vst_z.array() * Vst_z.array()).array() / (R.array()) +
// (-2 / lamda) *
// (Ast.x * Rst_x.array() + Ast.y * Rst_y.array() + Ast.z * Rst_z.array()).array() / (R.array()) -
// (-2 / lamda) *
// (Vst_x.array() * Rst_x.array() + Vst_y.array() * Rst_y.array() + Vst_z.array() * Rst_z.array()).array().pow(2) / (R.array().pow(3));
// 计算回波
Rx = R;// -(lamda / 2) * (fde * TRx + 0.5 * fr * TRx * TRx); // 斜距历程值
// 逐点计算 this->SigmaDatabasePtr->getAmp(covercls, localangle, polartype); // 后向散射系数 HH
for (long ii = 0; ii < dem_x.rows(); ii++) {
for (long jj = 0; jj < dem_y.cols(); jj++) {
sigam(ii, jj) = this->SigmaDatabasePtr->getAmp(dem_landcls_ptr[dem_x.cols() * ii + jj], localangle(ii, jj) * r2d, polartype);
}
}
if (sigam.maxCoeff() > 0) {}
else {
continue;
}
// projArea = 1 / std::cos(sloperAngle) * std::sin(localangle); // 投影面积系数，单位投影面积 1m x 1m --注意这里是假设，后期再补充
// echoAmp = projArea*TAntPattern * RAntPanttern * sigam / (4 * PI * R * R);
echoAmp = echoAmp.array() * sigam.array() * (1 / sloperAngle.array().cos() * localangle.array().sin()); // 反射强度
echoAmp = echoAmp.array() / (4 * PI * R.array().pow(2)); // 距离衰减
Rphi = -4 * PI / lamda * Rx.array();// 距离徙动相位
// 积分
TimeRange = ((2 * R.array() / LIGHTSPEED - TimgNearRange).array() * Fs).array();
double localAnglepoint = -1;
long prf_freq_id = 0;
for (long jj = 1; jj < dem_cols - 1; jj++) {
for (long ii = 1; ii < dem_rows - 1; ii++) {
prf_freq_id = std::floor(TimeRange(ii, jj));
if (prf_freq_id < 0 || prf_freq_id >= PlusePoint || localangle(ii, jj) < 0 || localangle(ii, jj) > PI / 2 || echoAmp(ii, jj) == 0) {
continue;
}
echofreq = echoAmp(ii, jj) * std::exp(Rphi(ii, jj) * Imag1);
echoPluse(prfidx, prf_freq_id) = echoPluse(prfidx, prf_freq_id) + echofreq;
}
}
#ifdef DEBUGSHOWDIALOG
ImageShowDialogClass* localangledialog = new ImageShowDialogClass(dialog);
localangledialog->show();
localangledialog->load_double_MatrixX_data(localangle.array() * r2d, "localangle");
ImageShowDialogClass* sigamdialog = new ImageShowDialogClass(dialog);
sigamdialog->show();
sigamdialog->load_double_MatrixX_data(TimeRange, "TimeRange");
ImageShowDialogClass* ampdialog = new ImageShowDialogClass(dialog);
ampdialog->show();
ampdialog->load_double_MatrixX_data(echoAmp, "echoAmp");
Eigen::MatrixXd echoPluseamp = echoPluse.array().abs().cast<double>().array();
ImageShowDialogClass* echoampdialog = new ImageShowDialogClass(dialog);
echoampdialog->show();
echoampdialog->load_double_MatrixX_data(echoPluseamp, "echoPluseamp");
dialog->exec();
#endif
//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << " end " << prfidx;
}
//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz")<<" step "<< prfcount_step;
omp_set_lock(&lock); // 回波整体赋值处理
for (long prfidx = 0; prfidx < prfcount_step; prfidx++) {
for (long freqidx = 0; freqidx < PlusePoint; freqidx++)
{
//qDebug() << prfidx << " " << freqidx << " " << echoPluse(prfidx, freqidx).real() << " + " << echoPluse(prfidx, freqidx).imag() << " j";
echo.get()[(prfidx + startprfidx) * PlusePoint + freqidx] = echo.get()[(prfidx + startprfidx) * PlusePoint + freqidx] + echoPluse(prfidx, freqidx);
}
}
//this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
omp_unset_lock(&lock); // 解锁
//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << " step 2" << prfcount_step;
}
omp_set_lock(&lock); // 保存文件
processNumber = processNumber + pluseStep;
this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
omp_unset_lock(&lock); // 解锁
qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << " \t " << start_ids << "\t--\t " << start_ids + line_invert << "\t/\t" << demxyz.height;
}
omp_destroy_lock(&lock); // 销毁锁
*/

36
ImageOperatorBase.h
View File

@ -287,7 +287,7 @@ std::shared_ptr<T> readDataArr(gdalImage& imgds, int start_row, int start_col, i
rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row; rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row;
cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col; cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col;
Eigen::MatrixXd datamatrix(rows_count, cols_count); //Eigen::MatrixXd datamatrix(rows_count, cols_count);
if (gdal_datatype == GDT_Byte) { if (gdal_datatype == GDT_Byte) {
char* temp = new char[rows_count * cols_count]; char* temp = new char[rows_count * cols_count];
@ -416,6 +416,40 @@ std::shared_ptr<T> readDataArr(gdalImage& imgds, int start_row, int start_col, i
} }
delete[] temp; delete[] temp;
} }
else if (gdal_datatype == GDT_CFloat32) {
float* temp = new float[rows_count * cols_count*2];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = std::complex<float>(temp[i * 2 ],
temp[i * 2 + 1]);
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_CFloat64) {
double* temp = new double[rows_count * cols_count*2];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
for (long i = 0; i < count; i++) {
result.get()[i] = std::complex<double>(temp[i * 2],
temp[i * 2 + 1]);
}
}
delete[] temp;
}
else { else {
} }
GDALClose((GDALDatasetH)rasterDataset); GDALClose((GDALDatasetH)rasterDataset);

8
SARSimulationImageL1.cpp
View File

@ -115,7 +115,7 @@ ErrorCode SARSimulationImageL1Dataset::OpenOrNew(QString folder, QString filenam
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 16, rowCount, 1, GDT_Float64, NULL)); std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 19, rowCount, 1, GDT_Float64, NULL));
GDALFlushCache((GDALDatasetH)poDstDS.get()); GDALFlushCache((GDALDatasetH)poDstDS.get());
poDstDS.reset(); poDstDS.reset();
omp_unset_lock(&lock); omp_unset_lock(&lock);
@ -563,8 +563,8 @@ std::shared_ptr<double> SARSimulationImageL1Dataset::getAntPos()
std::shared_ptr<double> temp = nullptr; std::shared_ptr<double> temp = nullptr;
if (gdal_datatype == GDT_Float64) { if (gdal_datatype == GDT_Float64) {
temp = std::shared_ptr<double>(new double[this->rowCount * 16], delArrPtr); temp = std::shared_ptr<double>(new double[this->rowCount * 19], delArrPtr);
demBand->RasterIO(GF_Read, 0, 0, 10, this->rowCount, temp.get(), 16, this->rowCount, gdal_datatype, 0, 0); demBand->RasterIO(GF_Read, 0, 0, 19, this->rowCount, temp.get(), 19, this->rowCount, gdal_datatype, 0, 0);
} }
else { else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)); qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR));
@ -591,7 +591,7 @@ ErrorCode SARSimulationImageL1Dataset::saveAntPos(std::shared_ptr<double> ptr)
long band_num = rasterDataset->GetRasterCount(); long band_num = rasterDataset->GetRasterCount();
if (gdal_datatype == GDT_Float64) { if (gdal_datatype == GDT_Float64) {
demBand->RasterIO(GF_Write, 0, 0, 16, this->rowCount, ptr.get(), 16, this->rowCount, gdal_datatype, 0, 0); demBand->RasterIO(GF_Write, 0, 0, 19, this->rowCount, ptr.get(), 19, this->rowCount, gdal_datatype, 0, 0);
} }
else { else {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)); qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR));

1
SARSimulationImageL1.h
View File

@ -14,7 +14,6 @@ enum RasterLevel {
RasterL2 RasterL2
}; };
class SARSimulationImageL1Dataset class SARSimulationImageL1Dataset
{ {
public: public: