修改DEM->XYZ的分块计算部分

Release-dev
陈增辉 2025-05-22 21:04:45 +08:00
parent fc4bdf7904
commit 70c4a62129
2 changed files with 121 additions and 84 deletions

View File

@ -1415,31 +1415,28 @@ void clipRaster(QString inRasterPath, QString outRasterPath, long minRow, long m
ErrorCode DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QString demsloperPath)
{
omp_lock_t lock;
omp_init_lock(&lock);
QString DEMPath = dempath;
QString XYZPath = demxyzpath;
QString SLOPERPath = demsloperPath;
gdalImage demds(DEMPath);
gdalImage demxyz = CreategdalImageDouble(XYZPath, demds.height, demds.width, 3, demds.gt, demds.projection, true, true, false);// X,Y,Z
// 分块计算并转换为XYZ
int64_t linecount = 0;
qDebug() << u8"start dem (lon,lat,ati) -> dem [X,Y,Z]";
int64_t line_invert= Memory1MB / 8.0 / demds.width * 200;
//Eigen::MatrixXd demArr = demds.getData(0, 0, demds.height, demds.width, 1);
//Eigen::MatrixXd demR = demArr;
#pragma omp parallel for
for (int64_t max_rows_ids = 0; max_rows_ids < demds.height; max_rows_ids = max_rows_ids + line_invert) {
double R = 0;
double dem_row = 0, dem_col = 0, dem_alt = 0;
long line_invert = Memory1MB / 8.0 / demds.width * 1000;
for (int max_rows_ids = 0; max_rows_ids < demds.height; max_rows_ids = max_rows_ids + line_invert) {
qDebug() << "dem -> XYZ [" << max_rows_ids*100.0/ demds.height << "] %" ;
Eigen::MatrixXd demdata = demds.getData(max_rows_ids, 0, line_invert, demds.width, 1);
long temp_line_invert = max_rows_ids + line_invert < demds.height ? line_invert : demds.height - max_rows_ids;
double R = 0;
double dem_row = 0, dem_col = 0, dem_alt = 0;
Eigen::MatrixXd demdata = demds.getData(max_rows_ids, 0, temp_line_invert, demds.width, 1);
Eigen::MatrixXd xyzdata_x = demdata.array() * 0;
Eigen::MatrixXd xyzdata_y = demdata.array() * 0;
Eigen::MatrixXd xyzdata_z = demdata.array() * 0;
@ -1447,7 +1444,7 @@ ErrorCode DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QStri
int datarows = demdata.rows();
int datacols = demdata.cols();
#pragma omp parallel for
for (int i = 0; i < datarows; i++) {
Landpoint LandP{ 0,0,0 };
Point3 GERpoint{ 0,0,0 };
@ -1463,45 +1460,48 @@ ErrorCode DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QStri
xyzdata_z(i, j) = GERpoint.z;
}
}
omp_set_lock(&lock);
demxyz.saveImage(xyzdata_x, max_rows_ids, 0, 1);
demxyz.saveImage(xyzdata_y, max_rows_ids, 0, 2);
demxyz.saveImage(xyzdata_z, max_rows_ids, 0, 3);
linecount = linecount + temp_line_invert;
qDebug() << "dem -> XYZ [" << linecount * 100.0 / demds.height << "] %";
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
}
// 计算坡向角
qDebug() << u8"start dem (lon,lat,ati) -> dem Sloper[X,Y,Z]";
gdalImage demsloperxyz = CreategdalImageDouble(SLOPERPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true, false);// X,Y,Z,cosangle
linecount = 0;
line_invert = Memory1MB / 8.0 / demds.width * 200;;
#pragma omp parallel for
for (int64_t start_ids = 1; start_ids < demds.height; start_ids = start_ids + line_invert) {
line_invert = 1000;
long start_ids = 0;
long dem_rows = 0, dem_cols = 0;
for (start_ids = 1; start_ids < demds.height; start_ids = start_ids + line_invert) {
Eigen::MatrixXd demdata = demds.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
long temp_line_invert = start_ids + line_invert < demds.height ? line_invert : demds.height - start_ids;
long dem_rows = 0, dem_cols = 0;
long startlineid = start_ids;
Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 2);
Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 3);
Eigen::MatrixXd demsloper_angle = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 4);
Eigen::MatrixXd demdata = demds.getData(start_ids - 1, 0, temp_line_invert + 2, demxyz.width, 1);
Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids , 0, temp_line_invert , demxyz.width, 1);
Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids , 0, temp_line_invert , demxyz.width, 2);
Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids , 0, temp_line_invert , demxyz.width, 3);
Eigen::MatrixXd demsloper_angle = demsloperxyz.getData(start_ids , 0, temp_line_invert , demxyz.width, 4);
dem_rows = demsloper_x.rows();
dem_cols = demsloper_x.cols();
#pragma omp parallel for
for (long i = 1; i < dem_rows - 1; i++) {
for (long i = 0; i < dem_rows ; i++) {
Landpoint p0, p1, p2, p3, p4, pslopeVector, pp;
Vector3D slopeVector;
double sloperAngle = 0;
Vector3D Zaxis = { 0,0,1 };
double rowidx = 0, colidx = 0;
for (long j = 1; j < dem_cols - 1; j++) {
rowidx = i + startlineid;
colidx = j;
for (long j = 0; j < dem_cols ; j++) {
rowidx = i + startlineid+1;
colidx = j+1;
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);
@ -1523,16 +1523,19 @@ ErrorCode DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QStri
}
}
omp_set_lock(&lock);
demsloperxyz.saveImage(demsloper_x, start_ids - 1, 0, 1);
demsloperxyz.saveImage(demsloper_y, start_ids - 1, 0, 2);
demsloperxyz.saveImage(demsloper_z, start_ids - 1, 0, 3);
demsloperxyz.saveImage(demsloper_angle, start_ids - 1, 0, 4);
linecount = linecount + temp_line_invert;
qDebug() << "dem -> Sloper [" << linecount * 100.0 / demds.height << "] %";
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
}
omp_destroy_lock(&lock); // 劫伙拷斤拷
return ErrorCode::SUCCESS;
}

View File

@ -488,6 +488,9 @@ ErrorCode GF3RDCreateLookTable(QString inxmlPath, QString indemPath, QString out
ErrorCode GF3OrthSLC( QString inRasterPath, QString inlooktablePath, QString outRasterPath)
{
gdalImage slcRaster(inRasterPath);//
gdalImage looktableRaster(inlooktablePath);//
gdalImage outRaster = CreategdalImage(outRasterPath, looktableRaster.height, looktableRaster.width, 1, looktableRaster.gt, looktableRaster.projection, true, true);// X,Y,Z
@ -497,59 +500,90 @@ ErrorCode GF3OrthSLC( QString inRasterPath, QString inlooktablePath, QString out
qDebug() << u8"look table size is not same as outRaster size"<< looktableRaster.height <<"!="<<outRaster.height<<" "<<looktableRaster.width<<"!="<<outRaster.width;
return ErrorCode::FAIL;
}
Eigen::MatrixXd slcImg = slcRaster.getData(0, 0, slcRaster.height, slcRaster.width, 1);
Eigen::MatrixXd sar_r = looktableRaster.getData(0, 0, looktableRaster.height, looktableRaster.width, 1);
Eigen::MatrixXd sar_c = looktableRaster.getData(0, 0, looktableRaster.height, looktableRaster.width, 2);
Eigen::MatrixXd outImg = outRaster.getData(0, 0, outRaster.height, outRaster.width, 1);
long lookRows = sar_r.rows();
long lookCols = sar_r.cols();
long slcRows = slcImg.rows();
long slcCols = slcImg.cols();
long lastr = 0;
long nextr = 0;
long lastc = 0;
long nextc = 0;
double Bileanervalue = 0;
// 插值
Landpoint p0{ 0,0,0 }, p11{ 0,0,0 }, p21{ 0,0,0 }, p12{ 0,0,0 }, p22{ 0,0,0 }, p{ 0,0,0 };
QProgressDialog progressDialog(u8"正射图像生成中", u8"终止", 0, lookRows);
omp_lock_t lock;
omp_init_lock(&lock);
QProgressDialog progressDialog(u8"正射图像生成中", u8"终止", 0, looktableRaster.height);
progressDialog.setWindowTitle(u8"ÕýÉäͼÏñÉú³ÉÖÐ");
progressDialog.setWindowModality(Qt::WindowModal);
progressDialog.setAutoClose(true);
progressDialog.setValue(0);
progressDialog.setMaximum(lookRows);
progressDialog.setMaximum(looktableRaster.height);
progressDialog.setMinimum(0);
progressDialog.show();
for (long i = 0; i < lookRows; i++) {
for (long j = 0; j < lookCols; j++) {
lastr = std::floor(sar_r(i, j));
nextr = std::ceil(sar_r(i, j));
lastc = std::floor(sar_c(i, j));
nextc = std::ceil(sar_c(i, j));
if (lastr < 0 || lastc < 0 || nextr >= slcRows || nextc >= slcCols) {
continue;
// 分块插值计算
long blockline = Memory1MB / looktableRaster.width / 8 * 300;//300M
blockline = blockline < 1 ? 1 : blockline;
#pragma omp parallel for
for (long startrowid = 0; startrowid < looktableRaster.height; startrowid = startrowid + blockline) {
long tempblockline = startrowid + blockline < looktableRaster.height ? blockline : looktableRaster.height - startrowid;
Eigen::MatrixXd sar_r = looktableRaster.getData(startrowid, 0, tempblockline, looktableRaster.width, 1);
Eigen::MatrixXd sar_c = looktableRaster.getData(startrowid, 0, tempblockline, looktableRaster.width, 2);
int64_t slc_min_rid = sar_r.array().minCoeff();
int64_t slc_max_rid = sar_r.array().maxCoeff();
int64_t slc_r_len = slc_max_rid - slc_min_rid+1;
slc_min_rid = slc_min_rid < 0 ? 0 : slc_min_rid-1;
slc_r_len = slc_r_len + 2;
Eigen::MatrixXd slcImg = slcRaster.getData(slc_min_rid, 0, slc_r_len, slcRaster.width, 1);
Eigen::MatrixXd outImg = outRaster.getData(startrowid, 0, tempblockline, outRaster.width, 1);
long lastr = 0;
long nextr = 0;
long lastc = 0;
long nextc = 0;
double Bileanervalue = 0;
// 插值
Landpoint p0{ 0,0,0 }, p11{ 0,0,0 }, p21{ 0,0,0 }, p12{ 0,0,0 }, p22{ 0,0,0 }, p{ 0,0,0 };
long lookRows = sar_r.rows();
long lookCols = sar_r.cols();
long slcRows = slcImg.rows();
long slcCols = slcImg.cols();
for (long i = 0; i < lookRows; i++) {
for (long j = 0; j < lookCols; j++) {
lastr = std::floor(sar_r(i, j));
nextr = std::ceil(sar_r(i, j));
lastc = std::floor(sar_c(i, j));
nextc = std::ceil(sar_c(i, j));
if (lastr < 0 || lastc < 0 || nextr >= slcRows || nextc >= slcCols) {
continue;
}
lastr = lastr - slc_min_rid;
nextr = nextr - slc_min_rid;
p0 = { sar_c(i,j) - lastc, sar_r(i,j) - lastr,0 };
p11 = Landpoint{ 0,0,slcImg(lastr,lastc) };
p21 = Landpoint{ 0,1,slcImg(nextr,lastc) };
p12 = Landpoint{ 1,0,slcImg(lastr,nextc) };
p22 = Landpoint{ 1,1,slcImg(nextr,nextc) };
Bileanervalue = Bilinear_interpolation(p0, p11, p21, p12, p22);
outImg(i, j) = Bileanervalue;
}
p0 = { sar_c(i,j)-lastc, sar_r(i,j)-lastr,0 };
p11 = Landpoint{ 0,0,slcImg(lastr,lastc) };
p21 = Landpoint{ 0,1,slcImg(nextr,lastc) };
p12 = Landpoint{ 1,0,slcImg(lastr,nextc) };
p22 = Landpoint{ 1,1,slcImg(nextr,nextc) };
Bileanervalue = Bilinear_interpolation(p0, p11, p21, p12, p22);
outImg(i, j) = Bileanervalue;
progressDialog.setValue(i);
}
progressDialog.setValue(i);
omp_set_lock(&lock);
outRaster.saveImage(outImg, startrowid, 0, 1);
int64_t processValue = progressDialog.value() + tempblockline;
processValue = processValue < progressDialog.maximum() ? processValue : progressDialog.maximum();
progressDialog.setValue(processValue);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
}
outRaster.saveImage(outImg, 0, 0, 1);
omp_destroy_lock(&lock); // 劫伙拷斤拷
progressDialog.close();
return ErrorCode::SUCCESS;
}
@ -757,6 +791,6 @@ ErrorCode GF3MainOrthProcess(QString inDEMPath, QString inTarFilepath, QString o
}
progressDialog.setValue(6);
progressDialog.close();
return ErrorCode::SUCCESS;
}