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

代码更新： 

1. 仿真部分代码修改为GPU
2. 删除了原CPU的仿真代码
3. GPU代码部分，由 SimulationSAR 文件，移动到 BaseLibraryCPP 文件中
4. 补充了成像平面的 matlab 代码
Release
陈增辉 2024-11-29 23:32:50 +08:00
parent e56b8f2d98
commit bee90080f2
18 changed files with 1446 additions and 612 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
BaseLibraryCPP

@ -1 +1 @@
Subproject commit 458be41d96524bbdad14aaa0b6184a5de387e7e5
Subproject commit ef94aa388ae2ed85e970ff4f857450a5f28cf0b2

237
README.md
View File

@ -1 +1,238 @@
# RasterProcessTool
extern "C" void AddCUDA(void* aPtr, void* bptr, void* cptr, long member, LAMPGPUDATETYPE datetype) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (member + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
}
// CUDA 核函数
__global__ void computeDistanceAndEchoID(float* antPx, float* antPy, float* antPz,
float* img_x, float* img_y, float* img_z,
complexfloat* echopluse, complexfloat* imgarr,
long rowcount, long colcount, long prfid,
float Rnear, float fs, float factorj) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
// 确保线程索引有效
if (idx < rowcount * colcount) {
// 计算距离
float dx = antPx[prfid] - img_x[idx];
float dy = antPy[prfid] - img_y[idx];
float dz = antPz[prfid] - img_z[idx];
float imgR = sqrt(dx * dx + dy * dy + dz * dz);
// 计算 EchoID
long echoID = floor(((imgR - Rnear) * 2 / LIGHTSPEED) * fs);//回波坐标
float Rftj = imgR * factorj; // 校正
//printf("%d:(%f,%f),%f,%f |||||, %f, %f, %f, %f, %f, %f, %f\n", idx, echopluse[echoID].real, Rftj, imgR, fs,
// antPx[prfid], antPy[prfid], antPz[prfid], img_x[idx], img_y[idx], img_z[idx], imgR);
if (echoID < 0 || echoID >= colcount) {
}
else {
complexfloat Rphi{ 0,Rftj };
Rphi = expComplex(Rphi);
imgarr[idx] = addComplex(imgarr[idx], mulComplex(echopluse[echoID], Rphi));
}
}
}
void TBPImageGPUBlock(float* antPx, float* antPy, float* antPz, float* img_x, float* img_y, float* img_z,
std::shared_ptr<std::complex<double>> echoArr, long prfcount, long plusecount,
std::shared_ptr<std::complex<double>> imageArr,
float freq, float fs, float Rnear, float Rfar, float factorj,
long startline, long stepline,
long rowcount, long colcount) {
long pixelcount = rowcount * colcount;
complexfloat* h_echopluse;
complexfloat* h_imgarr;
cudaMallocHost(&h_echopluse, sizeof(float) * 2 * plusecount); // 单个传感器的位置
cudaMallocHost(&h_imgarr, sizeof(float) * 2 * stepline * colcount);
for (long i = 0; i < stepline; i++) {
long rid = startline + i;
for (long j = 0; j < colcount; j++) {
h_imgarr[i * colcount + j].real = imageArr.get()[rid * colcount + j].real();
h_imgarr[i * colcount + j].imag = imageArr.get()[rid * colcount + j].imag();
}
}
std::cout << "h_imgarr init finished!!" << std::endl;
float* h_antPx, * h_antPy, * h_antPz;
cudaMallocHost(&h_antPx, sizeof(float) * prfcount); // 单个传感器的位置
cudaMallocHost(&h_antPy, sizeof(float) * prfcount);
cudaMallocHost(&h_antPz, sizeof(float) * prfcount);
// 初始化
for (long i = 0; i < prfcount; i++) {
h_antPx[i] = antPx[i];
h_antPy[i] = antPy[i];
h_antPz[i] = antPz[i];
}
float* h_img_x, * h_img_y, * h_img_z;
cudaMallocHost(&h_img_x, sizeof(float) * stepline * colcount);
cudaMallocHost(&h_img_y, sizeof(float) * stepline * colcount);
cudaMallocHost(&h_img_z, sizeof(float) * stepline * colcount);
// 初始化
long rid = 0;
for (long i = 0; i < stepline; i++) {
rid = startline + i;
for (long j = 0; j < colcount; j++) {
h_img_x[i * colcount + j] = img_x[rid * colcount + j];
h_img_y[i * colcount + j] = img_y[rid * colcount + j];
h_img_z[i * colcount + j] = img_z[rid * colcount + j];
}
}
std::cout << "h_img_x init finished!!" << std::endl;
// 分配设备内存
float* d_antPx, * d_antPy, * d_antPz, * d_img_x, * d_img_y, * d_img_z;
complexfloat* d_echopluse;
complexfloat* d_imgarr;
cudaMalloc(&d_echopluse, sizeof(float) * 2 * plusecount);
cudaMalloc(&d_imgarr, sizeof(float) * 2 * stepline * colcount);
cudaMalloc(&d_antPx, sizeof(float) * prfcount);
cudaMalloc(&d_antPy, sizeof(float) * prfcount);
cudaMalloc(&d_antPz, sizeof(float) * prfcount);
cudaMalloc(&d_img_x, sizeof(float) * rowcount * colcount);
cudaMalloc(&d_img_y, sizeof(float) * rowcount * colcount);
cudaMalloc(&d_img_z, sizeof(float) * rowcount * colcount);
// 将数据从主机拷贝到设备
cudaMemcpy(d_antPx, h_antPx, sizeof(float) * prfcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_antPx, h_antPx, sizeof(float) * prfcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_antPy, h_antPy, sizeof(float) * prfcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_antPz, h_antPz, sizeof(float) * prfcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_img_x, h_img_x, sizeof(float) * rowcount * colcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_img_y, h_img_y, sizeof(float) * rowcount * colcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_img_z, h_img_z, sizeof(float) * rowcount * colcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_imgarr, h_imgarr, sizeof(float) * 2 * stepline * colcount, cudaMemcpyHostToDevice);
int blockSize = 256; // 每个块的线程数
int numBlocks = (pixelcount + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
long eid = 0;
std::complex<double> Rphi;
for (long prfid = 0; prfid < prfcount; prfid++) {
for (long i = 0; i < plusecount; i++) {
h_echopluse[i].real = echoArr.get()[prfid * plusecount + i].real();
h_echopluse[i].imag = echoArr.get()[prfid * plusecount + i].imag();
}
cudaMemcpy(d_echopluse, h_echopluse, sizeof(float) * 2 * plusecount, cudaMemcpyHostToDevice);
computeDistanceAndEchoID << < numBlocks, blockSize >> > (d_antPx, d_antPy, d_antPz,
d_img_x, d_img_y, d_img_z,
d_echopluse, d_imgarr,
rowcount, colcount, prfid,
Rnear, fs, factorj);
cudaDeviceSynchronize();// 等待所有设备任务完成
if (prfid % 100 == 0) {
//std::cout << "\rprf " << prfid <<"/"<< prfcount << "\t\t\t";
}
//cudaMemcpy(h_echopluse, d_echopluse, sizeof(float) * 2 * stepline * colcount, cudaMemcpyDeviceToHost);
}
std::cout << std::endl;
// GPU -> CPU
cudaMemcpy(h_imgarr, d_imgarr, sizeof(float) * 2 * stepline * colcount, cudaMemcpyDeviceToHost);
for (long i = 0; i < stepline; i++) {
long rid = startline + i;
for (long j = 0; j < colcount; j++) {
imageArr.get()[rid * colcount + j] = std::complex<double>(h_imgarr[i * colcount + j].real, h_imgarr[i * colcount + j].imag);
}
}
// 清理资源
cudaFree(d_antPx);
cudaFree(d_antPy);
cudaFree(d_antPz);
cudaFree(d_img_x);
cudaFree(d_img_y);
cudaFree(d_img_z);
cudaFree(d_echopluse);
cudaFree(d_imgarr);
cudaFreeHost(h_antPx);
cudaFreeHost(h_antPy);
cudaFreeHost(h_antPz);
cudaFreeHost(h_img_x);
cudaFreeHost(h_img_y);
cudaFreeHost(h_img_z);
cudaFreeHost(h_echopluse);
cudaFreeHost(h_imgarr);
std::cout << "end GPU" << std::endl;
}
void RTPC(float* antx, float* anty, float* antz,
float* demx, float* demy, float* demz,
float* demslopex, float* demslopey, float* demslopez
) {}

6
RasterProcessTool.vcxproj
View File

@ -179,6 +179,7 @@
<QtMoc Include="SimulationSAR\QToolAbstract.h" />
<QtMoc Include="RegisterToolbox.h" />
<QtMoc Include="SimulationSAR\QSimulationBPImage.h" />
<ClInclude Include="BaseLibraryCPP\GPUTool.cuh" />
<ClInclude Include="SimulationSAR\TBPImageAlgCls.h" />
<QtMoc Include="QSimulationRTPCGUI.h" />
<QtMoc Include="GF3ProcessToolbox\QOrthSlrRaster.h" />
@ -206,8 +207,9 @@
<QtMoc Include="QMergeRasterProcessDialog.h" />
</ItemGroup>
<ItemGroup>
<CudaCompile Include="SimulationSAR\TBPGPU.cu">
<FileType>Document</FileType>
<CudaCompile Include="BaseLibraryCPP\GPUTool.cu">
<ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
</ExcludedFromBuild>
</CudaCompile>
</ItemGroup>
<Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />

7
RasterProcessTool.vcxproj.filters
View File

@ -219,6 +219,9 @@
<ClInclude Include="SimulationSAR\TBPImageAlgCls.h">
<Filter>SimulationSAR</Filter>
</ClInclude>
<ClInclude Include="BaseLibraryCPP\GPUTool.cuh">
<Filter>BaseLibraryCPP</Filter>
</ClInclude>
</ItemGroup>
<ItemGroup>
<QtMoc Include="QMergeRasterProcessDialog.h">
@ -297,8 +300,8 @@
</QtUic>
</ItemGroup>
<ItemGroup>
<CudaCompile Include="SimulationSAR\TBPGPU.cu">
<Filter>SimulationSAR</Filter>
<CudaCompile Include="BaseLibraryCPP\GPUTool.cu">
<Filter>BaseLibraryCPP</Filter>
</CudaCompile>
</ItemGroup>
</Project>

4
SimulationSAR/QSimulationBPImage.cpp
View File

@ -81,8 +81,10 @@ void QSimulationBPImage::onbtnaccepted()
TBPimag.setEchoL0(echoL0ds);
TBPimag.setImageL1(imagL1);
long cpucore_num = std::thread::hardware_concurrency();
TBPimag.setGPU(true);
TBPimag.Process(cpucore_num);
QMessageBox::information(this,u8"成像",u8"成像结束");
}
void QSimulationBPImage::onbtnrejected()

6
SimulationSAR/QSimulationBPImage.ui
View File

@ -44,6 +44,9 @@
<height>30</height>
</size>
</property>
<property name="text">
<string>D:\Programme\vs2022\RasterMergeTest\TestData\outData\GF3_Simulation.xml</string>
</property>
</widget>
</item>
<item>
@ -92,6 +95,9 @@
<height>30</height>
</size>
</property>
<property name="text">
<string>D:\Programme\vs2022\RasterMergeTest\TestData\outBPImage2\GF3_Simulation</string>
</property>
</widget>
</item>
<item>

718
SimulationSAR/RTPCProcessCls.cpp
View File

@ -13,25 +13,33 @@
#include <QDatetime>
#include <omp.h>
#include <QProgressDialog>
#include <QMessageBox>
#ifdef DEBUGSHOWDIALOG
#include "ImageShowDialogClass.h"
#endif
#ifdef __CUDANVCC___
#include "GPUTool.cuh"
#endif // __CUDANVCC___
RTPCProcessCls::RTPCProcessCls()
{
this->PluseCount = 0;
this->PlusePoint = 0;
this->TaskSetting = nullptr;
this->EchoSimulationData = nullptr;
this->DEMTiffPath="";
this->LandCoverPath="";
this->HHSigmaPath="";
this->HVSigmaPath="";
this->VHSigmaPath="";
this->DEMTiffPath = "";
this->LandCoverPath = "";
this->HHSigmaPath = "";
this->HVSigmaPath = "";
this->VHSigmaPath = "";
this->VVSigmaPath = "";
this->OutEchoPath = "";
@ -43,7 +51,7 @@ RTPCProcessCls::RTPCProcessCls()
this->VVSigmaPath.clear();
this->OutEchoPath.clear();
this->SigmaDatabasePtr=std::shared_ptr<SigmaDatabase>(new SigmaDatabase);
this->SigmaDatabasePtr = std::shared_ptr<SigmaDatabase>(new SigmaDatabase);
}
@ -53,13 +61,13 @@ RTPCProcessCls::~RTPCProcessCls()
void RTPCProcessCls::setTaskSetting(std::shared_ptr < AbstractSARSatelliteModel> TaskSetting)
{
this->TaskSetting= std::shared_ptr < AbstractSARSatelliteModel>(TaskSetting);
this->TaskSetting = std::shared_ptr < AbstractSARSatelliteModel>(TaskSetting);
qDebug() << "RTPCProcessCls::setTaskSetting";
}
void RTPCProcessCls::setEchoSimulationDataSetting(std::shared_ptr<EchoL0Dataset> EchoSimulationData)
{
this->EchoSimulationData =std::shared_ptr<EchoL0Dataset> (EchoSimulationData);
this->EchoSimulationData = std::shared_ptr<EchoL0Dataset>(EchoSimulationData);
qDebug() << "RTPCProcessCls::setEchoSimulationDataSetting";
}
@ -110,7 +118,7 @@ ErrorCode RTPCProcessCls::Process(long num_thread)
// RTPC 算法
qDebug() << u8"params init ....";
ErrorCode stateCode = this->InitParams();
if (stateCode != ErrorCode::SUCCESS){
if (stateCode != ErrorCode::SUCCESS) {
return stateCode;
}
else {}
@ -122,12 +130,12 @@ ErrorCode RTPCProcessCls::Process(long num_thread)
else {}
qDebug() << "RTPCMainProcess";
stateCode = this->RTPCMainProcess(num_thread);
if (stateCode != ErrorCode::SUCCESS) {
return stateCode;
}else{}
}
else {}
return ErrorCode::SUCCESS;
@ -135,7 +143,7 @@ ErrorCode RTPCProcessCls::Process(long num_thread)
ErrorCode RTPCProcessCls::InitParams()
{
if (nullptr==this->TaskSetting||this->DEMTiffPath.isEmpty() ||
if (nullptr == this->TaskSetting || this->DEMTiffPath.isEmpty() ||
this->LandCoverPath.isEmpty() || this->HHSigmaPath.isEmpty() ||
this->HVSigmaPath.isEmpty() || this->VHSigmaPath.isEmpty() ||
this->VVSigmaPath.isEmpty()) {
@ -157,7 +165,7 @@ ErrorCode RTPCProcessCls::InitParams()
// 初始化回波存放位置
qDebug() << "--------------Echo Data Setting ---------------------------------------";
this->EchoSimulationData=std::shared_ptr<EchoL0Dataset>(new EchoL0Dataset);
this->EchoSimulationData = std::shared_ptr<EchoL0Dataset>(new EchoL0Dataset);
this->EchoSimulationData->setCenterFreq(this->TaskSetting->getCenterFreq());
this->EchoSimulationData->setNearRange(this->TaskSetting->getNearRange());
this->EchoSimulationData->setFarRange(this->TaskSetting->getFarRange());
@ -222,7 +230,7 @@ ErrorCode RTPCProcessCls::DEMPreprocess()
// 计算坡向角
this->demsloperPath=QDir(tmpfolderPath).filePath("demsloper.tif");
this->demsloperPath = QDir(tmpfolderPath).filePath("demsloper.tif");
this->demmaskPath = QDir(tmpfolderPath).filePath("demmask.tif");
gdalImage demsloperxyz = CreategdalImage(this->demsloperPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true);// X,Y,Z,cosangle
@ -242,7 +250,7 @@ ErrorCode RTPCProcessCls::DEMPreprocess()
Eigen::MatrixXd demsloper_angle = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 4);
maskdata = maskdata.array() * 0;
Landpoint p0, p1, p2, p3, p4,pslopeVector,pp;
Landpoint p0, p1, p2, p3, p4, pslopeVector, pp;
Vector3D slopeVector;
dem_rows = maskdata.rows();
@ -257,10 +265,10 @@ ErrorCode RTPCProcessCls::DEMPreprocess()
rowidx = i + startlineid;
colidx = j;
demds.getLandPoint(rowidx, colidx, demdata(i, j), p0);
demds.getLandPoint(rowidx-1, colidx, demdata(i-1, j), p1);
demds.getLandPoint(rowidx, colidx-1, demdata(i, j-1), p2);
demds.getLandPoint(rowidx+1, colidx, demdata(i+1, j), p3);
demds.getLandPoint(rowidx, colidx+1, demdata(i, j+1), p4);
demds.getLandPoint(rowidx - 1, colidx, demdata(i - 1, j), p1);
demds.getLandPoint(rowidx, colidx - 1, demdata(i, j - 1), p2);
demds.getLandPoint(rowidx + 1, colidx, demdata(i + 1, j), p3);
demds.getLandPoint(rowidx, colidx + 1, demdata(i, j + 1), p4);
pslopeVector = getSlopeVector(p0, p1, p2, p3, p4); // 地面坡向矢量
slopeVector = { pslopeVector.lon,pslopeVector.lat,pslopeVector.ati };
@ -268,7 +276,7 @@ ErrorCode RTPCProcessCls::DEMPreprocess()
Zaxis.x = pp.lon;
Zaxis.y = pp.lat;
Zaxis.z = pp.ati;
sloperAngle = getCosAngle(slopeVector, Zaxis) ; // 地面坡向角
sloperAngle = getCosAngle(slopeVector, Zaxis); // 地面坡向角
demsloper_x(i, j) = slopeVector.x;
demsloper_y(i, j) = slopeVector.y;
@ -291,7 +299,7 @@ ErrorCode RTPCProcessCls::DEMPreprocess()
ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
{
omp_set_num_threads(num_thread);// 设置openmp 线程数量
double widthSpace = LIGHTSPEED/2/this->TaskSetting->getFs();
double widthSpace = LIGHTSPEED / 2 / this->TaskSetting->getFs();
double prf_time = 0;
double dt = 1 / this->TaskSetting->getPRF();// 获取每次脉冲的时间间隔
@ -299,10 +307,10 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
Landpoint LandP{ 0,0,0 };
Point3 GERpoint{ 0,0,0 };
double R = 0;
double dem_row = 0 , dem_col=0, dem_alt=0;
double dem_row = 0, dem_col = 0, dem_alt = 0;
long double imageStarttime = 0;
imageStarttime=this->TaskSetting->getSARImageStartTime();
imageStarttime = this->TaskSetting->getSARImageStartTime();
//std::vector<SatelliteOribtNode> sateOirbtNodes(this->PluseCount);
std::shared_ptr<SatelliteOribtNode[]> sateOirbtNodes(new SatelliteOribtNode[this->PluseCount], delArrPtr);
{ // 姿态计算不同
@ -310,7 +318,7 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
std::shared_ptr<double> antpos = this->EchoSimulationData->getAntPos();
double dAt = 1e-6;
double prf_time_dt = 0;
Landpoint InP{0,0,0}, outP{0,0,0};
Landpoint InP{ 0,0,0 }, outP{ 0,0,0 };
for (long prf_id = 0; prf_id < this->PluseCount; prf_id++) {
prf_time = dt * prf_id;
prf_time_dt = prf_time + dAt;
@ -327,8 +335,7 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
InP.lat = sateOirbtNode.Py;
InP.ati = sateOirbtNode.Pz;
outP=XYZ2LLA(InP);
outP = XYZ2LLA(InP);
antpos.get()[prf_id * 19 + 0] = prf_time + imageStarttime;
antpos.get()[prf_id * 19 + 1] = sateOirbtNode.Px;
@ -356,14 +363,10 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
qDebug() << "Ant position finished sucessfully !!!";
}
// 回波
long echoIdx = 0;
gdalImage demxyz (demxyzPath );// 地面点坐标
gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型
gdalImage demsloperxyz(this->demsloperPath);// 地面坡向
double NearRange = this->EchoSimulationData->getNearRange(); // 近斜据
double FarRange = this->EchoSimulationData->getFarRange();
@ -372,353 +375,390 @@ ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
double TimgFarRange = 2 * FarRange / LIGHTSPEED;
double Fs = this->TaskSetting->getFs(); // 距离向采样率
double Pt = this->TaskSetting->getPt()*this->TaskSetting->getGri();// 发射电压 1v
double Pt = this->TaskSetting->getPt() * this->TaskSetting->getGri();// 发射电压 1v
//double GainAntLen = -3;// -3dB 为天线半径
long pluseCount = this->PluseCount;
double lamda = this->TaskSetting->getCenterLamda(); // 波长
// 天线方向图
std::shared_ptr<AbstractRadiationPattern> TransformPattern = this->TaskSetting->getTransformRadiationPattern(); // 发射天线方向图
std::shared_ptr< AbstractRadiationPattern> ReceivePattern = this->TaskSetting->getReceiveRadiationPattern(); // 接收天线方向图
std::shared_ptr<AbstractRadiationPattern> ReceivePattern = this->TaskSetting->getReceiveRadiationPattern(); // 接收天线方向图
std::shared_ptr<std::complex<double>> echo = this->EchoSimulationData->getEchoArr();
long PlusePoint = this->EchoSimulationData->getPlusePoints();
// 初始化 为 0
for (long i = 0; i < pluseCount * PlusePoint; i++) {
echo.get()[i] = std::complex<double>(0,0);
echo.get()[i] = std::complex<double>(0, 0);
}
this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
POLARTYPEENUM polartype = this->TaskSetting->getPolarType();
long demrowStep = std::ceil(Memory1MB*10 / 8 / demxyz.width) ;
long line_invert = demrowStep > 3 ? demrowStep : 3;
#ifndef __CUDANVCC___
QMessageBox::information(this, u8"程序提示", u8"请确定安装了CUDA库");
#else
double lamda = this->TaskSetting->getCenterLamda(); // 波长
omp_lock_t lock; // 定义锁
omp_init_lock(&lock); // 初始化锁
// DEM计算
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;
QProgressDialog progressDialog(u8"RTPC回波生成中", u8"终止", 0, this->PluseCount);
progressDialog.setWindowTitle(u8"RTPC回波生成中");
progressDialog.setWindowModality(Qt::WindowModal);
progressDialog.setAutoClose(true);
progressDialog.setValue(0);
progressDialog.setMaximum(this->PluseCount);
progressDialog.setMinimum(0);
progressDialog.show();
#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;
}
// RTPC CUDA版本
if (pluseCount * 4 * 18 > Memory1MB * 100) {
long max = Memory1MB * 100 / 4 / 20 / PluseCount;
QMessageBox::warning(nullptr, u8"仿真场景太大了", u8"当前频点数下，脉冲数量最多为："+QString::number(max));
}
// 计算经过增益的能量
echoAmp = Pt * TransAnt.array() * ReciveAnt.array();
gdalImage demxyz(this->demxyzPath);// 地面点坐标
gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型
gdalImage demsloperxyz(this->demsloperPath);// 地面坡向
long demRow = demxyz.height;
long demCol = demxyz.width;
long blokline = 100;
// 每块 250MB*16 = 4GB
// dem 9
//
blokline = Memory1MB * 500 / 8 / demCol;
blokline = blokline < 1 ? 1 : blokline;
bool bloklineflag = false;
Eigen::MatrixXd dem_x = demxyz.getData(0, 0, blokline, demxyz.width, 1); // 地面坐标
long tempDemRows = dem_x.rows();
long tempDemCols = dem_x.cols();
Eigen::MatrixXd dem_y = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd dem_z = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd demsloper_x = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd demsloper_y = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd demsloper_z = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd sloperAngle = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
float* h_dem_x;
float* h_dem_y;
float* h_dem_z;
float* h_demsloper_x;
float* h_demsloper_y;
float* h_demsloper_z;
float* h_demsloper_angle;
float* d_dem_x;
float* d_dem_y;
float* d_dem_z;
float* d_demsloper_x;
float* d_demsloper_y;
float* d_demsloper_z;
float* d_demsloper_angle;
mallocCUDAHost((void*)h_dem_x, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_dem_y, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_dem_z, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_demsloper_x, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_demsloper_y, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_demsloper_z, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_demsloper_angle, sizeof(float) * blokline * tempDemCols);
mallocCUDADevice((void*)d_dem_x, sizeof(float) * blokline * tempDemCols); // 7
mallocCUDADevice((void*)d_dem_y, sizeof(float) * blokline * tempDemCols);
mallocCUDADevice((void*)d_dem_z, sizeof(float) * blokline * tempDemCols);
mallocCUDADevice((void*)d_demsloper_x, sizeof(float) * blokline * tempDemCols);
mallocCUDADevice((void*)d_demsloper_y, sizeof(float) * blokline * tempDemCols);
mallocCUDADevice((void*)d_demsloper_z, sizeof(float) * blokline * tempDemCols);
mallocCUDADevice((void*)d_demsloper_angle, sizeof(float) * blokline * tempDemCols);
float* h_dem_theta; // 天线方向图
float* h_dem_phi;
float* d_dem_theta;
float* d_dem_phi;
mallocCUDAHost((void*)h_dem_theta, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_dem_phi, sizeof(float) * blokline * tempDemCols);
mallocCUDADevice((void*)d_dem_theta, sizeof(float) * blokline * tempDemCols);// 9
mallocCUDADevice((void*)d_dem_phi, sizeof(float) * blokline * tempDemCols);
HostToDevice((void*)h_dem_theta, (void*)d_dem_theta, sizeof(float) * blokline * tempDemCols);
HostToDevice((void*)h_dem_phi, (void*)d_dem_phi, sizeof(float) * blokline * tempDemCols);
// 提前声明参数变量
float* h_R;// 辐射方向
float* h_localangle;//入射角
float* d_R;// 辐射方向
float* d_localangle;//入射角
mallocCUDAHost((void*)h_R, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_localangle, sizeof(float) * blokline * tempDemCols); // 11
mallocCUDADevice((void*)d_R, sizeof(float) * blokline * tempDemCols);
mallocCUDADevice((void*)d_localangle, sizeof(float) * blokline * tempDemCols);
float* h_RstX;
float* h_RstY;
float* h_RstZ;
float* d_RstX;
float* d_RstY;
float* d_RstZ;
mallocCUDAHost((void*)h_RstX, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_RstY, sizeof(float) * blokline * tempDemCols);
mallocCUDAHost((void*)h_RstZ, sizeof(float) * blokline * tempDemCols); // 14
mallocCUDADevice((void*)d_RstX, sizeof(float)* blokline* tempDemCols);
mallocCUDADevice((void*)d_RstY, sizeof(float)* blokline* tempDemCols);
mallocCUDADevice((void*)d_RstZ, sizeof(float)* blokline* tempDemCols);
float* h_sigma0;
float* h_TransAnt;
float* h_ReciveAnt;
float* d_sigma0;
float* d_TransAnt;
float* d_ReciveAnt;
mallocCUDAHost((void*)h_sigma0 , sizeof(float)* blokline* tempDemCols);
mallocCUDAHost((void*)h_TransAnt , sizeof(float)* blokline* tempDemCols);
mallocCUDAHost((void*)h_ReciveAnt, sizeof(float)* blokline* tempDemCols); // 17
mallocCUDADevice((void*)d_sigma0, sizeof(float)* blokline* tempDemCols);
mallocCUDADevice((void*)d_TransAnt, sizeof(float)* blokline* tempDemCols);
mallocCUDADevice((void*)d_ReciveAnt, sizeof(float)* blokline* tempDemCols);
// 回波
cuComplex* h_echoAmp;
cuComplex* d_echoAmp;
mallocCUDAHost((void*)h_echoAmp, sizeof(cuComplex)* blokline* tempDemCols);
mallocCUDADevice((void*)d_echoAmp, sizeof(cuComplex)* blokline* tempDemCols); //19
long* h_FreqID;
long* d_FreqID;
mallocCUDAHost((void*)h_FreqID, sizeof(long)* blokline* tempDemCols);
mallocCUDADevice((void*)d_FreqID, sizeof(long)* blokline* tempDemCols); //21
Eigen::MatrixXd landcover= Eigen::MatrixXd::Zero(blokline, tempDemCols);// 地面覆盖类型
long startline = 0;
for (startline = 0; startline < demRow; startline = startline + blokline) {
long newblokline = blokline;
if ((startline + blokline) < demRow) {
newblokline = demRow - startline;
bloklineflag = true;
}
dem_x = demxyz.getData(startline, 0, newblokline, demxyz.width, 1); // 地面坐标
dem_y = demxyz.getData(startline, 0, newblokline, demxyz.width, 2);
dem_z = demxyz.getData(startline, 0, newblokline, demxyz.width, 3);
demsloper_x = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 1);
demsloper_y = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 2);
demsloper_z = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 3);
sloperAngle = demsloperxyz.getData(startline, 0, newblokline, demsloperxyz.width, 4);
landcover = demlandcls.getData(startline, 0, newblokline, demlandcls.width, 1);
if (bloklineflag) {
FreeCUDAHost(h_dem_x); FreeCUDADevice(d_dem_x);
FreeCUDAHost(h_dem_y); FreeCUDADevice(d_dem_y);
FreeCUDAHost(h_dem_z); FreeCUDADevice(d_dem_z);
FreeCUDAHost(h_demsloper_x); FreeCUDADevice(d_demsloper_x);
FreeCUDAHost(h_demsloper_y); FreeCUDADevice(d_demsloper_y);
FreeCUDAHost(h_demsloper_z); FreeCUDADevice(d_demsloper_z); //6
FreeCUDAHost(h_demsloper_angle); FreeCUDADevice(d_demsloper_angle);//7
FreeCUDAHost(h_dem_theta); FreeCUDADevice(d_dem_theta);
FreeCUDAHost(h_dem_phi); FreeCUDADevice(d_dem_phi); //9
FreeCUDAHost(h_R); FreeCUDADevice(d_R);
FreeCUDAHost(h_localangle); FreeCUDADevice(d_localangle); //11
FreeCUDAHost(h_RstX); FreeCUDADevice(d_RstX);
FreeCUDAHost(h_RstY); FreeCUDADevice(d_RstY);
FreeCUDAHost(h_RstZ); FreeCUDADevice(d_RstZ); //14
FreeCUDAHost(h_sigma0 ); FreeCUDADevice(d_sigma0);
FreeCUDAHost(h_TransAnt ); FreeCUDADevice(d_TransAnt);
FreeCUDAHost(h_ReciveAnt); FreeCUDADevice(d_ReciveAnt); //17
FreeCUDAHost(h_echoAmp); FreeCUDADevice(d_echoAmp);//19
FreeCUDAHost(h_FreqID); FreeCUDADevice(d_FreqID);//20
mallocCUDAHost((void*)h_dem_x, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_dem_y, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_dem_z, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_demsloper_x, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_demsloper_y, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_demsloper_z, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_demsloper_angle, sizeof(float) * blokline * tempDemCols);//7
mallocCUDAHost((void*)h_dem_theta, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_dem_phi, sizeof(float) * newblokline * tempDemCols); //9
mallocCUDAHost((void*)h_R, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_localangle, sizeof(float) * newblokline * tempDemCols);//11
mallocCUDAHost((void*)h_RstX, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_RstY, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_RstZ, sizeof(float) * newblokline * tempDemCols);//14
mallocCUDAHost((void*)h_sigma0, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_TransAnt, sizeof(float) * newblokline * tempDemCols);
mallocCUDAHost((void*)h_ReciveAnt, sizeof(float) * newblokline * tempDemCols);//17
mallocCUDAHost((void*)h_echoAmp, sizeof(cuComplex) * newblokline * tempDemCols);//19
mallocCUDAHost((void*)h_FreqID, sizeof(long) * newblokline * tempDemCols);//20
mallocCUDADevice((void*)d_dem_x, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_dem_y, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_dem_z, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_demsloper_x, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_demsloper_y, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_demsloper_z, sizeof(float) * newblokline * tempDemCols);//6
mallocCUDADevice((void*)d_demsloper_angle, sizeof(float) * newblokline * tempDemCols);//7
mallocCUDADevice((void*)d_dem_theta, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_dem_phi, sizeof(float) * newblokline * tempDemCols);// 9
mallocCUDADevice((void*)d_R, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_localangle, sizeof(float) * newblokline * tempDemCols);//11
mallocCUDADevice((void*)d_RstX, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_RstY, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_RstZ, sizeof(float) * newblokline * tempDemCols);//14
mallocCUDADevice((void*)d_sigma0, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_TransAnt, sizeof(float) * newblokline * tempDemCols);
mallocCUDADevice((void*)d_ReciveAnt, sizeof(float) * newblokline * tempDemCols);//17
mallocCUDADevice((void*)d_echoAmp, sizeof(cuComplex) * newblokline * tempDemCols); //19
mallocCUDADevice((void*)d_FreqID, sizeof(long) * newblokline * tempDemCols); //20
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;
for (long i = 0; i < newblokline; i++) {
for (long j = 0; j < demxyz.width; j++) {
h_dem_x[i * demxyz.width + j] = dem_x(i, j);
h_dem_y[i * demxyz.width + j] = dem_y(i, j);
h_dem_z[i * demxyz.width + j] = dem_z(i, j);
h_demsloper_x[i * demxyz.width + j] = demsloper_x(i, j);
h_demsloper_y[i * demxyz.width + j] = demsloper_y(i, j);
h_demsloper_z[i * demxyz.width + j] = demsloper_z(i, j);
h_demsloper_angle[i * demxyz.width + j] = sloperAngle(i, j);
}
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();
HostToDevice((void*)h_dem_x, (void*)d_dem_x, sizeof(float) * newblokline * tempDemCols); // 复制 机器 -> GPU
HostToDevice((void*)h_dem_y, (void*)d_dem_y, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_z, (void*)d_dem_z, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_x, (void*)d_demsloper_x, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_y, (void*)d_demsloper_y, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_z, (void*)d_demsloper_z, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_angle, (void*)d_demsloper_angle, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_theta, (void*)d_dem_theta, sizeof(float) * blokline * tempDemCols);
HostToDevice((void*)h_dem_phi, (void*)d_dem_phi, sizeof(float) * blokline * tempDemCols);
long pixelcount = newblokline * tempDemCols;
for (long prfid = 0; prfid < pluseCount; prfid++) {
{// 计算
// 天线位置
float antpx = sateOirbtNodes[prfid].Px;
float antpy = sateOirbtNodes[prfid].Py;
float antpz = sateOirbtNodes[prfid].Pz;
float antvx = sateOirbtNodes[prfid].Vx;
float antvy = sateOirbtNodes[prfid].Vy;
float antvz = sateOirbtNodes[prfid].Vz; //6
float antdirectx = sateOirbtNodes[prfid].AntDirecX;
float antdirecty = sateOirbtNodes[prfid].AntDirecY;
float antdirectz = sateOirbtNodes[prfid].AntDirecZ; // 9 天线指向
float antXaxisX = sateOirbtNodes[prfid].AntXaxisX;
float antXaxisY = sateOirbtNodes[prfid].AntXaxisY;
float antXaxisZ = sateOirbtNodes[prfid].AntXaxisZ;//12 天线坐标系
float antYaxisX = sateOirbtNodes[prfid].AntYaxisX;
float antYaxisY = sateOirbtNodes[prfid].AntYaxisY;
float antYaxisZ = sateOirbtNodes[prfid].AntYaxisZ;//15
float antZaxisX = sateOirbtNodes[prfid].AntZaxisX;
float antZaxisY = sateOirbtNodes[prfid].AntZaxisY;
float antZaxisZ = sateOirbtNodes[prfid].AntZaxisZ;//18
make_VectorA_B(antpx, antpy, antpz, d_dem_x, d_dem_y, d_dem_z, d_RstX, d_RstY, d_RstZ, pixelcount); // Rst = Rs - Rt;
Norm_Vector(d_RstX, d_RstY, d_RstZ,d_R,pixelcount); // R
cosAngle_VA_AB(d_RstX, d_RstY, d_RstZ, d_demsloper_x, d_demsloper_y, d_demsloper_z,d_localangle ,pixelcount); // 局部入射角
SatelliteAntDirectNormal(d_RstX, d_RstY, d_RstZ,antXaxisX, antXaxisY, antXaxisZ,antYaxisX, antYaxisY, antYaxisZ,antZaxisX, antZaxisY, antZaxisZ,antdirectx, antdirecty, antdirectz,d_dem_theta, d_dem_phi , pixelcount);// 计算角度
DeviceToHost(h_dem_theta, d_dem_theta, sizeof(float)*pixelcount); // 从GPU -> 主机
DeviceToHost(h_dem_phi, d_dem_phi, sizeof(float)*pixelcount);
DeviceToHost(h_localangle, d_localangle, sizeof(float)*pixelcount);
for (long ii =0; ii <blokline; ii++) { // 计算发射天线方向图
for (long jj = 0; jj < tempDemCols; jj++) {
h_TransAnt[ii * tempDemCols + jj] =TransformPattern->getGainLearThetaPhi(h_dem_theta[ii*tempDemCols+ jj], h_dem_phi[ii * tempDemCols + jj] );
}
}
for (long ii = 0; ii < blokline; ii++) { // 计算接收天线方向图
for (long jj = 0; jj < tempDemCols; jj++) {
h_ReciveAnt[ii * tempDemCols + jj] = ReceivePattern->getGainLearThetaPhi(h_dem_theta[ii * tempDemCols + jj], h_dem_phi[ii * tempDemCols + jj]);
}
}
for (long ii = 0; ii < blokline; ii++) {// 后向散射图
for (long jj = 0; jj < tempDemCols; jj++) {
h_sigma0[ii * tempDemCols + jj] = this->SigmaDatabasePtr->getAmp(landcover(ii,jj), h_localangle[ii*tempDemCols+jj] * r2d, polartype);
}
}
HostToDevice((void*)h_sigma0, (void*)d_sigma0, sizeof(float)* pixelcount);
HostToDevice((void*)h_TransAnt, (void*)d_TransAnt, sizeof(float)* pixelcount);
HostToDevice((void*)h_ReciveAnt, (void*)d_ReciveAnt, sizeof(float)* pixelcount);
//// 计算多普勒中心频率 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); // 斜距历程值
calculationEcho(d_sigma0, d_TransAnt, d_ReciveAnt,d_localangle, d_R, d_demsloper_angle,NearRange, Fs, Pt, lamda, PlusePoint,d_echoAmp, d_FreqID, pixelcount);
DeviceToHost(h_echoAmp, d_echoAmp, sizeof(float)* pixelcount);
DeviceToHost(h_FreqID, d_FreqID, sizeof(long)* pixelcount);
// 回波存档
for (long i = 0; i < pixelcount; i++) {
echo.get()[i*PlusePoint+ h_FreqID[i]] =std::complex<double>(h_echoAmp[i].x, h_echoAmp[i].y);
}
// 逐点计算 this->SigmaDatabasePtr->getAmp(covercls, localangle, polartype); // 后向散射系数 HH
std::cout << "\r"<<"dem:\t"<<startline <<"\t-\t"<<startline+newblokline<<"\t:\t pluse :\t"<<prfid<<" / " <<pluseCount ;
}
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;
std::cout << std::endl;
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;
}
}
// 地面数据释放
FreeCUDAHost(h_dem_x); FreeCUDADevice(d_dem_x);
FreeCUDAHost(h_dem_y); FreeCUDADevice(d_dem_y);
FreeCUDAHost(h_dem_z); FreeCUDADevice(d_dem_z);
FreeCUDAHost(h_demsloper_x); FreeCUDADevice(d_demsloper_x);
FreeCUDAHost(h_demsloper_y); FreeCUDADevice(d_demsloper_y);
FreeCUDAHost(h_demsloper_z); FreeCUDADevice(d_demsloper_z); //6
FreeCUDAHost(h_demsloper_angle); FreeCUDADevice(d_demsloper_angle); //7
#ifdef DEBUGSHOWDIALOG
ImageShowDialogClass* localangledialog = new ImageShowDialogClass(dialog);
localangledialog->show();
localangledialog->load_double_MatrixX_data(localangle.array()*r2d, "localangle");
// 临时变量释放
FreeCUDAHost(h_dem_theta); FreeCUDADevice(d_dem_theta);
FreeCUDAHost(h_dem_phi); FreeCUDADevice(d_dem_phi);// 9
FreeCUDAHost(h_R); FreeCUDADevice(d_R);
FreeCUDAHost(h_localangle); FreeCUDADevice(h_localangle); //11
FreeCUDAHost(h_RstX); FreeCUDADevice(d_RstX);
FreeCUDAHost(h_RstY); FreeCUDADevice(d_RstY);
FreeCUDAHost(h_RstZ); FreeCUDADevice(d_RstZ); //14
FreeCUDAHost(h_sigma0); FreeCUDADevice(d_sigma0);
FreeCUDAHost(h_TransAnt); FreeCUDADevice(d_TransAnt);
FreeCUDAHost(h_ReciveAnt); FreeCUDADevice(d_ReciveAnt); //17
FreeCUDAHost(h_echoAmp); FreeCUDADevice(d_echoAmp);//19
FreeCUDAHost(h_FreqID); FreeCUDADevice(d_FreqID);//20
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;
processNumber = processNumber < progressDialog.maximum() ? processNumber : progressDialog.maximum();
progressDialog.setValue(processNumber);
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;
progressDialog.close();
}
omp_destroy_lock(&lock); // 销毁锁
this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
this->EchoSimulationData->saveToXml();
return ErrorCode::SUCCESS;
}
void RTPCProcessMain(long num_thread,QString TansformPatternFilePath, QString ReceivePatternFilePath, QString simulationtaskName,QString OutEchoPath,QString GPSXmlPath, QString TaskXmlPath, QString demTiffPath , QString LandCoverPath , QString HHSigmaPath , QString HVSigmaPath , QString VHSigmaPath , QString VVSigmaPath)
void RTPCProcessMain(long num_thread, QString TansformPatternFilePath, QString ReceivePatternFilePath, QString simulationtaskName, QString OutEchoPath, QString GPSXmlPath, QString TaskXmlPath, QString demTiffPath, QString LandCoverPath, QString HHSigmaPath, QString HVSigmaPath, QString VHSigmaPath, QString VVSigmaPath)
{
std::vector<RadiationPatternGainPoint> TansformPatternGainpoints = ReadGainFile(TansformPatternFilePath);
std::shared_ptr<AbstractRadiationPattern> TansformPatternGainPtr = CreateAbstractRadiationPattern(TansformPatternGainpoints);
@ -735,8 +775,8 @@ void RTPCProcessMain(long num_thread,QString TansformPatternFilePath, QString Re
else {
// 打印参数
qDebug() << "--------------Task Seting ---------------------------------------";
qDebug() << "SARImageStartTime: " << task->getSARImageStartTime() ;
qDebug() << "SARImageEndTime: " << task->getSARImageEndTime() ;
qDebug() << "SARImageStartTime: " << task->getSARImageStartTime();
qDebug() << "SARImageEndTime: " << task->getSARImageEndTime();
qDebug() << "BandWidth: " << task->getBandWidth();
qDebug() << "CenterFreq: " << task->getCenterFreq();
qDebug() << "PRF: " << task->getPRF();
@ -762,7 +802,7 @@ void RTPCProcessMain(long num_thread,QString TansformPatternFilePath, QString Re
}
else {}
std::shared_ptr<AbstractSatelliteOribtModel> SatelliteOribtModel = CreataPolyfitSatelliteOribtModel(nodes,task->getSARImageStartTime(),3); // 以成像开始时间作为 时间参考起点
std::shared_ptr<AbstractSatelliteOribtModel> SatelliteOribtModel = CreataPolyfitSatelliteOribtModel(nodes, task->getSARImageStartTime(), 3); // 以成像开始时间作为 时间参考起点
SatelliteOribtModel->setbeamAngle(task->getCenterLookAngle(), task->getIsRightLook()); // 设置天线方向图
if (nullptr == SatelliteOribtModel)
@ -778,11 +818,11 @@ void RTPCProcessMain(long num_thread,QString TansformPatternFilePath, QString Re
rtpc.setTaskSetting(task); //qDebug() << "setTaskSetting";
rtpc.setTaskFileName(simulationtaskName); //qDebug() << "setTaskFileName";
rtpc.setDEMTiffPath(demTiffPath); //qDebug() << "setDEMTiffPath";
rtpc.setLandCoverPath( LandCoverPath); //qDebug() << "setLandCoverPath";
rtpc.setHHSigmaPath( HHSigmaPath ); //qDebug() << "setHHSigmaPath";
rtpc.setHVSigmaPath( HVSigmaPath ); //qDebug() << "setHVSigmaPath";
rtpc.setVHSigmaPath( VHSigmaPath ); //qDebug() << "setVHSigmaPath";
rtpc.setVVSigmaPath( VVSigmaPath ); //qDebug() << "setVVSigmaPath";
rtpc.setLandCoverPath(LandCoverPath); //qDebug() << "setLandCoverPath";
rtpc.setHHSigmaPath(HHSigmaPath); //qDebug() << "setHHSigmaPath";
rtpc.setHVSigmaPath(HVSigmaPath); //qDebug() << "setHVSigmaPath";
rtpc.setVHSigmaPath(VHSigmaPath); //qDebug() << "setVHSigmaPath";
rtpc.setVVSigmaPath(VVSigmaPath); //qDebug() << "setVVSigmaPath";
rtpc.setOutEchoPath(OutEchoPath); //qDebug() << "setOutEchoPath";
qDebug() << "-------------- RTPC start---------------------------------------";
rtpc.Process(num_thread); // 处理程序

7
SimulationSAR/SARSatelliteSimulationAbstractCls.cpp
View File

@ -613,6 +613,13 @@ ErrorCode AbstractRadiationPattern::getGainLinear(double& theta, double& phi, do
}
double AbstractRadiationPattern::getGainLearThetaPhi(double theta, double phi)
{
double gainvlaue = 0;
this->getGainLinear(theta, phi, gainvlaue);
return gainvlaue;
}
ErrorCode AbstractRadiationPattern::setGain(double& theta, double& phi, double& GainValue)
{
this->GainMap.push_back(RadiationPatternGainPoint{ theta,phi,GainValue });

1
SimulationSAR/SARSatelliteSimulationAbstractCls.h
View File

@ -112,6 +112,7 @@ public:
virtual std::vector<RadiationPatternGainPoint> getGainList();
virtual ErrorCode getGain(double& theta, double& phi, double& GainValue);
virtual ErrorCode getGainLinear(double& theta, double& phi, double& GainValue);
double getGainLearThetaPhi(double theta, double phi);
virtual ErrorCode setGain(double& theta, double& phi, double& GainValue);
virtual ErrorCode RecontructGainMatrix(double threshold=-3);
virtual std::vector<double> getThetas();

20
SimulationSAR/SatelliteOribtModel.cpp
View File

@ -226,11 +226,17 @@ void PolyfitSatelliteOribtModel::setAntnnaAxisZ(double X, double Y, double Z)
///
ErrorCode PolyfitSatelliteOribtModel::getAntnnaDirection(SatelliteOribtNode& node)
{
bool flag = false;
double nexttime = node.time + 1e-6;
SatelliteOribtNode node1 = this->getSatelliteOribtNode(nexttime, flag);
// 1. 计算天线指向
Eigen::Vector3d axisZ0 = { -1*node.Px ,-1 * node.Py,-1 * node.Pz}; // z 轴 --波位角为0时，天线指向的反方向
Eigen::Vector3d axisX0 = { node.Vx ,node.Vy,node.Vz }; // x 轴 --飞行方向
Eigen::Vector3d axisZ0 = { -1 * node.Px ,-1 * node.Py,-1 * node.Pz }; // z 轴 --波位角为0时，天线指向的反方向
Eigen::Vector3d axisX0 = { (node1.Px - node.Px) , (node1.Py - node.Py), (node1.Pz - node.Pz) }; // x 轴 --飞行方向
Eigen::Vector3d axisY0 = axisZ0.cross(axisX0); // y 轴 --右手定则 -- 初始坐标系
double rotateAngle = this->RightLook ? -this->beamAngle : this->beamAngle; // 旋转角度 左（逆时针）：theta , 右（顺时针）： -theta
// 1.2. 根据波位角，确定卫星绕X轴-飞行轴
@ -240,7 +246,7 @@ ErrorCode PolyfitSatelliteOribtModel::getAntnnaDirection(SatelliteOribtNode& nod
axisX0=rotateMatrixBeam*axisX0;
// 1.3. 根据方位角，确定卫星绕Y轴旋转
double azangle = node.AzAngle;
Eigen::Matrix3d rotateMatrixAzAngle = rotationMatrix(axisY0, rotateAngle * d2r); // 旋转矩阵
Eigen::Matrix3d rotateMatrixAzAngle = rotationMatrix(axisY0, azangle * d2r); // 旋转矩阵
axisZ0 = rotateMatrixAzAngle * axisZ0; // 旋转矩阵
axisY0 = rotateMatrixAzAngle * axisY0;
axisX0 = rotateMatrixAzAngle * axisX0;
@ -270,9 +276,13 @@ ErrorCode PolyfitSatelliteOribtModel::getAntnnaDirection(SatelliteOribtNode& nod
ErrorCode PolyfitSatelliteOribtModel::getZeroDopplerAntDirect(SatelliteOribtNode& node)
{
bool flag = false;
double nexttime = node.time + 1e-6;
SatelliteOribtNode node1 = this->getSatelliteOribtNode(nexttime, flag);
// 1. 计算天线指向
Eigen::Vector3d axisZ0 = { -1 * node.Px ,-1 * node.Py,-1 * node.Pz }; // z 轴 --波位角为0时，天线指向的反方向
Eigen::Vector3d axisX0 = { node.Vx ,node.Vy,node.Vz }; // x 轴 --飞行方向
Eigen::Vector3d axisX0 = { (node1.Px-node.Px) , (node1.Py - node.Py) , (node1.Pz - node.Pz) }; // x 轴 --飞行方向
Eigen::Vector3d axisY0 = axisZ0.cross(axisX0); // y 轴 --右手定则 -- 初始坐标系
double rotateAngle = this->RightLook ? -this->beamAngle : this->beamAngle; // 旋转角度 左（逆时针）：theta , 右（顺时针）： -theta
@ -284,7 +294,7 @@ ErrorCode PolyfitSatelliteOribtModel::getZeroDopplerAntDirect(SatelliteOribtNode
axisX0 = rotateMatrixBeam * axisX0;
// 1.3. 根据方位角，确定卫星绕Y轴旋转
double azangle = 0;
Eigen::Matrix3d rotateMatrixAzAngle = rotationMatrix(axisY0, rotateAngle * d2r); // 旋转矩阵
Eigen::Matrix3d rotateMatrixAzAngle = rotationMatrix(axisY0, azangle * d2r); // 旋转矩阵
axisZ0 = rotateMatrixAzAngle * axisZ0; // 旋转矩阵
axisY0 = rotateMatrixAzAngle * axisY0;
axisX0 = rotateMatrixAzAngle * axisX0;

190
SimulationSAR/TBPGPU.cu
View File

@ -1,190 +0,0 @@
#include <iostream>
#include <memory>
#include <cmath>
#include <complex>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <device_launch_parameters.h>
#include "BaseConstVariable.h"
void checkCudaError(cudaError_t err, const char* msg) {
if (err != cudaSuccess) {
std::cerr << "CUDA error: " << msg << " (" << cudaGetErrorString(err) << ")" << std::endl;
exit(EXIT_FAILURE);
}
}
// CUDA 核函数
__global__ void computeDistanceAndEchoID(float* antPx, float* antPy, float* antPz,
float* img_x, float* img_y, float* img_z,
float* imgR, float* imgEchoID,
long rowcount, long colcount, long prfid,
float Rnear, float fs) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
// 确保线程索引有效
if (i < rowcount && j < colcount) {
int idx = i * colcount + j;
// 计算距离
float dx = antPx[prfid] - img_x[idx];
float dy = antPy[prfid] - img_y[idx];
float dz = antPz[prfid] - img_z[idx];
imgR[idx] = sqrt(dx * dx + dy * dy + dz * dz);
// 计算 EchoID
imgEchoID[idx] = (imgR[idx] - Rnear) / (2 * LIGHTSPEED) * fs;
}
}
/// <summary>
/// TBP GPU代码
/// </summary>
/// <param name="antpos_ptr">卫星轨道坐标</param>
/// <param name="echoArr">回波矩阵</param>
/// <param name="img_arr">图像矩阵</param>
void TBPImageGPUAlgBlock(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz, // 天线坐标
std::shared_ptr<float> antVx, std::shared_ptr<float> antVy, std::shared_ptr<float> antVz,
std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z, // 图像坐标
std::shared_ptr<std::complex<float>> echoArr,std::shared_ptr<std::complex<float>> img_arr,
float freq, float fs,float Rnear,float Rfar,
long rowcount,long colcount) {
float factorj = freq * 4 * PI / LIGHTSPEED;
std::shared_ptr<float> imgR(new float[rowcount * colcount]);
std::shared_ptr<long> imgEchoID(new long[rowcount * colcount]);
for (long prfid = 0; prfid < rowcount; prfid++) {
for (long i = 0; i < rowcount; i++) {
for (long j = 0; j < colcount; j++) {
imgR.get()[i * colcount + j] = std::sqrt(
std::pow(antPx.get()[prfid]-img_x.get()[i * colcount + j], 2)
+ std::pow(antPy.get()[prfid] - img_y.get()[i * colcount + j], 2)
+ std::pow(antPz.get()[prfid] - img_z.get()[i * colcount + j], 2)
);
imgEchoID.get()[i * colcount + j] = (imgR.get()[i * colcount + j] - Rnear) / 2 / LIGHTSPEED * fs;
}
}
}
// 分配主机内存
float* h_antPx, * h_antPy, * h_antPz, * h_img_x, * h_img_y, * h_img_z;
float* h_imgR, * h_imgEchoID;
cudaMallocHost(&h_antPx, sizeof(float) * 1); // 单个传感器的位置
cudaMallocHost(&h_antPy, sizeof(float) * 1);
cudaMallocHost(&h_antPz, sizeof(float) * 1);
cudaMallocHost(&h_img_x, sizeof(float) * rowcount * colcount);
cudaMallocHost(&h_img_y, sizeof(float) * rowcount * colcount);
cudaMallocHost(&h_img_z, sizeof(float) * rowcount * colcount);
cudaMallocHost(&h_imgR, sizeof(float) * rowcount * colcount);
cudaMallocHost(&h_imgEchoID, sizeof(float) * rowcount * colcount);
// 将数据初始化到主机内存（这里仅为示例，实际数据根据需求加载）
for (long i = 0; i < rowcount * colcount; i++) {
h_img_x[i] = 0.0f;
h_img_y[i] = 0.0f;
h_img_z[i] = 0.0f;
}
h_antPx[0] = 0.0f;
h_antPy[0] = 0.0f;
h_antPz[0] = 0.0f;
// 分配设备内存
float* d_antPx, * d_antPy, * d_antPz, * d_img_x, * d_img_y, * d_img_z;
float* d_imgR, * d_imgEchoID;
cudaMalloc(&d_antPx, sizeof(float) * 1);
cudaMalloc(&d_antPy, sizeof(float) * 1);
cudaMalloc(&d_antPz, sizeof(float) * 1);
cudaMalloc(&d_img_x, sizeof(float) * rowcount * colcount);
cudaMalloc(&d_img_y, sizeof(float) * rowcount * colcount);
cudaMalloc(&d_img_z, sizeof(float) * rowcount * colcount);
cudaMalloc(&d_imgR, sizeof(float) * rowcount * colcount);
cudaMalloc(&d_imgEchoID, sizeof(float) * rowcount * colcount);
// 将数据从主机拷贝到设备
cudaMemcpy(d_antPx, h_antPx, sizeof(float) * 1, cudaMemcpyHostToDevice);
cudaMemcpy(d_antPy, h_antPy, sizeof(float) * 1, cudaMemcpyHostToDevice);
cudaMemcpy(d_antPz, h_antPz, sizeof(float) * 1, cudaMemcpyHostToDevice);
cudaMemcpy(d_img_x, h_img_x, sizeof(float) * rowcount * colcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_img_y, h_img_y, sizeof(float) * rowcount * colcount, cudaMemcpyHostToDevice);
cudaMemcpy(d_img_z, h_img_z, sizeof(float) * rowcount * colcount, cudaMemcpyHostToDevice);
// 设置 CUDA 核函数的网格和块的尺寸
dim3 blockDim(16, 16); // 每个块处理 16x16 的像素点
dim3 gridDim((colcount + blockDim.x - 1) / blockDim.x,
(rowcount + blockDim.y - 1) / blockDim.y);
// 调用 CUDA 核函数
computeDistanceAndEchoID << <gridDim, blockDim >> > (d_antPx, d_antPy, d_antPz,
d_img_x, d_img_y, d_img_z,
d_imgR, d_imgEchoID,
rowcount, colcount, prfid,
RNEAR, LIGHTSPEED, FS);
// 等待所有设备任务完成
cudaDeviceSynchronize();
// 将结果从设备拷贝回主机
cudaMemcpy(h_imgR, d_imgR, sizeof(float) * rowcount * colcount, cudaMemcpyDeviceToHost);
cudaMemcpy(h_imgEchoID, d_imgEchoID, sizeof(float) * rowcount * colcount, cudaMemcpyDeviceToHost);
// 清理资源
cudaFree(d_antPx);
cudaFree(d_antPy);
cudaFree(d_antPz);
cudaFree(d_img_x);
cudaFree(d_img_y);
cudaFree(d_img_z);
cudaFree(d_imgR);
cudaFree(d_imgEchoID);
cudaFreeHost(h_antPx);
cudaFreeHost(h_antPy);
cudaFreeHost(h_antPz);
cudaFreeHost(h_img_x);
cudaFreeHost(h_img_y);
cudaFreeHost(h_img_z);
cudaFreeHost(h_imgR);
cudaFreeHost(h_imgEchoID);
}

88
SimulationSAR/TBPImageAlgCls.cpp
View File

@ -6,6 +6,8 @@
#include <cmath>
#include <QProgressDialog>
void TBPImageProcess(QString echofile, QString outImageFolder, QString imagePlanePath,long num_thread)
{
std::shared_ptr<EchoL0Dataset> echoL0ds(new EchoL0Dataset);
@ -63,7 +65,7 @@ std::shared_ptr<SARSimulationImageL1Dataset> TBPImageAlgCls::getImageL0()
ErrorCode TBPImageAlgCls::Process(long num_thread)
{
if (GPURUN) {
return this->ProcessGPU();
}
else {
return this->ProcessCPU(num_thread);
@ -337,36 +339,88 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
pixelZ.get()[i*colCount+ j] = Pz + AntDirectZ * R;
}
}
qDebug()<<"R: " << R;
this->L1ds->saveAntPos(antpos);
antpos.reset();
}
std::shared_ptr<std::complex<float>> imagarr(new std::complex<float>[rowCount * colCount]);
{
std::shared_ptr<std::complex<double>> Rasterarr = this->L1ds->getImageRaster();
for (long i = 0; i < pixelCount; i++) {
imagarr.get()[i] = Rasterarr.get()[i];
std::shared_ptr<std::complex<double>> echodataPtr = this->L0ds->getEchoArr();
for (long i = 0; i < rowCount; i++) {
for (long j = 0; j < colCount; j++) {
Rasterarr.get()[i * colCount + j] = Rasterarr.get()[i * colCount + j] * 0.0;
}
}
std::shared_ptr<std::complex<float>> echodata(new std::complex<float>[rowCount * colCount]);
{
std::shared_ptr<std::complex<double>> echodataPtr = this->L0ds->getEchoArr();
for (long i = 0; i < PRFCount; i++) {
for (long j = 0; j < PlusePoints; j++) {
echodata.get()[i * PlusePoints + j] = echodataPtr.get()[i * PlusePoints + j];
}
}
}
TBPImageGPUAlg(Pxs, Pys, Pzs, // ÌìÏß×ø±ê
Vxs, Vys, Vzs,
pixelX, pixelY, pixelZ, // ͼÏñ×ø±ê
echodata, imagarr, freq,fs,Rnear,Rfar,rowCount,colCount);
echodataPtr, Rasterarr,
freq, fs, Rnear, Rfar,
rowCount, colCount, this->L1ds);
qDebug() << "image writing:\t" << this->L1ds->getxmlFilePath();
return ErrorCode::SUCCESS;
}
void TBPImageAlgCls::setGPU(bool flag)
{
this->GPURUN = flag;
}
bool TBPImageAlgCls::getGPU( )
{
return this->GPURUN;
}
/// <summary>
/// TBP GPU代码
/// </summary>
/// <param name="antpos_ptr">卫星轨道坐标</param>
/// <param name="echoArr">回波矩阵</param>
/// <param name="img_arr">图像矩阵</param>
void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz, // 天线坐标
std::shared_ptr<float> antVx, std::shared_ptr<float> antVy, std::shared_ptr<float> antVz,
std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z, // 图像坐标
std::shared_ptr<std::complex<double>> echoArr, std::shared_ptr<std::complex<double>> img_arr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount, std::shared_ptr<SARSimulationImageL1Dataset> L1ds) {
float factorj = freq * 4 * PI / LIGHTSPEED;
qDebug() << "factorj:\t" << factorj;
qDebug() << "freq:\t" << freq;
qDebug() << "fs:\t" << fs;
qDebug() << "Rnear:\t" << Rnear;
qDebug() << "Rfar:\t" << Rfar;
qDebug() << "img_x:\t" << img_x.get()[0];
qDebug() << "img_y:\t" << img_y.get()[0];
qDebug() << "img_z:\t" << img_z.get()[0];
long blockline = Memory1MB * 1000 / sizeof(float) / colcount;
blockline = blockline < 10 ? 10 : blockline;
for (long startline = 0; startline < rowcount; startline = startline + blockline) {
long stepline = startline + blockline < rowcount ? blockline : rowcount - startline;
std::cout << startline << " \ " << rowcount << " "<< stepline << " start " << std::endl;
//TBPImageGPUBlock(antPx.get(), antPy.get(), antPz.get(), img_x.get(), img_y.get(), img_z.get(),
// echoArr, rowcount, colcount,
// img_arr,
// freq, fs, Rnear, Rfar, factorj, startline, stepline,
// stepline, colcount);
//std::cout << startline << " \ " << rowcount << " " << stepline << " end " << std::endl;
//L1ds->saveImageRaster(img_arr, 0, rowcount);
}
L1ds->saveImageRaster(img_arr, 0, rowcount);
L1ds->saveToXml();
}

11
SimulationSAR/TBPImageAlgCls.h
View File

@ -50,12 +50,21 @@ public:
public:
ErrorCode Process(long num_thread);
void setGPU(bool flag);
bool getGPU( );
private:
ErrorCode ProcessCPU(long num_thread);
ErrorCode ProcessGPU();
};
void TBPImageProcess(QString echofile,QString outImageFolder,QString imagePlanePath,long num_thread);
void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz, // ÌìÏß×ø±ê
std::shared_ptr<float> antVx, std::shared_ptr<float> antVy, std::shared_ptr<float> antVz,
std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z, // ͼÏñ×ø±ê
std::shared_ptr<std::complex<double>> echoArr, std::shared_ptr<std::complex<double>> img_arr,
float freq, float fs, float Rnear, float Rfar,
long rowcount, long colcount, std::shared_ptr<SARSimulationImageL1Dataset> L1ds);

BIN
script/DatasetDescription_PALSAR_Mosaic_FNF_revO.pdf Normal file
View File

Binary file not shown.

21
script/code分析/XYZ2LLA.m Normal file
View File

@ -0,0 +1,21 @@
%% WGS84直角坐标系转经纬度
% 东经正数，西经为负数
% 北纬为正数，南纬为负数
% 输入参数1：纬度；输入参数2：经度；输入参数3：高度
% 经纬度单位：度；高度单位：米
function [latitude,longitude,height]=XYZ2LLA(x,y,z)
a = 6378137.0;
b = 6356752.31424518;
c = sqrt(((a * a) - (b * b)) / (a * a));
d = sqrt(((a * a) - (b * b)) / (b * b));
p = sqrt((x * x) + (y * y));
q = atan((z * a)/ (p * b));
longitude = atan(y/x);
latitude = atan((z + (d * d) * b * sin(q)^3)/(p - (c * c) * a * cos(q)^3));
N = a / sqrt(1 - ((c * c) * sin(latitude)^2));
height = (p / cos(latitude)) - N;
longitude = longitude * 180.0 / pi+180;
latitude = latitude * 180.0 / pi;
end

226
script/code分析/mainAnayls.m Normal file
View File

@ -0,0 +1,226 @@
% n0= [-1944902.1250000000 4077521.2500000000 4487326.5000000000 ];
% n1= [-1944873.0000000000 4077511.0000000000 4487348.5000000000 ];
% n2= [-1944922.0000000000 4077512.0000000000 4487326.5000000000 ];
% n3= [-1944931.3750000000 4077531.5000000000 4487305.0000000000 ];
% n4= [-1944882.3750000000 4077530.7500000000 4487326.5000000000 ];
%
% n01 = n1 - n0, n02 = n2 - n0, n03 = n3 - n0, n04 = n4 - n0;
%
% VN=(cross(n01, n04) + cross(n04, n03) + cross(n03, n02) + cross(n02, n01)) * 0.25;
% VNorm=VN./norm(VN);
clear,clc,close all
% 打开二进制文件
fileID = fopen('D:\\Programme\\vs2022\\RasterMergeTest\\TestData\\outData\\GF3_Simulation.gpspos.data', 'rb'); % 假设二进制文件名为 'data.bin'
% 定义每个数据字段的数据类型
% 假设每个数据是双精度浮动数（8字节）
dataType = 'double';
% 读取数据：根据你的数据格式，数据长度为 16 * 8 字节（每个元素8字节）
% 例如：每16个数据项为一组（对应prf_id）
numFields = 19; % 每组数据的字段数
PRFnumRecords = 15382; % 假设有100个记录，你可以根据实际情况修改
% 读取所有记录的数据
rawData = fread(fileID, numFields * PRFnumRecords, dataType);
% 关闭文件
fclose(fileID);
% 解析数据为结构体数组
data = struct();
AntDirecX=zeros(PRFnumRecords,1);
AntDirecY=zeros(PRFnumRecords,1);
AntDirecZ=zeros(PRFnumRecords,1);
Px=zeros(PRFnumRecords,1);
Py=zeros(PRFnumRecords,1);
Pz=zeros(PRFnumRecords,1);
Vx=zeros(PRFnumRecords,1);
Vy=zeros(PRFnumRecords,1);
Vz=zeros(PRFnumRecords,1);
% 假设你的数据在 rawData 中按顺序排列
% 例如，解析 prf_id = 0, 1, 2...的相关数据
for prf_id = 0:PRFnumRecords-1
idx = prf_id * numFields + 1;
data(prf_id+1).prf_time = rawData(idx); % prf_time + imageStarttime
Px(prf_id+1) = rawData(idx+1); % sateOirbtNode.Px
Py(prf_id+1) = rawData(idx+2); % sateOirbtNode.Py
Pz(prf_id+1) = rawData(idx+3); % sateOirbtNode.Pz
Vx(prf_id+1) = rawData(idx+4); % sateOirbtNode.Vx
Vy(prf_id+1) = rawData(idx+5); % sateOirbtNode.Vy
Vz(prf_id+1) = rawData(idx+6); % sateOirbtNode.Vz
AntDirecX(prf_id+1) = rawData(idx+7); % sateOirbtNode.AntDirecX
AntDirecY(prf_id+1) = rawData(idx+8); % sateOirbtNode.AntDirecY
AntDirecZ(prf_id+1) = rawData(idx+9); % sateOirbtNode.AntDirecZ
data(prf_id+1).AVx = rawData(idx+10); % sateOirbtNode.AVx
data(prf_id+1).AVy = rawData(idx+11); % sateOirbtNode.AVy
data(prf_id+1).AVz = rawData(idx+12); % sateOirbtNode.AVz
data(prf_id+1).lon = rawData(idx+13); % outP.lon
data(prf_id+1).lat = rawData(idx+14); % outP.lat
data(prf_id+1).ati = rawData(idx+15); % outP.ati
end
% 查看解析后的数据
disp(data);
fileID = fopen('D:\\Programme\\vs2022\\RasterMergeTest\\TestData\\outData\\tmp\\ImageX.dat', 'rb'); % 假设二进制文件名为 'data.bin'
dataType = 'float';
numFields = 5400; % 每组数据的字段数
numRecords = 5400; % 假设有100个记录，你可以根据实际情况修改
demX = fread(fileID, numFields * numRecords, dataType);
demX=reshape(demX,[numRecords,numRecords]);
fclose(fileID);
fileID = fopen('D:\\Programme\\vs2022\\RasterMergeTest\\TestData\\outData\\tmp\\ImageY.dat', 'rb'); % 假设二进制文件名为 'data.bin'
dataType = 'float';
numFields = 5400; % 每组数据的字段数
numRecords = 5400; % 假设有100个记录，你可以根据实际情况修改
demY = fread(fileID, numFields * numRecords, dataType);
demY=reshape(demY,[numRecords,numRecords]);
fclose(fileID);
fileID = fopen('D:\\Programme\\vs2022\\RasterMergeTest\\TestData\\outData\\tmp\\ImageZ.dat', 'rb'); % 假设二进制文件名为 'data.bin'
dataType = 'float';
numFields = 5400; % 每组数据的字段数
numRecords = 5400; % 假设有100个记录，你可以根据实际情况修改
demZ = fread(fileID, numFields * numRecords, dataType);
demZ=reshape(demZ,[numRecords,numRecords]);
fclose(fileID);
%% 绘制成像平面
Rn=884053;
Rf=942069;
data_ant_norm=(AntDirecX.^2+AntDirecY.^2+AntDirecZ.^2).^0.5;
AntDirecX=AntDirecX./data_ant_norm;
AntDirecY=AntDirecY./data_ant_norm;
AntDirecZ=AntDirecZ./data_ant_norm;
Rlist=linspace(Rn,Rf,300);
ImageX=zeros(PRFnumRecords,300);
ImageY=zeros(PRFnumRecords,300);
ImageZ=zeros(PRFnumRecords,300);
for i=1:PRFnumRecords
ImageX(i,:)=Px(i,:) + Rlist.*AntDirecX(i,:);
ImageY(i,:)=Py(i,:) + Rlist.*AntDirecY(i,:);
ImageZ(i,:)=Pz(i,:) + Rlist.*AntDirecZ(i,:);
end
B=[AntDirecX(1,1),AntDirecY(1,1),AntDirecZ(1,1)];
A=[Px(1,1)-Px(2,1),Py(1,1)-Py(2,1),Pz(1,1)-Pz(2,1)] ;
V=[Vx(1,1),Vy(1,1),Vz(1,1)];
% 计算点积
dotProduct = dot(A,B);
% 计算模长
normA = norm(A);
normB = norm(B);
% 计算夹角的余弦值
cosTheta = acos(dotProduct / (normA * normB));
cosTheta*180/pi
tp=[-2023401.250000,4104205.500000,4428297.000000];
TP=Px.*0;
for i=1:PRFnumRecords
TP(i,:)=sqrt((Px(i,:) -tp(1)).^2+(Py(i,:) -tp(2)).^2+(Pz(i,:) -tp(3)).^2);
end
figure,
plot(TP )
PR=(Px.^2+Py.^2+Pz.^2).^0.5;
ImageR=(ImageX.^2+ImageY.^2+ImageZ.^2).^0.5;
figure,
plot(ImageR(:,1)),hold on
plot(ImageR(:,end)),hold on
plot(PR),hold on
legend(["Rn","Rf","Rs"])
a=[Px-mean(demX,"all"),Py-mean(demY,"all"),Pz-mean(demZ,"all")];
a=(a(:,1).^2+a(:,2).^2+a(:,3).^2).^0.5;
figure,
surfl(ImageX,ImageY,ImageZ),shading flat,hold on
scatter3(mean(demX,"all"),mean(demY,"all"),mean(demZ,"all"),'red','filled'),shading flat,hold on
plot3(Px,Py,Pz,'red')
xlabel("X")
ylabel("Y")
zlabel("Z")
view([-90.6000010944003 90]);
dl=sqrt((Px(1)-Px(end)).^2+(Py(1)-Py(end)).^2+(Pz(1)-Pz(end)).^2);
R=7;
ImageR=((Px(10000,1)-R-Px).^2+(Py(10000,1)-R-Py).^2+(Pz(10000,1)-R-Pz).^2).^0.5;
figure,
plot(ImageR)
P=zeros(100*100,3);
pp=1;
for i=1:100
for j=1:100
[latitude,longitude,height]=XYZ2LLA(ImageX(i,j),ImageY(i,j),ImageZ(i,j));
P(pp,:)=[latitude,longitude,height];
pp=pp+1;
end
end
% %% 计算
%
%
% Tx=-2001212.875000;
% Ty=4082207.000000;
% Tz=4471714.500000;
% R=((demX-Px(1)).^2+(demY-Py(1)).^2+(demZ-Pz(1)).^2).^0.5;
% figure,
% plot(R)
%
% imgX=ImageX(3000,150);
% imgY=ImageY(3000,150);
% imgZ=ImageZ(3000,150);
% R1= sqrt((ImageX(1,1)-Px(:,1)).^2+(ImageY(1,1)-Py(:,1)).^2+(ImageZ(1,1)-Pz(:,1)).^2);
% figure,
% plot(R1)
%
%
%
% Timgx=abs(ImageX-Tx);
% Timgy=abs(ImageY-Ty);
% Timgz=abs(ImageZ-Tz);
%
% [latitude,longitude,height]=XYZ2LLA(ImageX(3000,150),ImageY(3000,150),ImageZ(3000,150))
%
%
%

265
script/datasetCreate.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

141
script/gdalClipRaster.ipynb Normal file
View File

@ -0,0 +1,141 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 4,
"id": "7d84d0c7-8625-432c-acd7-89bdb0e66a43",
"metadata": {},
"outputs": [],
"source": [
"from osgeo import gdal\n",
"import os\n",
"\n",
"def crop_raster(input_raster_path, output_raster_path, bounds):\n",
" \"\"\"\n",
" 使用 GDAL 裁剪栅格影像到指定的四至范围。\n",
"\n",
" :param input_raster_path: 输入栅格影像路径\n",
" :param output_raster_path: 输出裁剪后影像的路径\n",
" :param bounds: 裁剪范围，格式为 (minx, miny, maxx, maxy)\n",
" \"\"\"\n",
" minx, miny, maxx, maxy = bounds\n",
"\n",
" # 打开输入影像\n",
" dataset = gdal.Open(input_raster_path)\n",
" if dataset is None:\n",
" print(f\"无法打开输入影像: {input_raster_path}\")\n",
" return\n",
"\n",
" # 获取影像的投影和变换信息\n",
" geo_transform = dataset.GetGeoTransform()\n",
" projection = dataset.GetProjection()\n",
"\n",
" # 计算裁剪区域在影像中的像素坐标\n",
" ulx = int((minx - geo_transform[0]) / geo_transform[1]) # 左上角 x\n",
" uly = int((maxy - geo_transform[3]) / geo_transform[5]) # 左上角 y\n",
" lrx = int((maxx - geo_transform[0]) / geo_transform[1]) # 右下角 x\n",
" lry = int((miny - geo_transform[3]) / geo_transform[5]) # 右下角 y\n",
"\n",
" # 检查像素坐标是否在影像范围内\n",
" if ulx < 0 or uly < 0 or lrx > dataset.RasterXSize or lry > dataset.RasterYSize:\n",
" print(\"裁剪范围超出影像范围！\")\n",
" return\n",
"\n",
" # 创建输出影像\n",
" driver = gdal.GetDriverByName('GTiff')\n",
" out_dataset = driver.Create(output_raster_path, lrx - ulx, lry - uly, dataset.RasterCount, gdal.GDT_Float32)\n",
"\n",
" # 设置输出影像的投影和变换信息\n",
" out_geo_transform = (minx, geo_transform[1], geo_transform[2], maxy, geo_transform[4], geo_transform[5])\n",
" out_dataset.SetGeoTransform(out_geo_transform)\n",
" out_dataset.SetProjection(projection)\n",
"\n",
" # 将裁剪后的数据写入输出影像\n",
" for i in range(1, dataset.RasterCount + 1):\n",
" band = dataset.GetRasterBand(i)\n",
" data = band.ReadAsArray(ulx, uly, lrx - ulx, lry - uly)\n",
" out_dataset.GetRasterBand(i).WriteArray(data)\n",
"\n",
" # 清理\n",
" out_dataset.FlushCache()\n",
" del out_dataset\n",
" del dataset\n",
"\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "3ac1f4a3-38f8-4c50-8bad-36852fcb18f6",
"metadata": {},
"outputs": [],
"source": [
"# 输入栅格影像路径\n",
"input_raster = \"115E39N_COP30.tif\"\n",
"\n",
"# 输出裁剪后影像的路径\n",
"output_raster = \"115E39N_COP30_clip.tif\"\n",
"\n",
"\n",
"# 定义裁剪范围 (minx, miny, maxx, maxy)\n",
"bounds =(115.5,43.5,117.0,45.0) # 示例范围\n",
"\n",
"# 裁剪影像\n",
"crop_raster(input_raster, output_raster, bounds)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "65dcf0f1-55e7-4dfa-8988-412ca0e6c41a",
"metadata": {},
"outputs": [],
"source": [
"(xmin, ymin, xmax, ymax)=(115.5,43.5,117.0,45.0)\n",
"# 定义裁剪范围 (minx, miny, maxx, maxy)\n",
"bounds = (xmin, ymin, xmax, ymax) # 例如 (100000, 200000, 110000, 210000)\n",
"# 裁剪影像\n",
"clip_raster(input_raster, output_raster, bounds)\n",
"print(f\"影像已裁剪并保存到: {output_raster}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ea8bcfa3-7612-477d-bf0b-b4e3537fcf33",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"id": "407416e7-729d-4d7e-ae09-318cdfce7739",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.19"
}
},
"nbformat": 4,
"nbformat_minor": 5
}