更新合并代码

main
陈增辉 2024-11-15 17:25:05 +08:00
parent 23277b911e
commit 26fd7c20ae
2 changed files with 154 additions and 89 deletions

View File

@ -21,8 +21,10 @@
#include <opencv2/opencv.hpp>
#include <QMessageBox>
#include <QDir>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <QProgressDialog>
/**
* ENVI
@ -569,28 +571,28 @@ Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count,
}
delete[] temp;
}
else if (gdal_datatype == GDT_UInt64) {
unsigned long* temp = new unsigned long[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
for (int j = 0; j <
cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_Int64) {
long* temp = new long[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
for (int j = 0; j <
cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
//else if (gdal_datatype == GDT_UInt64) {
// unsigned long* temp = new unsigned long[rows_count * cols_count];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
// rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
// for (int j = 0; j <
// cols_count; j++) {
// datamatrix(i, j) = temp[i * cols_count + j];
// }
// }
// delete[] temp;
//}
//else if (gdal_datatype == GDT_Int64) {
// long* temp = new long[rows_count * cols_count];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
// rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
// for (int j = 0; j <
// cols_count; j++) {
// datamatrix(i, j) = temp[i * cols_count + j];
// }
// }
// delete[] temp;
//}
else if(gdal_datatype == GDT_Float32) {
float* temp = new float[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
@ -695,28 +697,27 @@ Eigen::MatrixXi gdalImage::getDatai(int start_row, int start_col, int rows_count
}
delete[] temp;
}
else if (gdal_datatype == GDT_UInt64) {
unsigned long* temp = new unsigned long[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
for (int j = 0; j <
cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
else if (gdal_datatype == GDT_Int64) {
long* temp = new long[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
for (int j = 0; j <
cols_count; j++) {
datamatrix(i, j) = temp[i * cols_count + j];
}
}
delete[] temp;
}
//else if (gdal_datatype == GDT_UInt64) {
// unsigned long* temp = new unsigned long[rows_count * cols_count];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
// rows_count, gdal_datatype, 0, 0); for (int i = 0; i < rows_count; i++) {
// for (int j = 0; j <
// cols_count; j++) {
// datamatrix(i, j) = temp[i * cols_count + j];
// }
// }
// delete[] temp;
//}
//else if (gdal_datatype == GDT_Int64) {
// long* temp = new long[rows_count * cols_count];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
// for (int i = 0; i < rows_count; i++) {
// for (int j = 0; j < cols_count; j++) {
// datamatrix(i, j) = temp[i * cols_count + j];
// }
// }
// delete[] temp;
//}
else if (gdal_datatype == GDT_Float32) {
float* temp = new float[rows_count * cols_count];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count,
@ -773,6 +774,7 @@ GDALDataType gdalImage::getDataType()
void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0,
int band_ids = 1)
{
GDALDataType datetype = this->getDataType();
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
@ -789,7 +791,15 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
} else {
poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height,
this->band_num, GDT_Float32, NULL); // 斤拷锟斤拷
this->band_num, datetype, NULL); // 斤拷锟斤拷
if (nullptr == poDstDS) {
QString tip = u8"file path: " + this->img_path + " image size :( " + QString::number(this->height) + " , " + QString::number(this->width) + " ) " + " input size (" + QString::number(start_row + data.rows()) + ", " + QString::number(start_col + data.cols()) + ") ";
qDebug() << tip;
throw std::exception(tip.toUtf8().constData());
return;
}
poDstDS->SetProjection(this->projection.toUtf8().constData());
double gt_ptr[6];
@ -802,6 +812,9 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
//delete gt_ptr;
}
int datarows = data.rows();
int datacols = data.cols();
@ -817,7 +830,7 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
// poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
// datarows, GDT_Float32,band_ids, num,0,0,0);
poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
databuffer, datacols, datarows, GDT_Float32, 0, 0);
databuffer, datacols, datarows, datetype, 0, 0);
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
@ -829,6 +842,8 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
void gdalImage::saveImage(Eigen::MatrixXi data, int start_row, int start_col, int band_ids)
{
GDALDataType datetype=this->getDataType();
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
@ -872,7 +887,7 @@ void gdalImage::saveImage(Eigen::MatrixXi data, int start_row, int start_col, in
// poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols,
// datarows, GDT_Float32,band_ids, num,0,0,0);
poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
databuffer, datacols, datarows, GDT_Int64, 0, 0);
databuffer, datacols, datarows, datetype, 0, 0);
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
@ -1544,8 +1559,10 @@ int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path)
return -1;
}
ErrorCode MergeRasterProcess(QVector<QString> filepaths, QString outfileptah, QString mainString, MERGEMODE mergecode, bool isENVI)
ErrorCode MergeRasterProcess(QVector<QString> filepaths, QString outfileptah, QString mainString, MERGEMODE mergecode, bool isENVI, ShowProessAbstract* dia )
{
// 参数检查
if (!isExists(mainString)) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND) )<< "\t" << mainString;
@ -1588,8 +1605,6 @@ ErrorCode MergeRasterProcess(QVector<QString> filepaths, QString outfileptah, QS
}
}
Eigen::MatrixXd maingt = mainimg.getGeoTranslation();
Eigen::MatrixXd rgt = Eigen::MatrixXd::Zero(2,3);
RasterExtend mainExtend = mainimg.getExtend();
@ -1620,38 +1635,52 @@ ErrorCode MergeRasterProcess(QVector<QString> filepaths, QString outfileptah, QS
gdalImage resultImage = CreategdalImage(outfileptah, height, width, mainimg.band_num, rgt, EPSGCode, mainType, true, true, isENVI);
QString resultMaskString = addMaskToFileName(outfileptah, QString("_MASK"));
gdalImage maskImage = CreategdalImage(resultMaskString, height, width,1 , rgt, EPSGCode, GDT_UInt16, true, true, isENVI);
gdalImage maskImage = CreategdalImage(resultMaskString, height, width,1 , rgt, EPSGCode, GDT_Int32, true, true, isENVI);
// 初始化
// 分块计算
long resultline = Memory1MB * 100 / 8 / resultImage.width;
resultline = resultline < 100 ? resultline : 100; // 最多100行
long resultline = Memory1MB * 500 / 8 / resultImage.width;
resultline = resultline < 3000 ? resultline : 3000; // 最多100行
resultline = resultline > 0 ? resultline : 2;
long bandnum = resultImage.band_num + 1;
long starti = 0;
long rasterCount = imgdslist.count();
QProgressDialog progressDialog(u8"初始化影像", u8"终止", 0, resultImage.height);
progressDialog.setWindowTitle(u8"初始化影像");
progressDialog.setWindowModality(Qt::WindowModal);
progressDialog.setAutoClose(true);
progressDialog.setValue(0);
progressDialog.setMaximum(resultImage.height);
progressDialog.setMinimum(0);
progressDialog.show();
for (starti = 0; starti < resultImage.height; starti = starti + resultline) {
long blocklines = resultline;
blocklines = starti + blocklines < resultImage.height ? blocklines : resultImage.height - starti;
for (long b = 1; b < bandnum; b++) {
Eigen::MatrixXd data = resultImage.getData(starti, 0, blocklines, resultImage.width, b);
Eigen::MatrixXd maskdata = maskImage.getData(starti, 0, blocklines, resultImage.width, b);
Eigen::MatrixXi maskdata = maskImage.getDatai(starti, 0, blocklines, resultImage.width, b);
data = data.array() * 0;
maskdata = maskdata.array() * 0;
resultImage.saveImage(data, starti, 0, b);
maskImage.saveImage(maskdata, starti, 0, b);
}
if (nullptr != dia) {
dia->showProcess(starti * 1.0 / resultImage.height, u8"初始化影像数据");
}
progressDialog.setValue(starti+ blocklines);
}
progressDialog.close();
switch (mergecode)
{
case MERGE_GEOCODING:
return MergeRasterInGeoCoding(imgdslist, resultImage, maskImage);
return MergeRasterInGeoCoding(imgdslist, resultImage, maskImage, dia);
default:
break;
}
@ -1660,7 +1689,7 @@ ErrorCode MergeRasterProcess(QVector<QString> filepaths, QString outfileptah, QS
return ErrorCode::SUCCESS;
}
ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resultimg, gdalImage maskimg)
ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resultimg, gdalImage maskimg, ShowProessAbstract* dia)
{
omp_set_num_threads(Paral_num_thread);
// 逐点合并计算
@ -1680,12 +1709,22 @@ ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resulti
long processNumber = 0;
QProgressDialog progressDialog(u8"合并影像", u8"终止", 0, resultimg.height);
progressDialog.setWindowTitle(u8"合并影像");
progressDialog.setWindowModality(Qt::WindowModal);
progressDialog.setAutoClose(true);
progressDialog.setValue(0);
progressDialog.setMaximum(resultimg.height);
progressDialog.setMinimum(0);
progressDialog.show();
omp_lock_t lock;
omp_init_lock(&lock);
#pragma omp parallel for
//#pragma omp parallel for
for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
long blocklines = resultline;
blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
@ -1700,9 +1739,9 @@ ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resulti
long minRid = imgdslist[ir].height;
long maxRid = 0;
Eigen::MatrixXd ridlist = Eigen::MatrixXd::Zero(blocklines, resultimg.width);
Eigen::MatrixXd cidlist = Eigen::MatrixXd::Zero(blocklines, resultimg.width);
Eigen::MatrixXd ridlist = resultimg.getData(starti, 0, blocklines, resultimg.width, 1);
ridlist = ridlist.array() * 0;
Eigen::MatrixXd cidlist = ridlist.array() * 0;
for (long i = 0; i < blocklines; i++) {// 行号
rid = starti + i;
@ -1764,6 +1803,14 @@ ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resulti
r0 = ridtemp - std::floor(ridtemp);
c0 = cidtemp - std::floor(cidtemp);
lastr = std::floor(ridtemp);
nextr = std::ceil(ridtemp);
lastc = std::floor(cidtemp);
lastc = std::ceil(cidtemp);
lastr = lastr - minRid;
nextr = nextr - minRid;
p0 = Landpoint{ c0,r0,0 };
p11 = Landpoint{ 0,0,data(lastr,lastc) };
p21 = Landpoint{ 0,1,data(nextr,lastc) };
@ -1778,16 +1825,22 @@ ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resulti
maskimg.saveImage(resultmask, starti, 0, b);
}
}
omp_set_lock(&lock);
processNumber = processNumber + blocklines;
std::cout << "\rprocess bar:\t" << processNumber * 100.0 / resultimg.height << " % " << "\t\t\t";
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
if (nullptr != dia) {
dia->showProcess(processNumber * 1.0 / resultimg.height, u8"合并图像");
}
if (progressDialog.maximum() >= processNumber) {
processNumber = progressDialog.maximum() - 1;
}
progressDialog.setValue(processNumber);
omp_unset_lock(&lock);
}
omp_destroy_lock(&lock); // 劫伙拷斤拷
omp_destroy_lock(&lock);
std::cout << std::endl;
progressDialog.setWindowTitle(u8"影像掩膜");
// 全体处理
for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
long blocklines = resultline;
blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
@ -1798,6 +1851,7 @@ ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resulti
for (long i = 0; i < data.rows(); i++) {
for (long j = 0; j < data.cols(); j++) {
if (maskdata(i, j) == 0) {
data(i, j) = -9999;
continue;
}
data(i, j) = data(i, j) / maskdata(i, j);
@ -1807,9 +1861,10 @@ ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resulti
resultimg.saveImage(data, starti, 0, b);
maskimg.saveImage(maskdata, starti, 0, b);
}
progressDialog.setValue(starti + blocklines);
}
resultimg.setNoDataValue(-9999);
progressDialog.close();
return ErrorCode::SUCCESS;
}
@ -1862,6 +1917,9 @@ void gdalImageComplex::setData(Eigen::MatrixXcd data)
void gdalImageComplex::saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
@ -2183,6 +2241,8 @@ long GetEPSGFromRasterFile(QString filepath)
return -1;
}
std::cout << pszProjection << std::endl;
long epsgCode = atoi(oSRS.GetAuthorityCode(nullptr)); // 获取EPSG代码
if (epsgCode != 0) {
@ -2193,6 +2253,7 @@ long GetEPSGFromRasterFile(QString filepath)
else {
qDebug() << "EPSG code could not be determined from the spatial reference.";
GDALClose(poDataset);
return -1;
}
}
@ -2303,6 +2364,11 @@ CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long epsg_code)
}
}
void ShowProessAbstract::showProcess(double precent, QString tip)
{
}
void ShowProessAbstract::showToolInfo(QString tip)
{
}

View File

@ -74,6 +74,15 @@ enum GDALREADARRCOPYMETHOD {
class ShowProessAbstract {
public:
virtual void showProcess(double precent,QString tip);
virtual void showToolInfo( QString tip) ;
};
/// 根据经纬度获取
/// EPSG代码根据经纬度返回对应投影坐标系统其中如果在中华人民共和国境内默认使用
@ -99,13 +108,6 @@ std::shared_ptr<PointRaster> GetCenterPointInRaster(QString filepath);
CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long EPSGCODE);
// 文件打开 // 当指令销毁时调用GDALClose 销毁类型
std::shared_ptr<GDALDataset> OpenDataset(const QString& in_path, GDALAccess rwmode= GA_ReadOnly);
void CloseDataset(GDALDataset* ptr);
@ -248,13 +250,10 @@ enum MERGEMODE
};
ErrorCode MergeRasterProcess(QVector<QString> filepath, QString outfileptah, QString mainString, MERGEMODE mergecode = MERGEMODE::MERGE_GEOCODING, bool isENVI = false);
ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> inimgs, gdalImage resultimg,gdalImage maskimg);
ErrorCode MergeRasterProcess(QVector<QString> filepath, QString outfileptah, QString mainString, MERGEMODE mergecode = MERGEMODE::MERGE_GEOCODING, bool isENVI = false, ShowProessAbstract* dia=nullptr);
ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> inimgs, gdalImage resultimg,gdalImage maskimg, ShowProessAbstract* dia = nullptr);
template<typename T>
std::shared_ptr<T> readDataArr(gdalImage &img,int start_row, int start_col, int rows_count, int cols_count, int band_ids, GDALREADARRCOPYMETHOD method);