增加了极化转换代码

Release-dev
chenzenghui 2025-05-14 02:56:15 +08:00
parent 79859eb7e7
commit 2a466c16a0
5 changed files with 248 additions and 51 deletions

View File

@ -357,17 +357,20 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row;
cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col;
int64_t pixel_count64 = static_cast<int64_t>(rows_count)* static_cast<int64_t>(cols_count);
//Eigen::MatrixXd datamatrix(rows_count, cols_count);
if (gdal_datatype == GDT_Byte) {
char* temp = new char[rows_count * cols_count];
char* temp = new char[pixel_count64];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
@ -376,14 +379,14 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
delete[] temp;
}
else if (gdal_datatype == GDT_UInt16) {
unsigned short* temp = new unsigned short[rows_count * cols_count];
unsigned short* temp = new unsigned short[pixel_count64];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
@ -391,14 +394,14 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
delete[] temp;
}
else if (gdal_datatype == GDT_Int16) {
short* temp = new short[rows_count * cols_count];
short* temp = new short[pixel_count64];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
@ -406,14 +409,14 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
delete[] temp;
}
else if (gdal_datatype == GDT_UInt32) {
unsigned int* temp = new unsigned int[rows_count * cols_count];
unsigned int* temp = new unsigned int[pixel_count64];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
@ -421,14 +424,14 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
delete[] temp;
}
else if (gdal_datatype == GDT_Int32) {
int* temp = new int[rows_count * cols_count];
int* temp = new int[pixel_count64];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
@ -436,7 +439,7 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
delete[] temp;
}
// else if (gdal_datatype == GDT_UInt64) {
// unsigned long* temp = new unsigned long[rows_count * cols_count];
// unsigned long* temp = new unsigned long[pixel_count64];
// 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];
@ -445,7 +448,7 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
// delete[] temp;
// }
// else if (gdal_datatype == GDT_Int64) {
// long* temp = new long[rows_count * cols_count];
// long* temp = new long[pixel_count64];
// 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];
@ -454,15 +457,15 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
// delete[] temp;
// }
else if (gdal_datatype == GDT_Float32) {
float* temp = new float[rows_count * cols_count];
float* temp = new float[pixel_count64];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
@ -470,15 +473,15 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
delete[] temp;
}
else if (gdal_datatype == GDT_Float64) {
double* temp = new double[rows_count * cols_count];
double* temp = new double[pixel_count64];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i]);
}
@ -487,15 +490,15 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
}
//else if (gdal_datatype == GDT_CFloat32) {
// if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<double>>::value) {
// float* temp = new float[rows_count * cols_count * 2];
// float* temp = new float[pixel_count64 * 2];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
// result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
// result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
// if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
// std::memcpy(result.get(), temp, rows_count * cols_count);
// std::memcpy(result.get(), temp, pixel_count64);
// }
// else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
// long count = rows_count * cols_count;
// long count = pixel_count64;
// for (long i = 0; i < count; i++) {
// result.get()[i] = T(temp[i * 2],
// temp[i * 2 + 1]);
@ -509,15 +512,15 @@ inline std::shared_ptr<T> readDataArr(gdalImage& imgds, long start_row, long sta
//}
//else if (gdal_datatype == GDT_CFloat64 ) {
// if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<double>>::value) {
// double* temp = new double[rows_count * cols_count * 2];
// double* temp = new double[pixel_count64 * 2];
// demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
// result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
// result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
// if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
// std::memcpy(result.get(), temp, rows_count * cols_count);
// std::memcpy(result.get(), temp, pixel_count64);
// }
// else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
// long count = rows_count * cols_count;
// long count = pixel_count64;
// for (long i = 0; i < count; i++) {
// result.get()[i] = T(temp[i * 2],
// temp[i * 2 + 1]);
@ -556,19 +559,21 @@ inline std::shared_ptr<T> readDataArrComplex(gdalImageComplex& imgds, long start
rows_count = start_row + rows_count <= imgds.height ? rows_count : imgds.height - start_row;
cols_count = start_col + cols_count <= imgds.width ? cols_count : imgds.width - start_col;
int64_t pixel_count64 = static_cast<int64_t>(rows_count) * static_cast<int64_t>(cols_count);
//Eigen::MatrixXd datamatrix(rows_count, cols_count);
if (gdal_datatype == GDT_CFloat32) {
if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<double>>::value) {
float* temp = new float[rows_count * cols_count * 2];
float* temp = new float[pixel_count64 * 2];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i * 2],
temp[i * 2 + 1]);
@ -582,15 +587,15 @@ inline std::shared_ptr<T> readDataArrComplex(gdalImageComplex& imgds, long start
}
else if (gdal_datatype == GDT_CFloat64) {
if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<double>>::value) {
double* temp = new double[rows_count * cols_count * 2];
double* temp = new double[pixel_count64 * 2];
demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
result = std::shared_ptr<T>(new T[rows_count * cols_count], delArrPtr);
result = std::shared_ptr<T>(new T[pixel_count64], delArrPtr);
if (method == GDALREADARRCOPYMETHOD::MEMCPYMETHOD) {
std::memcpy(result.get(), temp, rows_count * cols_count);
std::memcpy(result.get(), temp, pixel_count64);
}
else if (method == GDALREADARRCOPYMETHOD::VARIABLEMETHOD) {
long count = rows_count * cols_count;
long count = pixel_count64;
for (long i = 0; i < count; i++) {
result.get()[i] = T(temp[i * 2],
temp[i * 2 + 1]);

View File

@ -1,2 +1,123 @@
#include "GF3CalibrationAndOrthLib.h"
#include "BaseConstVariable.h"
#include "RasterToolBase.h"
#include "ImageOperatorBase.h"
#include <QDebug>
#include <QFile>
#include <memory>
#include "FileOperator.h"
#include "GPUBaseTool.h"
#include "GPUTool.cuh"
GF3CALIBRATIONANDORTHLIB_EXPORT int GF3_Sigma0_HH2VV(QString in_SigmaHHRasterPath, QString in_IncidencAngleRasterPath, QString out_SigmaVVRasterPath)
{
// step 1 检查输入合法性
{
// 1. 检查输入合法性
if (in_SigmaHHRasterPath.isEmpty() || in_IncidencAngleRasterPath.isEmpty() || out_SigmaVVRasterPath.isEmpty())
{
qDebug() << "Input or output file path is empty.";
throw std::invalid_argument("Input or output file path is empty.");
return 1; // 代表输入或输出文件路径为空
}
// 2. 检查输入影像是否存在
if (!QFile::exists(in_SigmaHHRasterPath) || !QFile::exists(in_IncidencAngleRasterPath))
{
qDebug() << "Input file does not exist.";
throw std::runtime_error("Input file does not exist.");
return 2; // 代表输入文件不存在
}
gdalImage inSigmaHHRaster(in_SigmaHHRasterPath);
gdalImage inIncidencAngleRaster(in_IncidencAngleRasterPath);
// 3. 检查输入影像大小是否一致
if (inSigmaHHRaster.width != inIncidencAngleRaster.width || inSigmaHHRaster.height != inIncidencAngleRaster.height)
{
qDebug() << "Input images have different sizes.";
throw std::runtime_error("Input images have different sizes.");
return 3; // 代表输入影像大小不一致
}
// 复制输入文件到输出文件
copyFile(in_SigmaHHRasterPath, out_SigmaVVRasterPath);
// 4. 检查输出文件是否存在
if (!QFile::exists(out_SigmaVVRasterPath))
{
qDebug() << "Output file does not exist.";
throw std::runtime_error("Output file does not exist.");
return 4; // 代表输出文件不存在
}
}
// step 2 执行操作
{
gdalImage inSigmaHHRaster(in_SigmaHHRasterPath);
gdalImage inIncidencAngleRaster(in_IncidencAngleRasterPath);
gdalImage outSigmaVVRaster(out_SigmaVVRasterPath);
long width = inSigmaHHRaster.width;
long height = inSigmaHHRaster.height;
int64_t pixel_count64 = static_cast<int64_t>(height) * static_cast<int64_t>(width);
// 读取输入影像数据
std::shared_ptr<double> sigmaHHRasterData = readDataArr<double>(inSigmaHHRaster, 0, 0, height, width, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> incidenceAngleData = readDataArr<double>(inIncidencAngleRaster, 0, 0, height, width, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
std::shared_ptr<double> outSigmaVVRasterData = readDataArr<double>(outSigmaVVRaster, 0, 0, height, width, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
// 分块执行数据
int64_t block_count = Memory1MB / 8 * 500; // 500M内存 作为数据分块大小
// 创建交换内存
//double* sigmaHHRasterDataPtr_H = (double*)mallocCUDAHost(sizeof(double) * block_count); // 主机内存
double* sigmaHHRasterDataPtr_D = (double*)mallocCUDADevice(sizeof(double) * block_count);// 设备内存
//double* incidenceAngleDataPtr_H = (double*)mallocCUDAHost(sizeof(double) * block_count); // 主机内存
double* incidenceAngleDataPtr_D = (double*)mallocCUDADevice(sizeof(double) * block_count);// 设备内存
//double* outSigmaVVRasterDataPtr_H = (double*)mallocCUDAHost(sizeof(double) * block_count); // 主机内存
double* outSigmaVVRasterDataPtr_D = (double*)mallocCUDADevice(sizeof(double) * block_count);// 设备内存
for (int64_t i = 0; i < pixel_count64; i += block_count)
{
int64_t block_size = (i + block_count) > pixel_count64 ? pixel_count64 - i : block_count;
//采用内存拷贝 从主内存到分块内存
HostToDevice(sigmaHHRasterDataPtr_D, sigmaHHRasterData.get() + i, sizeof(double) * block_size);
HostToDevice(incidenceAngleDataPtr_D, incidenceAngleData.get() + i, sizeof(double) * block_size);
// 执行GPU计算
// 从设备内存到主机内存
DeviceToHost(outSigmaVVRasterDataPtr_D, outSigmaVVRasterData.get() + i, sizeof(double) * block_size);
}
// 强制释放内存
FreeCUDADevice(sigmaHHRasterDataPtr_D);
//FreeCUDAHost(sigmaHHRasterDataPtr_H);
FreeCUDADevice(incidenceAngleDataPtr_D);
//FreeCUDAHost(incidenceAngleDataPtr_H);
FreeCUDADevice(outSigmaVVRasterDataPtr_D);
//FreeCUDAHost(outSigmaVVRasterDataPtr_H);
sigmaHHRasterData.reset();
incidenceAngleData.reset();
outSigmaVVRasterData.reset();
}
return -1;// 代表执行成功
}

View File

@ -15,7 +15,14 @@
/// <summary>
/// 将HH极化的雷达散射系数转换为VV极化的雷达散射系数
/// </summary>
/// <param name="in_SigmaHHRasterPath"></param>
/// <param name="in_IncidencAngleRasterPath"></param>
/// <param name="out_SigmaVVRasterPath"></param>
/// <param name="isVVPolar"></param>
/// <returns></returns>
extern "C" GF3CALIBRATIONANDORTHLIB_EXPORT int GF3_Sigma0_HH2VV(QString in_SigmaHHRasterPath,QString in_IncidencAngleRasterPath,QString out_SigmaVVRasterPath);

View File

@ -13,7 +13,7 @@
#include "GPUTool.cuh"
__device__ __host__ float Computrer_polartionConver_rpolf(float inangle, float alpha)
__device__ __host__ float Computrer_polartionConver_rpol_f(float inangle, float alpha)
{
float tang = 0.0;
if (inangle <= 0.0) {
@ -31,7 +31,7 @@ __device__ __host__ float Computrer_polartionConver_rpolf(float inangle, float a
return rpol;
}
__device__ __host__ double Computrer_polartionConver_rpold(double inangle, double alpha)
__device__ __host__ double Computrer_polartionConver_rpol_d(double inangle, double alpha)
{
double tang = 0.0;
if (inangle <= 0.0) {
@ -48,3 +48,46 @@ __device__ __host__ double Computrer_polartionConver_rpold(double inangle, doubl
double rpol = pow((1.0 + 2.0 * tan_val * tan_val), 2) / pow((1.0 + alpha * tan_val * tan_val), 2);
return rpol;
}
__host__ __device__ float polartionConver_f(float insig, float inangle, float alpha = 1.0, bool isvv = true)
{
float rpol = Computrer_polartionConver_rpol_f(inangle, alpha);
// dB转线性
float insig_linear = powf(10.0, insig / 10.0);
float osig = 0.0;
if (isvv) { // C: VV->HH
osig = insig_linear * rpol;
// 返回线性值
return osig;
}
else { // L: HH->VV
osig = insig_linear / rpol;
// 返回dB值
return 10.0 * log10f(osig);
}
}
__host__ __device__ double polartionConver_d(double insig, double inangle, double alpha = 1.0, bool isvv = true)
{
double rpol = Computrer_polartionConver_rpol_d(inangle, alpha);
// dB转线性
double insig_linear = pow(10.0, insig / 10.0);
double osig = 0.0;
if (isvv) { // C: VV->HH
osig = insig_linear * rpol;
// 返回线性值
return osig;
} else { // L: HH->VV
osig = insig_linear / rpol;
// 返回dB值
return 10.0 * log10(osig);
}
}

View File

@ -16,7 +16,7 @@
/// <param name="inangle">入射角</param>
/// <param name="alpha">转换参数</param>
/// <returns></returns>
extern __device__ __host__ float Computrer_polartionConver_rpolf(float inangle, float alpha = 1);
extern __device__ __host__ float Computrer_polartionConver_rpol_f(float inangle, float alpha = 1);
/// <summary>
/// 计算极化转换系数 (double)
@ -24,7 +24,28 @@ extern __device__ __host__ float Computrer_polartionConver_rpolf(float inangle,
/// <param name="inangle">入射角</param>
/// <param name="alpha">转换参数</param>
/// <returns></returns>
extern __device__ __host__ double Computrer_polartionConver_rpold(double inangle, double alpha = 1);
extern __device__ __host__ double Computrer_polartionConver_rpol_d(double inangle, double alpha = 1);
/// <summary>
/// 后向散射系数变换函数
/// </summary>
/// <param name="insig">sigma0 (dB)</param>
/// <param name="inangle">入射角 (度)</param>
/// <param name="alpha">默认1基尔霍夫模型0.6为汤姆森</param>
/// <param name="isvv">true表示输入为VVfalse为HH</param>
/// <returns></returns>
extern __host__ __device__ double polartionConver_d(double insig, double inangle, double alpha = 1.0, bool isvv = true);
/// <summary>
/// 后向散射系数变换函数
/// </summary>
/// <param name="insig">sigma0 (dB)</param>
/// <param name="inangle">入射角 (度)</param>
/// <param name="alpha">默认1基尔霍夫模型0.6为汤姆森</param>
/// <param name="isvv">true表示输入为VVfalse为HH</param>
/// <returns></returns>
extern __host__ __device__ float polartionConver_f(float insig, float inangle, float alpha = 1.0, bool isvv = true);