LAMPCAE/src/PluginWBFZExchangePlugin/SARCalibration.cpp

385 lines
19 KiB
C++
Raw Normal View History

2024-04-07 16:19:33 +00:00
#include "SARCalibration.h"
#include "ImageOperatorBase.h"
#include "SARBaseTool.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <omp.h>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <gdal.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <proj.h>
#include <string.h>
#include <memory.h>
#include <memory>
using namespace std;
using namespace Eigen;
/**
* ENVI
*/
Eigen::Matrix2cd CalibrationMatrix(Eigen::MatrixXcd in_matrix, double calibrationValue)
2024-04-07 16:19:33 +00:00
{
return in_matrix.array() * calibrationValue;
}
int CalibrationSiglePolarSAR(QString out_path, QString in_sar_path, double calibrationValue)
2024-04-07 16:19:33 +00:00
{
return CalibrationComplex(out_path, in_sar_path, calibrationValue);
}
int CalibrationComplex(const QString& out_path, const QString& in_sar_path, double calibrationValue)
2024-04-07 16:19:33 +00:00
{
GDALAllRegister();
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(in_sar_path);
int width = rasterDataset->GetRasterXSize();
int height = rasterDataset->GetRasterYSize();
int band_num = rasterDataset->GetRasterCount();
double* gt = new double[6];
std::shared_ptr<double> gt_ptr(gt);
QString projectDef = rasterDataset->GetProjectionRef();
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
bool _nGt_flag = false;
if (projectDef == "") {
_nGt_flag = true;
}
else {
_nGt_flag = false;
}
CreateDataset(out_path, height, width, 1, gt_ptr.get(), projectDef, GDT_CFloat64, _nGt_flag);
// 计算大小
if (gdal_datatype == 0) {
return 1;
}
else if (gdal_datatype < 8) {
if (band_num != 2) { return 1; }
}
else if (gdal_datatype < 12) {
if (band_num != 1) { return 1; }
}
else {}
int block_num = block_num_pre_memory(width, height, gdal_datatype, 2e9);//
block_num = block_num > height ? height : block_num; // 行数
int line_num = block_num;
for (int i = 0; i < height; i= block_num+i) {
if (height - i < block_num) {
line_num = height - i;
}
else {}
// 构建矩阵块使用eigen 进行矩阵计算,加速计算
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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
// rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
// data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
// data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
// rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
// data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
// data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) =( imag_mat.array().cast<double>() * calibrationValue).array();
_flag = true;
}
else if (gdal_datatype == GDT_CInt16) {
Eigen::MatrixX<complex<short>> complex_short_mat(line_num*width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, i, 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>() * calibrationValue).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>() * calibrationValue).array();
_flag = true;
}
else if (gdal_datatype == GDT_CInt32) {
Eigen::MatrixX<complex<int>> complex_short_mat(line_num*width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, i, 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>() * calibrationValue).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>() * calibrationValue).array();
_flag = true;
}
else if (gdal_datatype == GDT_CFloat32) {
Eigen::MatrixX<complex<float>> complex_short_mat(line_num*width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, i, 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>() * calibrationValue).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>() * calibrationValue).array();
_flag = true;
}
else if (gdal_datatype == GDT_CFloat64) {
Eigen::MatrixX<complex<double>> complex_short_mat(line_num*width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, i, 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>() * calibrationValue).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>() * calibrationValue).array();
_flag = true;
}
else {}
// 保存数据
if (_flag) {
// 定义赋值矩阵
saveDataset(out_path, i,0,1, width, line_num, data_mat.data());
}
else {
}
}
return -1;
}
int CalibrationComplex2dB(const QString& out_path, const QString& in_sar_path, double calibrationValue)
2024-04-07 16:19:33 +00:00
{
GDALAllRegister();
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(in_sar_path);
int width = rasterDataset->GetRasterXSize();
int height = rasterDataset->GetRasterYSize();
int band_num = rasterDataset->GetRasterCount();
double* gt = new double[6];
std::shared_ptr<double> gt_ptr(gt);
QString projectDef = rasterDataset->GetProjectionRef();
GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
bool _nGt_flag = false;
if (projectDef == "") {
_nGt_flag = true;
}
else {
_nGt_flag = false;
}
CreateDataset(out_path, height, width, 1, gt_ptr.get(), projectDef, GDT_CFloat64, _nGt_flag);
// 计算大小
if (gdal_datatype == 0) {
return 1;
}
else if (gdal_datatype < 8) {
if (band_num != 2) { return 1; }
}
else if (gdal_datatype < 12) {
if (band_num != 1) { return 1; }
}
else {}
int block_num = block_num_pre_memory(width, height, gdal_datatype, 2e9);//
block_num = block_num > height ? height : block_num; // 行数
int line_num = block_num;
for (int i = 0; i < height; i = block_num + i) {
if (height - i < block_num) {
line_num = height - i;
}
else {}
// 构建矩阵块使用eigen 进行矩阵计算,加速计算
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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
// rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
// data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
// data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
// rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
// data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
// data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).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, i, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, i, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag
data_mat.col(0) = (real_mat.array().cast<double>() * calibrationValue).array();
data_mat.col(1) = (imag_mat.array().cast<double>() * calibrationValue).array();
_flag = true;
}
else if (gdal_datatype == GDT_CInt16) {
Eigen::MatrixX<complex<short>> complex_short_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, i, 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>() * calibrationValue).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>() * calibrationValue).array();
_flag = true;
}
else if (gdal_datatype == GDT_CInt32) {
Eigen::MatrixX<complex<int>> complex_short_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, i, 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>() * calibrationValue).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>() * calibrationValue).array();
_flag = true;
}
else if (gdal_datatype == GDT_CFloat32) {
Eigen::MatrixX<complex<float>> complex_short_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, i, 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>() * calibrationValue).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>() * calibrationValue).array();
_flag = true;
}
else if (gdal_datatype == GDT_CFloat64) {
Eigen::MatrixX<complex<double>> complex_short_mat(line_num * width, 1);
rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, i, 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>() * calibrationValue).array();
data_mat.col(1) = (complex_short_mat.imag().array().cast<double>() * calibrationValue).array();
_flag = true;
}
else {}
// 保存数据
if (_flag) {
// 定义赋值矩阵
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> data_out(line_num* width, 1);// 必须强制行优先
data_out.col(0) = Complex2dB(data_mat.col(0).array(), data_mat.col(1).array()).array();
saveDataset(out_path, i, 0, 1, width, line_num, data_out.data());
}
else {
}
}
return -1;
return 0;
}
int CalibrationComplex_DLL(const QString& out_path, const QString& in_sar_path, double calibrationValue)
2024-04-07 16:19:33 +00:00
{
return 0;
}
int CalibrationComplex2dB_DLL(const QString& out_path, const QString& in_sar_path, double calibrationValue)
2024-04-07 16:19:33 +00:00
{
return 0;
}