284 lines
12 KiB
C++
284 lines
12 KiB
C++
// SARBaseTool.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
|
||
//
|
||
|
||
#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>
|
||
#include <fftw3.h>
|
||
|
||
|
||
|
||
int LAMPTOOLAPI Complex2dB(QString in_tiff, QString out_dB_path)
|
||
{
|
||
GDALAllRegister();
|
||
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(in_tiff);
|
||
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_dB_path, height, width, 1, gt_ptr.get(), projectDef, GDT_Float64, _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; }
|
||
}
|
||
|
||
int block_num = block_num_pre_memory(width, height, gdal_datatype, 1e9);//
|
||
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, 1);// 必须强制行优先
|
||
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>()).array().pow(2) + (imag_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);
|
||
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>()).array().pow(2) + (imag_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);
|
||
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>()).array().pow(2) + (imag_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);
|
||
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>()).array().pow(2) + (imag_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);
|
||
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>()).array().pow(2) + (imag_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);
|
||
// 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>()).array().pow(2) + (imag_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);
|
||
// 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>()).array().pow(2) + (imag_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);
|
||
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>()).array().pow(2) + (imag_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);
|
||
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>()).array().pow(2) + (imag_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
|
||
_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>().array().pow(2) + complex_short_mat.imag().array().cast<double>().pow(2)).log10() * 10.0;
|
||
_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>().array().pow(2) + complex_short_mat.imag().array().cast<double>().pow(2)).log10() * 10.0;
|
||
_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>().array().pow(2) + complex_short_mat.imag().array().cast<double>().pow(2)).log10() * 10.0;
|
||
_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>().array().pow(2) + complex_short_mat.imag().array().cast<double>().pow(2)).log10() * 10.0;
|
||
_flag = true;
|
||
}
|
||
else {}
|
||
// 保存数据
|
||
if (_flag) {
|
||
// 定义赋值矩阵
|
||
saveDataset(out_dB_path, i, 0, 1, width, line_num, data_mat.data());
|
||
}
|
||
else {
|
||
}
|
||
}
|
||
return 0;
|
||
}
|
||
|
||
int LAMPTOOLAPI Amplitude2dB(QString in_tiff, QString out_dB_path)
|
||
{
|
||
|
||
GDALAllRegister();
|
||
std::shared_ptr<GDALDataset> rasterDataset = OpenDataset(in_tiff);
|
||
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_dB_path, height, width, 1, gt_ptr.get(), projectDef, GDT_Float64, _nGt_flag);
|
||
// 计算大小
|
||
if (gdal_datatype == 0) {
|
||
return 1;
|
||
}
|
||
else if (gdal_datatype < 8) {
|
||
if (band_num != 1) { return 1; }
|
||
}
|
||
else { return 1; }
|
||
|
||
|
||
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 = ReadMatrixDoubleData(i, width, line_num, rasterDataset, gdal_datatype, 1);
|
||
_flag = (data_mat.rows() > 0) && (data_mat.cols() > 0);
|
||
// 保存数据
|
||
if (_flag) {
|
||
// 定义赋值矩阵
|
||
saveDataset(out_dB_path, i, 0, 1, width, line_num, data_mat.data());
|
||
}
|
||
else {
|
||
}
|
||
}
|
||
return 0;
|
||
}
|
||
|
||
Eigen::MatrixXd LAMPTOOLAPI Complex2dB(Eigen::MatrixXcd in_matrix)
|
||
{
|
||
return Complex2dB(in_matrix.real().array(), in_matrix.imag().array());
|
||
}
|
||
|
||
Eigen::MatrixXd LAMPTOOLAPI Complex2dB(Eigen::MatrixXd in_matrix_real, Eigen::MatrixXd in_matrix_imag)
|
||
{
|
||
return (in_matrix_real.array().pow(2) + in_matrix_imag.array().pow(2)).log10() * 10;
|
||
}
|
||
|
||
double LAMPTOOLAPI Complex2Amplitude(std::complex<double> sign)
|
||
{
|
||
return sqrt(pow(sign.real(),2)+pow(sign.imag(),2));
|
||
}
|
||
|
||
double LAMPTOOLAPI Complex2phase(std::complex<double> sign)
|
||
{
|
||
//return (sign.real() == 0)*((sign.imag() < 0) * PI / 2)+ (sign.real() != 0)*(atan(sign.imag() / sign.real()) + (sign.real() < 0) * ((sign.imag() <= 0) * PI));
|
||
if (sign.real() != 0) {
|
||
return atan(sign.imag() / sign.real()) + (sign.real() < 0) * ((sign.imag() <= 0) * PI);
|
||
}
|
||
else {
|
||
return (sign.imag() < 0) * PI / 2;
|
||
}
|
||
}
|
||
|
||
Eigen::MatrixXd LAMPTOOLAPI Complex2Amplitude(Eigen::MatrixXcd in_matrix)
|
||
{
|
||
return in_matrix.array().abs().cast<double>().array();
|
||
}
|
||
|
||
Eigen::MatrixXd LAMPTOOLAPI Complex2phase(Eigen::MatrixXcd in_matrix)
|
||
{
|
||
// 复数转相位
|
||
Eigen::MatrixXd result = in_matrix.real().array();
|
||
int rows = in_matrix.rows();
|
||
int cols = in_matrix.cols();
|
||
for (int i = 0; i < rows; i++) { // 可能是性能瓶颈
|
||
for (int j = 0; j < cols; j++) {
|
||
result(i, j) = Complex2phase(in_matrix(i, j));
|
||
}
|
||
}
|
||
|
||
|
||
return result;
|
||
}
|
||
|
||
int LAMPTOOLAPI Complex2dB_DLL(QString out_path, QString in_sar_path)
|
||
{
|
||
return Complex2dB(in_sar_path, out_path);
|
||
}
|
||
|
||
int LAMPTOOLAPI Amplitude2dB_DLL(QString in_tiff, QString out_dB_path)
|
||
{
|
||
return Amplitude2dB(in_tiff, out_dB_path);
|
||
}
|
||
|
||
|