407 lines
16 KiB
C++
407 lines
16 KiB
C++
#include "ImageNetOperator.h"
|
||
#include "LogInfoCls.h"
|
||
#include "PrintMsgToQDebug.h"
|
||
#include <QDebug>
|
||
#include "ImageOperatorBase.h"
|
||
#include "GPUBaseTool.h"
|
||
#include "GPUBPImageNet.cuh"
|
||
#include "BaseTool.h"
|
||
|
||
|
||
void InitCreateImageXYZProcess(QString& outImageLLPath, QString& outImageXYZPath, QString& InEchoGPSDataPath,
|
||
double& NearRange, double& RangeResolution, int64_t& RangeNum)
|
||
{
|
||
qDebug() << "---------------------------------------------------------------------------------";
|
||
qDebug() << u8"创建粗成像平面斜距投影网格";
|
||
gdalImage antimg(InEchoGPSDataPath);
|
||
qDebug() << u8"1. 回波GPS坐标点文件参数：\t";
|
||
qDebug() << u8"文件路径：\t" << InEchoGPSDataPath;
|
||
qDebug() << u8"GPS 点数：\t" << antimg.height;
|
||
qDebug() << u8"文件列数：\t" << antimg.width;
|
||
qDebug() << u8"2.斜距网格参数：";
|
||
qDebug() << u8"近距离：\t" << NearRange;
|
||
qDebug() << u8"分辨率：\t" << RangeResolution;
|
||
qDebug() << u8"网格点数：\t" << RangeNum;
|
||
qDebug() << u8"3.输出文件参数：";
|
||
gdalImage outimgll = CreategdalImageDouble(outImageLLPath, antimg.height,RangeNum, 3,true,true);
|
||
gdalImage outimgxyz = CreategdalImageDouble(outImageXYZPath, antimg.height, RangeNum, 3, true, true);
|
||
qDebug() << u8"成像平面文件（经纬度）网格路径：\t" << outImageLLPath;
|
||
qDebug() << u8"成像平面文件（XYZ）网格路径：\t" << outImageXYZPath;
|
||
qDebug() << u8"4.开始创建成像网格XYZ";
|
||
long prfcount = antimg.height;
|
||
long rangeNum = RangeNum;
|
||
double Rnear = NearRange;
|
||
double dx = RangeResolution;
|
||
|
||
long blockRangeCount = Memory1GB / sizeof(double) / 4 / prfcount *6;
|
||
blockRangeCount = blockRangeCount < 1 ? 1 : blockRangeCount;
|
||
|
||
std::shared_ptr<double> Pxs((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> Pys((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> Pzs((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> AntDirectX((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> AntDirectY((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> AntDirectZ((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
|
||
{
|
||
long colnum = 19;
|
||
std::shared_ptr<double> antpos =readDataArr<double>(antimg,0,0,prfcount, colnum,1,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
|
||
double time = 0;
|
||
double Px = 0;
|
||
double Py = 0;
|
||
double Pz = 0;
|
||
for (long i = 0; i < prfcount; i++) {
|
||
|
||
Pxs.get()[i] = antpos.get()[i * 19 + 1]; // 卫星坐标
|
||
Pys.get()[i] = antpos.get()[i * 19 + 2];
|
||
Pzs.get()[i] = antpos.get()[i * 19 + 3];
|
||
AntDirectX.get()[i] = antpos.get()[i * 19 + 13];// zero doppler
|
||
AntDirectY.get()[i] = antpos.get()[i * 19 + 14];
|
||
AntDirectZ.get()[i] = antpos.get()[i * 19 + 15];
|
||
|
||
double NormAnt = std::sqrt(AntDirectX.get()[i] * AntDirectX.get()[i] +
|
||
AntDirectY.get()[i] * AntDirectY.get()[i] +
|
||
AntDirectZ.get()[i] * AntDirectZ.get()[i]);
|
||
AntDirectX.get()[i] = AntDirectX.get()[i] / NormAnt;
|
||
AntDirectY.get()[i] = AntDirectY.get()[i] / NormAnt;
|
||
AntDirectZ.get()[i] = AntDirectZ.get()[i] / NormAnt;// 归一化
|
||
}
|
||
antpos.reset();
|
||
}
|
||
|
||
std::shared_ptr<double> d_Pxs((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
|
||
std::shared_ptr<double> d_Pys((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
|
||
std::shared_ptr<double> d_Pzs((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
|
||
std::shared_ptr<double> d_AntDirectX((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
|
||
std::shared_ptr<double> d_AntDirectY((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
|
||
std::shared_ptr<double> d_AntDirectZ((double*)mallocCUDADevice(sizeof(double) * prfcount), FreeCUDADevice);
|
||
|
||
HostToDevice(Pxs.get(), d_Pxs.get(), sizeof(double) * prfcount);
|
||
HostToDevice(Pys.get(), d_Pys.get(), sizeof(double) * prfcount);
|
||
HostToDevice(Pzs.get(), d_Pzs.get(), sizeof(double) * prfcount);
|
||
HostToDevice(AntDirectX.get(), d_AntDirectX.get(), sizeof(double) * prfcount);
|
||
HostToDevice(AntDirectY.get(), d_AntDirectY.get(), sizeof(double) * prfcount);
|
||
HostToDevice(AntDirectZ.get(), d_AntDirectZ.get(), sizeof(double) * prfcount);
|
||
|
||
|
||
for (long startcolidx = 0; startcolidx < RangeNum; startcolidx = startcolidx + blockRangeCount) {
|
||
|
||
long tempechocol = blockRangeCount;
|
||
if (startcolidx + tempechocol >= RangeNum) {
|
||
tempechocol = RangeNum - startcolidx;
|
||
}
|
||
qDebug() << " imgxyz :\t" << startcolidx << "\t-\t" << startcolidx + tempechocol << " / " << RangeNum;
|
||
|
||
std::shared_ptr<double> demx = readDataArr<double>(outimgxyz, 0, startcolidx, prfcount, tempechocol, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
|
||
std::shared_ptr<double> demy = readDataArr<double>(outimgxyz, 0, startcolidx, prfcount, tempechocol, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
|
||
std::shared_ptr<double> demz = readDataArr<double>(outimgxyz, 0, startcolidx, prfcount, tempechocol, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
|
||
|
||
std::shared_ptr<double> h_demx((double*)mallocCUDAHost(sizeof(double) * prfcount * tempechocol), FreeCUDAHost);
|
||
std::shared_ptr<double> h_demy((double*)mallocCUDAHost(sizeof(double) * prfcount * tempechocol), FreeCUDAHost);
|
||
std::shared_ptr<double> h_demz((double*)mallocCUDAHost(sizeof(double) * prfcount * tempechocol), FreeCUDAHost);
|
||
|
||
#pragma omp parallel for
|
||
for (long ii = 0; ii < prfcount; ii++) {
|
||
for (long jj = 0; jj < tempechocol; jj++) {
|
||
h_demx.get()[ii * tempechocol + jj] = demx.get()[ii * tempechocol + jj];
|
||
h_demy.get()[ii * tempechocol + jj] = demy.get()[ii * tempechocol + jj];
|
||
h_demz.get()[ii * tempechocol + jj] = demz.get()[ii * tempechocol + jj];
|
||
}
|
||
}
|
||
|
||
std::shared_ptr<double> d_demx((double*)mallocCUDADevice(sizeof(double) * prfcount * tempechocol), FreeCUDADevice);
|
||
std::shared_ptr<double> d_demy((double*)mallocCUDADevice(sizeof(double) * prfcount * tempechocol), FreeCUDADevice);
|
||
std::shared_ptr<double> d_demz((double*)mallocCUDADevice(sizeof(double) * prfcount * tempechocol), FreeCUDADevice);
|
||
|
||
HostToDevice(h_demx.get(), d_demx.get(), sizeof(double) * prfcount * tempechocol);
|
||
HostToDevice(h_demy.get(), d_demy.get(), sizeof(double) * prfcount * tempechocol);
|
||
HostToDevice(h_demz.get(), d_demz.get(), sizeof(double) * prfcount * tempechocol);
|
||
|
||
|
||
TIMEBPCreateImageGrid(
|
||
d_Pxs.get(), d_Pys.get(), d_Pzs.get(),
|
||
d_AntDirectX.get(), d_AntDirectY.get(), d_AntDirectZ.get(),
|
||
d_demx.get(), d_demy.get(), d_demz.get(),
|
||
prfcount, tempechocol, 0,
|
||
Rnear + dx * startcolidx, dx // 更新最近修读
|
||
);
|
||
|
||
DeviceToHost(h_demx.get(), d_demx.get(), sizeof(double) * prfcount * tempechocol);
|
||
DeviceToHost(h_demy.get(), d_demy.get(), sizeof(double) * prfcount * tempechocol);
|
||
DeviceToHost(h_demz.get(), d_demz.get(), sizeof(double) * prfcount * tempechocol);
|
||
|
||
#pragma omp parallel for
|
||
for (long ii = 0; ii < prfcount; ii++) {
|
||
for (long jj = 0; jj < tempechocol; jj++) {
|
||
demx.get()[ii * tempechocol + jj] = h_demx.get()[ii * tempechocol + jj];
|
||
demy.get()[ii * tempechocol + jj] = h_demy.get()[ii * tempechocol + jj];
|
||
demz.get()[ii * tempechocol + jj] = h_demz.get()[ii * tempechocol + jj];
|
||
}
|
||
}
|
||
|
||
outimgxyz.saveImage(demx, 0, startcolidx, prfcount, tempechocol, 1);
|
||
outimgxyz.saveImage(demy, 0, startcolidx, prfcount, tempechocol, 2);
|
||
outimgxyz.saveImage(demz, 0, startcolidx, prfcount, tempechocol, 3);
|
||
|
||
// 将XYZ转换为经纬度
|
||
std::shared_ptr<double> demllx = readDataArr<double>(outimgll, 0, startcolidx, prfcount, tempechocol, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
|
||
std::shared_ptr<double> demlly = readDataArr<double>(outimgll, 0, startcolidx, prfcount, tempechocol, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
|
||
std::shared_ptr<double> demllz = readDataArr<double>(outimgll, 0, startcolidx, prfcount, tempechocol, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
|
||
|
||
#pragma omp parallel for
|
||
for (long ii = 0; ii < prfcount; ii++) {
|
||
for (long jj = 0; jj < tempechocol; jj++) {
|
||
double x = demx.get()[ii * tempechocol + jj];
|
||
double y = demy.get()[ii * tempechocol + jj];
|
||
double z = demz.get()[ii * tempechocol + jj];
|
||
Landpoint point;
|
||
XYZ2BLH_FixedHeight(x, y, z, 0, point);
|
||
demllx.get()[ii * tempechocol + jj] = point.lon;
|
||
demlly.get()[ii * tempechocol + jj] = point.lat;
|
||
demllz.get()[ii * tempechocol + jj] = point.ati;
|
||
}
|
||
}
|
||
|
||
|
||
outimgll.saveImage(demllx, 0, startcolidx, prfcount, tempechocol, 1);
|
||
outimgll.saveImage(demlly, 0, startcolidx, prfcount, tempechocol, 2);
|
||
outimgll.saveImage(demllz, 0, startcolidx, prfcount, tempechocol, 3);
|
||
|
||
}
|
||
|
||
qDebug() << u8"6.保存成像网格结果";
|
||
qDebug() << "---------------------------------------------------------------------------------";
|
||
}
|
||
|
||
bool OverlapCheck(QString& ImageLLPath, QString& ImageDEMPath)
|
||
{
|
||
// 检查DEM是否是WGS84坐标系
|
||
//long demEPSG = GetEPSGFromRasterFile(ImageDEMPath);
|
||
//if (demEPSG != 4326) {
|
||
// qDebug() << u8"DEM坐标系不是WGS84坐标系,ESPG:"<< demEPSG;
|
||
// return false;
|
||
//}
|
||
|
||
gdalImage demimg(ImageDEMPath);
|
||
gdalImage imgll(ImageLLPath);
|
||
|
||
long imgheight = imgll.height;
|
||
long imgwidth = imgll.width;
|
||
Eigen::MatrixXd imglonArr = imgll.getData(0, 0, imgheight, imgwidth, 1);
|
||
Eigen::MatrixXd imglatArr = imgll.getData(0, 0, imgheight, imgwidth, 2);
|
||
|
||
// 打印范围
|
||
qDebug() << u8"影像范围：";
|
||
qDebug() << u8"最小经度：\t" << imglonArr.minCoeff();
|
||
qDebug() << u8"最大经度：\t" << imglonArr.maxCoeff();
|
||
qDebug() << u8"最小纬度：\t" << imglatArr.minCoeff();
|
||
qDebug() << u8"最大纬度：\t" << imglatArr.maxCoeff();
|
||
qDebug() << u8"DEM范围：";
|
||
RasterExtend demextend = demimg.getExtend();
|
||
qDebug() << u8"最小经度：\t" << demextend.min_x;
|
||
qDebug() << u8"最大经度：\t" << demextend.max_x;
|
||
qDebug() << u8"最小纬度：\t" << demextend.min_y;
|
||
qDebug() << u8"最大纬度：\t" << demextend.max_y;
|
||
qDebug() << u8"影像大小：\t" << demimg.height << " * " << demimg.width;
|
||
|
||
|
||
|
||
for (long i = 0; i < imgheight; i++)
|
||
{
|
||
for (long j = 0; j < imgwidth; j++)
|
||
{
|
||
double lon = imglonArr(i, j); // X
|
||
double lat = imglatArr(i, j); // Y
|
||
Landpoint point = demimg.getRow_Col(lon, lat);
|
||
imglonArr(i, j) = point.lon;
|
||
imglatArr(i, j) = point.lat;
|
||
}
|
||
}
|
||
|
||
double minX = imglonArr.minCoeff();
|
||
double maxX = imglonArr.maxCoeff();
|
||
double minY = imglatArr.minCoeff();
|
||
double maxY = imglatArr.maxCoeff();
|
||
|
||
//打印范围
|
||
qDebug() << u8"dem 的范围：";
|
||
qDebug() << u8"minX:"<<minX<<"\t"<<demimg.width;
|
||
qDebug() << u8"maxX:"<<maxX << "\t" << demimg.width;
|
||
qDebug() << u8"minY:"<<minY << "\t" << demimg.height;
|
||
qDebug() << u8"maxY:"<<maxY << "\t" << demimg.height;
|
||
qDebug() << u8"图像行列：\t" << demimg.height << " , " << demimg.width;
|
||
|
||
if (minX<1 || maxX>demimg.width - 1 || minY<1 || maxY>demimg.height - 1) {
|
||
return false;
|
||
}
|
||
else {
|
||
return true;
|
||
}
|
||
|
||
|
||
|
||
}
|
||
|
||
bool GPSPointsNumberEqualCheck(QString& ImageLLPath, QString& InEchoGPSDataPath)
|
||
{
|
||
|
||
gdalImage antimg(InEchoGPSDataPath);
|
||
gdalImage imgll(ImageLLPath);
|
||
return antimg.height == imgll.height;
|
||
|
||
}
|
||
|
||
void InterploateAtiByRefDEM(QString& ImageLLPath, QString& ImageDEMPath, QString& outImageLLAPath, QString& InEchoGPSDataPath)
|
||
{
|
||
gdalImage demimg(ImageDEMPath);
|
||
gdalImage imgll(ImageLLPath);
|
||
gdalImage outimgll = CreategdalImageDouble(outImageLLAPath, imgll.height, imgll.width, 4, true, true); // 经度、纬度、高程、斜距
|
||
|
||
long imgheight = imgll.height;
|
||
long imgwidth = imgll.width;
|
||
|
||
Eigen::MatrixXd imglonArr = imgll.getData(0, 0, imgheight, imgwidth, 1);
|
||
Eigen::MatrixXd imglatArr = imgll.getData(0, 0, imgheight, imgwidth, 2);
|
||
Eigen::MatrixXd demArr = demimg.getData(0, 0, demimg.height, demimg.width, 1);
|
||
Eigen::MatrixXd imgatiArr = Eigen::MatrixXd::Zero(imgheight, imgwidth);
|
||
Eigen::MatrixXd imgRArr = Eigen::MatrixXd::Zero(imgheight, imgwidth);
|
||
|
||
outimgll.saveImage(imglonArr, 0, 0, 1);
|
||
outimgll.saveImage(imglatArr, 0, 0, 2);
|
||
|
||
|
||
double minX = imglonArr.minCoeff();
|
||
double maxX = imglonArr.maxCoeff();
|
||
double minY = imglatArr.minCoeff();
|
||
double maxY = imglatArr.maxCoeff();
|
||
//打印范围
|
||
qDebug() << u8"dem 的范围：";
|
||
qDebug() << u8"minX:" << minX << "\t" << demimg.width;
|
||
qDebug() << u8"maxX:" << maxX << "\t" << demimg.width;
|
||
qDebug() << u8"minY:" << minY << "\t" << demimg.height;
|
||
qDebug() << u8"maxY:" << maxY << "\t" << demimg.height;
|
||
qDebug() << u8"图像行列：\t" << demimg.height << " , " << demimg.width;
|
||
|
||
|
||
|
||
for (long i = 0; i < imgheight; i++) {
|
||
//printf("\rprocess:%f precent\t\t\t",i*100.0/imgheight);
|
||
for (long j = 0; j < imgwidth; j++) {
|
||
double lon = imglonArr(i, j);
|
||
double lat = imglatArr(i, j);
|
||
Landpoint point = demimg.getRow_Col(lon, lat);
|
||
Landpoint p0, p11, p21, p12, p22;
|
||
p0.lon = point.lon;
|
||
p0.lat = point.lat;
|
||
|
||
p11.lon = floor(p0.lon);
|
||
p11.lat = floor(p0.lat);
|
||
p11.ati = demArr(long(p11.lat), long(p11.lon));
|
||
|
||
p12.lon = ceil(p0.lon);
|
||
p12.lat = floor(p0.lat);
|
||
p12.ati = demArr(long(p12.lat), long(p12.lon));
|
||
|
||
p21.lon = floor(p0.lon);
|
||
p21.lat = ceil(p0.lat);
|
||
p21.ati = demArr(long(p21.lat), long(p21.lon));
|
||
|
||
p22.lon = ceil(p0.lon);
|
||
p22.lat = ceil(p0.lat);
|
||
p22.ati = demArr(long(p22.lat), long(p22.lon));
|
||
|
||
p0.lon = p0.lon - p11.lon;
|
||
p0.lat = p0.lat - p11.lat;
|
||
|
||
p12.lon = p12.lon - p11.lon;
|
||
p12.lat = p12.lat - p11.lat;
|
||
|
||
p21.lon = p21.lon - p11.lon;
|
||
p21.lat = p21.lat - p11.lat;
|
||
|
||
p22.lon = p22.lon - p11.lon;
|
||
p22.lat = p22.lat - p11.lat;
|
||
|
||
p11.lon = p11.lon - p11.lon;
|
||
p11.lat = p11.lat - p11.lat;
|
||
|
||
p0.ati=Bilinear_interpolation(p0, p11, p21, p12, p22);
|
||
imgatiArr(i, j) = p0.ati;
|
||
|
||
}
|
||
}
|
||
outimgll.saveImage(imgatiArr, 0, 0, 3);
|
||
qDebug() << u8"计算每个点的斜距值";
|
||
|
||
|
||
gdalImage antimg(InEchoGPSDataPath);
|
||
qDebug() << u8"1. 回波GPS坐标点文件参数：\t";
|
||
qDebug() << u8"文件路径：\t" << InEchoGPSDataPath;
|
||
qDebug() << u8"GPS 点数：\t" << antimg.height;
|
||
qDebug() << u8"文件列数：\t" << antimg.width;
|
||
|
||
long prfcount = antimg.height;
|
||
|
||
|
||
std::shared_ptr<double> Pxs((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> Pys((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> Pzs((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> AntDirectX((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> AntDirectY((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
std::shared_ptr<double> AntDirectZ((double*)mallocCUDAHost(sizeof(double) * prfcount), FreeCUDAHost);
|
||
|
||
{
|
||
long colnum = 19;
|
||
std::shared_ptr<double> antpos = readDataArr<double>(antimg, 0, 0, prfcount, colnum, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
|
||
double time = 0;
|
||
double Px = 0;
|
||
double Py = 0;
|
||
double Pz = 0;
|
||
for (long i = 0; i < prfcount; i++) {
|
||
|
||
Pxs.get()[i] = antpos.get()[i * 19 + 1]; // 卫星坐标
|
||
Pys.get()[i] = antpos.get()[i * 19 + 2];
|
||
Pzs.get()[i] = antpos.get()[i * 19 + 3];
|
||
AntDirectX.get()[i] = antpos.get()[i * 19 + 13];// zero doppler
|
||
AntDirectY.get()[i] = antpos.get()[i * 19 + 14];
|
||
AntDirectZ.get()[i] = antpos.get()[i * 19 + 15];
|
||
|
||
double NormAnt = std::sqrt(AntDirectX.get()[i] * AntDirectX.get()[i] +
|
||
AntDirectY.get()[i] * AntDirectY.get()[i] +
|
||
AntDirectZ.get()[i] * AntDirectZ.get()[i]);
|
||
AntDirectX.get()[i] = AntDirectX.get()[i] / NormAnt;
|
||
AntDirectY.get()[i] = AntDirectY.get()[i] / NormAnt;
|
||
AntDirectZ.get()[i] = AntDirectZ.get()[i] / NormAnt;// 归一化
|
||
}
|
||
antpos.reset();
|
||
}
|
||
|
||
|
||
|
||
#pragma omp parallel for
|
||
for (long prfid = 0; prfid < prfcount; prfid++) {
|
||
double Px = Pxs.get()[prfid];
|
||
double Py = Pys.get()[prfid];
|
||
double Pz = Pzs.get()[prfid];
|
||
double R = 0;
|
||
Landpoint LLA = {};
|
||
Point3 XYZ = {};
|
||
for (long j = 0; j < imgwidth; j++) {
|
||
LLA.lon = imglonArr(prfid, j);
|
||
LLA.lat = imglatArr(prfid, j);
|
||
LLA.ati = imgatiArr(prfid, j);
|
||
|
||
LLA2XYZ(LLA, XYZ);
|
||
|
||
R = sqrt(pow(Px - XYZ.x, 2) +
|
||
pow(Py - XYZ.y, 2) +
|
||
pow(Pz - XYZ.z, 2));
|
||
imgRArr(prfid, j) = R;
|
||
}
|
||
}
|
||
|
||
outimgll.saveImage(imgRArr, 0, 0, 4);
|
||
|
||
qDebug() << u8"插值完成";
|
||
}
|