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

1513 lines
53 KiB
C++
Raw Normal View History

代码初始化
2024-11-15 01:38:46 +00:00
#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.h>
#include <memory>
#include <iostream>
#include "FileOperator.h"
#include <opencv2/opencv.hpp>
#include <QMessageBox>
#include <QDir>
/**
* 输入数据是ENVI格式数据
*/
std::shared_ptr<GDALDataset> OpenDataset(const QString& in_path, GDALAccess rwmode)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* dataset_ptr = (GDALDataset*)(GDALOpen(in_path.toUtf8().constData(), rwmode));
std::shared_ptr<GDALDataset> rasterDataset(dataset_ptr, CloseDataset);
return rasterDataset;
}
void CloseDataset(GDALDataset* ptr)
{
GDALClose(ptr);
ptr = NULL;
}
int TIFF2ENVI(QString in_tiff_path, QString out_envi_path)
{
std::shared_ptr<GDALDataset> ds = OpenDataset(in_tiff_path);
const char* args[] = { "-of", "ENVI", NULL };
GDALTranslateOptions* psOptions = GDALTranslateOptionsNew((char**)args, NULL);
GDALClose(GDALTranslate(out_envi_path.toUtf8().constData(), ds.get(), psOptions, NULL));
GDALTranslateOptionsFree(psOptions);
return 0;
}
int ENVI2TIFF(QString in_envi_path, QString out_tiff_path)
{
std::shared_ptr<GDALDataset> ds = OpenDataset(in_envi_path);
const char* args[] = { "-of", "Gtiff", NULL };
GDALTranslateOptions* psOptions = GDALTranslateOptionsNew((char**)args, NULL);
GDALClose(GDALTranslate(out_tiff_path.toUtf8().constData(), ds.get(), psOptions, NULL));
GDALTranslateOptionsFree(psOptions);
return 0;
}
int CreateDataset(QString new_file_path, int height, int width, int band_num, double* gt,
QString projection, GDALDataType gdal_dtype, bool need_gt)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(new_file_path.toUtf8().constData(), width,
height, band_num, gdal_dtype, NULL));
if(need_gt) {
poDstDS->SetProjection(projection.toUtf8().constData());
poDstDS->SetGeoTransform(gt);
} else {
}
GDALFlushCache((GDALDatasetH)poDstDS.get());
return 0;
}
int saveDataset(QString new_file_path, int start_line, int start_cols, int band_ids, int datacols,
int datarows, void* databuffer)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
std::shared_ptr<GDALDataset> poDstDS = OpenDataset(new_file_path, GA_Update);
GDALDataType gdal_datatype = poDstDS->GetRasterBand(1)->GetRasterDataType();
poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_cols, start_line, datacols, datarows,
databuffer, datacols, datarows, gdal_datatype, 0, 0);
GDALFlushCache(poDstDS.get());
return 0;
}
int block_num_pre_memory(int block_width, int height, GDALDataType gdal_datatype, double memey_size)
{
// 计算大小
int size_meta = 0;
if(gdal_datatype == GDT_Byte) {
size_meta = 1;
} else if(gdal_datatype == GDT_UInt16) {
size_meta = 2; // 只有双通道才能构建 复数矩阵
} else if(gdal_datatype == GDT_UInt16) {
size_meta = 2;
} else if(gdal_datatype == GDT_Int16) {
size_meta = 2;
} else if(gdal_datatype == GDT_UInt32) {
size_meta = 4;
} else if(gdal_datatype == GDT_Int32) {
size_meta = 4;
}
// else if (gdal_datatype == GDT_UInt64) {
// size_meta = 8;
// }
// else if (gdal_datatype == GDT_Int64) {
// size_meta = 8;
// }
else if(gdal_datatype == GDT_Float32) {
size_meta = 4;
} else if(gdal_datatype == GDT_Float64) {
size_meta = 4;
} else if(gdal_datatype == GDT_CInt16) {
size_meta = 2;
} else if(gdal_datatype == GDT_CInt32) {
size_meta = 2;
} else if(gdal_datatype == GDT_CFloat32) {
size_meta = 4;
} else if(gdal_datatype == GDT_CFloat64) {
size_meta = 8;
} else {
}
int block_num = int(memey_size / (size_meta * block_width));
block_num = block_num > height ? height : block_num; // 行数
block_num = block_num < 1 ? 1 : block_num;
return block_num;
}
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
ReadComplexMatrixData(int start_line, int width, int line_num,
std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype)
{
int band_num = rasterDataset->GetRasterCount();
if(gdal_datatype == 0) {
return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
} else if(gdal_datatype < 8) {
if(band_num != 2) {
return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
}
} else if(gdal_datatype < 12) {
if(band_num != 1) {
return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
}
} else {
}
bool _flag = false;
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> data_mat(
line_num * width, 2); // 必须强制行优先
if(gdal_datatype == GDT_Byte) {
Eigen::MatrixX<char> real_mat(line_num * width, 1);
Eigen::MatrixX<char> imag_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num, gdal_datatype,
0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
imag_mat.data(), width, line_num, gdal_datatype,
0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_UInt16) {
Eigen::MatrixX<unsigned short> real_mat(line_num * width, 1);
Eigen::MatrixX<unsigned short> imag_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num, gdal_datatype,
0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
imag_mat.data(), width, line_num, gdal_datatype,
0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_Int16) {
Eigen::MatrixX<short> real_mat(line_num * width, 1);
Eigen::MatrixX<short> imag_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num, gdal_datatype,
0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
imag_mat.data(), width, line_num, gdal_datatype,
0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_UInt32) {
Eigen::MatrixX<unsigned int> real_mat(line_num * width, 1);
Eigen::MatrixX<unsigned int> imag_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num, gdal_datatype,
0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
imag_mat.data(), width, line_num, gdal_datatype,
0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_Int32) {
Eigen::MatrixX<int> real_mat(line_num * width, 1);
Eigen::MatrixX<int> imag_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num, gdal_datatype,
0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
imag_mat.data(), width, line_num, gdal_datatype,
0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
}
// else if (gdal_datatype == GDT_UInt64) {
// Eigen::MatrixX<unsigned long> real_mat(line_num * width, 1);
// Eigen::MatrixX<unsigned long> imag_mat(line_num * width, 1);
// rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
//real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
// rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
//imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) =
//(real_mat.array().cast<double>()).array(); data_mat.col(1) =
//(imag_mat.array().cast<double>()).array(); _flag = true;
// }
// else if (gdal_datatype == GDT_Int64) {
// Eigen::MatrixX<long> real_mat(line_num * width, 1);
// Eigen::MatrixX<long> imag_mat(line_num * width, 1);
// rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
//real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
// rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
//imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag data_mat.col(0) =
//(real_mat.array().cast<double>()).array(); data_mat.col(1) =
//(imag_mat.array().cast<double>()).array(); _flag = true;
// }
else if(gdal_datatype == GDT_Float32) {
Eigen::MatrixX<float> real_mat(line_num * width, 1);
Eigen::MatrixX<float> imag_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num, gdal_datatype,
0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
imag_mat.data(), width, line_num, gdal_datatype,
0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_Float64) {
Eigen::MatrixX<double> real_mat(line_num * width, 1);
Eigen::MatrixX<double> imag_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num, gdal_datatype,
0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
imag_mat.data(), width, line_num, gdal_datatype,
0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>()).array();
data_mat.col(1) = (imag_mat.array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_CInt16) {
Eigen::MatrixX<std::complex<short>> complex_short_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
complex_short_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = (complex_short_mat.real().array().cast<double>()).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_CInt32) {
Eigen::MatrixX<std::complex<int>> complex_short_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
complex_short_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = (complex_short_mat.real().array().cast<double>()).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_CFloat32) {
Eigen::MatrixX<std::complex<float>> complex_short_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
complex_short_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = (complex_short_mat.real().array().cast<double>()).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
_flag = true;
} else if(gdal_datatype == GDT_CFloat64) {
Eigen::MatrixX<std::complex<double>> complex_short_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
complex_short_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = (complex_short_mat.real().array().cast<double>()).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
_flag = true;
} else {
}
// 保存数据
if(_flag) {
return data_mat;
} else {
return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
0, 0); // 必须强制行优先;
}
}
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
ReadMatrixDoubleData(int start_line, int width, int line_num,
std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype,
int band_idx)
{
// 构建矩阵块，使用eigen 进行矩阵计算，加速计算
bool _flag = false;
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> data_mat(
line_num * width, 1); // 必须强制行优先
if(gdal_datatype == GDT_Byte) {
Eigen::MatrixX<char> real_mat(line_num * width, 1);
rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_UInt16) {
Eigen::MatrixX<unsigned short> real_mat(line_num * width, 1);
rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_Int16) {
Eigen::MatrixX<short> real_mat(line_num * width, 1);
rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_UInt32) {
Eigen::MatrixX<unsigned int> real_mat(line_num * width, 1);
rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_Int32) {
Eigen::MatrixX<int> real_mat(line_num * width, 1);
rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
}
// else if (gdal_datatype == GDT_UInt64) {
// Eigen::MatrixX<unsigned long> real_mat(line_num * width, 1);
// rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
//real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) =
//((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0; _flag = true;
// }
// else if (gdal_datatype == GDT_Int64) {
// Eigen::MatrixX<long> real_mat(line_num * width, 1);
// rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
//real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real data_mat.col(0) =
//((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0; _flag = true;
// }
else if(gdal_datatype == GDT_Float32) {
Eigen::MatrixX<float> real_mat(line_num * width, 1);
rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else if(gdal_datatype == GDT_Float64) {
Eigen::MatrixX<double> real_mat(line_num * width, 1);
rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
real_mat.data(), width, line_num,
gdal_datatype, 0, 0); // real
data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
_flag = true;
} else {
}
return data_mat;
}
Eigen::MatrixXd getGeoTranslationArray(QString in_path)
{
return Eigen::MatrixXd();
}
ImageGEOINFO getImageINFO(QString in_path)
{
std::shared_ptr<GDALDataset> df = OpenDataset(in_path);
int width = df->GetRasterXSize();
int heigh = df->GetRasterYSize();
int band_num = df->GetRasterCount();
ImageGEOINFO result;
result.width = width;
result.height = heigh;
result.bandnum = band_num;
return result;
}
GDALDataType getGDALDataType(QString fileptah)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
fileptah.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
return gdal_datatype;
}
gdalImage::gdalImage()
{
this->height = 0;
this->width = 0;
this->data_band_ids = 1;
this->start_row = 0;
this->start_col = 0;
}
/// <summary>
/// 斤拷图饺∮帮拷锟?1锟?7
/// </summary>
/// <param name="dem_path"></param>
gdalImage::gdalImage(const QString& raster_path)
{
omp_lock_t lock;
omp_init_lock(&lock); // 锟斤拷始斤拷斤拷
omp_set_lock(&lock); // 锟斤拷没斤拷锟?1锟?7
this->img_path = raster_path;
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
// 锟斤拷DEM影锟斤拷
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];
this->projection = rasterDataset->GetProjectionRef();
// 斤拷斤拷
// double* inv_gt = new double[6];;
// GDALInvGeoTransform(gt, inv_gt); // 斤拷斤拷
// 斤拷投影
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); // 劫伙拷斤拷
}
gdalImage::~gdalImage() {}
void gdalImage::setHeight(int height)
{
this->height = height;
}
void gdalImage::setWidth(int width)
{
this->width = width;
}
void gdalImage::setTranslationMatrix(Eigen::MatrixXd gt)
{
this->gt = gt;
}
void gdalImage::setData(Eigen::MatrixXd, int data_band_ids)
{
this->data = data;
this->data_band_ids = data_band_ids;
}
Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, int cols_count,
int band_ids = 1)
{
omp_lock_t lock;
omp_init_lock(&lock);
omp_set_lock(&lock);
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
this->img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids);
rows_count = start_row + rows_count <= this->height ? rows_count : this->height - start_row;
cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;
Eigen::MatrixXd datamatrix(rows_count, cols_count);
if(gdal_datatype == GDT_Byte) {
char* temp = new char[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_UInt16) {
unsigned short* temp = new unsigned short[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_Int16) {
short* temp = new short[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_UInt32) {
unsigned int* temp = new unsigned int[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_Int32) {
int* temp = new int[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,
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_Float64) {
double* temp = new double[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 {
}
GDALClose((GDALDatasetH)rasterDataset);
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return datamatrix;
}
Eigen::MatrixXd gdalImage::getGeoTranslation()
{
return this->gt;
}
GDALDataType gdalImage::getDataType()
{
GDALDataset* rasterDataset =
(GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
return gdal_datatype;
}
/// <summary>
///
/// </summary>
/// <param name="data"></param>
/// <param name="start_row"></param>
/// <param name="start_col"></param>
/// <param name="band_ids"></param>
void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0,
int band_ids = 1)
{
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+" 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());
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
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_Float32, 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();
float* databuffer =
new float[datarows * datacols]; // (float*)malloc(datarows * datacols * sizeof(float));
for(int i = 0; i < datarows; i++) {
for(int j = 0; j < datacols; j++) {
float temp = float(data(i, j));
databuffer[i * datacols + j] = temp;
}
}
// 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);
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
delete[] databuffer;
omp_unset_lock(&lock); // 锟酵放伙拷斤拷
omp_destroy_lock(&lock); // 劫伙拷斤拷
}
void gdalImage::saveImage()
{
this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
}
void gdalImage::setNoDataValue(double nodatavalue = -9999, int band_ids = 1)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
// GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue);
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
}
int gdalImage::InitInv_gt()
{
// 1 lon lat = x
// 1 lon lat = y
Eigen::MatrixXd temp = Eigen::MatrixXd::Zero(2, 3);
this->inv_gt = temp;
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);
double g = 1;
double det_gt = b * f - c * e;
if(det_gt == 0) {
return 0;
}
this->inv_gt(0, 0) = (c * d - a * f) / det_gt; // 2
this->inv_gt(0, 1) = f / det_gt; // lon
this->inv_gt(0, 2) = -c / det_gt; // lat
this->inv_gt(1, 0) = (a * e - b * d) / det_gt; // 1
this->inv_gt(1, 1) = -e / det_gt; // lon
this->inv_gt(1, 2) = b / det_gt; // lat
return 1;
}
Landpoint gdalImage::getRow_Col(double lon, double lat)
{
Landpoint p{ 0, 0, 0 };
p.lon = this->inv_gt(0, 0) + lon * this->inv_gt(0, 1) + lat * this->inv_gt(0, 2); // x
p.lat = this->inv_gt(1, 0) + lon * this->inv_gt(1, 1) + lat * this->inv_gt(1, 2); // y
return p;
}
Landpoint gdalImage::getLandPoint(double row, double col, double ati = 0)
{
Landpoint p{ 0, 0, 0 };
p.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); // x
p.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); // y
p.ati = ati;
return p;
}
void gdalImage::getLandPoint(double row, double col, double ati, Landpoint& Lp)
{
Lp.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); // x
Lp.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); // y
Lp.ati = ati;
}
double gdalImage::mean(int bandids)
{
double mean_value = 0;
double count = this->height * this->width;
int line_invert = 100;
int start_ids = 0;
do {
Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
mean_value = mean_value + sar_a.sum() / count;
start_ids = start_ids + line_invert;
} while(start_ids < this->height);
return mean_value;
}
double gdalImage::BandmaxValue(int bandids)
{
double max_value = 0;
bool state_max = true;
int line_invert = 100;
int start_ids = 0;
double temp_max = 0;
do {
Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
if(state_max) {
state_max = false;
max_value = sar_a.maxCoeff();
} else {
temp_max = sar_a.maxCoeff();
if(max_value < temp_max) {
max_value = temp_max;
}
}
start_ids = start_ids + line_invert;
} while(start_ids < this->height);
return max_value;
}
double gdalImage::BandminValue(int bandids)
{
double min_value = 0;
bool state_min = true;
int line_invert = 100;
int start_ids = 0;
double temp_min = 0;
do {
Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
if(state_min) {
state_min = false;
min_value = sar_a.minCoeff();
} else {
temp_min = sar_a.minCoeff();
if(min_value < temp_min) {
min_value = temp_min;
}
}
start_ids = start_ids + line_invert;
} while(start_ids < this->height);
return min_value;
}
GDALRPCInfo gdalImage::getRPC()
{
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
CPLSetConfigOption("GDAL_DATA", "./data");
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注斤拷锟斤拷
// 斤拷锟斤拷
GDALDataset* pDS = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
// 锟斤拷元斤拷锟叫伙拷取RPC锟斤拷息
char** papszRPC = pDS->GetMetadata("RPC");
// 斤拷取锟斤拷RPC锟斤拷息斤拷山峁癸拷锟?1锟?7
GDALRPCInfo oInfo;
GDALExtractRPCInfo(papszRPC, &oInfo);
GDALClose((GDALDatasetH)pDS);
return oInfo;
}
Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points)
{
if(points.cols() != 3) {
throw new std::exception("the size of points is equit 3!!!");
}
Eigen::MatrixXd result(points.rows(), 3);
result.col(2) = points.col(2); // 锟竭筹拷
points.col(2) = points.col(2).array() * 0 + 1;
points.col(0).swap(points.col(2)); // 斤拷
Eigen::MatrixXd gts(3, 2);
gts.col(0) = this->gt.row(0);
gts.col(1) = this->gt.row(1);
result.block(0, 0, points.rows(), 2) = points * gts;
return result;
}
Eigen::MatrixXd gdalImage::getHist(int bandids)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
// 锟斤拷DEM影锟斤拷
GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
this->img_path.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* xBand = rasterDataset->GetRasterBand(bandids);
double dfMin = this->BandminValue(bandids);
double dfMax = this->BandmaxValue(bandids);
int count = int((dfMax - dfMin) / 0.01);
count = count > 255 ? count : 255;
GUIntBig* panHistogram = new GUIntBig[count];
xBand->GetHistogram(dfMin, dfMax, count, panHistogram, TRUE, FALSE, NULL, NULL);
Eigen::MatrixXd result(count, 2);
double delta = (dfMax - dfMin) / count;
for(int i = 0; i < count; i++) {
result(i, 0) = dfMin + i * delta;
result(i, 1) = double(panHistogram[i]);
}
delete[] panHistogram;
GDALClose((GDALDatasetH)rasterDataset);
return result;
}
gdalImage CreategdalImage(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) {
gdalImage result_img(img_path);
return result_img;
} else {
throw "file has exist!!!";
exit(1);
}
}
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
GDT_Float32, 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(-9999);
}
GDALFlushCache((GDALDatasetH)poDstDS);
GDALClose((GDALDatasetH)poDstDS);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
gdalImage result_img(img_path);
return result_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 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;
}
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width,
int new_height, GDALResampleAlg eResample)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly);
if(pDSrc == NULL) {
return -1;
}
GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if(pDriver == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -2;
}
int width = pDSrc->GetRasterXSize();
int height = pDSrc->GetRasterYSize();
int nBandCount = pDSrc->GetRasterCount();
GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();
char* pszSrcWKT = NULL;
pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef());
// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
if(strlen(pszSrcWKT) <= 0) {
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
// oSRS.SetUTM(50, true); //锟斤拷锟斤拷锟斤拷 锟斤拷锟斤拷120锟斤拷
// oSRS.SetWellKnownGeogCS("WGS84");
oSRS.exportToWkt(&pszSrcWKT);
}
qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
void* hTransformArg;
hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,
FALSE, 0.0, 1);
qDebug() << "no proj " << Qt::endl;
//(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7)
if(hTransformArg == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -3;
}
qDebug() << "has proj " << Qt::endl;
double dGeoTrans[6] = { 0 };
int nNewWidth = 0, nNewHeight = 0;
if(GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg,
dGeoTrans, &nNewWidth, &nNewHeight)
!= CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -3;
}
// GDALDestroyGenImgProjTransformer(hTransformArg);
GDALDataset* pDDst =
pDriver->Create(pszOutFile, new_width, new_height, nBandCount, dataType, NULL);
if(pDDst == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -2;
}
pDDst->SetProjection(pszSrcWKT);
pDDst->SetGeoTransform(gt);
GDALWarpOptions* psWo = GDALCreateWarpOptions();
// psWo->papszWarpOptions = CSLDuplicate(NULL);
psWo->eWorkingDataType = dataType;
psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH)pDSrc;
psWo->hDstDS = (GDALDatasetH)pDDst;
qDebug() << "GDALCreateGenImgProjTransformer" << Qt::endl;
psWo->pfnTransformer = GDALGenImgProjTransform;
psWo->pTransformerArg = GDALCreateGenImgProjTransformer(
(GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
;
qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
psWo->nBandCount = nBandCount;
psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
for(int i = 0; i < nBandCount; i++) {
psWo->panSrcBands[i] = i + 1;
psWo->panDstBands[i] = i + 1;
}
GDALWarpOperation oWo;
if(oWo.Initialize(psWo) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
return -3;
}
qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
GDALFlushCache(pDDst);
qDebug() << "ChunkAndWarpImage over" << Qt::endl;
// GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
// GDALDestroyWarpOptions(psWo);
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
return 0;
}
int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_Update);
if(pDSrc == NULL) {
return -1;
}
GDALDataType gdal_datatype = pDSrc->GetRasterBand(1)->GetRasterDataType();
GDALRasterBand* demBand = pDSrc->GetRasterBand(band_ids);
int width = pDSrc->GetRasterXSize();
int height = pDSrc->GetRasterYSize();
int start_col = 0, start_row = 0, rows_count = 0, cols_count;
int row_delta = int(120000000 / width);
GDALRasterIOExtraArg psExtraArg;
INIT_RASTERIO_EXTRA_ARG(psExtraArg);
psExtraArg.eResampleAlg = eResample;
do {
rows_count = start_row + row_delta < height ? row_delta : height - start_row;
cols_count = width;
if(gdal_datatype == GDALDataType::GDT_UInt16) {
unsigned short* temp = new unsigned short[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);
demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp,
cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg);
delete[] temp;
} else if(gdal_datatype == GDALDataType::GDT_Int16) {
short* temp = new short[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);
demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp,
cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg);
delete[] temp;
} else if(gdal_datatype == GDALDataType::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, rows_count, gdal_datatype, 0, 0);
demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp,
cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg);
delete[] temp;
}
start_row = start_row + rows_count;
} while(start_row < height);
GDALClose((GDALDatasetH)pDSrc);
return 0;
}
int alignRaster(QString inputPath, QString referencePath, QString outputPath, GDALResampleAlg eResample)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
GDALDataset* pDSrc = (GDALDataset*)GDALOpen(inputPath.toUtf8().constData(), GA_ReadOnly);
if (pDSrc == NULL) {
return -1;
}
GDALDataset* tDSrc = (GDALDataset*)GDALOpen(referencePath.toUtf8().constData(), GA_ReadOnly);
if (tDSrc == NULL) {
return -1;
}
long new_width = tDSrc->GetRasterXSize();
long new_height = tDSrc->GetRasterYSize();
GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
if (pDriver == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -2;
}
int width = pDSrc->GetRasterXSize();
int height = pDSrc->GetRasterYSize();
int nBandCount = pDSrc->GetRasterCount();
GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();
char* pszSrcWKT = NULL;
pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef());
if (strlen(pszSrcWKT) <= 0) {
OGRSpatialReference oSRS;
oSRS.importFromEPSG(4326);
oSRS.exportToWkt(&pszSrcWKT);
}
qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
void* hTransformArg;
hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,FALSE, 0.0, 1);
qDebug() << "no proj " << Qt::endl;
if (hTransformArg == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -3;
}
qDebug() << "has proj " << Qt::endl;
std::shared_ptr<double> dGeoTrans(new double[6], delArrPtr);
int nNewWidth = 0,
nNewHeight = 0;
if (GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg,dGeoTrans.get(), &nNewWidth, &nNewHeight) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -3;
}
GDALDataset* pDDst = pDriver->Create(outputPath.toUtf8().constData(), new_width, new_height, pDSrc->GetRasterCount(), dataType, NULL);
if (pDDst == NULL) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
return -2;
}
std::shared_ptr<double> gt(new double[6], delArrPtr);
tDSrc->GetGeoTransform(gt.get());
pDDst->SetProjection(pszSrcWKT);
pDDst->SetGeoTransform(gt.get());
GDALWarpOptions* psWo = GDALCreateWarpOptions();
psWo->eWorkingDataType = dataType;
psWo->eResampleAlg = eResample;
psWo->hSrcDS = (GDALDatasetH)pDSrc;
psWo->hDstDS = (GDALDatasetH)pDDst;
psWo->pfnTransformer = GDALGenImgProjTransform;
psWo->pTransformerArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
psWo->nBandCount = nBandCount;
psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
for (int i = 0; i < nBandCount; i++) {
psWo->panSrcBands[i] = i + 1;
psWo->panDstBands[i] = i + 1;
}
GDALWarpOperation oWo;
if (oWo.Initialize(psWo) != CE_None) {
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)tDSrc);
return -3;
}
qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
GDALFlushCache(pDDst);
qDebug() << "ChunkAndWarpImage over" << Qt::endl;
GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
GDALClose((GDALDatasetH)(GDALDatasetH)tDSrc);
return 0;
}
int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path)
{
int rows = data.rows();
int cols = data.cols();
Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
gdalImage image_tiff =
CreategdalImage(out_tiff_path, rows, cols, 2, gt, "", false, true); // 注意这里保留仿真结果
// 保存二进制文件
Eigen::MatrixXd real_img = data.array().real();
Eigen::MatrixXd imag_img = data.array().imag();
image_tiff.saveImage(real_img, 0, 0, 1);
image_tiff.saveImage(imag_img, 0, 0, 2);
return -1;
}
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); //
}
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 = poBand->GetXSize();
int nYSize = 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;
}
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);
}