diff --git a/BaseConstVariable.h b/BaseConstVariable.h index bb099df..941cf24 100644 --- a/BaseConstVariable.h +++ b/BaseConstVariable.h @@ -3,7 +3,7 @@ #ifndef BASECONSTVARIABLE_H #define BASECONSTVARIABLE_H #define EIGEN_USE_MKL_ALL -//#define EIGEN_NO_DEBUG +#define EIGEN_NO_DEBUG #define EIGEN_USE_BLAS @@ -12,11 +12,17 @@ //#define DEBUGSHOWDIALOG +#define __CUDANVCC___ // 瀹氫箟CUDA鍑芥暟 + + + + //#include #include #include #include + #define MATPLOTDRAWIMAGE @@ -120,16 +126,19 @@ struct SatelliteAntPos { double Pz; double Vx; double Vy; - double Vz; + double Vz; //7 double AntDirectX; double AntDirectY; double AntDirectZ; double AVx; double AVy; - double AVz; + double AVz;//13 + double ZeroAntDiectX; + double ZeroAntDiectY; + double ZeroAntDiectZ; double lon; double lat; - double ati; // 15 + double ati; // 19 }; @@ -148,8 +157,6 @@ inline void delPointer(void* p) delete p; p = nullptr; } - - - + #endif \ No newline at end of file diff --git a/BaseLibraryCPP.rar b/BaseLibraryCPP.rar new file mode 100644 index 0000000..8bf6eaa Binary files /dev/null and b/BaseLibraryCPP.rar differ diff --git a/EchoDataFormat.cpp b/EchoDataFormat.cpp index 2e24a7b..b6ceecf 100644 --- a/EchoDataFormat.cpp +++ b/EchoDataFormat.cpp @@ -145,7 +145,7 @@ ErrorCode EchoL0Dataset::OpenOrNew(QString folder, QString filename, long PluseC CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); - std::shared_ptr poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 16, PluseCount, 1, GDT_Float64, NULL)); + std::shared_ptr poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 19, PluseCount, 1, GDT_Float64, NULL)); GDALFlushCache((GDALDatasetH)poDstDS.get()); poDstDS.reset(); omp_unset_lock(&lock); // @@ -296,22 +296,25 @@ SatelliteAntPos EchoL0Dataset::getSatelliteAntPos(long prf_id) { std::shared_ptr antpos = this->getAntPos(); SatelliteAntPos prfpos{}; - prfpos.time = antpos.get()[prf_id * 16 + 0]; - prfpos.Px = antpos.get()[prf_id * 16 + 1]; - prfpos.Py = antpos.get()[prf_id * 16 + 2]; - prfpos.Pz = antpos.get()[prf_id * 16 + 3]; - prfpos.Vx = antpos.get()[prf_id * 16 + 4]; - prfpos.Vy = antpos.get()[prf_id * 16 + 5]; - prfpos.Vz = antpos.get()[prf_id * 16 + 6]; - prfpos.AntDirectX = antpos.get()[prf_id * 16 + 7]; - prfpos.AntDirectY = antpos.get()[prf_id * 16 + 8]; - prfpos.AntDirectZ = antpos.get()[prf_id * 16 + 9]; - prfpos.AVx = antpos.get()[prf_id * 16 + 10]; - prfpos.AVy = antpos.get()[prf_id * 16 + 11]; - prfpos.AVz =antpos.get()[prf_id * 16 + 12]; - prfpos.lon =antpos.get()[prf_id * 16 + 13]; - prfpos.lat =antpos.get()[prf_id * 16 + 14]; - prfpos.ati =antpos.get()[prf_id * 16 + 15]; + prfpos.time = antpos.get()[prf_id *19 + 0]; + prfpos.Px = antpos.get()[prf_id *19 + 1]; + prfpos.Py = antpos.get()[prf_id *19 + 2]; + prfpos.Pz = antpos.get()[prf_id *19 + 3]; + prfpos.Vx = antpos.get()[prf_id *19 + 4]; + prfpos.Vy = antpos.get()[prf_id *19 + 5]; + prfpos.Vz = antpos.get()[prf_id *19 + 6]; + prfpos.AntDirectX = antpos.get()[prf_id *19 + 7]; + prfpos.AntDirectY = antpos.get()[prf_id *19 + 8]; + prfpos.AntDirectZ = antpos.get()[prf_id *19 + 9]; + prfpos.AVx = antpos.get()[prf_id *19 + 10]; + prfpos.AVy = antpos.get()[prf_id *19 + 11]; + prfpos.AVz =antpos.get()[prf_id *19 + 12]; + prfpos.ZeroAntDiectX = antpos.get()[prf_id *19 + 13]; + prfpos.ZeroAntDiectY = antpos.get()[prf_id *19 + 14]; + 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; } @@ -467,8 +470,8 @@ std::shared_ptr EchoL0Dataset::getAntPos() std::shared_ptr temp = nullptr; if (gdal_datatype == GDT_Float64) { - temp=std::shared_ptr(new double[this->PluseCount * 16],delArrPtr); - demBand->RasterIO(GF_Read, 0, 0, 16, this->PluseCount, temp.get(), 16, this->PluseCount, gdal_datatype, 0, 0); + temp=std::shared_ptr(new double[this->PluseCount * 19],delArrPtr); + demBand->RasterIO(GF_Read, 0, 0, 19, this->PluseCount, temp.get(), 19, this->PluseCount, gdal_datatype, 0, 0); } else { qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)) ; @@ -548,7 +551,7 @@ ErrorCode EchoL0Dataset::saveAntPos(std::shared_ptr ptr) long band_num = rasterDataset->GetRasterCount(); 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 { qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)); diff --git a/GPUTool.cu b/GPUTool.cu new file mode 100644 index 0000000..81aad62 --- /dev/null +++ b/GPUTool.cu @@ -0,0 +1,279 @@ + + +#include +#include +#include +#include +#include +#include +#include +#include + +#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 << > > ( 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 << > > (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 << > > (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 << > > (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 << > > (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 << > > ( 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 << > > ( sigma0, TransAnt,ReciveAnt, + localangle, R, slopeangle, + nearRange, Fs, pt, lamda, FreqIDmax, + echoAmp, FreqID, + len); + cudaDeviceSynchronize(); + +} + + + + + + +#endif \ No newline at end of file diff --git a/GPUTool.cuh b/GPUTool.cuh new file mode 100644 index 0000000..505fd09 --- /dev/null +++ b/GPUTool.cuh @@ -0,0 +1,382 @@ +#ifndef GPUTOOL_H +#define GPUTOOL_H +#ifdef __CUDANVCC___ +#include "BaseConstVariable.h" +#include +#include +#include +#include + +// 默认显存分布 + + +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(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 dem_landcls = readDataArr(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 echofreq; + std::complex 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 (maxRFarRange) { + 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().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); // 销毁锁 + + + + + + + + + + +*/ \ No newline at end of file diff --git a/ImageOperatorBase.h b/ImageOperatorBase.h index bc9a859..3e2d33a 100644 --- a/ImageOperatorBase.h +++ b/ImageOperatorBase.h @@ -287,7 +287,7 @@ std::shared_ptr readDataArr(gdalImage& imgds, int start_row, int start_col, i 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; - Eigen::MatrixXd datamatrix(rows_count, cols_count); + //Eigen::MatrixXd datamatrix(rows_count, cols_count); if (gdal_datatype == GDT_Byte) { char* temp = new char[rows_count * cols_count]; @@ -416,6 +416,40 @@ std::shared_ptr readDataArr(gdalImage& imgds, int start_row, int start_col, i } 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(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(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(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(temp[i * 2], + temp[i * 2 + 1]); + } + } + delete[] temp; + } else { } GDALClose((GDALDatasetH)rasterDataset); diff --git a/SARSimulationImageL1.cpp b/SARSimulationImageL1.cpp index 1e78618..4fae528 100644 --- a/SARSimulationImageL1.cpp +++ b/SARSimulationImageL1.cpp @@ -115,7 +115,7 @@ ErrorCode SARSimulationImageL1Dataset::OpenOrNew(QString folder, QString filenam CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); - std::shared_ptr poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 16, rowCount, 1, GDT_Float64, NULL)); + std::shared_ptr poDstDS(poDriver->Create(this->GPSPointFilePath.toUtf8().constData(), 19, rowCount, 1, GDT_Float64, NULL)); GDALFlushCache((GDALDatasetH)poDstDS.get()); poDstDS.reset(); omp_unset_lock(&lock); @@ -563,8 +563,8 @@ std::shared_ptr SARSimulationImageL1Dataset::getAntPos() std::shared_ptr temp = nullptr; if (gdal_datatype == GDT_Float64) { - temp = std::shared_ptr(new double[this->rowCount * 16], delArrPtr); - demBand->RasterIO(GF_Read, 0, 0, 10, this->rowCount, temp.get(), 16, this->rowCount, gdal_datatype, 0, 0); + temp = std::shared_ptr(new double[this->rowCount * 19], delArrPtr); + demBand->RasterIO(GF_Read, 0, 0, 19, this->rowCount, temp.get(), 19, this->rowCount, gdal_datatype, 0, 0); } else { qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)); @@ -591,7 +591,7 @@ ErrorCode SARSimulationImageL1Dataset::saveAntPos(std::shared_ptr ptr) long band_num = rasterDataset->GetRasterCount(); 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 { qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_GPSFILEFORMATERROR)); diff --git a/SARSimulationImageL1.h b/SARSimulationImageL1.h index c0ed6db..eb3cecf 100644 --- a/SARSimulationImageL1.h +++ b/SARSimulationImageL1.h @@ -14,7 +14,6 @@ enum RasterLevel { RasterL2 }; - class SARSimulationImageL1Dataset { public: