From 70c4a621293e4aced20debeee7d338d1ad7539c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E9=99=88=E5=A2=9E=E8=BE=89?= <3045316072@qq.com> Date: Thu, 22 May 2025 21:04:45 +0800 Subject: [PATCH] =?UTF-8?q?=E4=BF=AE=E6=94=B9DEM->XYZ=E7=9A=84=E5=88=86?= =?UTF-8?q?=E5=9D=97=E8=AE=A1=E7=AE=97=E9=83=A8=E5=88=86?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- BaseCommonLibrary/ImageOperatorFuntion.cpp | 85 +++++++------ .../GF3CalibrationAndGeocodingClass.cpp | 120 +++++++++++------- 2 files changed, 121 insertions(+), 84 deletions(-) diff --git a/BaseCommonLibrary/ImageOperatorFuntion.cpp b/BaseCommonLibrary/ImageOperatorFuntion.cpp index ee94390..ab4ae0a 100644 --- a/BaseCommonLibrary/ImageOperatorFuntion.cpp +++ b/BaseCommonLibrary/ImageOperatorFuntion.cpp @@ -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; } diff --git a/GF3CalibrationAndOrthLib/GF3CalibrationAndGeocodingClass.cpp b/GF3CalibrationAndOrthLib/GF3CalibrationAndGeocodingClass.cpp index 6eb4cea..e720544 100644 --- a/GF3CalibrationAndOrthLib/GF3CalibrationAndGeocodingClass.cpp +++ b/GF3CalibrationAndOrthLib/GF3CalibrationAndGeocodingClass.cpp @@ -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 <<"!="<= 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; }