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

拆分ImageOperator文件

pull/13/head
陈增辉 2025-03-25 16:41:05 +08:00
parent 66bed73702
commit 17aca87ef0
8 changed files with 2678 additions and 2446 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
BaseCommonLibrary/BaseCommonLibrary.vcxproj
View File

@ -287,14 +287,19 @@
<ClCompile Include="BaseTool\BaseTool.cpp" />
<ClCompile Include="BaseTool\EchoDataFormat.cpp" />
<ClCompile Include="BaseTool\FileOperator.cpp" />
<ClCompile Include="BaseTool\gdalImageComplexOperator.cpp" />
<ClCompile Include="BaseTool\gdalImageOperator.cpp" />
<ClCompile Include="BaseTool\GeoOperator.cpp" />
<ClCompile Include="BaseTool\ImageOperatorBase.cpp" />
<ClCompile Include="BaseTool\LogInfoCls.cpp" />
<ClCompile Include="BaseTool\MergeRasterOperator.cpp" />
<ClCompile Include="BaseTool\PrintMsgToQDebug.cpp" />
<ClCompile Include="BaseTool\QToolProcessBarDialog.cpp" />
<ClCompile Include="BaseTool\RasterToolBase.cpp" />
<ClCompile Include="BaseTool\SARSimulationImageL1.cpp" />
<ClCompile Include="BaseTool\ShowProessAbstract.cpp" />
<ClCompile Include="BaseTool\stdafx.cpp" />
<ClCompile Include="BaseTool\TestImageOperator.cpp" />
<ClCompile Include="dllmain.cpp" />
<ClCompile Include="pch.cpp">
<PrecompiledHeader Condition="'$(Configuration)|$(Platform)'=='Release|x64'">Create</PrecompiledHeader>

15
BaseCommonLibrary/BaseCommonLibrary.vcxproj.filters
View File

@ -104,6 +104,21 @@
<ClCompile Include="BaseTool\PrintMsgToQDebug.cpp">
<Filter>BaseTool</Filter>
</ClCompile>
<ClCompile Include="BaseTool\gdalImageOperator.cpp">
<Filter>BaseTool</Filter>
</ClCompile>
<ClCompile Include="BaseTool\gdalImageComplexOperator.cpp">
<Filter>BaseTool</Filter>
</ClCompile>
<ClCompile Include="BaseTool\ShowProessAbstract.cpp">
<Filter>BaseTool</Filter>
</ClCompile>
<ClCompile Include="BaseTool\TestImageOperator.cpp">
<Filter>源文件</Filter>
</ClCompile>
<ClCompile Include="BaseTool\MergeRasterOperator.cpp">
<Filter>BaseTool</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<QtMoc Include="BaseTool\QToolProcessBarDialog.h">

2453
BaseCommonLibrary/BaseTool/ImageOperatorBase.cpp
View File

File diff suppressed because it is too large Load Diff

488
BaseCommonLibrary/BaseTool/MergeRasterOperator.cpp Normal file
View File

@ -0,0 +1,488 @@
#include "stdafx.h"
#include "ImageOperatorBase.h"
#include "BaseTool.h"
#include "GeoOperator.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <omp.h>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <gdal.h>
#include <gdal_utils.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <proj.h>
#include <string.h>
#include <memory>
#include <iostream>
#include "FileOperator.h"
#include <opencv2/opencv.hpp>
#include <QMessageBox>
#include <QDir>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <QProgressDialog>
#include <gdal_priv.h>
#include <ogr_spatialref.h> // OGRSpatialReference 用于空间参考转换
#include <gdal_alg.h> // 用于 GDALWarp 操作
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;
return ErrorCode::FILENOFOUND;
}
else {}
gdalImage mainimg(mainString);
QVector<gdalImage> imgdslist(filepaths.count());
for (long i = 0; i < filepaths.count(); i++) {
if (!isExists(filepaths[i])) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND)) << "\t" << filepaths[i];
return ErrorCode::FILENOFOUND;
}
else {
imgdslist[i] = gdalImage(filepaths[i]);
if (imgdslist[i].band_num != mainimg.band_num) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTERBAND_NOTEQUAL)) << "\t" << imgdslist[i].band_num << " != " << mainimg.band_num;
return ErrorCode::RASTERBAND_NOTEQUAL;
}
}
}
// 检查坐标系是否统一
long EPSGCode = GetEPSGFromRasterFile(mainString);
long tempCode = 0;
for (long i = 0; i < filepaths.count(); i++) {
tempCode = GetEPSGFromRasterFile(filepaths[i]);
if (EPSGCode != tempCode) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSAME)) << "\t" << EPSGCode << "!=" << tempCode;
return ErrorCode::EPSGCODE_NOTSAME;
}
}
// 检查影像类型是否统一
GDALDataType mainType = mainimg.getDataType();
for (long i = 0; i < imgdslist.count(); i++) {
if (mainType != imgdslist[i].getDataType()) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTER_DATETYPE_NOTSAME)) << "\t" << mainType << "!=" << imgdslist[i].getDataType();
return ErrorCode::RASTER_DATETYPE_NOTSAME;
}
}
Eigen::MatrixXd maingt = mainimg.getGeoTranslation();
Eigen::MatrixXd rgt = Eigen::MatrixXd::Zero(2, 3);
RasterExtend mainExtend = mainimg.getExtend();
rgt(0, 1) = (mainExtend.max_x - mainExtend.min_x) / (mainimg.width - 1); //dx
rgt(1, 2) = -1 * std::abs(((mainExtend.max_y - mainExtend.min_y) / (mainimg.height - 1)));//dy
QVector<RasterExtend> extendlist(imgdslist.count());
for (long i = 0; i < imgdslist.count(); i++) {
extendlist[i] = imgdslist[i].getExtend();
mainExtend.min_x = mainExtend.min_x < extendlist[i].min_x ? mainExtend.min_x : extendlist[i].min_x;
mainExtend.max_x = mainExtend.max_x > extendlist[i].max_x ? mainExtend.max_x : extendlist[i].max_x;
mainExtend.min_y = mainExtend.min_y < extendlist[i].min_y ? mainExtend.min_y : extendlist[i].min_y;
mainExtend.max_y = mainExtend.max_y > extendlist[i].max_y ? mainExtend.max_y : extendlist[i].max_y;
}
rgt(0, 0) = mainExtend.min_x;
rgt(1, 0) = mainExtend.max_y;
// 计算数量
long width = std::ceil((mainExtend.max_x - mainExtend.min_x) / rgt(0, 1) + 1);
long height = std::ceil(std::abs((mainExtend.min_y - mainExtend.max_y) / rgt(1, 2)) + 1);
OGRSpatialReference oSRS;
if (oSRS.importFromEPSG(EPSGCode) != OGRERR_NONE) {
qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSUPPORT)) << "\t" << EPSGCode;
return ErrorCode::EPSGCODE_NOTSUPPORT;
}
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_Int32, true, true, isENVI);
// 初始化
long resultline = Memory1MB * 500 / 8 / resultImage.width;
resultline = resultline < 10000 ? resultline : 10000; // 最多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::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, dia);
default:
break;
}
return ErrorCode::SUCCESS;
}
ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resultimg, gdalImage maskimg, ShowProessAbstract* dia)
{
omp_set_num_threads(Paral_num_thread);
// 逐点合并计算
QVector<RasterExtend> extendlist(imgdslist.count());
for (long i = 0; i < imgdslist.count(); i++) {
extendlist[i] = imgdslist[i].getExtend();
imgdslist[i].InitInv_gt();
}
// 分块计算
long resultline = Memory1MB * 1000 / 8 / resultimg.width;
resultline = resultline < 300 ? resultline : 300; // 最多100行
long bandnum = resultimg.band_num + 1;
long starti = 0;
long rasterCount = imgdslist.count();
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
for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
long blocklines = resultline;
blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
long rid = starti;
long cid = 0;
Landpoint pp = { 0,0,0 };
Landpoint lp = { 0,0,0 };
for (long ir = 0; ir < rasterCount; ir++) {// 影像
long minRid = imgdslist[ir].height;
long maxRid = 0;
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;
for (long j = 0; j < resultimg.width; j++) {// 列号
cid = j;
resultimg.getLandPoint(rid, cid, 0, pp);
lp = imgdslist[ir].getRow_Col(pp.lon, pp.lat); // 获取点坐标
ridlist(i, j) = lp.lat;
cidlist(i, j) = lp.lon;
}
}
//ImageShowDialogClass* dialog = new ImageShowDialogClass;
//dialog->show();
//dialog->load_double_MatrixX_data(cidlist, u8"");
//dialog->exec();
if (ridlist.maxCoeff() < 0 || ridlist.minCoeff() >= imgdslist[ir].height) {
continue;
}
if (cidlist.maxCoeff() < 0 || cidlist.minCoeff() >= imgdslist[ir].width) {
continue;
}
minRid = std::floor(ridlist.minCoeff());
maxRid = std::ceil(ridlist.maxCoeff());
minRid = minRid < 0 ? 0 : minRid;
maxRid = maxRid < imgdslist[ir].height ? maxRid : imgdslist[ir].height - 1;
long rowlen = maxRid - minRid + 1;
if (rowlen <= 0) {
continue;
}
// 获取分配代码
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 rowcount = 0;
long colcount = 0;
double ridtemp = 0, cidtemp = 0;
long lastr = 0, nextr = 0;
long lastc = 0, nextc = 0;
double r0 = 0, c0 = 0;
for (long b = 1; b < bandnum; b++) {
Eigen::MatrixXd resultdata = resultimg.getData(starti, 0, blocklines, resultimg.width, b);
Eigen::MatrixXi resultmask = maskimg.getDatai(starti, 0, blocklines, resultimg.width, b);
Eigen::MatrixXd data = imgdslist[ir].getData(minRid, 0, rowlen, imgdslist[ir].width, b);
double nodata = imgdslist[ir].getNoDataValue(b);
for (long ii = 0; ii < data.rows(); ii++) {
for (long jj = 0; jj < data.cols(); jj++) {
if (std::abs(data(ii, jj) - nodata) < 1e-6) {
data(ii, jj) = 0;
}
}
}
rowcount = ridlist.rows();
colcount = ridlist.cols();
double Bileanervalue = 0;
for (long i = 0; i < rowcount; i++) {
for (long j = 0; j < colcount; j++) {
ridtemp = ridlist(i, j);
cidtemp = cidlist(i, j);
lastr = std::floor(ridtemp);
nextr = std::ceil(ridtemp);
lastc = std::floor(cidtemp);
nextc = std::ceil(cidtemp);
if (lastr < 0 || lastr >= imgdslist[ir].height
|| nextr < 0 || nextr >= imgdslist[ir].height
|| lastc < 0 || lastc >= imgdslist[ir].width
|| nextc < 0 || nextc >= imgdslist[ir].width) {
continue;
}
else {}
r0 = ridtemp - std::floor(ridtemp);
c0 = cidtemp - std::floor(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) };
p12 = Landpoint{ 1,0,data(lastr,nextc) };
p22 = Landpoint{ 1,1,data(nextr,nextc) };
Bileanervalue = Bilinear_interpolation(p0, p11, p21, p12, p22);
if (std::abs(Bileanervalue) < 1e-6 || resultmask(i, j) > 0) {
continue;
}
resultdata(i, j) = resultdata(i, j) + Bileanervalue;
resultmask(i, j) = resultmask(i, j) + 1;
}
}
resultimg.saveImage(resultdata, starti, 0, b);
maskimg.saveImage(resultmask, starti, 0, b);
}
}
omp_set_lock(&lock);
processNumber = processNumber + blocklines;
qDebug() << "\rprocess bar:\t" << processNumber * 100.0 / resultimg.height << " % " << "\t\t\t";
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);
progressDialog.setWindowTitle(u8"影像掩膜");
progressDialog.setLabelText(u8"影像掩膜");
for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
long blocklines = resultline;
blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
for (long b = 1; b < bandnum; b++) {
Eigen::MatrixXd data = resultimg.getData(starti, 0, blocklines, resultimg.width, b);
Eigen::MatrixXd maskdata = maskimg.getData(starti, 0, blocklines, maskimg.width, b);
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);
}
}
resultimg.saveImage(data, starti, 0, b);
maskimg.saveImage(maskdata, starti, 0, b);
}
if (nullptr != dia) {
dia->showProcess((starti + blocklines) * 1.0 / resultimg.height, u8"影像掩膜");
}
progressDialog.setValue(starti + blocklines);
}
progressDialog.close();
return ErrorCode::SUCCESS;
}
void MergeTiffs(QList<QString> inputFiles, QString outputFile) {
GDALAllRegister();
if (inputFiles.isEmpty()) {
fprintf(stderr, "No input files provided.\n");
return;
}
// Open the first file to determine the data type and coordinate system
GDALDataset* poFirstDS = (GDALDataset*)GDALOpen(inputFiles.first().toUtf8().constData(), GA_ReadOnly);
if (poFirstDS == nullptr) {
fprintf(stderr, "Failed to open the first file %s\n", inputFiles.first().toUtf8().constData());
return;
}
double adfGeoTransform[6];
CPLErr eErr = poFirstDS->GetGeoTransform(adfGeoTransform);
if (eErr != CE_None) {
fprintf(stderr, "Failed to get GeoTransform for the first file %s\n", inputFiles.first().toUtf8().constData());
GDALClose(poFirstDS);
return;
}
int nXSize = 0;
int nYSize = 0;
double minX = std::numeric_limits<double>::max();
double minY = std::numeric_limits<double>::max();
double maxX = std::numeric_limits<double>::lowest();
double maxY = std::numeric_limits<double>::lowest();
double pixelWidth = adfGeoTransform[1];
double pixelHeight = adfGeoTransform[5];
// Determine the bounding box and size of the output raster
for (const QString& inputFile : inputFiles) {
GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile.toUtf8().constData(), GA_ReadOnly);
if (poSrcDS == nullptr) {
fprintf(stderr, "Failed to open %s\n", inputFile.toUtf8().constData());
continue;
}
double adfThisTransform[6];
eErr = poSrcDS->GetGeoTransform(adfThisTransform);
if (eErr != CE_None) {
fprintf(stderr, "Failed to get GeoTransform for %s\n", inputFile.toUtf8().constData());
GDALClose(poSrcDS);
continue;
}
minX = std::min(minX, adfThisTransform[0]);
minY = std::min(minY, adfThisTransform[3] + adfThisTransform[5] * poSrcDS->GetRasterYSize());
maxX = std::max(maxX, adfThisTransform[0] + adfThisTransform[1] * poSrcDS->GetRasterXSize());
maxY = std::max(maxY, adfThisTransform[3]);
GDALClose(poSrcDS);
}
nXSize = static_cast<int>(std::ceil((maxX - minX) / pixelWidth));
nYSize = static_cast<int>(std::ceil((maxY - minY) / (-pixelHeight)));
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (poDriver == nullptr) {
fprintf(stderr, "GTiff driver not available.\n");
GDALClose(poFirstDS);
return;
}
char** papszOptions = nullptr;
GDALDataset* poDstDS = poDriver->Create(outputFile.toUtf8().constData(), nXSize, nYSize, 1, poFirstDS->GetRasterBand(1)->GetRasterDataType(), papszOptions);
if (poDstDS == nullptr) {
fprintf(stderr, "Creation of output file failed.\n");
GDALClose(poFirstDS);
return;
}
poDstDS->SetGeoTransform(adfGeoTransform);
const OGRSpatialReference* oSRS = poFirstDS->GetSpatialRef();
poDstDS->SetSpatialRef(oSRS);
float fillValue = std::numeric_limits<float>::quiet_NaN();
void* pafScanline = CPLMalloc(GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * nXSize);
memset(pafScanline, 0, GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * nXSize);
// Initialize all pixels to NaN
for (int iY = 0; iY < nYSize; ++iY) {
GDALRasterBand* poBand = poDstDS->GetRasterBand(1);
poBand->RasterIO(GF_Write, 0, iY, nXSize, 1, pafScanline, nXSize, 1, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
}
CPLFree(pafScanline);
// Read each source image and write into the destination image
for (const QString& inputFile : inputFiles) {
GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile.toUtf8().constData(), GA_ReadOnly);
if (poSrcDS == nullptr) {
fprintf(stderr, "Failed to open %s\n", inputFile.toUtf8().constData());
continue;
}
double adfThisTransform[6];
poSrcDS->GetGeoTransform(adfThisTransform);
int srcXSize = poSrcDS->GetRasterXSize();
int srcYSize = poSrcDS->GetRasterYSize();
int dstXOffset = static_cast<int>(std::round((adfThisTransform[0] - minX) / pixelWidth));
int dstYOffset = static_cast<int>(std::round((maxY - adfThisTransform[3]) / (-pixelHeight)));
GDALRasterBand* poSrcBand = poSrcDS->GetRasterBand(1);
GDALRasterBand* poDstBand = poDstDS->GetRasterBand(1);
void* pafBuffer = CPLMalloc(GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * srcXSize * srcYSize);
poSrcBand->RasterIO(GF_Read, 0, 0, srcXSize, srcYSize, pafBuffer, srcXSize, srcYSize, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
poDstBand->RasterIO(GF_Write, dstXOffset, dstYOffset, srcXSize, srcYSize, pafBuffer, srcXSize, srcYSize, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
CPLFree(pafBuffer);
GDALClose(poSrcDS);
}
GDALClose(poDstDS);
GDALClose(poFirstDS);
}

43
BaseCommonLibrary/BaseTool/ShowProessAbstract.cpp Normal file
View File

@ -0,0 +1,43 @@
#include "stdafx.h"
#include "ImageOperatorBase.h"
#include "BaseTool.h"
#include "GeoOperator.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <omp.h>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <gdal.h>
#include <gdal_utils.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <proj.h>
#include <string.h>
#include <memory>
#include <iostream>
#include "FileOperator.h"
#include <opencv2/opencv.hpp>
#include <QMessageBox>
#include <QDir>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <QProgressDialog>
#include <gdal_priv.h>
#include <ogr_spatialref.h> // OGRSpatialReference 用于空间参考转换
#include <gdal_alg.h> // 用于 GDALWarp 操作
void ShowProessAbstract::showProcess(double precent, QString tip)
{
}
void ShowProessAbstract::showToolInfo(QString tip)
{
}

146
BaseCommonLibrary/BaseTool/TestImageOperator.cpp Normal file
View File

@ -0,0 +1,146 @@
#include "stdafx.h"
#include "ImageOperatorBase.h"
#include "BaseTool.h"
#include "GeoOperator.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <omp.h>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <gdal.h>
#include <gdal_utils.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <proj.h>
#include <string.h>
#include <memory>
#include <iostream>
#include "FileOperator.h"
#include <opencv2/opencv.hpp>
#include <QMessageBox>
#include <QDir>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <QProgressDialog>
#include <gdal_priv.h>
#include <ogr_spatialref.h> // OGRSpatialReference 用于空间参考转换
#include <gdal_alg.h> // 用于 GDALWarp 操作
void testOutAmpArr(QString filename, float* amp, long rowcount, long colcount)
{
Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
for (long hii = 0; hii < rowcount; hii++) {
for (long hjj = 0; hjj < colcount; hjj++) {
h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
}
}
QString ampPath = getDebugDataPath(filename);
saveEigenMatrixXd2Bin(h_amp_img, ampPath);
qDebug() << filename.toLocal8Bit().constData();
qDebug() << "max:\t" << h_amp_img.maxCoeff();
qDebug() << "min:\t" << h_amp_img.minCoeff();
}
void testOutAmpArr(QString filename, double* amp, long rowcount, long colcount)
{
Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
for (long hii = 0; hii < rowcount; hii++) {
for (long hjj = 0; hjj < colcount; hjj++) {
h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
}
}
QString ampPath = getDebugDataPath(filename);
saveEigenMatrixXd2Bin(h_amp_img, ampPath);
qDebug() << filename.toLocal8Bit().constData();
qDebug() << "max:\t" << h_amp_img.maxCoeff();
qDebug() << "min:\t" << h_amp_img.minCoeff();
}
void testOutClsArr(QString filename, long* amp, long rowcount, long colcount) {
Eigen::MatrixXd h_amp_img = Eigen::MatrixXd::Zero(rowcount, colcount);
for (long hii = 0; hii < rowcount; hii++) {
for (long hjj = 0; hjj < colcount; hjj++) {
h_amp_img(hii, hjj) = amp[hii * colcount + hjj];
}
}
QString ampPath = getDebugDataPath(filename);
saveEigenMatrixXd2Bin(h_amp_img, ampPath);
qDebug() << filename.toLocal8Bit().constData();
qDebug() << "max:\t" << h_amp_img.maxCoeff();
qDebug() << "min:\t" << h_amp_img.minCoeff();
}
void testOutComplexDoubleArr(QString filename, std::complex<double>* data, long rowcount, long colcount)
{
QString ampPath = getDebugDataPath(filename);
gdalImageComplex compleximg = CreateEchoComplex(ampPath, rowcount, colcount, 1);
compleximg.saveImage(data, 0, 0, rowcount, colcount, 1);
return;
}
void testOutDataArr(QString filename, double* data, long rowcount, long colcount)
{
return testOutAmpArr(filename, data, rowcount, colcount);
}
void testOutDataArr(QString filename, float* data, long rowcount, long colcount)
{
return testOutAmpArr(filename, data, rowcount, colcount);
}
void testOutDataArr(QString filename, long* data, long rowcount, long colcount)
{
return testOutClsArr(filename, data, rowcount, colcount);
}
void testOutAntPatternTrans(QString antpatternfilename, double* antPatternArr,
double starttheta, double deltetheta,
double startphi, double deltaphi,
long thetanum, long phinum)
{
Eigen::MatrixXd antPatternMatrix(thetanum, phinum);
for (long t = 0; t < thetanum; ++t) {
for (long p = 0; p < phinum; ++p) {
long index = t * phinum + p;
if (index < thetanum * phinum) {
antPatternMatrix(t, p) = static_cast<double>(antPatternArr[index]); // Copy to Eigen matrix
}
}
}
Eigen::MatrixXd gt(2, 3);
gt(0, 0) = startphi;//x
gt(0, 1) = deltaphi;
gt(0, 2) = 0;
gt(1, 0) = starttheta;
gt(1, 1) = 0;
gt(1, 2) = deltetheta;
QString antpatternfilepath = getDebugDataPath(antpatternfilename);
gdalImage ds = CreategdalImageDouble(antpatternfilepath, thetanum, phinum, 1, gt, "", true, true, true);
ds.saveImage(antPatternMatrix, 0, 0, 1);
}

497
BaseCommonLibrary/BaseTool/gdalImageComplexOperator.cpp Normal file
View File

@ -0,0 +1,497 @@
#include "stdafx.h"
#include "ImageOperatorBase.h"
#include "BaseTool.h"
#include "GeoOperator.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <omp.h>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <gdal.h>
#include <gdal_utils.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <proj.h>
#include <string.h>
#include <memory>
#include <iostream>
#include "FileOperator.h"
#include <opencv2/opencv.hpp>
#include <QMessageBox>
#include <QDir>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <QProgressDialog>
#include <gdal_priv.h>
#include <ogr_spatialref.h> // OGRSpatialReference 用于空间参考转换
#include <gdal_alg.h> // 用于 GDALWarp 操作
gdalImageComplex::gdalImageComplex(const QString& raster_path)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
this->img_path = raster_path;
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
raster_path.toUtf8().constData(), GA_ReadOnly));
this->width = rasterDataset->GetRasterXSize();
this->height = rasterDataset->GetRasterYSize();
this->band_num = rasterDataset->GetRasterCount();
double* gt = new double[6];
rasterDataset->GetGeoTransform(gt);
this->gt = Eigen::MatrixXd(2, 3);
this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5];
double a = this->gt(0, 0);
double b = this->gt(0, 1);
double c = this->gt(0, 2);
double d = this->gt(1, 0);
double e = this->gt(1, 1);
double f = this->gt(1, 2);
this->projection = rasterDataset->GetProjectionRef();
// 释放投影
GDALFlushCache((GDALDatasetH)rasterDataset);
GDALClose((GDALDatasetH)rasterDataset);
rasterDataset = NULL; // 指矫匡拷
this->InitInv_gt();
delete[] gt;
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
}
gdalImageComplex::~gdalImageComplex() {}
void gdalImageComplex::setData(Eigen::MatrixXcd data)
{
this->data = 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);
if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
QString tip = u8"file path: " + this->img_path;
qDebug() << tip;
throw std::exception(tip.toUtf8().constData());
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
GDALDataset* poDstDS = nullptr;
if (exists_test(this->img_path)) {
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_CFloat64, NULL); // 斤拷锟斤拷
poDstDS->SetProjection(this->projection.toUtf8().constData());
double gt_ptr[6];
for (int i = 0; i < this->gt.rows(); i++) {
for (int j = 0; j < this->gt.cols(); j++) {
gt_ptr[i * 3 + j] = this->gt(i, j);
}
}
poDstDS->SetGeoTransform(gt_ptr);
//delete[] gt_ptr;
}
int datarows = data.rows();
int datacols = data.cols();
double* databuffer = new double[data.size() * 2];
for (int i = 0; i < data.rows(); i++) {
for (int j = 0; j < data.cols(); j++) {
databuffer[i * data.cols() * 2 + j * 2] = data(i, j).real();
databuffer[i * data.cols() * 2 + j * 2 + 1] = data(i, j).imag();
}
}
// 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_CFloat64, 0, 0);
GDALFlushCache(poDstDS);
delete databuffer;
GDALClose((GDALDatasetH)poDstDS);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
void gdalImageComplex::saveImage(std::shared_ptr<std::complex<double>> data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
if (start_row + rowCount > this->height || start_col + colCount > this->width) {
QString tip = u8"file path: " + this->img_path;
qDebug() << tip;
throw std::exception(tip.toUtf8().constData());
return;
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
GDALDataset* poDstDS = nullptr;
if (exists_test(this->img_path)) {
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_CFloat64, NULL); // 斤拷锟斤拷
poDstDS->SetProjection(this->projection.toUtf8().constData());
double gt_ptr[6];
for (int i = 0; i < this->gt.rows(); i++) {
for (int j = 0; j < this->gt.cols(); j++) {
gt_ptr[i * 3 + j] = this->gt(i, j);
}
}
poDstDS->SetGeoTransform(gt_ptr);
//delete[] gt_ptr;
}
double* databuffer = new double[rowCount * colCount * 2];
for (long i = 0; i < rowCount; i++) {
for (long j = 0; j < colCount; j++) {
databuffer[i * colCount * 2 + j * 2] = data.get()[i * colCount + j].real();
databuffer[i * colCount * 2 + j * 2 + 1] = data.get()[i * colCount + j].imag();
}
}
// 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, colCount, rowCount,
databuffer, colCount, rowCount, GDT_CFloat64, 0, 0);
GDALFlushCache(poDstDS);
delete databuffer;
GDALClose((GDALDatasetH)poDstDS);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
void gdalImageComplex::saveImage(std::complex<double>* data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
if (start_row + rowCount > this->height || start_col + colCount > this->width) {
QString tip = u8"file path: " + this->img_path;
qDebug() << tip;
throw std::exception(tip.toUtf8().constData());
return;
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
GDALDataset* poDstDS = nullptr;
if (exists_test(this->img_path)) {
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_CFloat64, NULL); // 斤拷锟斤拷
poDstDS->SetProjection(this->projection.toUtf8().constData());
double gt_ptr[6];
for (int i = 0; i < this->gt.rows(); i++) {
for (int j = 0; j < this->gt.cols(); j++) {
gt_ptr[i * 3 + j] = this->gt(i, j);
}
}
poDstDS->SetGeoTransform(gt_ptr);
//delete[] gt_ptr;
}
double* databuffer = new double[rowCount * colCount * 2];
for (long i = 0; i < rowCount; i++) {
for (long j = 0; j < colCount; j++) {
databuffer[i * colCount * 2 + j * 2] = data[i * colCount + j].real();
databuffer[i * colCount * 2 + j * 2 + 1] = data[i * colCount + j].imag();
}
}
// 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, colCount, rowCount,
databuffer, colCount, rowCount, GDT_CFloat64, 0, 0);
GDALFlushCache(poDstDS);
delete databuffer;
GDALClose((GDALDatasetH)poDstDS);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
}
Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col, int rows_count,
int cols_count, int band_ids)
{
GDALDataset* poDataset;
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// 打开TIFF文件
poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
if (poDataset == nullptr) {
QMessageBox::warning(nullptr, u8"错误", u8"无法打开文件：" + this->img_path);
qDebug() << u8"无法打开文件：" + this->img_path;
}
// 获取数据集的第一个波段
GDALRasterBand* poBand;
poBand = poDataset->GetRasterBand(1);
// 读取波段信息，假设是复数类型
int nXSize = cols_count; poBand->GetXSize();
int nYSize = rows_count; poBand->GetYSize();
double* databuffer = new double[nXSize * nYSize * 2];
poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count,
rows_count, GDT_CFloat64, 0, 0);
GDALClose((GDALDatasetH)poDataset);
Eigen::MatrixXcd rasterData(nYSize, nXSize); // 使用Eigen的MatrixXcd
for (size_t i = 0; i < nYSize; i++) {
for (size_t j = 0; j < nXSize; j++) {
rasterData(i, j) = std::complex<double>(databuffer[i * nXSize * 2 + j * 2],
databuffer[i * nXSize * 2 + j * 2 + 1]);
}
}
delete[] databuffer;
return rasterData;
}
std::shared_ptr<std::complex<double>> gdalImageComplex::getDataComplexSharePtr(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
{
GDALDataset* poDataset;
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// 打开TIFF文件
poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
if (poDataset == nullptr) {
QMessageBox::warning(nullptr, u8"错误", u8"无法打开文件：" + this->img_path);
qDebug() << u8"无法打开文件：" + this->img_path;
}
// 获取数据集的第一个波段
GDALRasterBand* poBand;
poBand = poDataset->GetRasterBand(1);
// 读取波段信息，假设是复数类型
double* databuffer = new double[rows_count * cols_count * 2];
poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count,
rows_count, GDT_CFloat64, 0, 0);
GDALClose((GDALDatasetH)poDataset);
std::shared_ptr<std::complex<double>> rasterData(new std::complex<double>[rows_count * cols_count], delArrPtr);
for (size_t i = 0; i < rows_count; i++) {
for (size_t j = 0; j < cols_count; j++) {
rasterData.get()[i * cols_count + j] = std::complex<double>(databuffer[i * cols_count * 2 + j * 2],
databuffer[i * cols_count * 2 + j * 2 + 1]);
}
}
delete[] databuffer;
return rasterData;
}
void gdalImageComplex::saveImage()
{
this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
}
void gdalImageComplex::savePreViewImage()
{
qDebug() << "void gdalImageComplex::savePreViewImage()";
Eigen::MatrixXd data_abs = Eigen::MatrixXd::Zero(this->height, this->width);
data_abs = (this->data.array().real().pow(2) + this->data.array().imag().pow(2))
.array()
.log10() * 10.0; // 计算振幅
double min_abs = data_abs.minCoeff(); // 最大值
double max_abs = data_abs.maxCoeff(); // 最小值
double delta = (max_abs - min_abs) / 1000; // 1000分位档
Eigen::MatrixX<size_t> data_idx =
((data_abs.array() - min_abs).array() / delta).array().floor().cast<size_t>();
// 初始化
double hist[1001];
for (size_t i = 0; i < 1001; i++) {
hist[i] = 0; // 初始化
}
for (size_t i = 0; i < this->height; i++) {
for (size_t j = 0; j < this->width; j++) {
hist[data_idx(i, j)]++;
}
}
// 统计
size_t count = this->height * this->width;
double precent = 0;
size_t curCount = 0;
double pre2 = 0;
bool findprec_2 = true;
double pre98 = 0;
bool findprec_98 = true;
for (size_t i = 0; i < 1001; i++) {
precent = precent + hist[i];
if (findprec_2 & precent / count > 0.02) {
pre2 = i * delta + min_abs;
findprec_2 = false;
}
if (findprec_98 & precent / count > 0.98) {
pre98 = (i - 1) * delta + min_abs;
findprec_98 = false;
}
}
// 拉伸
delta = (pre98 - pre2) / 200;
data_idx =
((data_abs.array() - pre2).array() / delta).array().floor().cast<size_t>();
for (size_t i = 0; i < this->height; i++) {
for (size_t j = 0; j < this->width; j++) {
if (data_idx(i, j) < 0) {
data_idx(i, j) = 0;
}
else if (data_idx(i, j) > 255) {
data_idx(i, j) = 255;
}
else {
}
}
}
// 赋值
QString filePath = this->img_path;
QFile file(filePath);
QFileInfo fileInfo(file);
QString directory = fileInfo.absolutePath();
qDebug() << "文件所在目录：" << directory;
QString baseName = fileInfo.completeBaseName();
qDebug() << "无后缀文件名：" << baseName;
// 创建文件路径
QString previewImagePath = JoinPath(directory, baseName + "_preview.png");
cv::Mat img(this->height, this->width, CV_8U, cv::Scalar(0));
for (size_t i = 0; i < this->height; i++) {
for (size_t j = 0; j < this->width; j++) {
img.at<uchar>(i, j) = (uchar)(data_idx(i, j));
}
}
//std::string outimgpath=previewImagePath.toUtf8().data();
cv::imwrite(previewImagePath.toUtf8().data(), img);
}
gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width,
int band_num, Eigen::MatrixXd gt, QString projection,
bool need_gt, bool overwrite)
{
if (exists_test(img_path.toUtf8().constData())) {
if (overwrite) {
gdalImageComplex result_img(img_path);
return result_img;
}
else {
throw "file has exist!!!";
exit(1);
}
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
GDT_CFloat64, NULL);
if (need_gt) {
poDstDS->SetProjection(projection.toUtf8().constData());
// 锟斤拷锟斤拷转锟斤拷锟斤拷锟斤拷
double gt_ptr[6] = { 0 };
for (int i = 0; i < gt.rows(); i++) {
for (int j = 0; j < gt.cols(); j++) {
gt_ptr[i * 3 + j] = gt(i, j);
}
}
poDstDS->SetGeoTransform(gt_ptr);
}
//for(int i = 1; i <= band_num; i++) {
// poDstDS->GetRasterBand(i)->SetNoDataValue(0); // 回波部分
//}
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImageComplex result_img(img_path);
return result_img;
}
gdalImageComplex CreategdalImageComplexNoProj(const QString& img_path, int height, int width, int band_num, bool overwrite)
{
// 创建新文件
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
GDALDataset* poDstDS = (poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_CFloat64, NULL));
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
//poDstDS.reset();
omp_unset_lock(&lock); //
omp_destroy_lock(&lock); //
gdalImageComplex result_img(img_path);
return result_img;
}
gdalImageComplex CreateEchoComplex(const QString& img_path, int height, int width, int band_num)
{
// 创建图像
Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
//Xgeo = GeoTransform[0] + Xpixel * GeoTransform[1] + Ypixel * GeoTransform[2]
//Ygeo = GeoTransform[3] + Xpixel * GeoTransform[4] + Ypixel * GeoTransform[5]
// X
gt(0, 0) = 0; gt(0, 2) = 1; gt(0, 2) = 0;
gt(1, 0) = 0; gt(1, 1) = 0; gt(1, 2) = 1;
// Y
QString projection = "";
gdalImageComplex echodata = CreategdalImageComplex(img_path, height, width, 1, gt, projection, false, true);
return echodata;
}

1477
BaseCommonLibrary/BaseTool/gdalImageOperator.cpp Normal file
View File

File diff suppressed because it is too large Load Diff