LAMPCAE/src/PluginWBFZExchangePlugin/SARCalibration.cpp

389 lines
20 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

#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)
{
return in_matrix.array() * calibrationValue;
}
int CalibrationSiglePolarSAR(QString out_path, QString in_sar_path, double calibrationValue)
{
return CalibrationComplex(out_path, in_sar_path, calibrationValue);
}
int CalibrationComplex(const QString& out_path, const QString& in_sar_path, double calibrationValue)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
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)
{
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
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)
{
return 0;
}
int CalibrationComplex2dB_DLL(const QString& out_path, const QString& in_sar_path, double calibrationValue)
{
return 0;
}